Home > atmlab > arts_usage > arts_radioocc_1D_slta.m

arts_radioocc_1D_slta

PURPOSE ^

ARTS_RADIOOCC_1D_SLTA Simulation of 1D radio occultation

SYNOPSIS ^

function [R,T] = arts_radioocc_1D_slta(Q,O,A)

DESCRIPTION ^

 ARTS_RADIOOCC_1D_SLTA   Simulation of 1D radio occultation

    This function is very similar to *arts_radioocc_1d*, but here more of
    the calculations are done inside ARTS. The ARTS calculations involve to
    determine exact paths between fixed transmitter and receiver positions,
    and this function is slower than *arts_radioocc_1d*. 

    See *arts_radioocc_1d* for overall nomenclature, the basic geometry
    assumed, comments around *A* and the purpose of the structure *T*. Note
    that the occultation here is specified by SLTA, while in *arts_radioocc_1d*
    impact height is used.

    The basic results are gathered into the structure R having fields
       z_impact    : The impact altitude, ie. the non-refracted tangent 
                     altitude. Returned as monotonously decreasing values.
       z_tan       : Tangent altitude.
       lat_tan     : Latitude of tangent point.
       lat_distance: The latitude distance between the satellites.
       bangle      : Bending angle.
       slta        : Straight line tangent altitude.
       l_geometric : The (straight line) geometrical distance between the 
                     satellites.
       l_optpath   : The optical path length.
       l_proppath  : The propagation path length.
       t           : Time relative to the occultation with max z_impact.
       tr_total    : Total transmission over the link. This transmission
                     is the product of the ones below.
       tr_space    : Transmission due to free space loss alone.
       tr_atmos    : Transmission due to atmospheric attenuation  alone.
       tr_space    : Transmission due to defocusing alone.
       faraday     : Faraday rotation [deg]. Only included if O.do_faraday
                     is true.
   The data in *R* are for an equidistant *R.lat_distance*.

   The occultation is specified by the structure O having mandatory fields
      rec_altitude : Altitude of receiver orbit.
      tra_altitude : Altitude of transmitter orbit.
      frequency    : Frequency of signal.
      lmax         : As ARTS WSV *lmax*. This setting affects here mainly
                     accuracy of tangent point.
      lraytrace    : As ARTS WSV *lraytrace*. Determines the overall
                     accuracy. A value in the order of 100 m should be OK.
      z_surface    : Surface altitude.
      slta_max     : Max SLTA of occultation.
      slta_min     : Min SLTA of occultation.
      slta_n       : Number of SLTA to include.
      z_impact4t0  : The impact height used as time reference (t=0 here).
      f_sampling   : Sampling frequency for data in T.
   And optional fields
      rec_movement : String describing movement. Options are 'none',
                     'approaching' or 'disappearing'. The last is default.
      tra_movement : String describing movement. Options are 'none',
                     'approaching' or 'disappearing'. The last is default.
      do_faraday   : Flag to also calculate Faraday rotation. 
                     Default is false.
      defoc_method : Method to calculated defocusing. See iyRadioLInk.
                     Default is 2.
      defoc_shift  : Angular shift for defocusing. See iyRadioLInk.
                     Default is 3e-3.
      r_planet     : Radius of reference ellipsoid. If set this replaces
                     the default value, that varies with planet.
      n_agenda     : Agenda for refractive index. If set this replaces
                     the default value, that varies with planet.

 FORMAT   R = arts_radioocc_1D_slta(Q,O,A)
        
 OUT   R   Result structure, see above.
 IN    Q   Initial settings of Q. If empty: Q = qarts;   
       O   Occultation structure, see above. If empty, some default 
           values set that work as demonstration.
       A   Atmospheric structure, see above. If empty, some planet specific
           settings are applied. For example, for Earth Fascode tropical is
           used. Free electrons are not included automatically for any planet.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

arts_radioocc_1D_slta.m

SOURCE CODE ^

