Home > atmlab > demos > qpack2_demo2.m

qpack2_demo2

PURPOSE ^

QPACK2_DEMO2 Another demonstration of the Qpack2 retrieval system

SYNOPSIS ^

function L2 = qpack2_demo2

DESCRIPTION ^

 QPACK2_DEMO2   Another demonstration of the Qpack2 retrieval system

    Similar to *qpack2_demo*, but illustrates rudimentary treatment
    tropospheric attenuation and shows how the polarisation response of the
    sensor is included.

 FORMAT   L2 = qpack2_demo2
        
 OUT      L2   L2 data output from *qpack2*.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

qpack2_demo2.m

SOURCE CODE ^

0001 % QPACK2_DEMO2   Another demonstration of the Qpack2 retrieval system
0002 %
0003 %    Similar to *qpack2_demo*, but illustrates rudimentary treatment
0004 %    tropospheric attenuation and shows how the polarisation response of the
0005 %    sensor is included.
0006 %
0007 % FORMAT   L2 = qpack2_demo2
0008 %
0009 % OUT      L2   L2 data output from *qpack2*.
0010 
0011 % 2013-06-17   Created by Patrick Eriksson.
0012 
0013 function L2 = qpack2_demo2
0014   
0015 errid = ['atmlab:' mfilename]; 
0016 
0017 %- Qarts settings
0018 %
0019 Q    = q_demo;           % Local file, found below
0020 
0021 
0022 %- Measurement data
0023 %
0024 Y = y_demo( Q );         % Local file, found below
0025 
0026 
0027 
0028 %- OEM variables
0029 %
0030 O = qp2_l2( Q );         % This ensures that OEM returns the variables needed
0031 %                          to fill the L2 structure, as defined in Q
0032 O.linear = false;
0033 %
0034 if ~O.linear 
0035   O.itermethod       = 'GN';
0036   O.stop_dx          = 0.01;
0037   O.maxiter          = 10;
0038 end
0039 
0040 
0041 %- Make inversion
0042 %
0043 L2 = qpack2( Q, O, Y );
0044 
0045 
0046 %- Convert from rel to vmr
0047 %
0048 L2 = qp2_rel2vmr( L2 );
0049 
0050 
0051 %- Plot, if no output argument
0052 %
0053 if ~nargout
0054   
0055   % Profiles
0056   figure(1),clf
0057   plot( L2(1).species1_x*1e6, p2z_simple(L2(1).species1_p)/1e3, 'b', ...
0058         L2(2).species1_x*1e6, p2z_simple(L2(2).species1_p)/1e3, 'r', ...
0059         L2(1).species1_xa*1e6, p2z_simple(L2(1).species1_p)/1e3, 'k-' );
0060   xlabel( 'Ozone [VMR]' );
0061   ylabel( 'Approximate altitude [km]' );
0062   axis( [ 0 10 0 90 ] )
0063   legend( 'Retrieval 1', 'Retrieval 2', 'True and a priori' );
0064 
0065   % Spectra
0066   figure(2),clf
0067   h=plot( L2(1).f/1e9, L2(1).y, 'k.', L2(2).f/1e9, L2(2).y, 'k.', ...
0068           L2(1).f/1e9, L2(1).yf, 'b-', L2(2).f/1e9, L2(2).yf, 'r-' );
0069   xlabel( 'Frequency [GHz]' );
0070   ylabel( 'Tb [K]' );
0071 %  axis( [ min(L2(1).f/1e9) max(L2(1).f/1e9) 0 18 ] )
0072   legend( h(3:end), 'Fitted 1', 'Fitted 2' );
0073 end
0074 
0075 return
0076 %-----------------------------------------------------------------------------
0077 %-----------------------------------------------------------------------------
0078 %-----------------------------------------------------------------------------
0079 
0080 
0081 
0082 
0083 
0084 function Q = q_demo
0085  
0086 errid = ['atmlab:' mfilename];
0087 %- Atmlab settings
0088 %
0089 arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
0090 arts_includes     = atmlab( 'ARTS_INCLUDES' );
0091 if isnan( arts_xmldata_path ) 
0092   error( errid,'You need to set ARTS_XMLDATA_PATH to run this exmaple.' );
0093 end
0094 if isnan( arts_includes )
0095   error( erird,'You need to ARTS_INCLUDES to run this example.' );
0096 end                                                      
0097 %
0098 fascod = fullfile( arts_xmldata_path, 'planets', 'Earth', 'Fascod' );
0099 
0100 
0101 %- Init Q
0102 %
0103 Q = qarts;
0104 %
0105 Q.INCLUDES            = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0106                           fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0107                           fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0108                           fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0109 Q.ATMOSPHERE_DIM      = 1;
0110 Q.STOKES_DIM          = 4;
0111 Q.J_DO                = true;
0112 Q.CLOUDBOX_DO         = false;
0113 
0114 
0115 %= Define agendas
0116 %
0117 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0118 Q.PPATH_STEP_AGENDA          = { 'ppath_step_agenda__GeometricPath'    };
0119 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0120 Q.IY_SPACE_AGENDA            = { 'iy_space_agenda__CosmicBackground'   };
0121 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0122 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0123 
0124 
0125 %- Radiative transfer
0126 %
0127 Q.IY_UNIT             = 'RJBT'; 
0128 Q.YCALC_WSMS          = { 'yCalc' };
0129 %
0130 Q.PPATH_LMAX          = 250;
0131 
0132 
0133 %- Surface
0134 %
0135 Q.Z_SURFACE           = 50;
0136 
0137 
0138 %- Absorption
0139 %
0140 Q.ABS_LINES           = fullfile( atmlab_example_data, 'o3line111ghz' );
0141 Q.ABS_LINES_FORMAT    = 'Arts';
0142 %
0143 Q.ABSORPTION          = 'OnTheFly';
0144 Q.ABS_NLS             = [];
0145 
0146 %- Pressure grid (a finer grid here, to better represent the troposphere)
0147 %
0148 z_toa                 = 95e3;
0149 %
0150 Q.P_GRID              = z2p_simple( Q.Z_SURFACE-500 : 500 : z_toa )';
0151 
0152 
0153 
0154 %- Frequency, polarisation, spectrometer and pencil beam antenna
0155 %
0156 % The hypothetical spectrometer has rectangular response functions
0157 %
0158 Q.F_GRID              = qarts_get( fullfile( atmlab_example_data , ...
0159                                                       'f_grid_111ghz.xml' ) );
0160 %
0161 H                     = qartsSensor;
0162 %
0163 H.SENSOR_NORM         = true;
0164 %
0165 H.BACKEND_DO          = true;
0166 df                    = 0.5e6;
0167 H.F_BACKEND           = ( min(Q.F_GRID)+df : df : max(Q.F_GRID)-df )';
0168 %
0169 B.name                = 'Spectrometer channel response function';
0170 B.gridnames           = { 'Frequency' };
0171 B.grids               = { [-df/2 df/2] };
0172 B.dataname            = 'Response';
0173 B.data                = [1 1];
0174 %
0175 H.BACKEND_CHANNEL_RESPONSE{1} = B;
0176 clear B
0177 %
0178 % We assume a H polarisation response. A number coding is applied, see:
0179 % arts -d sensor_pol
0180 H.SENSOR_POL          = { 6 };
0181 %
0182 Q.SENSOR_DO           = true;
0183 Q.SENSOR_RESPONSE     = H;
0184 %
0185 Q.ANTENNA_DIM         = 1;
0186 Q.MBLOCK_ZA_GRID      = 0;
0187 
0188 
0189 
0190 %- Correlation of thermal noise
0191 %
0192 f              = H.F_BACKEND;
0193 cl             = 1.4 * ( f(2) - f(1) );
0194 cfun           = 'gau';
0195 cco            = 0.05;
0196 %
0197 Q.TNOISE_C     = covmat1d_from_cfun( f, [], cfun, cl, cco );
0198 %
0199 clear H f
0200 
0201 
0202 %- Define L2 structure (beside retrieval quantities below)
0203 %
0204 Q.L2_EXTRA     = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', ...
0205                    'mresp', 'A', 'e', 'eo', 'es', 'date' };
0206 
0207 
0208 %- Temperature
0209 %
0210 Q.T.RETRIEVE   = false;
0211 Q.T.ATMDATA    = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
0212                         'msis90', 'msis90.t.xml' ), 'Temperature', 't_field' );
0213 
0214 %- Determine altitudes through HSE
0215 %
0216 Q.HSE.ON       = true;
0217 Q.HSE.P        = Q.P_GRID(1);
0218 Q.HSE.ACCURACY = 0.1;
0219 
0220 
0221 %- Species
0222 
0223 % Ozone, only species is retrieved here
0224 Q.ABS_SPECIES(1).TAG      = { 'O3' };
0225 Q.ABS_SPECIES(1).RETRIEVE = true;
0226 Q.ABS_SPECIES(1).L2       = true;
0227 % Note that grid below is a subset of Q.P_GRID
0228 Q.ABS_SPECIES(1).GRIDS    = { Q.P_GRID(1:4:end), [], [] };
0229 Q.ABS_SPECIES(1).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0230       'midlatitude-winter', 'midlatitude-winter.O3.xml' ), 'O3', 'vmr_field' );
0231 Q.ABS_SPECIES(1).UNIT     = 'rel';
0232 Q.ABS_SPECIES(1).SX       = ...
0233                    covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0234                                              'lin', 0.2, 0.00, @log10 ) + ...
0235                    covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0236                                              'lin', 0.5, 0.00, @log10 );
0237 Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0238   
0239   
0240 %- Water
0241 %
0242 Q.ABS_SPECIES(2).TAG      = { 'H2O-PWR98' };
0243 Q.ABS_SPECIES(2).RETRIEVE = true;
0244 Q.ABS_SPECIES(2).L2       = false;
0245 % Note that grid below is the first part of Q.P_GRID
0246 Q.ABS_SPECIES(2).GRIDS    = { Q.P_GRID(1:21), [], [] };
0247 Q.ABS_SPECIES(2).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0248     'midlatitude-winter', 'midlatitude-winter.H2O.xml' ), 'H2O', 'vmr_field' );
0249 Q.ABS_SPECIES(2).UNIT     = 'rel';
0250 Q.ABS_SPECIES(2).SX       = ...
0251                    covmat1d_from_cfun( Q.ABS_SPECIES(2).GRIDS{1}, 0.5, ...
0252                                              'lin', 0.5, 0.00, @log10 );
0253 Q.ABS_SPECIES(2).MINMAX  = [1e-6 10]';
0254 
0255 
0256 %- Liquid clouds (non-precip)
0257 %
0258 % There is no possibility to distinguish vapour and liquid water in this
0259 % example, and this part is included only for demonstration purpose.
0260 % If you activate liquid cloud retrieval, don't expect stable retrievals in
0261 % this demo case.
0262 % Note that the Jacobian calculation not works if any "vmr" is zero, and we
0263 % must set LWC to a small value if retrived.
0264 % Unit for LWC in ARTS is kg/m3, but is flagged as 'vmr' in the retrieval part.
0265 %
0266 Q.ABS_SPECIES(3).TAG      = { 'liquidcloud-MPM93' };
0267 Q.ABS_SPECIES(3).RETRIEVE = false;
0268 Q.ABS_SPECIES(3).L2       = true;
0269 % Note that grid below is the first part of Q.P_GRID
0270 Q.ABS_SPECIES(3).GRIDS    = { Q.P_GRID(1:11), [], [] };
0271 % Set LWC to zero at all altitudes
0272 Q.ABS_SPECIES(3).ATMDATA      = atmdata_empty(0);
0273 if Q.ABS_SPECIES(3).RETRIEVE
0274   Q.ABS_SPECIES(3).ATMDATA.DATA = 1e-10;
0275 else
0276   Q.ABS_SPECIES(3).ATMDATA.DATA = 0;
0277 end
0278 %
0279 Q.ABS_SPECIES(3).UNIT     = 'vmr';  % Yes 'vmr', even if unit is kg/m3
0280 Q.ABS_SPECIES(3).SX       = ...
0281                    covmat1d_from_cfun( Q.ABS_SPECIES(3).GRIDS{1}, 0.02e-3, ...
0282                                              'drc', 0, 0.00, @log10 );
0283 Q.ABS_SPECIES(3).MINMAX  = [1e-10 0.1e-3]';
0284 
0285 
0286 %- Oxygen
0287 %
0288 Q.ABS_SPECIES(4).TAG      = { 'O2-PWR93' };
0289 Q.ABS_SPECIES(4).RETRIEVE = false;
0290 Q.ABS_SPECIES(4).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0291       'midlatitude-winter', 'midlatitude-winter.O2.xml' ), 'O2', 'vmr_field' );
0292 
0293 %- Nitrogen
0294 %
0295 Q.ABS_SPECIES(5).TAG      = { 'N2-SelfContStandardType' };
0296 Q.ABS_SPECIES(5).RETRIEVE = false;
0297 Q.ABS_SPECIES(5).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0298       'midlatitude-winter', 'midlatitude-winter.N2.xml' ), 'N2', 'vmr_field' ); 
0299 
0300 
0301 
0302 %- Polyfit
0303 %
0304 % Not active, but if you want to test ...
0305 %
0306 Q.POLYFIT.RETRIEVE        = false;
0307 Q.POLYFIT.ORDER           = 1;
0308 Q.POLYFIT.L2              = true;
0309 Q.POLYFIT.SX0             = 1^2; 
0310 Q.POLYFIT.SX1             = 0.5^2; 
0311 
0312 return
0313 %-----------------------------------------------------------------------------
0314 %-----------------------------------------------------------------------------
0315 %-----------------------------------------------------------------------------
0316 
0317 
0318 
0319 
0320 function Y = y_demo(Q)
0321 
0322 % The data should be loaded from one or several files, but are here genereted
0323 % by a forward model call to show how qpack2 can also be used to generate
0324 % simulaled measurements (matching a priori assumptions).
0325 
0326 % The simulated data model airborn measurements at two different zenith
0327 % angles, from two nearby positions.
0328   
0329 % Init Y
0330 %
0331 Y = qp2_y;
0332 
0333 % Set a date
0334 %
0335 Y.YEAR  = 2008;
0336 Y.MONTH = 2;
0337 Y.DAY   = 25;
0338 
0339 % Lat / lon
0340 %
0341 Y.LATITUDE  = 45;
0342 Y.LONGITUDE = exp(1);
0343 
0344 % An airborn measurement assumed here
0345 %
0346 Y.Z_PLATFORM = 50;
0347 Y.ZA         = 0;
0348 
0349 % Reference point for hydrostatic equilibrium
0350 %
0351 Y.HSE_P = 1000e2;
0352 Y.HSE_Z = 0;
0353 
0354 % Set backend frequencies
0355 %
0356 Y.F = Q.SENSOR_RESPONSE.F_BACKEND;
0357 
0358 % Thermal noise standard deviation
0359 %
0360 Y.TNOISE = 0.1;
0361 % To test varying noise
0362 %Y.TNOISE = linspace( 0.03, 0.07, length(Y.F) )';
0363 
0364 % Simulate a measurement
0365 %
0366 Y0   = Y;
0367 Y0.Y = [];                          % A flag to tell qpack2 to calculate the
0368 %                                     spectrum ({} signifies undefined!).
0369 Y = qpack2( Q, oem, Y0 );           % Dummy oem structure OK here
0370 
0371 
0372 
0373 % Add a second measurement with H2O disturbed, and higher za
0374 %
0375 fascod = fullfile( atmlab( 'ARTS_XMLDATA_PATH' ), ...
0376                                                 'planets', 'Earth', 'Fascod' );
0377 Q.ABS_SPECIES(2).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0378     'midlatitude-summer', 'midlatitude-summer.H2O.xml' ), 'H2O', 'vmr_field' );
0379 %
0380 Y2 = qpack2( Q, oem, Y0 );
0381 
0382 
0383 % Append and add thermal noise
0384 %
0385 Y(2) = Y2;
0386 %
0387 % The correlation specified in Q is included
0388 %
0389 for i = 1 : length(Y)
0390   Y(i).Y = Y(i).Y + Y(i).TNOISE .* make_noise(1,Q.TNOISE_C);
0391 end
0392 
0393

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