Home > atmlab > tests > test_arts_ppath.m

test_arts_ppath

PURPOSE ^

TEST_ARTS_PPATH Testing of ARTS ppath calculations

SYNOPSIS ^

function test_arts_ppath( atmosphere_dim, n, varargin )

DESCRIPTION ^

 TEST_ARTS_PPATH   Testing of ARTS ppath calculations

    Random test cases are generated. The following variables are
    randomised:
       length of p_grid
       length of lat_grid
       length of lon_grid
       z_surface
       z_field
       sensor_pos
       sensor_los

    A ppath calculation for ARTS is performed for each test case. A basic check
    of the ARTS calculations performed, by using atmlab functions. If any
    error is found, "keyboard" is called. The work folder is then left for
    debugging purposes.

 FORMAT   test_arts_ppath( atmosphere_dim, n, [unsphere, do_refr] )
        
 IN    atmosphere_dim   Atmospheric dimensionality
       n                Run n cases
 OPT   unsphere         Deviation from spherical geometry. Default is 0.
                        Possible choices:
                          0: 1D geometry
                          1: z_field disturbed
                          2: 1 + non-spherical ellipsoid
       do_refr          Run with refraction ray tracing WSM (but with n=1).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

test_arts_ppath.m

SOURCE CODE ^

0001 % TEST_ARTS_PPATH   Testing of ARTS ppath calculations
0002 %
0003 %    Random test cases are generated. The following variables are
0004 %    randomised:
0005 %       length of p_grid
0006 %       length of lat_grid
0007 %       length of lon_grid
0008 %       z_surface
0009 %       z_field
0010 %       sensor_pos
0011 %       sensor_los
0012 %
0013 %    A ppath calculation for ARTS is performed for each test case. A basic check
0014 %    of the ARTS calculations performed, by using atmlab functions. If any
0015 %    error is found, "keyboard" is called. The work folder is then left for
0016 %    debugging purposes.
0017 %
0018 % FORMAT   test_arts_ppath( atmosphere_dim, n, [unsphere, do_refr] )
0019 %
0020 % IN    atmosphere_dim   Atmospheric dimensionality
0021 %       n                Run n cases
0022 % OPT   unsphere         Deviation from spherical geometry. Default is 0.
0023 %                        Possible choices:
0024 %                          0: 1D geometry
0025 %                          1: z_field disturbed
0026 %                          2: 1 + non-spherical ellipsoid
0027 %       do_refr          Run with refraction ray tracing WSM (but with n=1).
0028 
0029 % 2012-02-22   Created by Patrick.
0030 
0031 function test_arts_ppath( atmosphere_dim, n, varargin  )
0032 %
0033 [unsphere,do_refr] = optargs( varargin, { 0, false } );
0034 
0035 do_check = 1;
0036   
0037 %- Fixed settings
0038 %
0039 Q = qarts;
0040 %
0041 Q.INCLUDES              = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0042                             fullfile( 'ARTS_INCLUDES', 'agendas.arts' ) };
0043 Q.ATMOSPHERE_DIM        = atmosphere_dim;
0044 Q.CLOUDBOX_DO           = false;
0045 Q.F_GRID                = 100e9;
0046 Q.ABS_SPECIES(1).TAG{1} = 'H2O';
0047 %
0048 Q.PPATH_AGENDA          = { 'ppath_agenda__FollowSensorLosPath'   };
0049 %
0050 if do_refr
0051   Q.PPATH_LRAYTRACE     = 1e3;
0052   Q.PPATH_STEP_AGENDA   = { 'ppath_step_agenda__RefractedPath'    };
0053   Q.REFR_INDEX_AIR_AGENDA = { ...
0054                             'Ignore(f_grid)', ...
0055                             'Ignore(rtp_pressure)', ...
0056                             'Ignore(rtp_temperature)', ...
0057                             'Ignore(rtp_vmr)', ...
0058                             'NumericSet(refr_index_air,1.0)', ...
0059                             'NumericSet(refr_index_air_group,1.0)' };
0060 else
0061   Q.PPATH_STEP_AGENDA   = { 'ppath_step_agenda__GeometricPath'    };
0062 end
0063 
0064 Q.PPATH_LMAX            = 10e3;
0065 %
0066 if atmosphere_dim == 1  ||  unsphere < 2
0067   refell                = ellipsoidmodels( 'SphericalEarth' );
0068   re_string             = 'refellipsoidEarth( refellipsoid, "Sphere" )';
0069 else
0070   refell                = ellipsoidmodels( 'WGS84' );
0071   re_string             = 'refellipsoidEarth( refellipsoid, "WGS84" )'; 
0072 end
0073 
0074 
0075 %- Load a test atmosphere
0076 %
0077 atmtest =  fullfile( atmlab('ARTS_XMLDATA_PATH'), 'planets', 'Earth', ...
0078                         'Fascod', 'midlatitude-winter', 'midlatitude-winter' );
0079 Az      = xmlLoad( [ atmtest, '.z.xml' ] ); 
0080 At      = xmlLoad( [ atmtest, '.t.xml' ] ); 
0081 Ah2o    = xmlLoad( [ atmtest, '.H2O.xml' ] ); 
0082 
0083 
0084 
0085 %- Create a workfolder
0086 %
0087 workfolder = create_tmpfolder;
0088 
0089 
0090 %- To avoid getting same random sequence between sessions:
0091 %
0092 RandStream.setGlobalStream(RandStream('mt19937ar','seed',sum(100*clock)));
0093 
0094 
0095 %- Run through n random cases
0096 %
0097 emaxatm  = 0;
0098 emaxsurf = 0;
0099 %
0100 for c = 1 : n
0101 
0102   fprintf(' > Doing %d/%d\r',c,n)
0103 
0104   %- Pressure grid
0105   %
0106   np                    = 30 + round( 170*rand );
0107   Q.P_GRID              = z2p_simple( linspace( 0, 80e3, np )' );
0108 
0109   %- Number of latitude and longitude grid points
0110   %
0111   nlat                = 1;
0112   nlon                = 1;
0113   %
0114   if atmosphere_dim == 2
0115     nlat                = 10 + round( 531*rand );
0116     Q.LAT_GRID          = linspace( -135, 135, nlat )';
0117     nlon                = 1;
0118   elseif atmosphere_dim == 3
0119     nlat                = 10 + round( 190*rand );
0120     nlon                = 10 + round( 190*rand );
0121     Q.LAT_GRID          = linspace( -80, 80, nlat )';
0122     Q.LON_GRID          = linspace( -90, 90, nlon )';
0123   end
0124 
0125   % Atmospheric fields
0126   Q.Z_FIELD = repmat( interpp( Az.grids{1}, Az.data, Q.P_GRID ), ...
0127                                                            [ 1, nlat, nlon ] ); 
0128   if unsphere
0129     Q.Z_FIELD = Q.Z_FIELD + 5*randn( size(Q.Z_FIELD) );
0130   end
0131   Q.T_FIELD = repmat( interpp( At.grids{1}, At.data, Q.P_GRID ), ...
0132                                                            [ 1, nlat, nlon ] ); 
0133   Q.VMR_FIELD = repmat( interpp( Ah2o.grids{1}, Ah2o.data, Q.P_GRID ), ...
0134                                                            [ 1, nlat, nlon ] ); 
0135   Q.VMR_FIELD = reshape( Q.VMR_FIELD, [1 np nlat nlon ] );
0136   %
0137   if atmosphere_dim == 3
0138     % Fix border conditions
0139     Q.Z_FIELD(:,:,end) = Q.Z_FIELD(:,:,1);
0140     Q.Z_FIELD(:,1,:)   = repmat(Q.Z_FIELD(:,1,1),[1 1 nlon]);
0141     Q.Z_FIELD(:,end,:) = repmat(Q.Z_FIELD(:,end,1),[1 1 nlon]);
0142   end
0143   
0144   %- Surface altitude
0145   %
0146   if ~unsphere
0147     Q.Z_SURFACE           = repmat( 500, nlat, nlon );
0148   else
0149     Q.Z_SURFACE           = 500*rand + 5e3*rand( nlat, nlon );
0150     if atmosphere_dim == 3
0151       % Avoid strong topography at the poles:
0152       Q.Z_SURFACE        = Q.Z_SURFACE .* repmat(cosd(Q.LAT_GRID),1,nlon);
0153       % Fix border conditions
0154       Q.Z_SURFACE(:,end) = Q.Z_SURFACE(:,1);
0155     end
0156     Q.Z_SURFACE = Q.Z_SURFACE + 500;  
0157   end
0158   
0159   %- Sensor pos/los
0160   %
0161   lat0                  = 0;
0162   lon0                  = 0;
0163   aa0                   = 0;
0164   %
0165   if rand < 0.33
0166     Q.SENSOR_POS        = 7e3 + 40e3*rand;
0167     Q.SENSOR_LOS        = 180*rand;    
0168   else
0169     Q.SENSOR_POS        = 700e3 + 50e3*randn;
0170     if rand < 0.5
0171       Q.SENSOR_LOS      = 140  + 40*rand;    
0172     else
0173       Q.SENSOR_LOS      = 111 + 5*rand;    
0174     end
0175   end
0176   %
0177   if atmosphere_dim == 2
0178     lat0                = -90 + 180*rand;
0179     Q.SENSOR_POS(2)     = lat0;
0180     if rand < 0.5
0181       Q.SENSOR_LOS      = -Q.SENSOR_LOS;
0182       aa0               = 180;
0183     end
0184   elseif atmosphere_dim == 3
0185     lat0                = -45 + 90*rand;    
0186     lon0                = -15 + 30*rand;    
0187     Q.SENSOR_POS(2)     = lat0;
0188     Q.SENSOR_POS(3)     = lon0;
0189     aa0                 = -180 + 360*rand;
0190 %    aa0 = 0
0191     Q.SENSOR_LOS(2)     = aa0;
0192   end
0193 
0194   %- Define calculation part
0195   %
0196   Q.WSMS_AT_START{end+1} = 'IndexSet(stokes_dim,1)';
0197   Q.WSMS_AT_START{end+1} = re_string;
0198   Q.WSMS_AT_END         = { 
0199     'VectorExtractFromMatrix(rte_pos,sensor_pos,0,"row")', ...
0200     'VectorExtractFromMatrix(rte_los,sensor_los,0,"row")', ...
0201     'VectorSet(rte_pos2,[])', ...
0202     'ppathCalc', ...
0203     sprintf('ppathWriteXMLPartial("ascii",ppath,"%s/ppath.xml")',workfolder) };
0204 
0205 
0206   %- Run ARTS and load result
0207   %
0208   S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0209                         'RteSet', 'CloseF' }, workfolder );
0210   cfile = fullfile( workfolder, 'cfile.arts' );
0211   strs2file( cfile, S );
0212   notok = arts( cfile, true );
0213   %
0214   if notok
0215     fprintf('\n!!! Error while running ARTS !!!\n');
0216     keyboard
0217   end
0218 
0219   if do_check
0220   
0221     ppath = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );
0222 
0223 
0224     %- Check result
0225     %
0226     if ppath.np > 1
0227       
0228       if atmosphere_dim == 1  ||  ~unsphere
0229         re = refell(1);
0230       else
0231         re = interp1( Q.LAT_GRID, ellipsoidradii(refell,Q.LAT_GRID), lat0 );
0232       end
0233       
0234       [x0,y0,z0,dx,dy,dz] = geocentricposlos2cart( re+Q.SENSOR_POS(1), ...
0235                                          lat0, lon0, abs(Q.SENSOR_LOS(1)), aa0 );
0236 
0237       % Check ppath points
0238       lt    = ppath.end_lstep + [ 0; cumsum(ppath.lstep) ];
0239       e     = zeros( size(ppath.pos,1), 1 );
0240       for i = 1 : size(ppath.pos,1)
0241               
0242         if atmosphere_dim == 3
0243           lons = ppath.pos(i,3);
0244         else
0245           lons = 0;
0246         end
0247         
0248         [xt,yt,zt] = geocentric2cart( ppath.r(i), ppath.pos(i,2), lons );
0249         dd = norm( [ xt-(x0+lt(i)*dx), yt-(y0+lt(i)*dy), zt-(z0+lt(i)*dz) ] );
0250         e(i) = dd;
0251         if dd > 10
0252           fprintf('\n!!! Deviating ppath point found !!!\n');
0253           [rt,latt,lont,zat,aat] = cartposlos2geocentric( x0+lt(i)*dx, ...
0254                            y0+lt(i)*dy, z0+lt(i)*dz, dx, dy, dz, ...
0255                            (re+Q.SENSOR_POS(1))*sind(abs(Q.SENSOR_LOS(1))), ...
0256                            lat0, lon0, abs(Q.SENSOR_LOS(1)), aa0);
0257           keyboard
0258         end
0259       end
0260       
0261       if strcmp( ppath.background, 'surface' )
0262         if max(e(end)) > emaxsurf, emaxsurf = e(end); end
0263         e(end) = 0;
0264       end
0265       %
0266       if max(e) > emaxatm, emaxatm = max(e); end
0267     end
0268   end
0269 end
0270 %
0271 if do_check
0272   fprintf( '\n' );
0273   fprintf( 'Max atmospheric deviation: %.1e m\n', emaxatm );
0274   fprintf( '    Max surface deviation: %.1e m\n', emaxsurf );
0275 else
0276   fprintf( '\n' );
0277 end
0278 
0279 
0280 %- Remove workfolder
0281 %
0282 delete_tmpfolder( workfolder );
0283 
0284 
0285 return
0286 for i = 1:122
0287   [x0,y0,z0] = geocentric2cart( ppath0.r(i), ppath0.pos(i,2), ppath0.pos(i,3) );
0288   [x,y,z] = geocentric2cart( ppath.r(i), ppath.pos(i,2), ppath.pos(i,3) );
0289   s = norm( [ x0-x, y0-y, z0-z ] );
0290   fprintf( '%3d: %.2e\n', i, s );
0291 end
0292 
0293

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