Home > atmlab > arts_usage > arts_radioocc_1D.m

arts_radioocc_1D

PURPOSE ^

ARTS_RADIOOCC_1D Simulation of 1D radio occultation, no power

SYNOPSIS ^

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

DESCRIPTION ^

 ARTS_RADIOOCC_1D   Simulation of 1D radio occultation, no power

    The function simulates an occultation measurement between two
    satellites, denoted as transmitter and receiver, respectively.  
    
    For NWP radio occultation the transmitter is the GPS satellite and the
    receiver the LEO satellite. The assumed geometry is explained using this
    case. The occultation is considered to take place exactly in the LEO orbit
    plane. Angular distances are denoted as "latitude" following the ARTS
    nomenclature. The latitude gives here the angular distance from the GPS
    satellite. The flight direction of both satellites can be adjusted.

    If either receiver or transmitter is not in orbit around the planet,
    but for example placed on Earth, use then the option for 'none'
    movement. 

    Defocusing is calculated inside the function (ie. not inside ARTS),
    where in the ARTS nomenclature "method 2" is used (note the field
    O.defoc_shift, that influences the result obtained). 

    Settings and output are first described assuming that attenuation
    variables are not considered.

    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.
   The data in *R* are for an equidistant *R.z_impact*.

   The results are also returned as a function of time in the structure T. 
   The field T.t is the time grid used. The spacing of this grid follows
   O.f_sampling. The remaining fields of T are matrices (in contrast to R 
   only having vector fields), where the columns represent "multi-pathing". 
   That is, each column represents a path (not to be confused with the 
   pencil beam paths used i ARTS). The limit between these paths correspond 
   to discontinuities in the arrival time (for pencil beam paths).
   This structure T contains the same fields as R. 

   The occultation is specified by the structure O having the fields
      tra_altitude : Altitude of transmitter orbit.
      rec_altitude : Altitude of receiver 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.
      z_impact_max : Start altitude of occultation.
      z_impact_dz  : Vertical spacing between path calculations.
      z_impact4t0  : The impact height used as time reference (t=0 here).
      f_sampling   : Sampling frequency for data in T.
   And optional fields
      z_impact_min : End altitude of occultation. Default is to use z_surface.
      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.
      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.

   If O.do_atten is set to true, R also has these fields
      tr_atmos : Transmission due to attenuation of atmospheric
                 constituents (absorption and scattering).
      tr_defoc : Transmission due to defocusing.
      tr_space : Transmission due to free space loss (defined as 1/(4*pi*l*l))
      tr_total : Total transmission, the product of the terms above
   If also O.do_absspecies is, R includes
      tr_absspecies : The atmospheric attenuation for each absorption
                      species separately. This field is ignored for T.
   If also O.do_faraday is true, R includes
      faraday  : Total Faraday rotation [deg].
   Default for both O.do_atten and O.do_faraday is false. O.do_absspecies
   and O.do_faraday only considered if O.do_atten is true. 
   Additional settings when O.do_atten is true
      defoc_shift  : Angular shift for defocusing. See iyRadioLInk.
                     Default is 3e-3.

   Finally, the atmosphere is specified by the structure A having at least
   these fields
      planet  : The planet to model, such as 'earth'.
      atmfunc : Handle to function inserting an atmosphere, such as
                @arts_include_fascode.
   The structure A is passed on to the atmosphere function. That is, A
   likely also includes fields required by the "atmfunc". See
   *arts_include_fascode* for how such a function is assumed to work

 FORMAT   R = arts_radioocc_1D(Q,O,A[,workfolder])
        
 OUT   R   Result structure, see above.
 IN    Q   Initial settings of Q. If empty: Q = qarts;   
           Note that not all fields are considered, this including
           WSMS_BEFORE_RTE and YCALC_WSMS.
       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.
 OPT   workfolder   Path to a folder where to place all files generated.
                    Default is to create a temporary folder, following the
                    general atmlab settings.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

