Home > atmlab > arts > arts_x2QR.m

arts_x2QR

PURPOSE ^

ARTS_x2QR Maps an x state to Q and R fields

SYNOPSIS ^

function [Q,R] = arts_x2QR(Q,R,x)

DESCRIPTION ^

 ARTS_x2QR   Maps an x state to Q and R fields

    The function updates *Q* and *R* to match the present state of *x*.

    No forward model calculations are performed.

 FORMAT   [Q,R] = arts_x2QR( Q, R, x )
        
 OUT   Q      Modified Q structure.
       R      Modified R structure.
 IN    Q      Qarts structure.
       R      Retrieval data structure. 
       x      State vector.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

arts_x2QR.m

SOURCE CODE ^

0001 % ARTS_x2QR   Maps an x state to Q and R fields
0002 %
0003 %    The function updates *Q* and *R* to match the present state of *x*.
0004 %
0005 %    No forward model calculations are performed.
0006 %
0007 % FORMAT   [Q,R] = arts_x2QR( Q, R, x )
0008 %
0009 % OUT   Q      Modified Q structure.
0010 %       R      Modified R structure.
0011 % IN    Q      Qarts structure.
0012 %       R      Retrieval data structure.
0013 %       x      State vector.
0014 
0015 % 2010-01-11   Started by Patrick Eriksson.
0016 
0017 
0018 function [Q,R] = arts_x2QR(Q,R,x)
0019 
0020 
0021 %- General variables
0022 %
0023 nq       = length( R.jq );                  % Number of retrieval quantities.
0024 %
0025 i_asj    = find([Q.ABS_SPECIES.RETRIEVE]);  % Index of retrieval absorption
0026 %                                             species.
0027 any_gas  = false;                           % Any gas species among rq
0028 %
0029 bl       = 0;                               % Sum of all baseline terms
0030 %
0031 j_loaded = false;                           % To avoid repeated loading of J
0032 %
0033 do_sensor = false;                          % Recalculation of sensor needed
0034 
0035 
0036 %- Limit to set min and max values
0037 %
0038 ind    = find( x < R.minmax(:,1) );
0039 x(ind) = R.minmax(ind,1);
0040 ind    = find( x > R.minmax(:,2) );
0041 x(ind) = R.minmax(ind,2);
0042   
0043 
0044 %- Loop retrieval quantities
0045 %
0046 for i = 1 : nq
0047 
0048   ind = R.ji{i}{1} : R.ji{i}{2};
0049 
0050   switch R.jq{i}.maintag
0051 
0052    case 'Absorption species'   %----------------------------------------------
0053     %
0054     if ~any_gas
0055       vmr = R.vmr_field0;    % Will be filled with latest VMR for all species
0056     end                      % that are retrieved
0057     any_gas   = true;
0058     ig        = i_asj(i);    % Gas species index
0059     %
0060     X = x2field( Q.ATMOSPHERE_DIM, R.jq{i}.grids, x(ind) );
0061     X = arts_regrid( Q.ATMOSPHERE_DIM, R.jq{i}.grids, X, Q );
0062     %
0063     if strcmp( R.jq{i}.mode, 'rel' )
0064       vmr(ig,:,:,:) = getdims( vmr(ig,:,:,:), 2:Q.ATMOSPHERE_DIM+1 ) .* X;
0065     elseif strcmp( R.jq{i}.mode, 'vmr' )
0066       vmr(ig,:,:,:) = X;
0067     elseif strcmp( R.jq{i}.mode, 'nd' )
0068       vmr(ig,:,:,:) = nd2vmr( X, repmat(R.p_grid,[1 size2(X,2:3)]), ...
0069                              qarts_get(fullfile(R.workfolder,'t_field.xml')) );
0070     elseif strcmp( R.jq{i}.mode, 'logrel' )
0071       vmr(ig,:,:,:) = getdims( vmr(ig,:,:,:), 2:Q.ATMOSPHERE_DIM+1 ) .* exp(X);
0072     else 
0073       assert( false );
0074     end
0075     %
0076     clear X
0077     
0078    case 'Atmospheric temperatures'   %----------------------------------------
0079     %
0080     X = x2field( Q.ATMOSPHERE_DIM, R.jq{i}.grids, x(ind) );
0081     X = arts_regrid( Q.ATMOSPHERE_DIM, R.jq{i}.grids, X, Q );
0082     xmlStore( fullfile(R.workfolder,'t_field.xml'), X, 'Tensor3' );
0083     %
0084     clear X
0085 
0086    case 'Wind'   %------------------------------------------------------------
0087     %
0088     X = x2field( Q.ATMOSPHERE_DIM, R.jq{i}.grids, x(ind) );
0089     X = arts_regrid( Q.ATMOSPHERE_DIM, R.jq{i}.grids, X, Q );
0090     if strcmp( R.jq{i}.subtag, 'u' )
0091       xmlStore( fullfile(R.workfolder,'wind_u_field.xml'), X, 'Tensor3' );
0092     elseif strcmp( R.jq{i}.subtag, 'v' )
0093       xmlStore( fullfile(R.workfolder,'wind_v_field.xml'), X, 'Tensor3' );
0094     elseif strcmp( R.jq{i}.subtag, 'w' )
0095       xmlStore( fullfile(R.workfolder,'wind_w_field.xml'), X, 'Tensor3' );
0096     else
0097       assert( false );
0098     end
0099     %
0100     clear X
0101     
0102    case 'Sensor pointing'   %-------------------------------------------------
0103     %
0104     assert( strcmp( R.jq{i}.subtag, 'Zenith angle off-set' ) );
0105     %
0106     sensor_los      = R.los;
0107     sensor_los(:,1) = R.los(:,1) + arts_polybasis_func( Q.SENSOR_TIME, ...
0108                                                     R.jq{i}.grids{1} )* x(ind); 
0109     xmlStore( fullfile(R.workfolder,'sensor_los.xml'), sensor_los, 'Matrix' );
0110    
0111    case 'Frequency'   %-------------------------------------------------------
0112     %
0113     assert( length(ind) == 1 );
0114     if ~do_sensor
0115       Q.SENSOR_RESPONSE = R.sensor_response;
0116     end    
0117     %
0118     fb = qarts_get(Q.SENSOR_RESPONSE.F_BACKEND);
0119     %
0120     if strcmp( R.jq{i}.subtag, 'Shift' )
0121       df = x(ind(1));
0122     elseif strcmp( R.jq{i}.subtag, 'Stretch' )
0123       df = arts_polybasis_func( fb, 1 ) * x(ind(1));
0124     else
0125       assert( false );
0126     end
0127     %
0128     Q.SENSOR_RESPONSE.F_BACKEND = fb + df;
0129     do_sensor                   = true;
0130     clear   fb;
0131       
0132    case 'Polynomial baseline fit'   %-----------------------------------------
0133     %
0134     if ~j_loaded
0135       J        = xmlLoad( fullfile( R.workfolder, 'jacobian.xml' ) );
0136       j_loaded = 1;
0137     end
0138     %
0139     bl = bl + J(:,ind) * x(ind);
0140 
0141    case 'Sinusoidal baseline fit'   %-----------------------------------------
0142     %
0143     if ~j_loaded
0144       J        = xmlLoad( fullfile( R.workfolder, 'jacobian.xml' ) );
0145       j_loaded = 1;
0146     end
0147     %
0148     bl = bl + J(:,ind) * x(ind);
0149     
0150     otherwise   %--------------------------------------------------------------
0151       error('Unknown retrieval quantitity.'); 
0152   end 
0153 end 
0154 
0155 
0156 %- Data to save ?
0157 %
0158 if any_gas
0159   xmlStore( fullfile(R.workfolder,'vmr_field.xml'), vmr, 'Tensor4' );
0160 end
0161 
0162 
0163 %- Update sensor characteristics
0164 %
0165 if do_sensor
0166   arts_sensor( Q, R.workfolder );
0167 end
0168 
0169 
0170 %- Total baseline
0171 %
0172 R.bl = bl;
0173 
0174 return
0175 
0176 
0177 
0178 %--- Internal sub-functions ----------------------------------------------------
0179 
0180 function X = x2field( dim, rgrids, x )
0181   %
0182   map = ones( 1, 3 );
0183   for i = 1 : dim
0184     map(i) = length( rgrids{i} );
0185   end
0186   X = reshape( x, map );
0187 return

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