Home > atmlab > retrieval > qpack2 > qpack2.m

qpack2

PURPOSE ^

QPACK2 Main function of Qpack2

SYNOPSIS ^

function [L2,X] = qpack2(Q,O,Y)

DESCRIPTION ^

 QPACK2   Main function of Qpack2

    Read qpack2.pdf for an introduction to Qpack2 and user instructions.

    The main task to make OEM type inversions. The measurement data in *Y* are
    inverted by combining the function oem.m and Qarts (interface to arts-2).
    Forward model and retrieval variables are packed into the Qarts structure
    Q. OEM specific parameters are found in O.

    The data fields of the output structure L2 depends on selected retrieval
    quantities and settings in Q. See further *qp2_l2*.

 FORMAT   [L2,X] = qpack2(Q,O,Y)

 IN    Q        Qarts setting structure. See *qarts*.
       O        OEM settings. See *oem*.
       Y        Measurement data. See *qp2_y*.
 OUT   L2       Retrieved data. See *qp2_l2*.
       X        As the output of *oem*, providing additional information.

    The measurement matching the apriori conditions can be calculated through
    this function. This can be useful for test retrievals, or comparison with
    som measurement data. Such calculations are triggered by setting all
    measurement vectors in Y (that is, Y.Y) to be empty ([]). Note that the
    field shall be set to [], and not to {} that inside qpack2 and qarts
    identifies an undefined field.

    The calculated spctre are inserted into Y, and the format is:

 FORMAT   Y = qpack2(Q,O,Y)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

qpack2.m

SOURCE CODE ^

