Home > atmlab > demos > qarts_scattering_demo.m

qarts_scattering_demo

PURPOSE ^

QARTS_SCATTERING_DEMO Demonstration of scattering calculations by Qarts

SYNOPSIS ^

function [Q,f,ztan,y_c,y,dy] = qarts_scattering_demo( varargin )

DESCRIPTION ^

 QARTS_SCATTERING_DEMO  Demonstration of scattering calculations by Qarts

   ARTS has three modules for performing scattering calculations: DOIT, MC,
   and TR (iyTransmissionStandard). The DOIT method is implemented for both 1D
   and 3D, but is only recommended for 1D. MC is only implemented for 3D. TR
   is defined for all atmospheric dimensionalities. This function runs a
   simple scenario with DOIT-1D, MC (3D) or TR (1/2/3D). 

   The 1D scenario used for DOIT is expanded to 3D for MC. The cloud is then
   given a rectangular shape (in lat and lon). The 1D cloud values are applied
   inside a range of [-15,15] degrees, in both latitude and longitude. This
   extended by a linear transition to clear sky conditions over 3 degrees.
   That is, a perfect match between DOIT and MC can not always be expected as
   the "3D cloud" has a smaller extension than the 1D one, this especially if
   a very low tangent altitude is selected (your path can then even be
   totally outside of the cloud box).

   The cloud is expanded in same way for TR, when 2D or 3D is selected.

 FORMAT [Q,f,ztan,y_c,y,dy] = 
                  qarts_scattering_demo([ztan,method,m_arg,iwpfac,do_refr])

  OUT   Q          Qarts setting structure.
        f          Frequency grid.
        ztan       Set of tangent altitudes used.
        y_c        Calcualted clearsky radiances [RJ BT].
        y          Calculated cloudy randiances  [RJ BT].
        dy         Error estimate. Only valid for MC.
 OPT    ztan       Tangent altitudes. Default is [7 10 13] km.
        method     Scattering method to use. Default is 'mc'. Other options
                   are 'doit', 'and 'tr'. (Upper or lower case letters 
                   do not matter).
        m_arg      Method argument. Not used for DOIT. For MC it is the
                   Target calculation accuracy for MC, in K. Default is 5K.
                   For TR this is atmospheric dimensionality to use.
        iwcfac     Scaling factor of the IWC used. Default is 1. Set to 0
                   to obtain "clear-sky" also inside cloudbox.
        do_refr    Flag to run with refraction. Default is false.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

qarts_scattering_demo.m

SOURCE CODE ^

