Home > atmlab > tests > test_arts_refraction.m

test_arts_refraction

PURPOSE ^

TEST_ARTS_REFRACTION Testing of ARTS ppaths with refraction

SYNOPSIS ^

function test_arts_refraction( varargin )

DESCRIPTION ^

 TEST_ARTS_REFRACTION   Testing of ARTS ppaths with refraction

   Makes 1D, 2D and 3D calculations with refraction and gives estimate
   of the error. The correct value is estimated by making calculations 
   with a 100 m step length, which makes the function slow.

 FORMAT   test_arts_refraction( [lrt,geomztan,switchaa] )
        
 OPT   lrt        Ray tracing step length
       geomztan   Geometrical tangent altitude
       switchaa   Make calculations with negative za for 2D and with a 180 
                  deg azimuth angle for 3D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

test_arts_refraction.m

SOURCE CODE ^

0001 % TEST_ARTS_REFRACTION   Testing of ARTS ppaths with refraction
0002 %
0003 %   Makes 1D, 2D and 3D calculations with refraction and gives estimate
0004 %   of the error. The correct value is estimated by making calculations
0005 %   with a 100 m step length, which makes the function slow.
0006 %
0007 % FORMAT   test_arts_refraction( [lrt,geomztan,switchaa] )
0008 %
0009 % OPT   lrt        Ray tracing step length
0010 %       geomztan   Geometrical tangent altitude
0011 %       switchaa   Make calculations with negative za for 2D and with a 180
0012 %                  deg azimuth angle for 3D
0013 
0014 % 2012-03-15   Created by Patrick.
0015 
0016 function test_arts_refraction( varargin )
0017 %
0018 [lrt,geomztan,switchaa] = optargs( varargin, { 2e3, 10e3, false } );
0019   
0020   
0021 %- Fixed settings
0022 %
0023 Q = qarts;
0024 Q.INCLUDES              = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0025                             fullfile( 'ARTS_INCLUDES', 'agendas.arts' ) };
0026 %
0027 Q.PPATH_AGENDA          = { 'ppath_agenda__FollowSensorLosPath'   };
0028 Q.PPATH_STEP_AGENDA     = { 'ppath_step_agenda__RefractedPath'    };
0029 Q.REFR_INDEX_AIR_AGENDA = { 'refr_index_air_agenda__GasThayer'    };
0030 %
0031 Q.F_GRID                = 100e9;
0032 Q.CLOUDBOX_DO           = false;
0033 Q.ABS_SPECIES(1).TAG{1} = 'H2O';
0034 refell                  = ellipsoidmodels( 'SphericalEarth' );
0035 re_string               = 'refellipsoidEarth( refellipsoid, "Sphere" )';
0036 
0037 
0038 %- Load a test atmosphere
0039 %
0040 atmtest =  fullfile( atmlab('ARTS_XMLDATA_PATH'), 'planets', 'Earth', ...
0041                                        'Fascod', 'tropical', 'tropical' );
0042 Az      = xmlLoad( [ atmtest, '.z.xml' ] ); 
0043 At      = xmlLoad( [ atmtest, '.t.xml' ] ); 
0044 Ah2o    = xmlLoad( [ atmtest, '.H2O.xml' ] ); 
0045 
0046 
0047 %- 1D atmosphere
0048 %
0049 Q.P_GRID              = z2p_simple( linspace( 0, 80e3, 320 )' );
0050 Q.Z_FIELD             = interpp( Az.grids{1}, Az.data, Q.P_GRID );
0051 Q.T_FIELD             = interpp( At.grids{1}, At.data, Q.P_GRID );
0052 Q.VMR_FIELD           = interpp( Ah2o.grids{1}, Ah2o.data, Q.P_GRID );
0053 Q.VMR_FIELD           = reshape( Q.VMR_FIELD, [1 length(Q.P_GRID) 1 1 ] );
0054   
0055 %- Surface altitude
0056 %
0057 Q.Z_SURFACE           = 500;
0058   
0059 %- Sensor pos/los
0060 Q.SENSOR_POS          = 600e3;
0061 Q.SENSOR_LOS          = geomztan2za( refell(1), Q.SENSOR_POS, geomztan );
0062 
0063 
0064 %- Determine expected tangent altitude
0065 %
0066 eh2o = Q.P_GRID .* Q.VMR_FIELD(1,:)';
0067 n    = 1 + ( 77.6e-8 * ( Q.P_GRID - eh2o ) + ...
0068                     ( 64.8e-8 + 3.776e-3 ./ Q.T_FIELD ) .* eh2o ) ./ Q.T_FIELD;
0069 x    = n .* (refell(1)+Q.Z_FIELD);
0070 ztan = interp1( x, Q.Z_FIELD, (refell(1)+Q.SENSOR_POS)*sind(Q.SENSOR_LOS), ...
0071                 'cubic' );
0072 % Cubic interpolation was found best in comparison to ppath0 result below,
0073 % closely followed by spline.
0074 
0075 
0076 %- Create a workfolder
0077 %
0078 workfolder = create_tmpfolder;
0079 
0080 
0081 %- Define calculation part
0082 %
0083 Q.WSMS_AT_START       = { 'IndexSet(stokes_dim,1)', re_string };
0084 Q.WSMS_AT_END         = { 
0085     'VectorExtractFromMatrix(rte_pos,sensor_pos,0,"row")', ...
0086     'VectorExtractFromMatrix(rte_los,sensor_los,0,"row")', ...
0087     'VectorSet(rte_pos2,[])', ...
0088     'ppathCalc', ...
0089     sprintf('ppathWriteXMLPartial("ascii",ppath,"%s/ppath.xml")',workfolder) };
0090 
0091 
0092 %- Run ARTS for 1D, with short step length
0093 %
0094 Q.ATMOSPHERE_DIM      = 1;
0095 Q.PPATH_LRAYTRACE     = 5;
0096 Q.PPATH_LMAX          = Q.PPATH_LRAYTRACE;
0097 %
0098 S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0099                       'RteSet', 'CloseF' }, workfolder );
0100 cfile = fullfile( workfolder, 'cfile.arts' );
0101 strs2file( cfile, S );
0102 notok = arts( cfile, true );
0103 %
0104 if notok
0105   fprintf('\n!!! Error while running ARTS for 1D !!!\n');
0106   keyboard
0107 end
0108 %
0109 ppath0 = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );
0110 %
0111 [ztan0,imin] = min(ppath0.pos(:,1));
0112 lat0         = ppath0.pos(imin,2);
0113 if switchaa
0114   lat0 = -lat0;
0115 end
0116 fprintf( '\nGeometrical tangent altitude : %8.3f km\n', geomztan/1e3 );
0117 fprintf( '  Refracted tangent altitude : %8.3f km\n', ztan0/1e3 );
0118 fprintf( 'Rough error estimate of ztan : %8.3f m\n',  ztan0-ztan );
0119 fprintf( '     Ray tracing step length : %8.3f km\n', lrt/1e3 );
0120   
0121 
0122 
0123 %- Run ARTS for 1D
0124 %
0125 Q.ATMOSPHERE_DIM      = 1;
0126 Q.PPATH_LRAYTRACE     = lrt;
0127 Q.PPATH_LMAX          = Q.PPATH_LRAYTRACE;
0128 %
0129 S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0130                       'RteSet', 'CloseF' }, workfolder );
0131 cfile = fullfile( workfolder, 'cfile.arts' );
0132 strs2file( cfile, S );
0133 notok = arts( cfile, true );
0134 %
0135 if notok
0136   fprintf('\n!!! Error while running ARTS for 1D !!!\n');
0137   keyboard
0138 end
0139 %
0140 ppath1 = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );
0141 %
0142 if switchaa
0143   ppath1.pos(:,2) = -ppath1.pos(:,2);
0144 end
0145 %
0146 [zmin1,imin] = min(ppath1.pos(:,1));
0147 lat1         = ppath1.pos(imin,2);
0148 %
0149 fprintf( '       Altitude error for 1D : %8.3f m\n', zmin1-ztan0 );
0150 fprintf( '       Latitude error for 1D : %8.3f m\n', (lat1-lat0)*111e3 );
0151   
0152 
0153 %- 2D
0154 %
0155 Q.ATMOSPHERE_DIM      = 2;
0156 %
0157 Q.SENSOR_POS(2)       = 0; 
0158 Q.LAT_GRID            = [ 10 : 35 ]';
0159 if switchaa
0160   Q.SENSOR_LOS(1)     = -Q.SENSOR_LOS(1);
0161   Q.LAT_GRID          = sort( -Q.LAT_GRID );
0162 end
0163 nlat                  = length( Q.LAT_GRID );
0164 Q.Z_FIELD             = repmat( Q.Z_FIELD, 1 , nlat );
0165 Q.T_FIELD             = repmat( Q.T_FIELD, 1 , nlat );
0166 Q.VMR_FIELD           = repmat( Q.VMR_FIELD, [ 1, 1 , nlat] );
0167 Q.Z_SURFACE           = repmat( Q.Z_SURFACE, nlat, 1 );
0168 %
0169 S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0170                       'RteSet', 'CloseF' }, workfolder );
0171 cfile = fullfile( workfolder, 'cfile.arts' );
0172 strs2file( cfile, S );
0173 notok = arts( cfile, true );
0174 %
0175 if notok
0176   fprintf('\n!!! Error while running ARTS for 2D !!!\n');
0177   keyboard
0178 end
0179 %
0180 ppath2 = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );
0181 %
0182 [zmin2,imin] = min(ppath2.pos(:,1));
0183 lat2         = ppath2.pos(imin,2);
0184 fprintf( '       Altitude error for 2D : %8.3f m\n', zmin2-ztan0 );
0185 fprintf( '       Latitude error for 2D : %8.3f m\n', (lat2-lat0)*111e3 );
0186 
0187 
0188 
0189 %- 3D
0190 %
0191 Q.ATMOSPHERE_DIM      = 3;
0192 %
0193 Q.SENSOR_POS(3)       = 10; 
0194 Q.SENSOR_LOS(2)       = 0; 
0195 if switchaa
0196   Q.SENSOR_LOS(1)     = -Q.SENSOR_LOS(1);
0197   Q.SENSOR_LOS(2)     = 180;   
0198 end
0199 Q.LON_GRID            = Q.SENSOR_POS(3) + [ -5 : 2 : 5 ]';
0200 nlon                  = length( Q.LON_GRID );
0201 Q.Z_FIELD             = repmat( Q.Z_FIELD, [ 1 , 1, nlon ] );
0202 Q.T_FIELD             = repmat( Q.T_FIELD, [ 1 , 1, nlon ] );
0203 Q.VMR_FIELD           = repmat( Q.VMR_FIELD, [ 1, 1 , 1, nlon] );
0204 Q.Z_SURFACE           = repmat( Q.Z_SURFACE, [ 1, nlon ] );
0205 
0206 %
0207 S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0208                       'RteSet', 'CloseF' }, workfolder );
0209 cfile = fullfile( workfolder, 'cfile.arts' );
0210 strs2file( cfile, S );
0211 notok = arts( cfile, true );
0212 %
0213 if notok
0214   fprintf('\n!!! Error while running ARTS for 3D !!!\n');
0215   keyboard
0216 end
0217 %
0218 ppath3 = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );
0219 %
0220 [zmin3,imin] = min(ppath3.pos(:,1));
0221 lat3         = ppath3.pos(imin,2);
0222 lon3         = ppath3.pos(imin,3)-Q.SENSOR_POS(3);
0223 fprintf( '       Altitude error for 3D : %8.3f m\n', zmin3-ztan0 );
0224 fprintf( '       Latitude error for 3D : %8.3f m\n', (lat3-lat0)*111e3 );
0225 fprintf( '      Longitude error for 3D : %8.3f m\n\n', lon3*111e3 );
0226   
0227 
0228 
0229 %- Remove workfolder
0230 %
0231 delete_tmpfolder( workfolder );
0232

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