0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 function L2 = qpack2_demo
0023
0024 errid = ['atmlab:' mfilename];
0025
0026
0027
0028 Q = q_demo;
0029
0030
0031
0032
0033 Y = y_demo( Q );
0034
0035
0036
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
0045
0046 O = qp2_l2( Q );
0047
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
0058
0059 L2 = qpack2( Q, O, Y );
0060
0061
0062
0063
0064 if ~nargout
0065
0066
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
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
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
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
0128
0129
0130
0131
0132
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
0143
0144 Q.IY_UNIT = 'RJBT';
0145 Q.YCALC_WSMS = { 'yCalc' };
0146
0147 Q.PPATH_LMAX = 250;
0148
0149
0150
0151
0152 Q.Z_SURFACE = 1e3;
0153
0154
0155
0156
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
0165
0166 z_toa = 95e3;
0167
0168 Q.P_GRID = z2p_simple( Q.Z_SURFACE-1e3 : 2e3 : z_toa )';
0169
0170
0171
0172
0173
0174
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
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
0217
0218 Q.L2_EXTRA = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', ...
0219 'mresp', 'A', 'e', 'eo', 'es', 'date' };
0220
0221
0222
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
0229
0230 Q.HSE.ON = true;
0231 Q.HSE.P = Q.P_GRID(1);
0232 Q.HSE.ACCURACY = 0.1;
0233
0234
0235
0236
0237
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
0246
0247
0248
0249 switch 1
0250 case 1
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
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
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
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
0291
0292
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
0302
0303
0304
0305
0306
0307
0308 Q.WIND_V_FIELD = [];
0309
0310 Q.WIND_V.RETRIEVE = false;
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
0322
0323
0324
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
0335
0336 Q.FSHIFTFIT.RETRIEVE = true;
0337 Q.FSHIFTFIT.DF = 25e3;
0338 Q.FSHIFTFIT.SX = 50e3^2;
0339 Q.FSHIFTFIT.L2 = true;
0340
0341
0342
0343 Q.FSTRETCHFIT.RETRIEVE = false;
0344 Q.FSTRETCHFIT.DF = 25e3;
0345 Q.FSTRETCHFIT.SX = 50e3^2;
0346 Q.FSTRETCHFIT.L2 = true;
0347
0348
0349
0350
0351
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
0363
0364
0365
0366 Q.SINEFIT.RETRIEVE = false;
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
0384
0385
0386
0387
0388
0389
0390
0391
0392 Y = qp2_y;
0393
0394
0395
0396 Y.YEAR = 2008;
0397 Y.MONTH = 2;
0398 Y.DAY = 25;
0399
0400
0401
0402 Y.LATITUDE = 45;
0403 Y.LONGITUDE = exp(1);
0404
0405
0406
0407 Y.Z_PLATFORM = 10.5e3;
0408 Y.ZA = 50;
0409
0410
0411
0412 Y.HSE_P = 100e2;
0413 Y.HSE_Z = 16e3;
0414
0415
0416
0417 Y.F = Q.SENSOR_RESPONSE.F_BACKEND;
0418
0419
0420
0421 Y.TNOISE = 0.1;
0422
0423
0424
0425
0426
0427 Y.Y = [];
0428
0429
0430
0431
0432 Y(2) = Y(1);
0433
0434 Y(2).LONGITUDE = pi;
0435 Y(2).ZA = 45;
0436
0437
0438
0439
0440 Y = qpack2( Q, oem, Y );
0441
0442
0443
0444
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
0451
0452 Y(2).Y = Y(2).Y + 1;
0453