0001 % QPACK2   Main function of Qpack2
0002 %
0003 %    Read qpack2.pdf for an introduction to Qpack2 and user instructions.
0004 %
0005 %    The main task to make OEM type inversions. The measurement data in *Y* are
0006 %    inverted by combining the function oem.m and Qarts (interface to arts-2).
0007 %    Forward model and retrieval variables are packed into the Qarts structure
0008 %    Q. OEM specific parameters are found in O.
0009 %
0010 %    The data fields of the output structure L2 depends on selected retrieval
0011 %    quantities and settings in Q. See further *qp2_l2*.
0012 %
0013 % FORMAT   [L2,X] = qpack2(Q,O,Y)
0014 %
0015 % IN    Q        Qarts setting structure. See *qarts*.
0016 %       O        OEM settings. See *oem*.
0017 %       Y        Measurement data. See *qp2_y*.
0018 % OUT   L2       Retrieved data. See *qp2_l2*.
0019 %       X        As the output of *oem*, providing additional information.
0020 %
0021 %    The measurement matching the apriori conditions can be calculated through
0022 %    this function. This can be useful for test retrievals, or comparison with
0023 %    som measurement data. Such calculations are triggered by setting all
0024 %    measurement vectors in Y (that is, Y.Y) to be empty ([]). Note that the
0025 %    field shall be set to [], and not to {} that inside qpack2 and qarts
0026 %    identifies an undefined field.
0027 %
0028 %    The calculated spctre are inserted into Y, and the format is:
0029 %
0030 % FORMAT   Y = qpack2(Q,O,Y)
0031 
0032 % 2010-05-10   Created by Patrick Eriksson.
0033 
0034 
0035 function [L2,X] = qpack2(Q,O,Y)
0036 
0037 
0038 
0039 %------------------------------------------------------------------------------
0040 %- Check input
0041 %------------------------------------------------------------------------------
0042 if atmlab( 'STRICT_ASSERT' )
0043   rqre_nargin( 3, nargin );
0044   rqre_datatype( Q, @isstruct );
0045   rqre_datatype( O, @isstruct );
0046   rqre_datatype( Y, @isstruct );
0047   qcheck( @qarts, Q );
0048   qcheck( @oem, O );
0049   qcheck( @qp2_y, Y );
0050 end
0051 %
0052 empty_y = checkY( Y );
0053 
0054 
0055 
0056 %------------------------------------------------------------------------------
0057 %- Loop measurements
0058 %------------------------------------------------------------------------------
0059 %
0060 L2 = [];
0061 %
0062 workfolder = create_tmpfolder;
0063 cu = onCleanup( @()delete_tmpfolder( workfolder ) );
0064 %
0065 ny = length(Y);
0066 %
0067 for m = 1 : ny
0068   
0069   %- Simulate a (noise free) measurement if y is empty
0070   if empty_y
0071     %
0072     Qm = qp2_y2Q( Q, Y, m );
0073     %
0074     if m == 1
0075       L2    = Y(m);
0076     else
0077       L2(m) = Y(m);
0078     end
0079     %
0080     if out(1)
0081       fprintf( 'Simulating spectrum %d/%d\n', m, ny );
0082     end
0083     %
0084     L2(m).Y = arts_y( Qm, workfolder );
0085   
0086     
0087   %- Inversion
0088   else
0089     
0090     %- Init retrieval
0091     %
0092     if ny > 1
0093       O.msg = sprintf( 'Inversion case %d (of %d)', m, ny );
0094     end
0095     %
0096     [Qm,Se]     = qp2_y2Q( Q, Y, m );
0097     [Qm,O,R,xa] = arts_oem_init( Qm, O, workfolder );
0098     
0099     %- Sx
0100     if m == 1
0101       %
0102       Sx = arts_sx( Qm, R );
0103       %
0104       si = sqrt( diag( Sx ) );
0105       if max(si) / min(si) <= 1e3
0106         Sxinv    = Sx \ speye(size(Sx));
0107         O.sxnorm = 0;
0108       else
0109         Sxinv    = [];
0110         O.sxnorm = 1;
0111       end
0112     end
0113 
0114     %- Perform retrieval
0115     %
0116     [X,R] = oem( O, Qm, R, @arts_oem, Sx, Se, Sxinv, [], xa, Y(m).Y );
0117     
0118     %- Fill return structure
0119     %
0120     if m == 1
0121       L2    = qp2_l2( Qm, R, xa, O, X, Y, m );
0122     else
0123       L2(m) = qp2_l2( Qm, R, xa, O, X, Y, m );
0124     end
0125   
0126   end % empty_y
0127 end % for m
0128 %
0129 return
0130 %
0131 %------------------------------------------------------------------------------
0132 %------------------------------------------------------------------------------
0133 %------------------------------------------------------------------------------
0134 
0135 
0136 
0137 % Checks if all Y.Y have the same length
0138 %
0139 function empty_y = checkY( Y )
0140 
0141   errid = ['atmlab:' mfilename];
0142   empty_y = isempty( Y(1).Y );
0143   
0144   if atmlab( 'STRICT_ASSERT' )
0145     if empty_y
0146       for i = 2 : length(Y)                                          
0147         if ~isempty(Y(i).Y)                                        
0148           error( errid, ['First spectrum in Y is empty. Then all other '... 
0149                          'spectra must also be empty.'] );                  
0150         end                                                          
0151       end                                                            
0152       
0153     else            
0154       rqre_datatype( Y(1).Y, @istensor1, 'Y(1).Y' );                 
0155       n = length( Y(1).Y );                                          
0156       for i = 2 : length(Y)                                          
0157         rqre_datatype( Y(i).Y, @istensor1, sprintf('Y(%d).Y',i) );   
0158         if length(Y(i).Y) ~= n                                       
0159           error( errid, 'All spectra must have the same length.' );         
0160         end                                                          
0161       end                                                            
0162     end
0163   end
0164 return      
0165 
0166 
0167 
0168 
0169 
0170 function  [Q,Se] = qp2_y2Q( Q, Y, m )
0171 
0172   % If output does not include Se, it assumed that no retrieval will
0173   % be performed (for all m).
0174     
0175   %- Stuff only done for first measurement
0176   %
0177   errid = ['atmlab:' mfilename];
0178   do_checks = atmlab( 'STRICT_ASSERT' );
0179   %
0180   if m == 1
0181 
0182     Q.ATMOSPHERE_DIM = qarts_get( Q.ATMOSPHERE_DIM );  
0183 
0184     if do_checks
0185       if ~qarts_isset(Q.ATMOSPHERE_DIM)                                       
0186         error( errid, 'Qpack2 requires that Q.ATMOSPHERE_DIM is set.' );
0187       end                                                                     
0188       rqre_datatype( Q.ATMOSPHERE_DIM, @istensor0, Q.ATMOSPHERE_DIM );         
0189      %if Q.ATMOSPHERE_DIM ~= 1
0190      %   error( errid, 'Qpack2 requires that Q.ATMOSPHERE_DIM is set to 1.' );
0191      % end
0192       %
0193       if ~qarts_isset(Q.ABSORPTION)                                           
0194         error( errid, 'Qpack2 requires that Q.ABSORPTION is set.' ); 
0195       end                                                                     
0196       rqre_datatype( Q.ABSORPTION, @ischar, Q.ABSORPTION );                   
0197       if strcmp( Q.ABSORPTION, 'CalcTable' )                                  
0198         error( errid, ...
0199                      'Qpack2 is not handling Q.ABSORPTION == ''CalcTable''.' );
0200       end
0201       %
0202       if ~qarts_isset(Q.HSE)                                                  
0203         error( errid, 'Qpack2 requires that Q.HSE is set.' ); 
0204       end                                                                     
0205       rqre_datatype( Q.HSE.ON, @isboolean, Q.HSE.ON );                        
0206       %
0207       if qarts_isset(Q.RAW_ATMOSPHERE)                                        
0208         error( errid, 'Qpack2 requires that Q.RAW_ATMOSPHERE is unset.' ); 
0209       end                                                                     
0210       if ~( isfield(Q.ABS_SPECIES,'ATMDATA') | qarts_isset(Q.VMR_FIELD) )
0211         error( errid, ['Qpack2 requires that Q.VMR_FIELD is set or ',...
0212                        'Q.ABS_SPECIES has the field ATMDATA.'] );
0213       end                                                                     
0214       if ~( isfield(Q.T,'ATMDATA') | qarts_isset(Q.T_FIELD) )
0215         error( errid, ['Qpack2 requires that Q.T_FIELD is set or ',...
0216                        'Q.T has the field ATMDATA.'] );
0217       end                                                                     
0218       if ~( Q.HSE.ON | isfield(Q.Z,'ATMDATA') | qarts_isset(Q.Z_FIELD) )
0219         error( errid, ['Qpack2 requires that either Q.Z_FIELD is set, ',...
0220                        'Q.HSE.ON is true or Q.Z has the field ATMDATA.'] );
0221       end
0222     end
0223     
0224     if nargout < 2
0225       Q.J_DO = 0;
0226     else
0227       Q.J_DO = 1;
0228       Q.TNOISE_C = qarts_get( Q.TNOISE_C ); 
0229       if do_checks
0230         if ~qarts_isset(Q.TNOISE_C)                                          
0231           error( errid, 'Qpack2 requires that Q.TNOISE_C is set.' ); 
0232         end                                                                  
0233         if ~issparse(Q.TNOISE_C)                                             
0234           error( errid, ...
0235                      'Q.TNOISE_C is required to be sparse (to save memory).' );
0236         elseif dimens(Q.TNOISE_C)~=2 || size(Q.TNOISE_C,1)~=size(Q.TNOISE_C,2)
0237           error( errid, 'Q.TNOISE_C must be square.' );
0238         end                                                                 
0239       end  
0240     end % nargout < 2
0241   end % if m==1
0242 
0243   
0244   % Checks of Y(m)
0245   %
0246   if do_checks
0247     rqre_datatype( Y(m).LATITUDE, @istensor0, sprintf('Y(%d).LATITUDE',m) );  
0248     rqre_datatype( Y(m).LONGITUDE, @istensor0, sprintf('Y(%d).LONGITUDE',m) );
0249     rqre_datatype( Y(m).HSE_P, @istensor0, sprintf('Y(%d).HSE_P',m) );        
0250     rqre_datatype( Y(m).HSE_Z, @istensor0, sprintf('Y(%d).HSE_Z',m) );        
0251     rqre_datatype( Y(m).YEAR, @istensor0, sprintf('Y(%d).YEAR',m) );          
0252     rqre_datatype( Y(m).MONTH, @istensor0, sprintf('Y(%d).MONTH',m) );        
0253     rqre_datatype( Y(m).DAY, @istensor0, sprintf('Y(%d).DAY',m) );            
0254     rqre_datatype( Y(m).Z_PLATFORM, @istensor1, ...                           
0255                                              sprintf('Y(%d).Z_PLATFORM',m) ); 
0256     rqre_datatype( Y(m).ZA, @istensor1, sprintf('Y(%d).ZA',m) );
0257     if length(Y(m).ZA) ~= length(Y(m).Z_PLATFORM)
0258       error( errid, ...
0259              'Y(%d).ZA and Y(%d).Z_PLATFORM must have the same length', m, m );
0260     end
0261     rqre_datatype( Y(m).F, @istensor1, sprintf('Y(%d).F',m) ); 
0262     if Q.ATMOSPHERE_DIM == 3
0263       rqre_datatype( Y(m).AA, @istensor1, sprintf('Y(%d).AA',m) );
0264       if length(Y(m).ZA) ~= length(Y(m).AA) 
0265         error( errid, ...
0266                      'Y(%d).AA and Y(%d).ZA must have the same length', m, m );
0267       end 
0268     end
0269   end
0270 
0271   %- Map ATMDATA (if set) to VMR_FIELD, Z_FIELD and T_FIELD, and apply HSE
0272   %
0273   Q.LAT_TRUE  = Y(m).LATITUDE;
0274   Q.LON_TRUE  = Y(m).LONGITUDE;
0275   %
0276   if qarts_isset( Y(m).HOUR )
0277     ho = Y(m).HOUR;
0278   else
0279     ho = 12;
0280   end
0281   if qarts_isset( Y(m).MINUTE )
0282     mi = Y(m).MINUTE;
0283   else
0284     mi = 0;
0285   end
0286   if qarts_isset( Y(m).SECOND )
0287     se = Y(m).SECOND;
0288   else
0289     se = 0;
0290   end
0291   %
0292   if do_checks
0293     rqre_datatype( ho, @istensor0, sprintf('Y(%d).HOUR',m) );
0294     rqre_datatype( mi, @istensor0, sprintf('Y(%d).MINUTE',m) );
0295     rqre_datatype( se, @istensor0, sprintf('Y(%d).SECOND',m) );
0296   end
0297   %
0298   mjd = date2mjd( Y(m).YEAR, Y(m).MONTH, Y(m).DAY, ho, mi, se );
0299   %
0300   if isfield( Q.ABS_SPECIES, 'ATMDATA' )
0301     Q.VMR_FIELD = qarts_vmr_field( Q, mjd, ho );
0302   end
0303   if isfield( Q.T, 'ATMDATA' )
0304     Q.T_FIELD   = qarts_atm_field( Q, 't', mjd, ho );
0305   end
0306   if isfield( Q.Z, 'ATMDATA' )
0307     Q.Z_FIELD   = qarts_atm_field( Q, 'z', mjd, ho );
0308   end
0309   if Q.HSE.ON
0310     Q.Z_FIELD   = qarts_hse( Q, Y(m).HSE_P, Y(m).HSE_Z );
0311   end
0312   %
0313   Q.HSE.P = Y(m).HSE_P;
0314 
0315   
0316   %- Wind components
0317   %
0318   if isfield( Q.WIND_U, 'ATMDATA' )
0319     Q.WIND_U_FIELD   = qarts_atm_field( Q, 'wind_u', mjd, ho );
0320   end
0321   if isfield( Q.WIND_V, 'ATMDATA' )
0322     Q.WIND_V_FIELD   = qarts_atm_field( Q, 'wind_v', mjd, ho );
0323   end
0324   if isfield( Q.WIND_W, 'ATMDATA' )
0325     Q.WIND_W_FIELD   = qarts_atm_field( Q, 'wind_w', mjd, ho );
0326   end
0327 
0328 
0329   %- Magnetic field
0330   %
0331   if isfield( Q.MAG_U, 'ATMDATA' )
0332     Q.MAG_U_FIELD   = qarts_atm_field( Q, 'mag_u', mjd, ho );
0333   end
0334   if isfield( Q.MAG_V, 'ATMDATA' )
0335     Q.MAG_V_FIELD   = qarts_atm_field( Q, 'mag_v', mjd, ho );
0336   end
0337   if isfield( Q.MAG_W, 'ATMDATA' )
0338     Q.MAG_W_FIELD   = qarts_atm_field( Q, 'mag_w', mjd, ho );
0339   end
0340 
0341   
0342   %- Sensor pos and los
0343   %
0344   nza = length( Y(m).ZA );
0345   %
0346   if Q.ATMOSPHERE_DIM == 1
0347     Q.SENSOR_POS  = Y(m).Z_PLATFORM;
0348     Q.SENSOR_LOS  = Y(m).ZA;
0349   elseif Q.ATMOSPHERE_DIM == 2
0350     Q.SENSOR_POS  = [ Y(m).Z_PLATFORM, repmat( [Y(m).LATITUDE], nza, 1 ) ];
0351     Q.SENSOR_LOS  = [ Y(m).ZA ];
0352   elseif Q.ATMOSPHERE_DIM == 3
0353     Q.SENSOR_POS  = [ Y(m).Z_PLATFORM, ...
0354                       repmat( [Y(m).LATITUDE,Y(m).LONGITUDE], nza, 1 ) ];
0355     Q.SENSOR_LOS  = [ Y(m).ZA, Y(m).AA ];
0356   end
0357 
0358 
0359   %- Set a dummy sensor_time, if not already defined
0360   %
0361   if ~qarts_isset( Q.SENSOR_TIME )
0362     Q.SENSOR_TIME = [ 1 : nza ]';
0363   end
0364   
0365   %- Se
0366   %
0367   if nargout == 2
0368     
0369     nnv = size( Y(m).TNOISE, 1 );  % Number of noise values
0370     nf  = length( Y(m).F );
0371     
0372     if do_checks
0373       rqre_datatype( Y(m).TNOISE, @istensor2, sprintf('Y(%d).TNOISE',m) );
0374       if nnv ~= 1  &&  nnv ~= nf
0375         error( errid, ...
0376          'Row size of Y(%d).TNOISE must be 1 or equal to length of Y.F.', m );
0377       end
0378       if size( Y(m).TNOISE, 2 ) ~= nza
0379         error( errid, ...
0380          'Column size of Y(%d).TNOISE must match length of Y(%d).ZA.', m, m );
0381       end
0382     end
0383     
0384     [i,j,s] = find( Q.TNOISE_C );
0385     ni      = length( i );
0386     [i2,j2,s2] = deal( zeros( ni*nza, 1 ) );
0387 
0388     for z = 1 : nza
0389       ind     = (z-1)*ni + (1:ni);
0390       i2(ind) = (z-1)*nf + i;
0391       j2(ind) = (z-1)*nf + j;
0392       %
0393       if nnv == 1
0394         s2(ind) = Y(m).TNOISE(1,z)^2 * s; 
0395       else
0396         for k = 1 : ni
0397           s2(ind(1)-1+k) = prod( Y(m).TNOISE([i(k) j(k)],z) ) * s(k);
0398         end
0399       end
0400     end
0401     %
0402     Se = sparse(i2,j2,s2);  
0403   end 
0404 return
0405

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