Home > atmlab > demos > qpack2_wind3d_demo.m

qpack2_wind3d_demo

PURPOSE ^

QPACK2_WIND3D_DEMO Demonstration of Qpack2 with scanning in azimuth

SYNOPSIS ^

function L2 = qpack2_wind3d_demo

DESCRIPTION ^

 QPACK2_WIND3D_DEMO   Demonstration of Qpack2 with scanning in azimuth

    This file exemplifies how to set Qpack2 for an instrument scanning in
    azimuth (and azimuth angle matters). The example treats a measurement
    aimimg at retrieving both horisontal wind components.

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

qpack2_wind3d_demo.m

SOURCE CODE ^

0001 % QPACK2_WIND3D_DEMO   Demonstration of Qpack2 with scanning in azimuth
0002 %
0003 %    This file exemplifies how to set Qpack2 for an instrument scanning in
0004 %    azimuth (and azimuth angle matters). The example treats a measurement
0005 %    aimimg at retrieving both horisontal wind components.
0006 %
0007 % FORMAT   L2 = qpack2_wind3d_demo
0008 %
0009 % OUT      L2   L2 data output from *qpack2*.
0010 
0011 % 2010-05-12   Created by Patrick Eriksson.
0012 
0013 function L2 = qpack2_wind3d_demo
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 u_wind = 60;
0025 v_wind = -40;
0026 %
0027 Y = y_demo( Q, u_wind, v_wind );
0028 
0029 
0030 %- Check that all frequencies are OK
0031 %
0032 if ~qp2_check_f( Q, Y, 1e3 );
0033   error( errid, ...
0034             'Some mismatch between Q.F_BACKEND and frequencies of spectra.' );
0035 end
0036 
0037 
0038 %- OEM variables
0039 %
0040 O = qp2_l2( Q );         % This ensures that OEM returns the variables needed
0041 %                          to fill the L2 structure, as defined in Q
0042 O.linear = true;
0043 %
0044 if ~O.linear 
0045   O.itermethod = 'GN';
0046   O.stop_dx    = 0.01;
0047   O.maxiter    = 5;
0048 end
0049 
0050 
0051 %- Make inversion
0052 %
0053 L2 = qpack2( Q, O, Y );
0054 
0055 
0056 %- Plot, if no output argument
0057 %
0058 if ~nargout
0059 
0060   % Wind profiles
0061   figure(1),clf  
0062   plot( [0 0], [10 90], 'k-', ...
0063         L2(1).wind_u_x, p2z_simple(L2(1).wind_u_p)/1e3, 'b', ...
0064         u_wind*[1 1], [10 90], 'm', ...
0065         L2(1).wind_v_x, p2z_simple(L2(1).wind_v_p)/1e3, 'r', ...
0066         v_wind*[1 1], [10 90], 'c' );
0067   xlabel( 'Wind components [m/s]' );
0068   ylabel( 'Approximate altitude [km]' );
0069   axis( [ -80 80 10 90 ] )
0070   legend( 'A priori, U and V', 'Retrieved U', 'True U', ...
0071                                'Retrieved V', 'True V' );
0072 
0073   % Ozone profile
0074   figure(2),clf
0075   plot( L2(1).species1_x*1e2, p2z_simple(L2(1).species1_p)/1e3, 'b', ...
0076         L2(1).species1_xa*1e2, p2z_simple(L2(1).species1_p)/1e3, 'k-' );
0077   xlabel( 'Ozone compared to a priori [%]' );
0078   ylabel( 'Approximate altitude [km]' );
0079   axis( [ 0 130 10 90 ] )
0080   legend( 'Retrieval', 'True and a priori' );
0081 
0082   % Spectra
0083   figure(3),clf
0084   plot( L2(1).y, 'k.' );
0085   hold on
0086   h(1) = plot( L2(1).yf, 'b-' );
0087   h(2) = plot( L2(1).bl, 'b-.');
0088   xlabel( 'Data index' );
0089   ylabel( 'Tb [K]' );
0090 %  axis( [ min(L2(1).f/1e9) max(L2(1).f/1e9) 0 18 ] )
0091   legend( h, 'Fitted', 'Baseline' );
0092 end
0093 
0094 
0095 return
0096 %-----------------------------------------------------------------------------
0097 %-----------------------------------------------------------------------------
0098 %-----------------------------------------------------------------------------
0099 
0100 
0101 
0102 
0103 
0104 function Q = q_demo
0105  
0106 errid = ['atmlab:' mfilename];
0107 %- Atmlab settings
0108 %
0109 arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
0110 arts_includes     = atmlab( 'ARTS_INCLUDES' );
0111 if isnan( arts_xmldata_path ) 
0112   error( errid,'You need to set ARTS_XMLDATA_PATH to run this exmaple.' );
0113 end
0114 if isnan( arts_includes )
0115   error( erird,'You need to ARTS_INCLUDES to run this example.' );
0116 end                                                      
0117 %
0118 fascod = fullfile( arts_xmldata_path, 'planets', 'Earth', 'Fascod' );
0119 
0120 
0121 %- Init Q
0122 %
0123 Q = qarts;
0124 %
0125 Q.INCLUDES            = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0126                           fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0127                           fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0128                           fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0129 Q.ATMOSPHERE_DIM      = 3;
0130 Q.STOKES_DIM          = 1;
0131 Q.J_DO                = true;
0132 Q.CLOUDBOX_DO         = false;
0133 
0134 
0135 %= Define agendas
0136 %
0137 % Here we do it by using the predefined agenda templates
0138 %   (found in arts/controlfiles/general/agendas.arts)
0139 % This works only if the pre-defined agenda is names following the pattern:
0140 %   name_of_agenda__(Something)
0141 %
0142 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0143 Q.PPATH_STEP_AGENDA          = { 'ppath_step_agenda__GeometricPath'    };
0144 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0145 Q.IY_SPACE_AGENDA            = { 'iy_space_agenda__CosmicBackground'   };
0146 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0147 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0148 
0149 
0150 %- Radiative transfer
0151 %
0152 Q.IY_UNIT             = 'RJBT'; 
0153 Q.YCALC_WSMS          = { 'yCalc' };
0154 %
0155 Q.PPATH_LMAX          = 250;
0156 
0157 
0158 %- Absorption
0159 %
0160 Q.ABS_LINES           = fullfile( atmlab_example_data, 'o3line111ghz' );
0161 Q.ABS_LINES_FORMAT    = 'Arts';
0162 %
0163 Q.ABSORPTION          = 'OnTheFly';
0164 Q.ABS_NLS             = [];
0165 
0166 %- Atmospheric grids
0167 %
0168 z_toa                 = 95e3;
0169 z_surf                = 1e3;
0170 %
0171 Q.P_GRID              = z2p_simple( z_surf-1e3 : 2e3 : z_toa )';
0172 Q.LAT_GRID            = [ 30 : 10 : 60 ]';
0173 Q.LON_GRID            = [ -20 : 20 : 20 ]';
0174 
0175 
0176 %- Surface
0177 %
0178 Q.Z_SURFACE           = repmat( z_surf, length(Q.LAT_GRID), ...
0179                                         length(Q.LON_GRID) );
0180 
0181 
0182 %- Frequency, spectrometer and pencil beam antenna
0183 %
0184 % The hypothetical spectrometer has rectangular response functions
0185 %
0186 Q.F_GRID              = qarts_get( fullfile( atmlab_example_data , ...
0187                                                       'f_grid_111ghz.xml' ) );
0188 %
0189 H                     = qartsSensor;
0190 %
0191 H.SENSOR_NORM         = true;
0192 %
0193 H.BACKEND_DO          = true;
0194 df                    = 0.05e6;
0195 H.F_BACKEND           = ( min(Q.F_GRID)+df : df : max(Q.F_GRID)-df )';
0196 %
0197 B.name                = 'Spectrometer channel response function';
0198 B.gridnames           = { 'Frequency' };
0199 B.grids               = { [-df/2 df/2] };
0200 B.dataname            = 'Response';
0201 B.data                = [1 1];
0202 %
0203 H.BACKEND_CHANNEL_RESPONSE{1} = B;
0204 clear B
0205 %
0206 Q.SENSOR_DO           = true;
0207 Q.SENSOR_RESPONSE     = H;
0208 %
0209 Q.ANTENNA_DIM         = 1;
0210 Q.MBLOCK_ZA_GRID      = 0;
0211 
0212 
0213 
0214 %- Correlation of thermal noise
0215 %
0216 Q.TNOISE_C = speye( length(H.F_BACKEND) );
0217 %
0218 clear H 
0219 
0220 
0221 %- Define L2 structure (beside retrieval quantities below)
0222 %
0223 Q.L2_EXTRA     = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', ...
0224                    'mresp', 'A', 'e', 'eo', 'es', 'date' };
0225 
0226 
0227 %- Temperature
0228 %
0229 Q.T.RETRIEVE   = false;
0230 Q.T.ATMDATA    = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
0231                         'msis90', 'msis90.t.xml' ), 'Temperature', 't_field' );
0232 
0233 %- Determine altitudes through HSE
0234 %
0235 Q.HSE.ON       = true;
0236 Q.HSE.P        = Q.P_GRID(1);
0237 Q.HSE.ACCURACY = 0.1;
0238 
0239 
0240 %- Species
0241 
0242 % Ozone, only species is retrieved here
0243 Q.ABS_SPECIES(1).TAG      = { 'O3' };
0244 Q.ABS_SPECIES(1).RETRIEVE = true;
0245 Q.ABS_SPECIES(1).L2       = true;
0246 Q.ABS_SPECIES(1).GRIDS    = { Q.P_GRID, mean(Q.LAT_GRID), mean(Q.LON_GRID)};
0247 Q.ABS_SPECIES(1).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0248       'midlatitude-winter', 'midlatitude-winter.O3.xml' ), 'O3', 'vmr_field' );
0249 %
0250 Q.ABS_SPECIES(1).UNIT     = 'rel';
0251 Q.ABS_SPECIES(1).SX       = ...
0252                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0253                                                'lin', 0.2, 0.00, @log10 ) + ...
0254                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0255                                                'lin', 0.5, 0.00, @log10 );
0256 Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0257   
0258   
0259 %- Water
0260 %
0261 % This generates no absorption, as linefile has no H2O lines
0262 %
0263 Q.ABS_SPECIES(2).TAG      = { 'H2O' };
0264 Q.ABS_SPECIES(2).RETRIEVE = false;
0265 Q.ABS_SPECIES(2).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0266     'midlatitude-winter', 'midlatitude-winter.H2O.xml' ), 'H2O', 'vmr_field' );
0267 
0268 
0269 
0270 %- Horsontal wind components vind
0271 %
0272 Q.WIND_V_FIELD    = [];      % This is a short-cut for zero wind.
0273 %
0274 Q.WIND_V.RETRIEVE = true;   % Set to true to activate retrieval
0275 Q.WIND_V.L2       = true;
0276 Q.WIND_V.GRIDS    = { Q.P_GRID(1:2:end), mean(Q.LAT_GRID), mean(Q.LON_GRID) };
0277 Q.WIND_V.SX       = covmat1d_from_cfun( Q.WIND_V(1).GRIDS{1}, 50, ...
0278                                                     'lin', 0.8, 0.00, @log10 );
0279 if Q.WIND_V.RETRIEVE
0280   Q.WSMS_AT_START{end+1} = 'IndexSet(abs_f_interp_order,3)';
0281 end
0282 %
0283 % Exactly same for u-component
0284 Q.WIND_U_FIELD    = Q.WIND_V_FIELD;
0285 Q.WIND_U          = Q.WIND_V;
0286 
0287 
0288 %- Frequency shift
0289 %
0290 Q.FSHIFTFIT.RETRIEVE      = true;
0291 Q.FSHIFTFIT.DF            = 25e3; 
0292 Q.FSHIFTFIT.SX            = 50e3^2;  
0293 Q.FSHIFTFIT.L2            = true;
0294 
0295 
0296 %- Polyfit
0297 %
0298 Q.POLYFIT.RETRIEVE        = true;
0299 Q.POLYFIT.ORDER           = 0;
0300 Q.POLYFIT.L2              = true;
0301 Q.POLYFIT.SX0             = 1.0^2 * speye(5); 
0302 
0303 
0304 return
0305 %-----------------------------------------------------------------------------
0306 %-----------------------------------------------------------------------------
0307 %-----------------------------------------------------------------------------
0308 
0309 
0310 
0311 
0312 
0313 function Y = y_demo(Q,u_wind,v_wind)
0314 
0315 % The data should be loaded from one or several files, but are here genereted
0316 % by a forward model call to show how qpack2 can also be used to generate
0317 % simulaled measurements (matching a priori assumptions).
0318 
0319 % The simulated data model airborn measurements at two different zenith
0320 % angles, from two nearby positions.
0321   
0322 % Init Y
0323 %
0324 Y = qp2_y;
0325 
0326 % Set a date
0327 %
0328 Y.YEAR  = 2008;
0329 Y.MONTH = 2;
0330 Y.DAY   = 25;
0331 
0332 % Lat / lon
0333 %
0334 Y.LATITUDE  = 45;
0335 Y.LONGITUDE = exp(1);
0336 
0337 % An airborn measurement assumed here
0338 %
0339 Y.AA         = [ 0 -135 -45 45  135 ]';
0340 naa          = length( Y.AA );
0341 Y.Z_PLATFORM = repmat( 10.5e3, naa, 1 );
0342 Y.ZA         = [ 0; repmat( 70, naa-1, 1 ) ];
0343 
0344 % Reference point for hydrostatic equilibrium
0345 %
0346 Y.HSE_P = 100e2;
0347 Y.HSE_Z = 16e3;
0348 
0349 % Set backend frequencies
0350 %
0351 Y.F = qarts_get( Q.SENSOR_RESPONSE.F_BACKEND );
0352 nf  = length( Q.SENSOR_RESPONSE.F_BACKEND );
0353 
0354 % Thermal noise standard deviation
0355 %
0356 Y.TNOISE = repmat( 0.05, 1, naa );
0357 % To test varying noise
0358 %Y.TNOISE = repmat( linspace( 0.03, 0.07, length(Y.F) )', 1, naa );
0359 
0360 
0361 % Include winds, constant with altitude
0362 %
0363 Q.WIND_U_FIELD = repmat( u_wind, [ length(Q.P_GRID), length(Q.LAT_GRID), ...
0364                                                      length(Q.LON_GRID) ] );
0365 Q.WIND_V_FIELD = repmat( v_wind, [ length(Q.P_GRID), length(Q.LAT_GRID), ...
0366                                                      length(Q.LON_GRID) ] );
0367 
0368 % Simulate a measurement
0369 %
0370 Y.Y = [];                           % A flag to tell qpack2 to calculate the
0371 %                                     spectrum ({} signifies undefined!).
0372 
0373 % Calculate simulated spectra
0374 %
0375 Y = qpack2( Q, oem, Y );            % Dummy oem structure OK here
0376 
0377 
0378 % Add thermal noise and baseline off-set
0379 %
0380 % The correlation specified in Q is included
0381 %
0382 for i = 1 : length(Y)
0383   for j = 1 : length( Y.ZA )
0384     ind = (j-1)*nf + [1:nf];
0385     Y(i).Y(ind) = Y(i).Y(ind) + ...   % Noise free spectrum
0386                   2*randn(1)  + ...   % A constant off-set
0387                   Y(i).TNOISE(:,j) .* make_noise(1,Q.TNOISE_C); % Thermal noise
0388   end
0389 end
0390

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