Description of qpack2

0001 % QPACK2_DEMO   Demonstration of the Qpack2 retrieval system
0002 %
0003 %    The main features of Qpack2 are demonstrated. The example case is airborne
0004 %    measurements of ozone at 110.8 GHz. Synthetic measurement data are
0005 %    generated internally. See the code and internal comments for details.
0006 %
0007 %    Everything is here put into a single file. For practical retrievals it
0008 %    is probably better to put the definitions of Q (together with O?) in a
0009 %    separate function (i.e. [Q,O] = q_mycase). A function to import
0010 %    measurement data into the "Y format" (see *qp2_y*) is needed. The
0011 %    retrieval result is returned in the L2 format produced by *qp2_l2*.
0012 %
0013 %    This script focuses on giving an introduction, indicating different
0014 %    retrieval units and retrieval of other variables. See also *qpack2_demo2*.
0015 %
0016 % FORMAT   L2 = qpack2_demo
0017 %
0018 % OUT      L2   L2 data output from *qpack2*.
0019 
0020 % 2010-05-12   Created by Patrick Eriksson.
0021 
0022 function L2 = qpack2_demo
0023   
0024 errid = ['atmlab:' mfilename]; 
0025 
0026 %- Qarts settings
0027 %
0028 Q    = q_demo;           % Local file, found below
0029 
0030 
0031 %- Measurement data
0032 %
0033 Y = y_demo( Q );         % Local file, found below
0034 
0035 
0036 %- Check that all frequencies are OK
0037 %
0038 if ~qp2_check_f( Q, Y, 1e3 );
0039   error( errid, ...
0040             'Some mismatch between Q.F_BACKEND and frequencies of spectra.' );
0041 end
0042 
0043 
0044 %- OEM variables
0045 %
0046 O = qp2_l2( Q );         % This ensures that OEM returns the variables needed
0047 %                          to fill the L2 structure, as defined in Q
0048 O.linear = false;
0049 %
0050 if ~O.linear 
0051   O.itermethod = 'GN';
0052   O.stop_dx    = 0.01;
0053   O.maxiter    = 5;
0054 end
0055 
0056 
0057 %- Make inversion
0058 %
0059 L2 = qpack2( Q, O, Y );
0060 
0061 
0062 %- Plot, if no output argument
0063 %
0064 if ~nargout
0065   
0066   % Profiles
0067   figure(1),clf
0068   plot( L2(1).species1_x*1e6, p2z_simple(L2(1).species1_p)/1e3, 'b', ...
0069         L2(2).species1_x*1e6, p2z_simple(L2(2).species1_p)/1e3, 'r', ...
0070         L2(1).species1_xa*1e6, p2z_simple(L2(1).species1_p)/1e3, 'k-' );
0071   xlabel( 'Ozone [VMR]' );
0072   ylabel( 'Approximate altitude [km]' );
0073   axis( [ 0 8 10 90 ] )
0074   legend( 'Retrieval 1', 'Retrieval 2', 'True and a priori' );
0075 
0076   % Spectra
0077   figure(2),clf
0078   h=plot( L2(1).f/1e9, L2(1).y, 'k.', L2(2).f/1e9, L2(2).y, 'k.', ...
0079           L2(1).f/1e9, L2(1).yf, 'b-', L2(2).f/1e9, L2(2).yf, 'r-', ...
0080           L2(1).f/1e9, L2(1).bl, 'b-.', L2(2).f/1e9, L2(2).bl, 'r-.' );
0081   xlabel( 'Frequency [GHz]' );
0082   ylabel( 'Tb [K]' );
0083   axis( [ min(L2(1).f/1e9) max(L2(1).f/1e9) 0 18 ] )
0084   legend( h(3:end), 'Fitted 1', 'Fitted 2', 'Baseline 1', 'Baseline 2' );
0085 end
0086 
0087 return
0088 %-----------------------------------------------------------------------------
0089 %-----------------------------------------------------------------------------
0090 %-----------------------------------------------------------------------------
0091 
0092 
0093 
0094 
0095 
0096 function Q = q_demo
0097  
0098 errid = ['atmlab:' mfilename];
0099 %- Atmlab settings
0100 %
0101 arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
0102 arts_includes     = atmlab( 'ARTS_INCLUDES' );
0103 if isnan( arts_xmldata_path ) 
0104   error( errid,'You need to set ARTS_XMLDATA_PATH to run this exmaple.' );
0105 end
0106 if isnan( arts_includes )
0107   error( erird,'You need to ARTS_INCLUDES to run this example.' );
0108 end                                                      
0109 %
0110 fascod = fullfile( arts_xmldata_path, 'planets', 'Earth', 'Fascod' );
0111 
0112 
0113 %- Init Q
0114 %
0115 Q = qarts;
0116 %
0117 Q.INCLUDES            = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0118                           fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0119                           fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0120                           fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0121 Q.ATMOSPHERE_DIM      = 1;
0122 Q.STOKES_DIM          = 1;
0123 Q.J_DO                = true;
0124 Q.CLOUDBOX_DO         = false;
0125 
0126 
0127 %= Define agendas
0128 %
0129 % Here we do it by using the predefined agenda templates
0130 %   (found in arts/controlfiles/general/agendas.arts)
0131 % This works only if the pre-defined agenda is names following the pattern:
0132 %   name_of_agenda__(Something)
0133 %
0134 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0135 Q.PPATH_STEP_AGENDA          = { 'ppath_step_agenda__GeometricPath'    };
0136 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0137 Q.IY_SPACE_AGENDA            = { 'iy_space_agenda__CosmicBackground'   };
0138 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0139 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0140 
0141 
0142 %- Radiative transfer
0143 %
0144 Q.IY_UNIT             = 'RJBT'; 
0145 Q.YCALC_WSMS          = { 'yCalc' };
0146 %
0147 Q.PPATH_LMAX          = 250;
0148 
0149 
0150 %- Surface
0151 %
0152 Q.Z_SURFACE           = 1e3;          % Just a dummy value. A 10 km
0153                                       % observation altitude is assumed here
0154 
0155 
0156 %- Absorption
0157 %
0158 Q.ABS_LINES           = fullfile( atmlab_example_data, 'o3line111ghz' );
0159 Q.ABS_LINES_FORMAT    = 'Arts';
0160 %
0161 Q.ABSORPTION          = 'OnTheFly';
0162 Q.ABS_NLS             = [];
0163 
0164 %- Pressure grid
0165 %
0166 z_toa                 = 95e3;
0167 %
0168 Q.P_GRID              = z2p_simple( Q.Z_SURFACE-1e3 : 2e3 : z_toa )';
0169 
0170 
0171 
0172 %- Frequency, spectrometer and pencil beam antenna
0173 %
0174 % The hypothetical spectrometer has rectangular response functions
0175 %
0176 Q.F_GRID              = qarts_get( fullfile( atmlab_example_data , ...
0177                                                       'f_grid_111ghz.xml' ) );
0178 %
0179 H                     = qartsSensor;
0180 %
0181 H.SENSOR_NORM         = true;
0182 %
0183 H.BACKEND_DO          = true;
0184 df                    = 0.5e6;
0185 H.F_BACKEND           = ( min(Q.F_GRID)+df : df : max(Q.F_GRID)-df )';
0186 %
0187 B.name                = 'Spectrometer channel response function';
0188 B.gridnames           = { 'Frequency' };
0189 B.grids               = { [-df/2 df/2] };
0190 B.dataname            = 'Response';
0191 B.data                = [1 1];
0192 %
0193 H.BACKEND_CHANNEL_RESPONSE{1} = B;
0194 clear B
0195 %
0196 Q.SENSOR_DO           = true;
0197 Q.SENSOR_RESPONSE     = H;
0198 %
0199 Q.ANTENNA_DIM         = 1;
0200 Q.MBLOCK_ZA_GRID      = 0;
0201 
0202 
0203 
0204 %- Correlation of thermal noise
0205 %
0206 f              = H.F_BACKEND;
0207 cl             = 1.4 * ( f(2) - f(1) );
0208 cfun           = 'gau';
0209 cco            = 0.05;
0210 %
0211 Q.TNOISE_C     = covmat1d_from_cfun( f, [], cfun, cl, cco );
0212 %
0213 clear H f
0214 
0215 
0216 %- Define L2 structure (beside retrieval quantities below)
0217 %
0218 Q.L2_EXTRA     = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', ...
0219                    'mresp', 'A', 'e', 'eo', 'es', 'date' };
0220 
0221 
0222 %- Temperature
0223 %
0224 Q.T.RETRIEVE   = false;
0225 Q.T.ATMDATA    = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
0226                         'msis90', 'msis90.t.xml' ), 'Temperature', 't_field' );
0227 
0228 %- Determine altitudes through HSE
0229 %
0230 Q.HSE.ON       = true;
0231 Q.HSE.P        = Q.P_GRID(1);
0232 Q.HSE.ACCURACY = 0.1;
0233 
0234 
0235 %- Species
0236 
0237 % Ozone, only species is retrieved here
0238 Q.ABS_SPECIES(1).TAG      = { 'O3' };
0239 Q.ABS_SPECIES(1).RETRIEVE = true;
0240 Q.ABS_SPECIES(1).L2       = true;
0241 Q.ABS_SPECIES(1).GRIDS    = { Q.P_GRID, [], [] };
0242 Q.ABS_SPECIES(1).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0243       'midlatitude-winter', 'midlatitude-winter.O3.xml' ), 'O3', 'vmr_field' );
0244 %
0245 % If you don't apply a min value (by MINMAX), you could need to active this:
0246 %Q.VMR_NEGATIVE_OK = true;
0247 %
0248 % For demonstration, setting for several units are provided:
0249 switch 1
0250  case 1  % Constant VMR
0251   Q.ABS_SPECIES(1).UNIT     = 'vmr';
0252   Q.ABS_SPECIES(1).SX       = ...
0253                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 1.5e-6, ...
0254                                                'lin', 0.2, 0.00, @log10 ) + ...
0255                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.3e-6, ...
0256                                                'lin', 0.5, 0.00, @log10 );
0257    Q.ABS_SPECIES(1).MINMAX  = 1e-12;
0258  case 2 % Constant rel
0259   Q.ABS_SPECIES(1).UNIT     = 'rel';
0260   Q.ABS_SPECIES(1).SX       = ...
0261                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0262                                                'lin', 0.2, 0.00, @log10 ) + ...
0263                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0264                                                'lin', 0.5, 0.00, @log10 );
0265   Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0266  case 3 % Mimic case 2 in vmr
0267   Q.ABS_SPECIES(1).UNIT     = 'vmr';
0268   Q.ABS_SPECIES(1).SX       = ...
0269                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, ...
0270                                        [ Q.ABS_SPECIES(1).ATMDATA.GRID1,...
0271                                    0.5 * Q.ABS_SPECIES(1).ATMDATA.DATA ],...
0272                                                'lin', 0.2, 0.00, @log10 ) + ...
0273                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, ...
0274                                        [ Q.ABS_SPECIES(1).ATMDATA.GRID1, ...
0275                                    0.1 * Q.ABS_SPECIES(1).ATMDATA.DATA ],...
0276                                                'lin', 0.5, 0.00, @log10 );
0277   Q.ABS_SPECIES(1).MINMAX  = 1e-12;
0278  case 4 % Constant logrel
0279   Q.ABS_SPECIES(1).UNIT     = 'logrel';
0280   Q.ABS_SPECIES(1).SX       = ...
0281                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0282                                                'lin', 0.2, 0.00, @log10 ) + ...
0283                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0284                                                'lin', 0.5, 0.00, @log10 );
0285   Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0286 end  
0287   
0288   
0289   
0290 %- Water
0291 %
0292 % This generates no absorption, as linefile has no H2O lines
0293 %
0294 Q.ABS_SPECIES(2).TAG      = { 'H2O' };
0295 Q.ABS_SPECIES(2).RETRIEVE = false;
0296 Q.ABS_SPECIES(2).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0297     'midlatitude-winter', 'midlatitude-winter.H2O.xml' ), 'H2O', 'vmr_field' );
0298 
0299 
0300 
0301 %- Wind
0302 %
0303 % Here demonstrated for v-component, that should the component of main
0304 % concern for ground-based instruments. This component can be seen as the
0305 % part of the horizontal wind going along the viewing direction, where
0306 % positive values mean a movement away from the sensor.
0307 %
0308 Q.WIND_V_FIELD    = [];      % This is a short-cut for zero wind.
0309 %
0310 Q.WIND_V.RETRIEVE = false;   % Set to true to activate retrieval
0311 Q.WIND_V.L2       = true;
0312 Q.WIND_V.GRIDS    = { Q.P_GRID(1:2:end), [], [] };
0313 Q.WIND_V.SX       = covmat1d_from_cfun( Q.WIND_V(1).GRIDS{1}, 40, ...
0314                                                     'lin', 0.5, 0.00, @log10 );
0315 if Q.WIND_V.RETRIEVE
0316   Q.WSMS_AT_START{end+1} = 'IndexSet(abs_f_interp_order,3)';
0317 end
0318 
0319 
0320 
0321 %- Pointing
0322 %
0323 % Here just included for testing purposes, of little interest for ground-based
0324 % spectrometers.
0325 %
0326 Q.POINTING.RETRIEVE       = false;
0327 Q.POINTING.DZA            = 0.1; 
0328 Q.POINTING.POLY_ORDER     = 0;
0329 Q.POINTING.CALCMODE       = 'recalc';
0330 Q.POINTING.SX             = 1;  
0331 Q.POINTING.L2             = true;
0332 
0333 
0334 %- Frequency shift
0335 %
0336 Q.FSHIFTFIT.RETRIEVE      = true;
0337 Q.FSHIFTFIT.DF            = 25e3; 
0338 Q.FSHIFTFIT.SX            = 50e3^2;  
0339 Q.FSHIFTFIT.L2            = true;
0340 
0341 %- Frequency stretch
0342 %
0343 Q.FSTRETCHFIT.RETRIEVE    = false;   % Set to true to activate retrieval
0344 Q.FSTRETCHFIT.DF          = 25e3; 
0345 Q.FSTRETCHFIT.SX          = 50e3^2;  
0346 Q.FSTRETCHFIT.L2          = true;
0347 
0348 
0349 %- Polyfit
0350 %
0351 % A polynomial of order 3 is used for "baseline fit".
0352 %
0353 Q.POLYFIT.RETRIEVE        = true;
0354 Q.POLYFIT.ORDER           = 3;
0355 Q.POLYFIT.L2              = true;
0356 Q.POLYFIT.SX0             = 1^2; 
0357 Q.POLYFIT.SX1             = 0.5^2; 
0358 Q.POLYFIT.SX2             = 0.2^2;
0359 Q.POLYFIT.SX3             = 0.1^2; 
0360 
0361 
0362 %- Sinefit
0363 %
0364 % Here demonstrated with two period lengths
0365 %
0366 Q.SINEFIT.RETRIEVE        = false;   % Set to true to activate retrieval
0367 Q.SINEFIT.PERIODS         = [ 75e3 200e3 ]';
0368 Q.SINEFIT.L2              = true;
0369 Q.SINEFIT.SX1             = 0.2^2; 
0370 Q.SINEFIT.SX2             = 0.4^2; 
0371 
0372 return
0373 %-----------------------------------------------------------------------------
0374 %-----------------------------------------------------------------------------
0375 %-----------------------------------------------------------------------------
0376 
0377 
0378 
0379 
0380 
0381 function Y = y_demo(Q)
0382 
0383 % The data should be loaded from one or several files, but are here genereted
0384 % by a forward model call to show how qpack2 can also be used to generate
0385 % simulaled measurements (matching a priori assumptions).
0386 
0387 % The simulated data model airborn measurements at two different zenith
0388 % angles, from two nearby positions.
0389   
0390 % Init Y
0391 %
0392 Y = qp2_y;
0393 
0394 % Set a date
0395 %
0396 Y.YEAR  = 2008;
0397 Y.MONTH = 2;
0398 Y.DAY   = 25;
0399 
0400 % Lat / lon
0401 %
0402 Y.LATITUDE  = 45;
0403 Y.LONGITUDE = exp(1);
0404 
0405 % An airborn measurement assumed here
0406 %
0407 Y.Z_PLATFORM = 10.5e3;
0408 Y.ZA         = 50;
0409 
0410 % Reference point for hydrostatic equilibrium
0411 %
0412 Y.HSE_P = 100e2;
0413 Y.HSE_Z = 16e3;
0414 
0415 % Set backend frequencies
0416 %
0417 Y.F = Q.SENSOR_RESPONSE.F_BACKEND;
0418 
0419 % Thermal noise standard deviation
0420 %
0421 Y.TNOISE = 0.1;
0422 % To test varying noise
0423 %Y.TNOISE = linspace( 0.03, 0.07, length(Y.F) )';
0424 
0425 % Simulate a measurement
0426 %
0427 Y.Y = [];                           % A flag to tell qpack2 to calculate the
0428 %                                     spectrum ({} signifies undefined!).
0429 
0430 % Add a second measurement
0431 %
0432 Y(2) = Y(1);
0433 %
0434 Y(2).LONGITUDE = pi;
0435 Y(2).ZA        = 45;
0436 
0437 
0438 % Calculate simulated spectra
0439 %
0440 Y = qpack2( Q, oem, Y );            % Dummy oem structure OK here
0441 
0442 % Add thermal noise
0443 %
0444 % The correlation specified in Q is included
0445 %
0446 for i = 1 : length(Y)
0447   Y(i).Y = Y(i).Y + Y(i).TNOISE .* make_noise(1,Q.TNOISE_C);
0448 end
0449 
0450 % Add a constant "baseline shift" for measurement 2
0451 %
0452 Y(2).Y = Y(2).Y + 1;
0453
qpack2_demo

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

DOWNLOAD

SOURCE CODE