0001 % QARTS_SCATTERING_DEMO  Demonstration of scattering calculations by Qarts
0002 %
0003 %   ARTS has three modules for performing scattering calculations: DOIT, MC,
0004 %   and TR (iyTransmissionStandard). The DOIT method is implemented for both 1D
0005 %   and 3D, but is only recommended for 1D. MC is only implemented for 3D. TR
0006 %   is defined for all atmospheric dimensionalities. This function runs a
0007 %   simple scenario with DOIT-1D, MC (3D) or TR (1/2/3D).
0008 %
0009 %   The 1D scenario used for DOIT is expanded to 3D for MC. The cloud is then
0010 %   given a rectangular shape (in lat and lon). The 1D cloud values are applied
0011 %   inside a range of [-15,15] degrees, in both latitude and longitude. This
0012 %   extended by a linear transition to clear sky conditions over 3 degrees.
0013 %   That is, a perfect match between DOIT and MC can not always be expected as
0014 %   the "3D cloud" has a smaller extension than the 1D one, this especially if
0015 %   a very low tangent altitude is selected (your path can then even be
0016 %   totally outside of the cloud box).
0017 %
0018 %   The cloud is expanded in same way for TR, when 2D or 3D is selected.
0019 %
0020 % FORMAT [Q,f,ztan,y_c,y,dy] =
0021 %                  qarts_scattering_demo([ztan,method,m_arg,iwpfac,do_refr])
0022 %
0023 %  OUT   Q          Qarts setting structure.
0024 %        f          Frequency grid.
0025 %        ztan       Set of tangent altitudes used.
0026 %        y_c        Calcualted clearsky radiances [RJ BT].
0027 %        y          Calculated cloudy randiances  [RJ BT].
0028 %        dy         Error estimate. Only valid for MC.
0029 % OPT    ztan       Tangent altitudes. Default is [7 10 13] km.
0030 %        method     Scattering method to use. Default is 'mc'. Other options
0031 %                   are 'doit', 'and 'tr'. (Upper or lower case letters
0032 %                   do not matter).
0033 %        m_arg      Method argument. Not used for DOIT. For MC it is the
0034 %                   Target calculation accuracy for MC, in K. Default is 5K.
0035 %                   For TR this is atmospheric dimensionality to use.
0036 %        iwcfac     Scaling factor of the IWC used. Default is 1. Set to 0
0037 %                   to obtain "clear-sky" also inside cloudbox.
0038 %        do_refr    Flag to run with refraction. Default is false.
0039 
0040 % 2010-12-01  Extended to also handle TR by Patrick Eriksson.
0041 % 2010-02-12  Extended to also handle FOS by Patrick Eriksson.
0042 % 2007-07-27  Extended to also handle MC by Patrick Eriksson.
0043 % 2005-06-13  Created by Claudia Emde
0044 
0045 
0046 function [Q,f,ztan,y_c,y,dy] = qarts_scattering_demo( varargin )
0047 %
0048 [ztan,method,m_arg,iwcfac,do_refr] = ...
0049                       optargs( varargin, { [7 10 13]*1e3, 'mc', 5, 1, false });
0050 
0051 
0052 %= Init Q structure
0053 %
0054 Q  = qarts;
0055 
0056 
0057 %= Overall settings
0058 %
0059 Q.CLOUDBOX_DO       = false;
0060 Q.J_DO              = false;
0061 Q.SENSOR_DO         = false;
0062 %
0063 Q.INCLUDES          = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0064                         fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0065                         fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0066                         fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0067 
0068 
0069 %= Set standard clear-sky agendas
0070 %
0071 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0072 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0073 Q.IY_SPACE_AGENDA            = { 'iy_space_agenda__CosmicBackground'   };
0074 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0075 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0076 %
0077 if do_refr
0078   Q.PPATH_STEP_AGENDA        = { 'ppath_step_agenda__RefractedPath'    };
0079   Q.REFR_INDEX_AIR_AGENDA    = { 'refr_index_air_agenda__GasThayer'    };
0080 else
0081   Q.PPATH_STEP_AGENDA        = { 'ppath_step_agenda__GeometricPath'    };
0082 end
0083 
0084 %- Blackbody surface
0085 Q.SURFACE_RTPROP_AGENDA      = ...
0086                  { 'surface_rtprop_agenda__Blackbody_SurfTFromt_field' };
0087 
0088 
0089 %= General part
0090 %
0091 f                   = linspace( 501.18e9, 501.58e9, 3 )';
0092 %f                   = linspace( 299e9, 301e9, 2 )';
0093 Q.F_GRID            = f;
0094 Q.STOKES_DIM        = 2;
0095 %
0096 Q.ABS_SPECIES(1).TAG{1} = 'ClO';
0097 Q.ABS_SPECIES(2).TAG{1} = 'O3';
0098 Q.ABS_SPECIES(3).TAG{1} = 'N2O';
0099 Q.ABS_SPECIES(4).TAG{1} = 'H2O-*-490e9-510e9';  % Some local lines not in PWR98
0100 Q.ABS_SPECIES(4).TAG{2} = 'H2O-PWR98';
0101 Q.ABS_SPECIES(5).TAG{1} = 'N2-SelfContStandardType';
0102 
0103 
0104 %= Define atmosphere and surface
0105 %
0106 Q.ATMOSPHERE_DIM      = 1;
0107 %
0108 Q.P_GRID              = z2p_simple( [0:500:45e3 46e3:1e3:80e3] )';
0109 %
0110 arts_xmldata_path     = atmlab( 'ARTS_XMLDATA_PATH' );
0111 if isnan( arts_xmldata_path )
0112   error('You need to ARTS_XMLDATA_PATH to run this example.');
0113 end
0114 %
0115 Q.RAW_ATMOSPHERE      = fullfile( arts_xmldata_path, 'planets', 'Earth', ...
0116                                             'Fascod', 'tropical', 'tropical' );
0117 Q.RAW_ATM_EXPAND_1D   = false;
0118 %
0119 Q.REFELLIPSOID        = ellipsoidmodels( 'SphericalEarth' );
0120 Q.Z_SURFACE           = 500;
0121 
0122 
0123 %= Absorption
0124 %
0125 Q.ABS_LINES_FORMAT    = 'Arts';
0126 Q.ABS_LINES           = fullfile( atmlab_example_data , 'lines501.4' );
0127 Q.ABS_NLS             = [];
0128 %
0129 Q                     = qarts_abstable( Q );
0130 
0131 
0132 %= Set RTE variables
0133 %
0134 Q.YCALC_WSMS          = { 'sensor_checkedCalc','yCalc' };
0135 %
0136 Q.PPATH_LMAX          = 5e3/max([1,6*iwcfac]);
0137 Q.PPATH_LRAYTRACE     = Q.PPATH_LMAX;
0138 %
0139 zplat                 = 600e3;
0140 sensor_pos            = zplat;
0141 Q.SENSOR_POS          = repmat( sensor_pos, length(ztan), 1 );
0142 Q.SENSOR_LOS          = geomztan2za( Q.REFELLIPSOID(1), zplat, ztan )';
0143 %
0144 if strcmp( lower(method), 'tr' )
0145   Q.IY_MAIN_AGENDA        = { 'iyTransmissionStandard' };
0146   Q.IY_TRANSMITTER_AGENDA = { 'Ignore(rtp_pos)', 'Ignore(rtp_los)', ...
0147                               'MatrixUnitIntensity(iy,stokes_dim,f_grid)' };
0148   Q.IY_UNIT               = '1';
0149 else
0150   Q.IY_UNIT               = 'RJBT';
0151 end
0152   
0153 
0154 
0155 % Calculate clearsky
0156 %
0157 y_c = arts_y( Q );
0158 
0159 
0160 
0161 %= Setting structure for cloudbox and scattering solution method
0162 %
0163 C                      = qartsCloudbox;
0164 Q.WSMS_AT_START{end+1} = 'FlagOff( use_mean_scat_data )';
0165 
0166 
0167 %- Define 1D cloudbox and particles
0168 
0169 C.LIMITS = [ 6e3 16e3 ];
0170 
0171 % Temperature for scattering data (if only one is given, no
0172 % temperature interpolation is performed in ARTS calculation.
0173 T_grid = [200 270];   % MC does not accept a single temperature!
0174 
0175 
0176 % Create scattering properties for a simple cloud assumption
0177 % (mono disperese particle distribution)
0178 
0179 % Calculate refractive indices
0180 rfr_index = zeros( length(Q.F_GRID), length(T_grid) );
0181 %
0182 for i = 1 : length(T_grid)     % Note .' below. We do not want to conjugate
0183   rfr_index(:,i) = sqrt(eps_ice_matzler06( Q.F_GRID, T_grid(i) ) ).';
0184 end
0185 
0186 
0187 % Scattering angles
0188 theta = 0:10:180;
0189 % Particle size [m]
0190 r = 50e-6;
0191 % ice mass content
0192 imc = iwcfac * 0.005*1e-3;
0193 % cloud altitude
0194 alt = [11e3 13e3];
0195 
0196 % Calculate scattering data using Mie
0197 C.SCAT_DATA{1} = mie_arts_scat_data( Q.F_GRID,  T_grid, rfr_index, theta, r );
0198 
0199 % Calculate a pnd field for a homogeneous cloud
0200 C.PND_FIELD{1} = box_pnd_mono_size_1d( alt, imc, r );
0201 
0202 
0203 % Set up scattering method
0204 %
0205 if strcmp( lower(method), 'mc' )
0206   %
0207   C.METHOD                        = 'MC';
0208   Q.IY_MAIN_AGENDA                = { 'Ignore( rte_pos2 )', 'iyMC', ...
0209                                       'Touch(ppath)' };
0210   Q.IY_AUX_VARS                   = { 'Error (uncorrelated)' };
0211   %
0212   C.METHOD_PRMTRS.STD_ERR         = m_arg;
0213   C.METHOD_PRMTRS.MAX_TIME        = -1;
0214   C.METHOD_PRMTRS.MAX_ITER        = -1;
0215   C.METHOD_PRMTRS.MIN_ITER        = 100;
0216   %
0217 elseif strcmp( lower(method), 'doit' )
0218   %
0219   C.METHOD                        = 'DOIT';
0220   % Need a trick here for setting RJ
0221   Q.IY_UNIT                       = '1';
0222   Q.WSMS_BEFORE_RTE{end+1}        = 'StringSet(iy_unit,"RJBT")';  
0223   % Angular grids
0224   C.METHOD_PRMTRS.N_ZA_GRID       =  19;
0225   C.METHOD_PRMTRS.N_AA_GRID       =  10;
0226   C.METHOD_PRMTRS.ZA_GRID_OPT_FILE = fullfile( atmlab_example_data, ...
0227                                                            'doit_za_grid.xml');
0228   C.METHOD_PRMTRS.EPSILON         = [ 0.1 0.01 0.01 0.01 ];
0229   C.METHOD_PRMTRS.SCAT_ZA_INTERP  = 'polynomial';
0230   C.METHOD_PRMTRS.ALL_F           = false;
0231   %
0232   C.OPT_PROP_PART_AGENDA = { 'ext_matInit', 'abs_vecInit', 'ext_matAddPart', ...
0233                                                            'abs_vecAddPart' };
0234   C.SPT_CALC_AGENDA = { 'opt_prop_sptFromMonoData' };
0235 
0236 %elseif strcmp( lower(method), 'fos' )
0237 %  %
0238 %  C.METHOD                        = 'FOS';
0239 %  Q.IY_MAIN_AGENDA                = { 'iyFOS(fos_n=0)' };
0240 
0241 elseif strcmp( lower(method), 'tr' )
0242   %
0243   C.METHOD = 'none';
0244 
0245 else
0246   error( ...
0247     'Allowed options for *method* are ''doit'', ''mc'' and ''tr''.' );
0248 end
0249 
0250 
0251 %- Map to 2D or 3D?
0252 %
0253 latlon_grid = [-45 -20:5:20 45]';
0254 csize       = 10;
0255 %
0256 if any( strcmp( lower(method), {'fos','tr'} ) )  &  m_arg == 2
0257   %
0258   Q.ATMOSPHERE_DIM                = 2;
0259   Q.RAW_ATM_EXPAND_1D             = true;
0260   Q.LAT_GRID                      = latlon_grid;
0261   Q.Z_SURFACE                     = repmat( Q.Z_SURFACE, length(latlon_grid),1);
0262   %
0263   Q.SENSOR_POS                    = repmat( Q.SENSOR_POS, 1, 2 );
0264   Q.SENSOR_POS(:,2)               = -23;
0265   %
0266   C.LIMITS                        = [ C.LIMITS csize*[-2 2] ];
0267   C.PND_FIELD{1}.grids{2}         = [ -90 csize*[-2 -1.5 1.5 2] 90 ];
0268   p                               = C.PND_FIELD{1}.data;
0269   C.PND_FIELD{1}.data             = zeros( length( C.PND_FIELD{1}.grids{1} ),...
0270                                            length( C.PND_FIELD{1}.grids{2} ) );
0271   C.PND_FIELD{1}.data(:,3:4)      = repmat(p,[1 2]);
0272   %
0273 elseif strcmp( lower(method), 'mc' )  |  ...
0274      ( any( strcmp( lower(method), {'fos','tr'} ) )  &  m_arg == 3 )  
0275   %
0276   Q.ATMOSPHERE_DIM                = 3;
0277   Q.RAW_ATM_EXPAND_1D             = true;
0278   Q.LAT_GRID                      = latlon_grid;
0279   Q.LON_GRID                      = latlon_grid;
0280   Q.Z_SURFACE                     = repmat( Q.Z_SURFACE, length(latlon_grid),...
0281                                                          length(latlon_grid) );
0282   %
0283   Q.SENSOR_POS                    = repmat( Q.SENSOR_POS, 1, 3 );
0284   Q.SENSOR_POS(:,2)               = -20;
0285   Q.SENSOR_POS(:,3)               = 0;
0286   Q.SENSOR_LOS                    = repmat( Q.SENSOR_LOS, 1, 2 );
0287   Q.SENSOR_LOS(:,2)               = 0;
0288   %
0289   C.LIMITS                        = [ C.LIMITS csize*[-2 2 -2 2] ];
0290   C.PND_FIELD{1}.grids{2}         = [ -90 csize*[-2 -1.5 1.5 2] 90 ];
0291   C.PND_FIELD{1}.grids{3}         = C.PND_FIELD{1}.grids{2};
0292   p                               = C.PND_FIELD{1}.data;
0293   C.PND_FIELD{1}.data             = zeros( length( C.PND_FIELD{1}.grids{1} ),...
0294                                            length( C.PND_FIELD{1}.grids{2} ),...
0295                                            length( C.PND_FIELD{1}.grids{3} ) );
0296   C.PND_FIELD{1}.data(:,3:4,3:4)  = repmat(p,[1 2 2]);
0297   %
0298 end
0299 
0300 
0301 
0302 %- Activate cloudbox
0303 %
0304 Q.CLOUDBOX_DO = true;
0305 Q.CLOUDBOX    = C;
0306 
0307 
0308 
0309 % Calculate radiances with scattering
0310 %
0311 if strcmp( lower(method), 'mc' )
0312   [y,y_aux] = arts_y( Q );
0313   dy = y_aux{1};
0314 else
0315   y  = arts_y( Q );
0316   dy = NaN;
0317 end
0318 
0319 return
0320 
0321 
0322

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