Home > atmlab > demos > arts_oem_demo.m

arts_oem_demo

PURPOSE ^

ARTS_OEM_DEMO Demonstration of inversions using ARTS and OEM

SYNOPSIS ^

function [X,R] = arts_oem_demo

DESCRIPTION ^

 ARTS_OEM_DEMO   Demonstration of inversions using ARTS and OEM

    The example treats ground-based measurements of ozone at 110.8 GHz.

    A figure with the inversion result is generated if the function is
    called with no output arguments. Otherwise, the output from *oem* is
    returned (work folder is deleted).

 FORMAT   [X,R] = arts_oem_demo

 OUT   X   As same output from *oem*.
       R   As same output from *oem*.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

arts_oem_demo.m

SOURCE CODE ^

0001 % ARTS_OEM_DEMO   Demonstration of inversions using ARTS and OEM
0002 %
0003 %    The example treats ground-based measurements of ozone at 110.8 GHz.
0004 %
0005 %    A figure with the inversion result is generated if the function is
0006 %    called with no output arguments. Otherwise, the output from *oem* is
0007 %    returned (work folder is deleted).
0008 %
0009 % FORMAT   [X,R] = arts_oem_demo
0010 %
0011 % OUT   X   As same output from *oem*.
0012 %       R   As same output from *oem*.
0013 
0014 % 2006-09-07   Created by Patrick Eriksson.
0015 
0016 function [X,R] = arts_oem_demo
0017 
0018 
0019 %- Check if needed atmlab settings are found
0020 %
0021 arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
0022 if isnan( arts_xmldata_path )
0023   error('You need to ARTS_XMLDATA_PATH to run this example.');
0024 end
0025 %
0026 arts_includes = atmlab( 'ARTS_INCLUDES' );
0027 if isnan( arts_includes )
0028   error('You need to ARTS_INCLUDES to run this example.');
0029 end
0030 
0031 
0032 
0033 %----------------------------------------------------------------------------
0034 %- Init control structures
0035 %----------------------------------------------------------------------------
0036 %
0037 Q      = qarts;
0038 Q.J_DO = true;
0039 O      = oem;
0040 
0041 
0042 %----------------------------------------------------------------------------
0043 %- OEM settings
0044 %----------------------------------------------------------------------------
0045 %
0046 O.itermethod = 'GN';
0047 O.stop_dx    = 0.01;
0048 O.maxiter    = 20;
0049 %
0050 O.cost       = true;
0051 O.e          = true; 
0052 O.eo         = true;
0053 O.es         = true;
0054 O.yf         = true;
0055 
0056 
0057 
0058 %----------------------------------------------------------------------------
0059 %- Forward model parameters
0060 %----------------------------------------------------------------------------
0061 
0062 %- General
0063 %
0064 Q.INCLUDES            = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0065                           fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0066                           fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0067                           fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0068 Q.ATMOSPHERE_DIM      = 1;
0069 Q.STOKES_DIM          = 1;
0070 
0071 %- Atmosphere
0072 %
0073 Q.Z_SURFACE           = 0;
0074 Q.P_GRID              = z2p_simple( -250:250:90e3 )';
0075 %
0076 Q.RAW_ATMOSPHERE      = fullfile( arts_xmldata_path, 'planets', 'Earth', ...
0077                         'Fascod', 'midlatitude-winter', 'midlatitude-winter' );
0078 Q.CLOUDBOX_DO         = false;
0079 
0080 
0081 %= Define agendas
0082 %
0083 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0084 Q.PPATH_STEP_AGENDA          = { 'ppath_step_agenda__GeometricPath'    };
0085 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0086 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0087 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0088 
0089 
0090 %- Exclude cosmic background radiation. A relative change of the spectrum
0091 %- corresponds then to a relative change of the ozone profile.
0092 %
0093 Q.IY_SPACE_AGENDA     = { 'Ignore(rtp_los)',  'Ignore(rtp_pos)', ...
0094                           'nelemGet(nrows,f_grid)', ...
0095                           'Copy(ncols,stokes_dim)', ...
0096                           'MatrixSetConstant(iy,nrows,ncols,0)' };
0097 
0098 
0099 %- RTE
0100 %
0101 Q.YCALC_WSMS          = { 'yCalc' };
0102 %
0103 Q.PPATH_LMAX          = 250;
0104 Q.IY_UNIT              = 'RJBT';
0105 %
0106 Q.SENSOR_LOS          = 70;
0107 Q.SENSOR_POS          = Q.Z_SURFACE;
0108 
0109 
0110 %- Frequency, spectrometer and pencil beam antenna
0111 %
0112 % The hypothetical spectrometer has rectangular response functions
0113 %
0114 Q.F_GRID              = qarts_get( fullfile( atmlab_example_data , ...
0115                                               'f_grid_111ghz.xml' ) );
0116 %
0117 H                     = qartsSensor;
0118 %
0119 H.SENSOR_NORM         = true;
0120 %
0121 df                    = 0.5e6;
0122 H.F_BACKEND           = [ min(Q.F_GRID)+df : df : max(Q.F_GRID)-df ]';
0123 %
0124 B.name                = 'Spectrometer channel response function';
0125 B.gridnames           = { 'Frequency' };
0126 B.grids               = { [-df/2 df/2] };
0127 B.dataname            = 'Response';
0128 B.data                = [1 1];
0129 %
0130 H.BACKEND_CHANNEL_RESPONSE{1} = B;
0131 H.BACKEND_DO = true;
0132 clear B
0133 %
0134 Q.SENSOR_DO           = true;
0135 Q.SENSOR_RESPONSE     = H;
0136 %
0137 Q.ANTENNA_DIM         = 1;
0138 Q.MBLOCK_ZA_GRID      = 0;
0139 
0140 
0141 
0142 %----------------------------------------------------------------------------
0143 %- Species
0144 %----------------------------------------------------------------------------
0145 
0146 % Water vapour has to be inluded (needed by absorption part)
0147 %
0148 Q.ABS_SPECIES(1).TAG      = { 'H2O' };
0149 Q.ABS_SPECIES(1).RETRIEVE = false;
0150 
0151 % Ozone, retrieved
0152 %
0153 Q.ABS_SPECIES(2).TAG      = { 'O3' };
0154 Q.ABS_SPECIES(2).RETRIEVE = true;
0155 Q.ABS_SPECIES(2).GRIDS    = { z2p_simple( 0e3:1e3:90e3 )', [], [] };
0156 Q.ABS_SPECIES(2).UNIT     = 'vmr';
0157 Q.ABS_SPECIES(2).MINMAX   = 1e-12;
0158 Q.ABS_SPECIES(2).SX       = covmat1d_from_cfun( Q.ABS_SPECIES(2).GRIDS{1}, ...
0159                                     2e-6, 'lin', 0.2, 0.00, @log10 ) + ...
0160                             covmat1d_from_cfun( Q.ABS_SPECIES(2).GRIDS{1}, ...
0161                                   0.5e-6, 'lin', 0.5, 0.00, @log10 );
0162 
0163 
0164 %----------------------------------------------------------------------------
0165 %- Set-up/create absorption look-up table
0166 %----------------------------------------------------------------------------
0167 
0168 Q.ABS_LINES_FORMAT         = 'Arts';
0169 Q.ABS_LINES                = fullfile( atmlab_example_data , 'o3line111ghz' );
0170 Q.ABS_NLS                  = [];
0171 %
0172 Q                          = qarts_abstable( Q ); 
0173 % The absorption table will with the settings above be calculated
0174 % automatically. If you want to pre-calculate the table, for later re-usage,
0175 % do:
0176 %Q.ABS_LOOKUP               = arts_abstable( Q );
0177 
0178 
0179 
0180 %----------------------------------------------------------------------------
0181 %- Init retrieval variables (such as xa) and create Sx
0182 %----------------------------------------------------------------------------
0183   
0184 workfolder = create_tmpfolder;
0185 cu = onCleanup( @()delete_tmpfolder( workfolder ) );
0186 %
0187 [Qoem,O,R,xa] = arts_oem_init( Q, O, workfolder );
0188 %
0189 Sx = arts_sx( Q, R );
0190 
0191 
0192 
0193 %----------------------------------------------------------------------------
0194 %- Create a measurement
0195 %----------------------------------------------------------------------------
0196 
0197 %- Noise statistics
0198 %
0199 si = 0.01;
0200 %
0201 SE.FORMAT    = 'param';
0202 SE.SI        = si;
0203 SE.CCO       = 0;  
0204 SE.CFUN1     = 'drc';
0205 SE.CL1       = 0; 
0206 
0207 %- Create Se
0208 %
0209 Se = covmat3d( 1, SE, Q.SENSOR_RESPONSE.F_BACKEND );
0210 
0211 %- Calculate a spectrum matching a priori
0212 %
0213 y = arts_y( Q );
0214 
0215 %- Crate a test measurement
0216 %
0217 y = 1.1 * y + si*randn(size(Se,1),1);
0218 
0219 
0220 
0221 
0222 %----------------------------------------------------------------------------
0223 %- Calculations
0224 %----------------------------------------------------------------------------
0225 
0226 %- Perform retrieval
0227 %
0228 X = oem( O, Qoem, R, @arts_oem, Sx, Se, [], [], xa, y );
0229 %
0230 clear Qoem 
0231 
0232 
0233 %- Plot
0234 %
0235 if ~nargout
0236   %
0237   z = p2z_simple( Q.ABS_SPECIES(2).GRIDS{1} );
0238   %
0239   if O.e & O.eo
0240     h= plot( 1e6*X.x, z/1e3, 'k-', ...
0241              1e6*(repmat(X.x,1,2)+[-X.e X.e]), z/1e3, 'b:', ...
0242              1e6*(repmat(X.x,1,2)+[-X.eo X.eo]), z/1e3, 'r--' );
0243     legend( h([1 4 2]), 'Retrieved profile', 'Observation error', ...
0244                                                                'Total error' );
0245   else
0246     plot( 1e6*X.x, z/1e3 )
0247   end
0248   xlabel( 'Ozone [ppm]' )
0249   ylabel( 'Approximative altitude [km]' );
0250 end

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