0001 % ARTS_RADIOOCC_1D_SLTA   Simulation of 1D radio occultation
0002 %
0003 %    This function is very similar to *arts_radioocc_1d*, but here more of
0004 %    the calculations are done inside ARTS. The ARTS calculations involve to
0005 %    determine exact paths between fixed transmitter and receiver positions,
0006 %    and this function is slower than *arts_radioocc_1d*.
0007 %
0008 %    See *arts_radioocc_1d* for overall nomenclature, the basic geometry
0009 %    assumed, comments around *A* and the purpose of the structure *T*. Note
0010 %    that the occultation here is specified by SLTA, while in *arts_radioocc_1d*
0011 %    impact height is used.
0012 %
0013 %    The basic results are gathered into the structure R having fields
0014 %       z_impact    : The impact altitude, ie. the non-refracted tangent
0015 %                     altitude. Returned as monotonously decreasing values.
0016 %       z_tan       : Tangent altitude.
0017 %       lat_tan     : Latitude of tangent point.
0018 %       lat_distance: The latitude distance between the satellites.
0019 %       bangle      : Bending angle.
0020 %       slta        : Straight line tangent altitude.
0021 %       l_geometric : The (straight line) geometrical distance between the
0022 %                     satellites.
0023 %       l_optpath   : The optical path length.
0024 %       l_proppath  : The propagation path length.
0025 %       t           : Time relative to the occultation with max z_impact.
0026 %       tr_total    : Total transmission over the link. This transmission
0027 %                     is the product of the ones below.
0028 %       tr_space    : Transmission due to free space loss alone.
0029 %       tr_atmos    : Transmission due to atmospheric attenuation  alone.
0030 %       tr_space    : Transmission due to defocusing alone.
0031 %       faraday     : Faraday rotation [deg]. Only included if O.do_faraday
0032 %                     is true.
0033 %   The data in *R* are for an equidistant *R.lat_distance*.
0034 %
0035 %   The occultation is specified by the structure O having mandatory fields
0036 %      rec_altitude : Altitude of receiver orbit.
0037 %      tra_altitude : Altitude of transmitter orbit.
0038 %      frequency    : Frequency of signal.
0039 %      lmax         : As ARTS WSV *lmax*. This setting affects here mainly
0040 %                     accuracy of tangent point.
0041 %      lraytrace    : As ARTS WSV *lraytrace*. Determines the overall
0042 %                     accuracy. A value in the order of 100 m should be OK.
0043 %      z_surface    : Surface altitude.
0044 %      slta_max     : Max SLTA of occultation.
0045 %      slta_min     : Min SLTA of occultation.
0046 %      slta_n       : Number of SLTA to include.
0047 %      z_impact4t0  : The impact height used as time reference (t=0 here).
0048 %      f_sampling   : Sampling frequency for data in T.
0049 %   And optional fields
0050 %      rec_movement : String describing movement. Options are 'none',
0051 %                     'approaching' or 'disappearing'. The last is default.
0052 %      tra_movement : String describing movement. Options are 'none',
0053 %                     'approaching' or 'disappearing'. The last is default.
0054 %      do_faraday   : Flag to also calculate Faraday rotation.
0055 %                     Default is false.
0056 %      defoc_method : Method to calculated defocusing. See iyRadioLInk.
0057 %                     Default is 2.
0058 %      defoc_shift  : Angular shift for defocusing. See iyRadioLInk.
0059 %                     Default is 3e-3.
0060 %      r_planet     : Radius of reference ellipsoid. If set this replaces
0061 %                     the default value, that varies with planet.
0062 %      n_agenda     : Agenda for refractive index. If set this replaces
0063 %                     the default value, that varies with planet.
0064 %
0065 % FORMAT   R = arts_radioocc_1D_slta(Q,O,A)
0066 %
0067 % OUT   R   Result structure, see above.
0068 % IN    Q   Initial settings of Q. If empty: Q = qarts;
0069 %       O   Occultation structure, see above. If empty, some default
0070 %           values set that work as demonstration.
0071 %       A   Atmospheric structure, see above. If empty, some planet specific
0072 %           settings are applied. For example, for Earth Fascode tropical is
0073 %           used. Free electrons are not included automatically for any planet.
0074 
0075 % 2013-09-29   Created by Patrick Eriksson.
0076 
0077 function [R,T] = arts_radioocc_1D_slta(Q,O,A)
0078 %
0079 if isempty(Q)
0080   Q = qarts;
0081 end
0082 %
0083 if isempty(O)
0084   O.rec_altitude = 820e3;
0085   O.tra_altitude = 20200e3;
0086   O.frequency    = 1575.42e6;
0087   O.lmax         = 2e3;
0088   O.lraytrace    = 100;
0089   O.z_surface    = 0;
0090 
0091   O.slta_max     = 30e3;
0092   O.slta_min     = -100e3;
0093   O.slta_n       = 101;
0094 
0095   O.z_impact4t0  = 20e3;
0096   O.f_sampling   = 4;  
0097 end
0098 %
0099 if isempty(A)
0100   A.planet       = 'earth';
0101   A.atmfunc      = @qarts_add_fascode;
0102   A.fascode_atm  = 'tropical';
0103 end
0104 
0105 
0106 % Check movements
0107 %
0108 rec_movement = safegetfield( O, 'rec_movement', 'disappearing' );
0109 tra_movement = safegetfield( O, 'tra_movement', 'disappearing' );
0110 %
0111 if strcmp( rec_movement, 'none' )  &  strcmp( tra_movement, 'none' )
0112   error( 'Both O.rec_movement and O.tra_movement can not be ''none'',' );
0113 end
0114 
0115 
0116 %- Create a temporary workfolder
0117 %
0118 workfolder = create_tmpfolder;
0119 cu = onCleanup( @()delete_tmpfolder( workfolder ) );
0120 %
0121 cfile = fullfile( workfolder, 'cfile.arts' );
0122 
0123 
0124 %- Set basic and observation related fields of Q
0125 %
0126 Q = q_basic_local( Q, O, workfolder );
0127 
0128 
0129 %- Set basic and ppath parts of Q
0130 %
0131 [Q,M] = q_planet_local( Q, A, workfolder );
0132 %
0133 if isfield(O,'r_planet') & ~isempty(O.r_planet)
0134   Q.REFELLIPSOID = [ O.r_planet, 0 ]';
0135 end
0136 %
0137 if isfield(O,'n_agenda') & ~isempty(O.n_agenda)
0138   Q.REFR_INDEX_AIR_AGENDA = O.n_agenda;
0139 end
0140 
0141 
0142 %- Radiii
0143 %
0144 r_rec     = Q.REFELLIPSOID(1) + O.rec_altitude;
0145 r_tra     = Q.REFELLIPSOID(1) + O.tra_altitude;
0146 
0147 
0148 %- Determine latitude distances to use and set-up R
0149 %
0150 dlat_min = acosd( (Q.REFELLIPSOID(1)+O.slta_max) / r_tra ) + ...
0151            acosd( (Q.REFELLIPSOID(1)+O.slta_max) / r_rec );
0152 dlat_max = acosd( (Q.REFELLIPSOID(1)+O.slta_min) / r_tra ) + ...
0153            acosd( (Q.REFELLIPSOID(1)+O.slta_min) / r_rec );
0154 %
0155 R.lat_distance = linspace( dlat_min, dlat_max, O.slta_n );
0156 n0             = length( R.lat_distance );
0157 %
0158 [R.z_tan,R.lat_tan,R.bangle,R.z_impact]         = deal( repmat( NaN, n0, 1 ) );
0159 [R.slta,R.l_geometric,R.l_optpath,R.l_proppath] = deal( repmat( NaN, n0, 1 ) );
0160 [R.tr_total,R.tr_space,R.tr_atmos,R.tr_defoc]   = deal( repmat( NaN, n0, 1 ) );
0161 %
0162 do_faraday = safegetfield( O, 'do_faraday', false );
0163 %
0164 if do_faraday
0165   R.faraday = deal( repmat( NaN, n0, 1 ) );
0166 end
0167 
0168 
0169 %- Run ARTS
0170 %
0171 for i = 1 : n0
0172 
0173   Q.TRANSMITTER_POS = [ O.rec_altitude, R.lat_distance(i) ];
0174 
0175   % Run ARTS
0176   %
0177   S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0178                         'RteSet', 'CloseF' }, workfolder );
0179   strs2file( cfile, S );  
0180   notok = arts( cfile, true );
0181   %
0182   if notok
0183     error('\n!!! Error while running ARTS !!!\n');
0184   end
0185 
0186   ppath  = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );   
0187   iy     = xmlLoad( fullfile( workfolder, 'iy.xml' ) );   
0188   iy_aux = xmlLoad( fullfile( workfolder, 'iy_aux.xml' ) );   
0189 
0190   if strcmp( ppath.background, 'transmitter' )
0191   
0192     % Extract tangent point
0193     [u,j]        = min( ppath.pos(:,1) );
0194     R.z_tan(i)   = ppath.pos(j,1);
0195     R.lat_tan(i) = ppath.pos(j,2);
0196 
0197     % Impact height
0198     R.z_impact(i) = ppath.constant - Q.REFELLIPSOID(1);
0199     
0200     % Bending angle
0201     R.bangle(i) = ppath.start_los(1) - ppath.end_los(1) + ppath.start_pos(2);
0202 
0203     % Path lenghts
0204     lppath = ppath.start_lstep + ppath.end_lstep + sum( ppath.lstep );
0205     lopt   = ppath.start_lstep + ppath.end_lstep + sum( ppath.lstep .* ...
0206                        ( 0.5 * (ppath.ngroup(1:end-1)+ppath.ngroup(2:end)) ) );
0207 
0208     % SLTA and geometrical distance
0209     lgeom        = sqrt( r_rec^2 + r_tra^2 - 2*r_rec*r_tra*...
0210                                                      cosd(R.lat_distance(i)) );
0211     R.slta(i)    = r_rec * r_tra * sind(R.lat_distance(i)) / lgeom - ...
0212                                                              Q.REFELLIPSOID(1);
0213     R.l_proppath(i)  = lppath;
0214     R.l_optpath(i)   = lopt;
0215     R.l_geometric(i) = lgeom;
0216 
0217     R.tr_total(i)  = iy(1);
0218     R.tr_space(i)  = iy_aux{1};
0219     R.tr_atmos(i)  = iy_aux{2};
0220     R.tr_defoc(i)  = iy_aux{3}; 
0221 
0222     if do_faraday
0223       R.faraday(i) = iy_aux{4}; 
0224     end
0225   end
0226 end
0227 
0228 
0229 % Speed of increase in angular distance
0230 %
0231 G         = constants( 'GRAVITATIONAL_CONSTANT' );
0232 t_rec     = sqrt( (4*pi^2*r_rec^3) / (G*M) );
0233 %
0234 switch lower(rec_movement)
0235  case 'none'
0236   degpersec = 0;
0237  case 'approaching'
0238   degpersec = - 360 / t_rec;
0239  case 'disappearing'
0240   degpersec = + 360 / t_rec;
0241  otherwise
0242   error( 'Unrecognised choice for *O.rec_movement* (%s).', rec_movement );
0243 end
0244 %
0245 t_tra     = sqrt( (4*pi^2*r_tra^3) / (G*M) );
0246 %
0247 switch lower(tra_movement)
0248  case 'none'
0249   ;
0250  case 'approaching'
0251   degpersec = degpersec - 360 / t_tra;
0252  case 'disappearing'
0253   degpersec = degpersec + 360 / t_tra;
0254  otherwise
0255   error( 'Unrecognised choice for *O.tra_movement* (%s).', tra_movement );
0256 end
0257 
0258 
0259 % Relative arrival time
0260 %
0261 ind = find( ~isnan( R.z_impact ) );
0262 %
0263 R.t = ( R.lat_distance - interp1( R.z_impact(ind), R.lat_distance(ind), ...
0264                                                 O.z_impact4t0 ) ) / degpersec;
0265 
0266 if nargout == 1
0267   return;   % --->
0268 end
0269 
0270 % Create T
0271 %
0272 % Time grid to use
0273 T.t   = min(R.t) : 1/O.f_sampling : max(R.t);
0274 %
0275 % Find start and end of each time chunck
0276 istart = [];
0277 iend   = [];
0278 i1     = NaN; % Start of present chunck%
0279 
0280 for i = 1 : n0
0281   % If Nan, then we are inbetween chuncks
0282   if isnan( R.t(i) )
0283     i1 = NaN;
0284 
0285   else
0286     % If i1==NaN, this is the start of a new chunck
0287     if isnan(i1)
0288       i1 = i;
0289     end
0290     
0291     % End point if last point or next is NaN
0292     if i == n0 | isnan(R.t(i+1))
0293       if i > i1 % Ignore length-1 chuncks
0294         istart = [ istart i1 ];
0295         iend   = [ iend   i  ];
0296       end
0297       i1     = NaN;
0298     
0299     % Also end point if sign change of t-derivate
0300     elseif i>1  &  ~isnan(R.t(i-1))  &  ...
0301           sign(R.t(i+1)-R.t(i)) ~= sign(R.t(i)-R.t(i-1))
0302       istart = [ istart i1 ];
0303       iend   = [ iend   i  ];
0304       i1     = i;      
0305     end    
0306   end
0307 end
0308 %
0309 np    = length( istart );
0310 nt    = length( T.t );
0311 %
0312 [T.z_tan,T.lat_tan,T.bangle,T.lat_distance]   = deal( repmat( NaN, nt, np ) );
0313 [T.slta,T.l_geometric,T.l_proppath]           = deal( repmat( NaN, nt, np ) );
0314 [T.l_optpath,T.defocus,T.z_impact]            = deal( repmat( NaN, nt, np ) );
0315 [T.tr_total,T.tr_space,T.tr_atmos,T.tr_defoc] = deal( repmat( NaN, nt, np ) );
0316 if do_faraday
0317   T.faraday = deal( repmat( NaN, nt, np ) );
0318 end
0319 %
0320 for i = 1 : np
0321 
0322   % Interpolate over time of chunk
0323   ir = istart(i) : iend(i);
0324   tmin = min( [ R.t(ir([1,end])) ] );
0325   tmax = max( [ R.t(ir([1,end])) ] );
0326   it = find( T.t >= tmin  &  T.t < tmax );  
0327   
0328   T.z_tan(it,i)        = interp1( R.t(ir), R.z_tan(ir),        T.t(it) );
0329   T.lat_tan(it,i)      = interp1( R.t(ir), R.lat_tan(ir),      T.t(it) );
0330   T.lat_distance(it,i) = interp1( R.t(ir), R.lat_distance(ir), T.t(it) );
0331   T.bangle(it,i)       = interp1( R.t(ir), R.bangle(ir),       T.t(it) );
0332   T.slta(it,i)         = interp1( R.t(ir), R.slta(ir),         T.t(it) );
0333   T.l_geometric(it,i)  = interp1( R.t(ir), R.l_geometric(ir),  T.t(it) );
0334   T.l_proppath(it,i)   = interp1( R.t(ir), R.l_proppath(ir),   T.t(it) );
0335   T.l_optpath(it,i)    = interp1( R.t(ir), R.l_optpath(ir),    T.t(it) );
0336   T.z_impact(it,i)     = interp1( R.t(ir), R.z_impact(ir),     T.t(it) );
0337 
0338   T.tr_total(it,i)     = interp1( R.t(ir), R.tr_total(ir),     T.t(it) );
0339   T.tr_space(it,i)     = interp1( R.t(ir), R.tr_space(ir),     T.t(it) );
0340   T.tr_atmos(it,i)     = interp1( R.t(ir), R.tr_atmos(ir),     T.t(it) );
0341   T.tr_defoc(it,i)     = interp1( R.t(ir), R.tr_defoc(ir),     T.t(it) );
0342 
0343   if do_faraday
0344     T.faraday(it,i)    = interp1( R.t(ir), R.faraday(ir),      T.t(it) );
0345   end
0346 end
0347 
0348 
0349 % Defocusing
0350 %
0351 for ip = 1 : np
0352   for it = 1 : nt
0353     
0354     if ~isnan( T.bangle(it,ip) )
0355       
0356       r_tan2 = ( Q.REFELLIPSOID(1) + T.z_tan(it,ip) )^2;
0357       lr     = sqrt( r_rec^2 - r_tan2 ); 
0358       lt     = sqrt( r_tra^2 - r_tan2 ); 
0359       
0360       fac = (pi/180) * (lr*lt) / (lr+lt);
0361 
0362       if it == 1  |  isnan(T.bangle(it-1,ip))
0363         if it < nt   % If at upper end, we have to leave a NaN
0364           T.defocus(it,ip) = 1 / ( 1 - fac * ...
0365                             (  T.bangle(it+1,ip)  -  T.bangle(it,ip)  ) /  ...
0366                             ( T.z_impact(it+1,ip) - T.z_impact(it,ip) ) );
0367         end
0368       elseif it == nt  |  isnan(T.bangle(it+1,ip))
0369         if it > 1  % If at lower end, we have to leave a NaN
0370           T.defocus(it,ip) = 1 / ( 1 - fac * ...
0371                             (  T.bangle(it,ip)  -  T.bangle(it-1,ip)  ) /  ...
0372                             ( T.z_impact(it,ip) - T.z_impact(it-1,ip) ) );
0373         end
0374       else
0375         T.defocus(it,ip) = 1 / ( 1 - fac * ...
0376                           (  T.bangle(it+1,ip)  -  T.bangle(it-1,ip)  ) /  ...
0377                           ( T.z_impact(it+1,ip) - T.z_impact(it-1,ip) ) );
0378       end
0379     end
0380   end
0381 end
0382 
0383 
0384 R = orderfields( R );
0385 T = orderfields( T );
0386 
0387 return
0388 
0389 
0390 %----------------------------------------------------------------------------
0391 function Q = q_basic_local(Q,O,workfolder)
0392   %
0393   do_faraday = safegetfield( O, 'do_faraday', false );
0394   %
0395   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'general.arts' );
0396   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'agendas.arts' ); 
0397   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'continua.arts' );
0398   %
0399   Q.ATMOSPHERE_DIM        = 1;
0400   if do_faraday
0401     Q.STOKES_DIM          = 3;
0402   else
0403     Q.STOKES_DIM          = 1;
0404   end
0405   Q.F_GRID                = O.frequency;
0406   Q.PPATH_LMAX            = O.lmax;
0407   Q.PPATH_LRAYTRACE       = O.lraytrace;
0408   %
0409   Q.CLOUDBOX_DO           = false;
0410   Q.J_DO                  = false;
0411   %
0412   Q.PPATH_AGENDA          = { 'ppath_agenda__TransmitterReceiverPath'   };
0413   Q.PPATH_STEP_AGENDA     = { 'ppath_step_agenda__RefractedPath'    };
0414   %
0415   Q.SENSOR_POS            = O.tra_altitude;
0416   Q.SENSOR_LOS            = [];
0417   %
0418   Q.Z_SURFACE             = O.z_surface;
0419   %
0420   Q.WSMS_AT_END{end+1} = 'VectorSet(rte_los,[])';
0421   Q.WSMS_AT_END{end+1} = 'VectorExtractFromMatrix(rte_pos,sensor_pos,0,"row")';
0422   Q.WSMS_AT_END{end+1} = ...
0423                    'VectorExtractFromMatrix(rte_pos2,transmitter_pos,0,"row")';
0424   %
0425   Q.ABSORPTION         = 'OnTheFly';
0426   Q.ABS_LINES_FORMAT   = 'None';
0427   %
0428   if do_faraday
0429     Q.PROPMAT_CLEARSKY_AGENDA = { ...
0430                                  'propmat_clearsky_agenda__OnTheFly_Faraday' };
0431   else
0432     Q.PROPMAT_CLEARSKY_AGENDA = { 'propmat_clearsky_agenda__OnTheFly' };
0433   end
0434   %
0435   Q.IY_TRANSMITTER_AGENDA = { 'iy_transmitter_agenda__UnitUnpolIntensity' };
0436   %
0437   Q.WSMS_AT_END{end+1} = 'AgendaSet( iy_main_agenda ){';
0438   dmethod = safegetfield( O, 'defoc_method', 2 );
0439   dshift  = safegetfield( O, 'defoc_shift', 3e-3 );
0440   Q.WSMS_AT_END{end+1} = sprintf(...
0441       '  iyRadioLink( defocus_method = %d, defocus_shift = %.3e )',...
0442                                                              dmethod, dshift );
0443   Q.WSMS_AT_END{end+1} = '}';
0444   %
0445   saux = '"Free space loss","Atmospheric loss","Defocusing loss"';
0446   %
0447   if do_faraday
0448     saux = [ saux, ',"Faraday rotation"' ];
0449   end
0450   Q.WSMS_AT_END{end+1} = [ 'ArrayOfStringSet(iy_aux_vars,[',saux,'])' ];
0451   %
0452   Q.WSMS_AT_END{end+1} = 'iyCalc';
0453   Q.WSMS_AT_END{end+1} = ...
0454                       sprintf('WriteXML("ascii",iy,"%s/iy.xml")',workfolder);
0455   Q.WSMS_AT_END{end+1} = ...
0456                 sprintf('WriteXML("ascii",iy_aux,"%s/iy_aux.xml")',workfolder);
0457   Q.WSMS_AT_END{end+1} = ...
0458       sprintf('ppathWriteXMLPartial("ascii",ppath,"%s/ppath.xml")',workfolder);
0459   %
0460 return
0461 
0462 
0463 %----------------------------------------------------------------------------
0464 function [Q,M] = q_planet_local(Q,A,workfolder)
0465   %
0466   rqre_field( A, 'planet' );
0467   rqre_field( A, 'atmfunc' );
0468   %
0469   if ~isa( A.atmfunc, 'function_handle' )
0470     error( 'The field *A.atmfunc* must be a function handle' ); 
0471   end
0472   if ~exist( func2str(A.atmfunc), 'file' )
0473     error( 'The function *%s* is not found (selected by A.atmfunc).', ...
0474                                                         func2str(A.atmfunc) ); 
0475   end
0476 
0477   % Set basics and default for planet:
0478   %
0479   switch upper( A.planet )
0480 
0481    case 'VENUS'
0482     %
0483     M = 4.8676e24;
0484     %
0485     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_venus.arts' );
0486     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalVenus' );
0487    
0488    case 'EARTH'
0489     %
0490     M = constants( 'EARTH_MASS' );
0491     %
0492     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' );
0493     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalEarth' );
0494     %
0495     if ~isfield( A, 'ABS_SPECIES' )
0496       A.ABS_SPECIES(1).TAG{1} = 'N2-SelfContStandardType';
0497       A.ABS_SPECIES(2).TAG{1} = 'O2-PWR98';
0498       A.ABS_SPECIES(3).TAG{1} = 'H2O-PWR98';
0499       A.ABS_SPECIES(4).TAG{1} = 'free_electrons';
0500     end
0501 
0502    case 'MARS'
0503     %
0504     M = 6.4185e23;
0505     %
0506     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_mars.arts' );
0507     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalEarth' );
0508     
0509    case 'JUPITER'
0510     %
0511     M = 1.8986e27;
0512     %
0513     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_jupiter.arts' );
0514     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalJupiter' );
0515    
0516    otherwise
0517     error( 'Unrecognised choice for A.planet (%s).', A.planet );
0518   end
0519 
0520   switch upper( A.planet )
0521    case 'EARTH'
0522     %
0523     Q.REFR_INDEX_AIR_AGENDA = { 
0524         'NumericSet( refr_index_air, 1.0 )', ...
0525         'NumericSet( refr_index_air_group, 1.0 )',...
0526         'refr_index_airThayer', ...
0527         'refr_index_airFreeElectrons( demand_vmr_value = 0 )' };
0528     %
0529    otherwise 
0530     %
0531     Q.REFR_INDEX_AIR_AGENDA = { 
0532         'NumericSet( refr_index_air, 1.0 )', ...
0533         'NumericSet( refr_index_air_group, 1.0 )',...
0534         'refr_index_airMWgeneral', ...
0535         'refr_index_airFreeElectrons( demand_vmr_value = 0 )' };
0536     %
0537   end
0538   if 0   % Set to 1 if you want to test without refraction
0539     Q.REFR_INDEX_AIR_AGENDA = { 'refr_index_air_agenda__NoRefrac' };
0540   end
0541   
0542   % Include atmospheric specific data
0543   Q = A.atmfunc( A, Q, workfolder );
0544   
0545     
0546 return

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