arts_radioocc_1D.m

SOURCE CODE ^

0001 % ARTS_RADIOOCC_1D   Simulation of 1D radio occultation, no power
0002 %
0003 %    The function simulates an occultation measurement between two
0004 %    satellites, denoted as transmitter and receiver, respectively.
0005 %
0006 %    For NWP radio occultation the transmitter is the GPS satellite and the
0007 %    receiver the LEO satellite. The assumed geometry is explained using this
0008 %    case. The occultation is considered to take place exactly in the LEO orbit
0009 %    plane. Angular distances are denoted as "latitude" following the ARTS
0010 %    nomenclature. The latitude gives here the angular distance from the GPS
0011 %    satellite. The flight direction of both satellites can be adjusted.
0012 %
0013 %    If either receiver or transmitter is not in orbit around the planet,
0014 %    but for example placed on Earth, use then the option for 'none'
0015 %    movement.
0016 %
0017 %    Defocusing is calculated inside the function (ie. not inside ARTS),
0018 %    where in the ARTS nomenclature "method 2" is used (note the field
0019 %    O.defoc_shift, that influences the result obtained).
0020 %
0021 %    Settings and output are first described assuming that attenuation
0022 %    variables are not considered.
0023 %
0024 %    The basic results are gathered into the structure R having fields
0025 %       z_impact    : The impact altitude, ie. the non-refracted tangent
0026 %                     altitude. Returned as monotonously decreasing values.
0027 %       z_tan       : Tangent altitude.
0028 %       lat_tan     : Latitude of tangent point.
0029 %       lat_distance: The latitude distance between the satellites.
0030 %       bangle      : Bending angle.
0031 %       slta        : Straight line tangent altitude.
0032 %       l_geometric : The (straight line) geometrical distance between the
0033 %                     satellites.
0034 %       l_optpath   : The optical path length.
0035 %       l_proppath  : The propagation path length.
0036 %       t           : Time relative to the occultation with max z_impact.
0037 %   The data in *R* are for an equidistant *R.z_impact*.
0038 %
0039 %   The results are also returned as a function of time in the structure T.
0040 %   The field T.t is the time grid used. The spacing of this grid follows
0041 %   O.f_sampling. The remaining fields of T are matrices (in contrast to R
0042 %   only having vector fields), where the columns represent "multi-pathing".
0043 %   That is, each column represents a path (not to be confused with the
0044 %   pencil beam paths used i ARTS). The limit between these paths correspond
0045 %   to discontinuities in the arrival time (for pencil beam paths).
0046 %   This structure T contains the same fields as R.
0047 %
0048 %   The occultation is specified by the structure O having the fields
0049 %      tra_altitude : Altitude of transmitter orbit.
0050 %      rec_altitude : Altitude of receiver orbit.
0051 %      frequency    : Frequency of signal.
0052 %      lmax         : As ARTS WSV *lmax*. This setting affects here mainly
0053 %                     accuracy of tangent point.
0054 %      lraytrace    : As ARTS WSV *lraytrace*. Determines the overall
0055 %                     accuracy. A value in the order of 100 m should be OK.
0056 %      z_surface    : Surface altitude.
0057 %      z_impact_max : Start altitude of occultation.
0058 %      z_impact_dz  : Vertical spacing between path calculations.
0059 %      z_impact4t0  : The impact height used as time reference (t=0 here).
0060 %      f_sampling   : Sampling frequency for data in T.
0061 %   And optional fields
0062 %      z_impact_min : End altitude of occultation. Default is to use z_surface.
0063 %      rec_movement : String describing movement. Options are 'none',
0064 %                     'approaching' or 'disappearing'. The last is default.
0065 %      tra_movement : String describing movement. Options are 'none',
0066 %                     'approaching' or 'disappearing'. The last is default.
0067 %      r_planet     : Radius of reference ellipsoid. If set this replaces
0068 %                     the default value, that varies with planet.
0069 %      n_agenda     : Agenda for refractive index. If set this replaces
0070 %                     the default value, that varies with planet.
0071 %
0072 %   If O.do_atten is set to true, R also has these fields
0073 %      tr_atmos : Transmission due to attenuation of atmospheric
0074 %                 constituents (absorption and scattering).
0075 %      tr_defoc : Transmission due to defocusing.
0076 %      tr_space : Transmission due to free space loss (defined as 1/(4*pi*l*l))
0077 %      tr_total : Total transmission, the product of the terms above
0078 %   If also O.do_absspecies is, R includes
0079 %      tr_absspecies : The atmospheric attenuation for each absorption
0080 %                      species separately. This field is ignored for T.
0081 %   If also O.do_faraday is true, R includes
0082 %      faraday  : Total Faraday rotation [deg].
0083 %   Default for both O.do_atten and O.do_faraday is false. O.do_absspecies
0084 %   and O.do_faraday only considered if O.do_atten is true.
0085 %   Additional settings when O.do_atten is true
0086 %      defoc_shift  : Angular shift for defocusing. See iyRadioLInk.
0087 %                     Default is 3e-3.
0088 %
0089 %   Finally, the atmosphere is specified by the structure A having at least
0090 %   these fields
0091 %      planet  : The planet to model, such as 'earth'.
0092 %      atmfunc : Handle to function inserting an atmosphere, such as
0093 %                @arts_include_fascode.
0094 %   The structure A is passed on to the atmosphere function. That is, A
0095 %   likely also includes fields required by the "atmfunc". See
0096 %   *arts_include_fascode* for how such a function is assumed to work
0097 %
0098 % FORMAT   R = arts_radioocc_1D(Q,O,A[,workfolder])
0099 %
0100 % OUT   R   Result structure, see above.
0101 % IN    Q   Initial settings of Q. If empty: Q = qarts;
0102 %           Note that not all fields are considered, this including
0103 %           WSMS_BEFORE_RTE and YCALC_WSMS.
0104 %       O   Occultation structure, see above. If empty, some default
0105 %           values set that work as demonstration.
0106 %       A   Atmospheric structure, see above. If empty, some planet specific
0107 %           settings are applied. For example, for Earth Fascode tropical is
0108 %           used. Free electrons are not included automatically for any planet.
0109 % OPT   workfolder   Path to a folder where to place all files generated.
0110 %                    Default is to create a temporary folder, following the
0111 %                    general atmlab settings.
0112 
0113 % 2013-10-16   Created by Patrick Eriksson.
0114 
0115 function [R,T] = arts_radioocc_1D(Q,O,A,workfolder)
0116 %
0117 if isempty(Q)
0118   Q = qarts;
0119 end
0120 %
0121 if isempty(O)
0122   O.rec_altitude = 820e3;
0123   O.tra_altitude = 20200e3;
0124   O.tra_movement = 'disappearing';
0125   O.frequency    = 1575.42e6;
0126   O.lmax         = 2e3;
0127   O.lraytrace    = 100;
0128   O.z_surface    = 0;
0129 
0130   O.z_impact_max = 30e3;
0131   O.z_impact_dz  = 200;
0132   O.z_impact4t0  = 20e3;
0133   O.f_sampling   = 4;
0134 
0135 end
0136 %
0137 if isempty(A)
0138   A.planet       = 'earth';
0139   A.atmfunc      = @qarts_add_fascode;
0140   A.fascode_atm  = 'tropical';
0141 end
0142 %
0143 if nargin < 4, workfolder = []; end
0144 
0145 
0146 % Check movements
0147 %
0148 rec_movement = safegetfield( O, 'rec_movement', 'disappearing' );
0149 tra_movement = safegetfield( O, 'tra_movement', 'disappearing' );
0150 %
0151 if strcmp( rec_movement, 'none' )  &  strcmp( tra_movement, 'none' )
0152   error( 'Both O.rec_movement and O.tra_movement can not be ''none'',' );
0153 end
0154 
0155 % Get basic calculation options
0156 %
0157 do_atten      = safegetfield( O, 'do_atten', false );
0158 do_absspecies = safegetfield( O, 'do_absspecies', [] );
0159 do_faraday    = safegetfield( O, 'do_faraday', false );
0160 
0161 
0162 
0163 %- Create a temporary workfolder
0164 %
0165 if isempty( workfolder )
0166   workfolder = create_tmpfolder;
0167   cu = onCleanup( @()delete_tmpfolder( workfolder ) );
0168 end
0169 %
0170 cfile = fullfile( workfolder, 'cfile.arts' );
0171 
0172 
0173 %- Set basic and observation related fields of Q
0174 %
0175 Q = q_basic_local( Q, O, workfolder, do_atten, do_faraday );
0176 
0177 
0178 %- Set basic and ppath parts of Q, and get mass of planet
0179 %
0180 [Q,M] = q_planet_local( Q, A, workfolder, do_atten, do_faraday, do_absspecies );
0181 %
0182 if isfield(O,'r_planet') & ~isempty(O.r_planet)
0183   Q.REFELLIPSOID = [ O.r_planet, 0 ]';
0184 end
0185 %
0186 if isfield(O,'n_agenda') & ~isempty(O.n_agenda)
0187   Q.REFR_INDEX_AIR_AGENDA = O.n_agenda;
0188 end
0189 
0190 
0191 %- Set up z_impact grid and R
0192 %
0193 R.z_impact = [ O.z_impact_max : -O.z_impact_dz : ...
0194                              safegetfield( O, 'z_impact_min', O.z_surface ) ]';
0195 n0         = length( R.z_impact );
0196 %
0197 [R.z_tan,R.lat_tan,R.bangle,R.lat_distance]     = deal( repmat( NaN, n0, 1 ) );
0198 [R.slta,R.l_geometric,R.l_optpath,R.l_proppath] = deal( repmat( NaN, n0, 1 ) );
0199 %
0200 if do_atten
0201   [R.tr_defoc,R.tr_space,R.tr_atmos] = deal( repmat( NaN, n0, 1 ) );
0202   if do_faraday
0203     R.faraday                        = deal( repmat( NaN, n0, 1 ) );
0204   end
0205 end
0206 
0207 
0208 %- Radiii
0209 %
0210 r_tra = Q.REFELLIPSOID(1) + O.tra_altitude;
0211 r_rec = Q.REFELLIPSOID(1) + O.rec_altitude;
0212 
0213 
0214 %- Run ARTS
0215 %
0216 [za_tra] = deal( repmat( NaN, n0, 1 ) );
0217 %
0218 for i = 1 : n0
0219 
0220   % Zenith angle to use
0221   %
0222   za_tra(i) = geomztan2za( Q.REFELLIPSOID(1), O.tra_altitude, R.z_impact(i) );
0223   %
0224   Q.SENSOR_LOS = za_tra(i);
0225 
0226   % Run ARTS
0227   %
0228   S = qarts2cfile( Q, { 'Generl', 'AtmSrf', 'AbsrptSave', 'CldBox', ...
0229                         'RteSet', 'CloseF' }, workfolder );
0230   strs2file( cfile, S );  
0231   notok = arts( cfile, true );
0232   %
0233   if notok
0234     error('\n!!! Error while running ARTS !!!\n');
0235   end
0236 
0237   ppath = xmlLoad( fullfile( workfolder, 'ppath.xml' ) );   
0238 
0239   if strcmp( ppath.background, 'space' )  &  ppath.np > 1
0240   
0241     % Extract tangent point
0242     [u,j]        = min( ppath.pos(:,1) );
0243     R.z_tan(i)   = ppath.pos(j,1);
0244     R.lat_tan(i) = ppath.pos(j,2);
0245 
0246     % Bending angle
0247     R.bangle(i) = ppath.start_los(1) - ppath.end_los(1) + ppath.start_pos(2);
0248 
0249     % Calculate angular distance between transmitter and receiver
0250     c                 = ppath.r(end) * sind(ppath.start_los);
0251     za_rec            = asind(c/r_rec); 
0252     dlat              = ppath.start_los - za_rec;
0253     R.lat_distance(i) = ppath.start_pos(2) + dlat;
0254 
0255     % Path lenghts
0256     % Part covered by ppath
0257     lppath = ppath.end_lstep + sum( ppath.lstep );
0258     lopt   = ppath.end_lstep + sum( ppath.lstep .* ( 0.5 * ...
0259                                (ppath.nreal(1:end-1)+ppath.nreal(2:end)) ) );
0260     % Add remaining part
0261     lext   = r_rec * sind(dlat) / sind(ppath.start_los); 
0262     lppath = lppath + lext;
0263     lopt   = lopt   + lext;
0264     
0265     % SLTA and geometrical distance
0266     lgeom        = sqrt( r_tra^2 + r_rec^2 - 2*r_tra*r_rec*...
0267                                                      cosd(R.lat_distance(i)) );
0268     R.slta(i)    = r_tra * r_rec * sind(R.lat_distance(i)) / lgeom - ...
0269                                                              Q.REFELLIPSOID(1);
0270     R.l_proppath(i)  = lppath;
0271     R.l_optpath(i)   = lopt;
0272     R.l_geometric(i) = lgeom;
0273 
0274     if do_atten
0275       
0276       iy     = xmlLoad( fullfile( workfolder, 'iy.xml' ) );   
0277       iy_aux = xmlLoad( fullfile( workfolder, 'iy_aux.xml' ) );   
0278 
0279       R.tr_space(i) = 1 / ( 4 * pi * lppath^2 );
0280       R.tr_atmos(i) = iy(1);            
0281 
0282       if do_faraday
0283         R.faraday(i) = iy_aux{1}; 
0284       end
0285 
0286       % Transmission for each abs_species
0287       if ~isempty(do_absspecies)
0288         nas = length(iy_aux) - do_faraday;
0289         if ~isfield( R, 'tr_absspecies' )
0290           R.tr_absspecies = deal( repmat( NaN, n0, nas ) );
0291         end
0292         %
0293         for j = 1 : nas
0294           ga = squeeze( iy_aux{j+do_faraday}(1,1,1,:) );
0295           R.tr_absspecies(i,j) = exp( -0.5 * sum( ppath.lstep .* ...
0296                                                ( ga(1:end-1) + ga(2:end) ) ) );
0297         end
0298       end
0299       
0300       % Remaining attenuation calculated as post-processing
0301     end
0302   end
0303 end
0304 
0305 
0306 if do_atten
0307   
0308   % Defocusing
0309   %
0310   defoc_shift = safegetfield( O, 'defoc_shift', 3e-3 );
0311   ind         = find( ~isnan( R.bangle ) );
0312   %
0313   for i = 1 : n0
0314     if ~isnan( R.bangle(i) )
0315    
0316       r_tan = Q.REFELLIPSOID(1) + R.z_tan(i);
0317       l1    = sqrt( r_tra^2 - r_tan^2 );
0318       l2    = sqrt( r_rec^2 - r_tan^2 );
0319       fac   = pi/180 * (l1*l2)/(l1+l2); 
0320    
0321       za = za_tra(i) - defoc_shift;
0322       a1 = r_tra * sind( za );
0323       b1 = interp1( za_tra(ind), R.bangle(ind), za, 'linear', 'extrap' );
0324       za = za_tra(i) + defoc_shift;
0325       a2 = r_tra * sind( za );
0326       b2 = interp1( za_tra(ind), R.bangle(ind), za, 'linear', 'extrap' );
0327 
0328       R.tr_defoc(i) = 1 / ( 1 - fac*(b2-b1)/(a2-a1) );
0329     end
0330   end
0331 
0332 
0333   % Total transmission
0334   R.tr_total = R.tr_space .* R.tr_atmos .* R.tr_defoc;
0335 end
0336 
0337 % Speed of increase in angular distance
0338 %
0339 G         = constants( 'GRAVITATIONAL_CONSTANT' );
0340 t_rec     = sqrt( (4*pi^2*r_rec^3) / (G*M) );
0341 %
0342 switch lower(rec_movement)
0343  case 'none'
0344   degpersec = 0;
0345  case 'approaching'
0346   degpersec = - 360 / t_rec;
0347  case 'disappearing'
0348   degpersec = + 360 / t_rec;
0349  otherwise
0350   error( 'Unrecognised choice for *O.rec_movement* (%s).', rec_movement );
0351 end
0352 %
0353 t_tra     = sqrt( (4*pi^2*r_tra^3) / (G*M) );
0354 %
0355 switch lower(tra_movement)
0356  case 'none'
0357   ;
0358  case 'approaching'
0359   degpersec = degpersec - 360 / t_tra;
0360  case 'disappearing'
0361   degpersec = degpersec + 360 / t_tra;
0362  otherwise
0363   error( 'Unrecognised choice for *O.tra_movement* (%s).', tra_movement );
0364 end
0365 
0366 
0367 % Relative arrival time
0368 %
0369 ind = find( ~isnan( R.lat_distance ) );
0370 %
0371 R.t = ( R.lat_distance - interp1( R.z_impact(ind), R.lat_distance(ind), ...
0372                                                 O.z_impact4t0 ) ) / degpersec;
0373 
0374 if nargout == 1
0375   return;   % --->
0376 end
0377 
0378 % Create T
0379 %
0380 % Time grid to use
0381 T.t   = [ min(R.t) : 1/O.f_sampling : max(R.t) ]';
0382 %
0383 % Find start and end of each time chunck
0384 istart = [];
0385 iend   = [];
0386 i1     = NaN; % Start of present chunck%
0387 
0388 for i = 1 : n0
0389   % If Nan, then we are inbetween chuncks
0390   if isnan( R.t(i) )
0391     i1 = NaN;
0392 
0393   else
0394     % If i1==NaN, this is the start of a new chunck
0395     if isnan(i1)
0396       i1 = i;
0397     end
0398     
0399     % End point if last point or next is NaN
0400     if i == n0 | isnan(R.t(i+1))
0401       if i > i1 % Ignore length-1 chuncks
0402         istart = [ istart i1 ];
0403         iend   = [ iend   i  ];
0404       end
0405       i1     = NaN;
0406     
0407     % Also end point if sign change of t-derivate
0408     elseif i>1  &  ~isnan(R.t(i-1))  &  ...
0409           sign(R.t(i+1)-R.t(i)) ~= sign(R.t(i)-R.t(i-1))
0410       istart = [ istart i1 ];
0411       iend   = [ iend   i  ];
0412       i1     = i;      
0413     end    
0414   end
0415 end
0416 
0417 
0418 % Allocate T
0419 %
0420 np    = length( istart );
0421 nt    = length( T.t );
0422 %
0423 rfields = fieldnames( R );
0424 %
0425 for j = 1 : length( rfields )
0426   fname = rfields{j};
0427   if ~strcmp( fname, 't' )
0428     T.(fname) = repmat( NaN, nt, np );
0429   end
0430 end  
0431 
0432 
0433 % Interpolate in T.t
0434 for i = 1 : np
0435 
0436   % Interpolate over time of chunk
0437   ir = istart(i) : iend(i);
0438   tmin = min( [ R.t(ir([1,end])) ] );
0439   tmax = max( [ R.t(ir([1,end])) ] );
0440   it = find( T.t >= tmin  &  T.t < tmax );  
0441 
0442   for j = 1 : length( rfields )
0443     fname = rfields{j};
0444     if ~strcmp( fname, { 't', 'tr_absspecies'} )
0445       T.(fname)(it,i) = interp1( R.t(ir), R.(fname)(ir), T.t(it) );
0446     end
0447   end
0448 end
0449 
0450 
0451 R = orderfields( R );
0452 T = orderfields( T );
0453 
0454 return
0455 
0456 
0457 %----------------------------------------------------------------------------
0458 function Q = q_basic_local(Q,O,workfolder,do_atten,do_faraday)
0459   %
0460   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'general.arts' );
0461   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'agendas.arts' ); 
0462   Q.INCLUDES{end+1}       = fullfile( 'ARTS_INCLUDES', 'continua.arts' );
0463   %
0464   Q.ATMOSPHERE_DIM        = 1;
0465   Q.STOKES_DIM            = 1;   % Can be changed below
0466   Q.F_GRID                = O.frequency;
0467   Q.PPATH_LMAX            = O.lmax;
0468   Q.PPATH_LRAYTRACE       = O.lraytrace;
0469   %
0470   Q.CLOUDBOX_DO           = false;
0471   %
0472   Q.PPATH_AGENDA          = { 'ppath_agenda__FollowSensorLosPath'   };
0473   Q.PPATH_STEP_AGENDA     = { 'ppath_step_agenda__RefractedPath'    };
0474   %
0475   Q.SENSOR_POS            = O.tra_altitude;
0476   %
0477   Q.Z_SURFACE             = O.z_surface;
0478   %
0479   Q.WSMS_AT_END{end+1} = 'VectorExtractFromMatrix(rte_pos,sensor_pos,0,"row")';
0480   Q.WSMS_AT_END{end+1} = 'VectorExtractFromMatrix(rte_los,sensor_los,0,"row")';
0481   Q.WSMS_AT_END{end+1} = 'VectorSet(rte_pos2,[])';
0482   %
0483   if ~do_atten
0484     Q.WSMS_AT_END{end+1} = 'ppathCalc';
0485   else
0486     %
0487     Q.WSMS_AT_END{end+1} = 'jacobianOff';
0488     %
0489     Q.ABSORPTION         = 'OnTheFly';
0490     Q.ABS_LINES_FORMAT   = 'None';
0491     %
0492     if do_faraday
0493       Q.STOKES_DIM              = 3; 
0494       Q.PROPMAT_CLEARSKY_AGENDA = { ...
0495                                  'propmat_clearsky_agenda__OnTheFly_Faraday' };
0496     else
0497       Q.PROPMAT_CLEARSKY_AGENDA = { 'propmat_clearsky_agenda__OnTheFly' };
0498     end
0499     %
0500     Q.IY_TRANSMITTER_AGENDA = { 'iy_transmitter_agenda__UnitUnpolIntensity' };
0501     %
0502     Q.WSMS_AT_END{end+1} = 'AgendaSet(iy_main_agenda){iyTransmissionStandard}';                    
0503     %
0504     Q.WSMS_AT_END{end+1} = 'iyCalc';
0505     Q.WSMS_AT_END{end+1} = ...
0506                       sprintf('WriteXML("ascii",iy,"%s/iy.xml")',workfolder);
0507     Q.WSMS_AT_END{end+1} = ...
0508                 sprintf('WriteXML("ascii",iy_aux,"%s/iy_aux.xml")',workfolder);
0509   end
0510   %
0511   Q.WSMS_AT_END{end+1} = ...
0512       sprintf('ppathWriteXMLPartial("ascii",ppath,"%s/ppath.xml")',workfolder);
0513   %
0514 return
0515 
0516 
0517 %----------------------------------------------------------------------------
0518 function [Q,M] = q_planet_local(Q,A,workfolder,do_atten,do_faraday, ...
0519                                                         do_absspecies)
0520   %
0521   rqre_field( A, 'planet' );
0522   rqre_field( A, 'atmfunc' );
0523   %
0524   if ~isa( A.atmfunc, 'function_handle' )
0525     error( 'The field *A.atmfunc* must be a function handle' ); 
0526   end
0527   if ~exist( func2str(A.atmfunc), 'file' )
0528     error( 'The function *%s* is not found (selected by A.atmfunc).', ...
0529                                                         func2str(A.atmfunc) ); 
0530   end
0531 
0532   % Set basics and default for planet:
0533   %
0534   switch upper( A.planet )
0535 
0536    case 'VENUS'
0537     %
0538     M = 4.8676e24;
0539     %
0540     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_venus.arts' );
0541     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalVenus' );
0542    
0543    case 'EARTH'
0544     %
0545     M = constants( 'EARTH_MASS' );
0546     %
0547     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' );
0548     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalEarth' );
0549     %
0550     if ~isfield( A, 'ABS_SPECIES' )
0551       A.ABS_SPECIES(1).TAG{1} = 'N2-SelfContStandardType';
0552       A.ABS_SPECIES(2).TAG{1} = 'O2-PWR98';
0553       A.ABS_SPECIES(3).TAG{1} = 'H2O-PWR98';
0554       A.ABS_SPECIES(4).TAG{1} = 'free_electrons';
0555     end
0556 
0557    case 'MARS'
0558     %
0559     M = 6.4185e23;
0560     %
0561     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_mars.arts' );
0562     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalMars' );
0563     
0564    case 'JUPITER'
0565     %
0566     M = 1.8986e27;
0567     %
0568     Q.INCLUDES{end+1} = fullfile( 'ARTS_INCLUDES', 'planet_jupiter.arts' );
0569     Q.REFELLIPSOID    = ellipsoidmodels( 'SphericalJupiter' );
0570    
0571    otherwise
0572     error( 'Unrecognised choice for A.planet (%s).', A.planet );
0573   end
0574 
0575   switch upper( A.planet )
0576    case 'EARTH'
0577     %
0578     Q.REFR_INDEX_AIR_AGENDA = { 
0579         'NumericSet( refr_index_air, 1.0 )', ...
0580         'NumericSet( refr_index_air_group, 1.0 )',...
0581         'refr_index_airThayer', ...
0582         'refr_index_airFreeElectrons( demand_vmr_value = 0 )' };
0583     %
0584    otherwise 
0585     %
0586     Q.REFR_INDEX_AIR_AGENDA = { 
0587         'NumericSet( refr_index_air, 1.0 )', ...
0588         'NumericSet( refr_index_air_group, 1.0 )',...
0589         'refr_index_airMWgeneral', ...
0590         'refr_index_airFreeElectrons( demand_vmr_value = 0 )' };
0591     %
0592   end
0593   
0594   % Include atmospheric specific data
0595   Q = A.atmfunc( A, Q, workfolder );
0596   
0597   if 0   % Set to 1 if you want to test without refraction
0598     Q.REFR_INDEX_AIR_AGENDA = { 'refr_index_air_agenda__NoRefrac' };
0599   end
0600 
0601   % Set iy_aux_vars
0602   saux = [];
0603   if do_atten
0604     if do_faraday
0605       if isempty(do_absspecies)
0606         saux = '"Faraday rotation"';
0607       else
0608         saux = '"Faraday rotation",';        
0609       end
0610     end
0611     if ~isempty(do_absspecies)
0612       n = length( do_absspecies );
0613       for i = 1 : n
0614         k = do_absspecies(i) - 1;
0615         if i < n
0616           saux = [ saux, sprintf('"Absorption, species %d",', k ) ];
0617         else
0618           saux = [ saux, sprintf('"Absorption, species %d"', k ) ];
0619         end
0620       end
0621     end
0622     Q.WSMS_AT_START{end+1} = [ 'ArrayOfStringSet(iy_aux_vars,[',saux,'])' ];    
0623   end
0624   
0625 return

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