Home > atmlab > arts > arts_oem.m

arts_oem

PURPOSE ^

ARTS_OEM Interface between *oem* and the ARTS forward model

SYNOPSIS ^

function [R,y,J] = arts_oem( Q, R, x, iter )

DESCRIPTION ^

 ARTS_OEM   Interface between *oem* and the ARTS forward model

    See *arts_oem_init* for usage of this function.

 FORMAT   [R,y,J] = arts_oem( Q, R, x, iter )
        
 OUT   R      Retrieval data structure. See above.
       y      Calculated measurement vector.
       J      Jacobian.
 IN    Q      Qarts structure, as returned by *arts_oem_init*.
       R      Retrieval data structure. See above.
       x      State vector.
       iter   Iteration number. The case of x = xa is treated as iter=1. If 
              another initial linearisation point is used, iter must be set
              to >= 2.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

arts_oem.m

SOURCE CODE ^

0001 % ARTS_OEM   Interface between *oem* and the ARTS forward model
0002 %
0003 %    See *arts_oem_init* for usage of this function.
0004 %
0005 % FORMAT   [R,y,J] = arts_oem( Q, R, x, iter )
0006 %
0007 % OUT   R      Retrieval data structure. See above.
0008 %       y      Calculated measurement vector.
0009 %       J      Jacobian.
0010 % IN    Q      Qarts structure, as returned by *arts_oem_init*.
0011 %       R      Retrieval data structure. See above.
0012 %       x      State vector.
0013 %       iter   Iteration number. The case of x = xa is treated as iter=1. If
0014 %              another initial linearisation point is used, iter must be set
0015 %              to >= 2.
0016 
0017 % 2010-01-11   Started by Patrick Eriksson.
0018 
0019 
0020 function [R,y,J] = arts_oem( Q, R, x, iter )
0021 
0022 
0023 % Checks done *arts_oem_init* are not repeated here
0024 
0025 
0026 %--- Map x to variables ----------------------------------------------
0027 
0028 if iter > 1
0029   [Q,R] = arts_x2QR( Q, R, x );
0030 end
0031 
0032 
0033 %--- Run ARTS --------------------------------------------------------
0034 
0035 if nargout == 3
0036   do_j  = 1;
0037   cfile = R.cfile_yj;
0038 else
0039   do_j  = 0;
0040   cfile = R.cfile_y;
0041 end
0042 %
0043 arts( cfile );
0044 %
0045 y = xmlLoad( fullfile( R.workfolder, 'y.xml' ) );
0046 %
0047 if do_j  
0048   
0049   % Load Jacobian
0050   J = xmlLoad( fullfile( R.workfolder, 'jacobian.xml' ) );
0051 
0052   % Jacobian calculated for x, but for "rel" it should be with respect to xa:
0053   % (as arts takes x as xa, no scaling needed for "logrel", and no scaling
0054   %  needed for first calculation)
0055   if iter > 1 & ~isempty( R.i_rel )
0056     % The solution below sets an upper limit on the size of J at roughly 33%
0057     % of free memory.
0058     % J(:,R.i_rel) =  J(:,R.i_rel) ./ repmat( x(R.i_rel)', size(J,1), 1 );
0059     % To avoid hitting the memory roof, we do this instead column-by-column.
0060     % And this turned out to be a faster option, at least for large J!
0061     for i = vec2row(R.i_rel)
0062       J(:,i) = J(:,i) / x(i);
0063     end
0064   end  
0065 end
0066 
0067 
0068 %- Add baseline
0069 %
0070 if iter > 1
0071   y = y + R.bl;
0072 end
0073

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005