Home > atmlab > demos > qarts_mcdoit_demo.m

qarts_mcdoit_demo

PURPOSE ^

QARTS_MCDOIT_DEMO

SYNOPSIS ^

function [y0,yd,ym,ymerror] = qarts_mcdoit_demo( varargin )

DESCRIPTION ^

 QARTS_MCDOIT_DEMO

   A demenstration / test of setting up MC and DOIT 

   A satellite measurement is simulated, from 830 km and a zenith angle of
   135 deg.

   The frequency is set by the user. Cloud extends between zcloud +- 1km. 
   Particles are set to be (solid) ice spheres, all with the same 
   size (*dpart*).
 
   The surface reflectivi matrix is set to:
     rsurf*[  1  0.1 0 0 
             0.1  1  0 0   
              0   0  1 0   
              0   0  0 1 ]

 FORMAT [y0,yd,ym,ymerror] = qarts_mcdoit_demo( 
                                    [CMdu,CMmu,iwp,zcloud,freq,dpart,rsurf] )

 OUT   y0           Clear-sky brightness temperature [K]
       yd           DOIT result
       ym           MC result
       ymerror      Error estimate for *ym*
 OPT   CMdu         User DOIT settings for C.METHOD_PRMTRS. To see defaults
                    look into the function, the defaulrs are set a few
                    lines this text. CMd is the variable to check out.
       CMdm         As above but refers to MC.
       iwp          IWP. Default is 0.5 kg/m2
       zcloud       Centre altitude of cloud. Default is 4 km.
       freq         Frequency. Default is 89 GHz
       dpart        Particle diameter. Default is 1 mm
       rsurf        Surface reflectivity. Default is 0.1

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

qarts_mcdoit_demo.m

SOURCE CODE ^

0001 % QARTS_MCDOIT_DEMO
0002 %
0003 %   A demenstration / test of setting up MC and DOIT
0004 %
0005 %   A satellite measurement is simulated, from 830 km and a zenith angle of
0006 %   135 deg.
0007 %
0008 %   The frequency is set by the user. Cloud extends between zcloud +- 1km.
0009 %   Particles are set to be (solid) ice spheres, all with the same
0010 %   size (*dpart*).
0011 %
0012 %   The surface reflectivi matrix is set to:
0013 %     rsurf*[  1  0.1 0 0
0014 %             0.1  1  0 0
0015 %              0   0  1 0
0016 %              0   0  0 1 ]
0017 %
0018 % FORMAT [y0,yd,ym,ymerror] = qarts_mcdoit_demo(
0019 %                                    [CMdu,CMmu,iwp,zcloud,freq,dpart,rsurf] )
0020 %
0021 % OUT   y0           Clear-sky brightness temperature [K]
0022 %       yd           DOIT result
0023 %       ym           MC result
0024 %       ymerror      Error estimate for *ym*
0025 % OPT   CMdu         User DOIT settings for C.METHOD_PRMTRS. To see defaults
0026 %                    look into the function, the defaulrs are set a few
0027 %                    lines this text. CMd is the variable to check out.
0028 %       CMdm         As above but refers to MC.
0029 %       iwp          IWP. Default is 0.5 kg/m2
0030 %       zcloud       Centre altitude of cloud. Default is 4 km.
0031 %       freq         Frequency. Default is 89 GHz
0032 %       dpart        Particle diameter. Default is 1 mm
0033 %       rsurf        Surface reflectivity. Default is 0.1
0034 
0035 function [y0,yd,ym,ymerror] = qarts_mcdoit_demo( varargin )
0036 %
0037 % Defaults for DOIT and MC
0038 CMd.N_ZA_GRID            =  19;
0039 CMd.N_AA_GRID            =  10;
0040 CMd.ZA_GRID_OPT_FILE     = '';
0041 CMd.EPSILON              = [ 0.1 0.01 0.01 0.01 ];
0042 CMd.SCAT_ZA_INTERP       = 'polynomial';
0043 CMd.ALL_F                = true;
0044 CMd.NORMALIZE            = true;
0045 CMd.NORM_ERROR_THRESHOLD = 0.5;
0046 %
0047 CMm.STD_ERR              = 1;
0048 CMm.MAX_TIME             = 20;
0049 CMm.MAX_ITER             = -1;
0050 CMm.MIN_ITER             = 100;
0051 %
0052 [CMdu,CMmu,iwp,zcloud,freq,dpart,rsurf] = optargs( varargin, ...
0053                                        { [], [], 0.5, 4e3, 89e9, 1e-3, 0.1 } );
0054 %
0055 if zcloud<2e3, error( '*zcloud must be >= 2 km.' ); end
0056 
0057 % Move settings in CMdu and CMmu to Cmd and CMm,respectively
0058 %
0059 CMd = optargs_struct( CMdu, CMd );
0060 CMm = optargs_struct( CMmu, CMm );
0061 
0062 
0063 %= Create a temporary workfolder
0064 %
0065 workfolder = create_tmpfolder;
0066 cu = onCleanup( @()delete_tmpfolder( workfolder ) );
0067   
0068   
0069 %= Init Q structure
0070 %
0071 Q  = qarts;
0072 
0073 
0074 %= Overall settings
0075 %
0076 Q.CLOUDBOX_DO       = false;
0077 Q.J_DO              = false;
0078 Q.SENSOR_DO         = false;
0079 %
0080 Q.INCLUDES          = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0081                         fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0082                         fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0083                         fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0084 
0085 
0086 %= Set standard clear-sky agendas
0087 %
0088 Q.WSMS_AT_START{1} = 'Copy(iy_main_agenda,iy_main_agenda__Emission)';
0089 Q.WSMS_AT_START{2} = 'Copy(ppath_agenda,ppath_agenda__FollowSensorLosPath)';
0090 Q.WSMS_AT_START{3} = ...
0091                     'Copy(ppath_step_agenda,ppath_step_agenda__GeometricPath)';
0092 Q.WSMS_AT_START{4} = ...
0093          'Copy(blackbody_radiation_agenda,blackbody_radiation_agenda__Planck)';
0094 Q.WSMS_AT_START{5} = 'Copy(iy_space_agenda,iy_space_agenda__CosmicBackground)';
0095 Q.WSMS_AT_START{6} = ...
0096                  'Copy(iy_surface_agenda,iy_surface_agenda__UseSurfaceRtprop)';
0097 
0098 
0099 %= General part
0100 %
0101 Q.F_GRID            = freq;
0102 Q.STOKES_DIM        = 4;
0103 %
0104 Q.ABS_SPECIES(1).TAG{1} = 'H2O-PWR98';
0105 Q.ABS_SPECIES(2).TAG{1} = 'N2-SelfContStandardType';
0106 Q.ABS_SPECIES(3).TAG{1} = 'O2-PWR93';
0107 
0108 
0109 %= Define atmosphere
0110 %
0111 Q.ATMOSPHERE_DIM      = 1;
0112 %
0113 Q.P_GRID              = z2p_simple( [0:200:10e3 11e3:1e3:25e3] )';
0114 %
0115 arts_xmldata_path     = atmlab( 'ARTS_XMLDATA_PATH' );
0116 if isnan( arts_xmldata_path )
0117   error('You need to ARTS_XMLDATA_PATH to run this example.');
0118 end
0119 %
0120 Q.RAW_ATMOSPHERE        = fullfile( arts_xmldata_path, 'planets', 'Earth', ...
0121                      'Fascod', 'midlatitude-winter', 'midlatitude-winter' );
0122 Q.RAW_ATM_EXPAND_1D   = false;
0123 
0124 
0125 %= Surface
0126 %
0127 Q.REFELLIPSOID        = ellipsoidmodels( 'SphericalEarth' );
0128 Q.Z_SURFACE           = text2cfile( 'Extract( z_surface, z_field, 0 )' );
0129 %
0130 R        = zeros( 1, Q.STOKES_DIM, Q.STOKES_DIM ); 
0131 R(1,:,:) = rsurf*eye( Q.STOKES_DIM );
0132 
0133 R(1,1,2) = rsurf*0.1; 
0134 R(1,2,1) = rsurf*0.1; 
0135 
0136 fname = fullfile(workfolder,'surface_reflectivity.xml');
0137 xmlStore( fname, R, 'Tensor3' );
0138 Q.WSMS_AT_START{end+1} = sprintf('ReadXML(surface_reflectivity,"%s")',fname);
0139 %
0140 Q.SURFACE_RTPROP_AGENDA{1} = 'specular_losCalc';
0141 Q.SURFACE_RTPROP_AGENDA{2} = 'InterpAtmFieldToPosition( out=surface_skin_t';
0142 Q.SURFACE_RTPROP_AGENDA{3} = '  , field=t_field )';
0143 Q.SURFACE_RTPROP_AGENDA{4} = 'surfaceFlatReflectivity';
0144 
0145 
0146 %= Create an absorption table
0147 %
0148 Q.ABS_LINES_FORMAT    = 'none';
0149 %
0150 Q.ABS_NLS             = [];
0151 Q                     = qarts_abstable( Q );
0152 arts_abstable( Q, workfolder );
0153 
0154 
0155 %= Set RTE variables
0156 %
0157 Q.IY_UNIT             = 'RJBT';
0158 Q.YCALC_WSMS          = { 'yCalc' };
0159 %
0160 Q.PPATH_LMAX          = 5e3/max([1,iwp/0.5]);
0161 %
0162 Q.SENSOR_POS          = 830e3;
0163 Q.SENSOR_LOS          = 135;
0164 
0165 
0166 % Calculate clearsky
0167 %
0168 y0 = arts_y( Q, workfolder );
0169 %
0170 if nargout < 2, return, end
0171 
0172 
0173 %= Setting structure for cloudbox and scattering solution method
0174 %
0175 C               = qartsCloudbox;
0176 
0177 %- Define 1D cloudbox and particles
0178 
0179 C.LIMITS = [ -1e3 zcloud+2e3 ];
0180 
0181 % Temperature for scattering data (if only one is given, no
0182 % temperature interpolation is performed in ARTS calculation.
0183 T_grid = [200 235 270];   % MC does not accept a single temperature!
0184 
0185 
0186 % Create scattering properties for a simple cloud assumption
0187 % (mono disperese particle distribution)
0188 
0189 % Calculate refractive indices
0190 fgrid     = Q.F_GRID + [-1e6,1e6]';  
0191 rfr_index = zeros( length(fgrid), length(T_grid) );
0192 %
0193 for i = 1 : length(T_grid)     % Note .' below. We do not want to conjugate
0194   rfr_index(:,i) = sqrt(eps_ice_matzler06(fgrid, T_grid(i) ) ).';
0195 end
0196 
0197 % Scattering angles
0198 theta = 0:2:180;
0199 % Particle size [m]
0200 r = dpart/2;
0201 % cloud altitude
0202 alt = zcloud + [-1e3 1e3];
0203 
0204 % Calculate scattering data using Mie
0205 C.SCAT_DATA{1} = mie_arts_scat_data( fgrid,  T_grid, rfr_index, theta, r );
0206 
0207 % Calculate a pnd field for a homogeneous cloud
0208 C.PND_FIELD{1} = box_pnd_mono_size_1d( alt, iwp/diff(alt), r );
0209 
0210 % Set up scattering method
0211 %
0212 C.METHOD                        = 'DOIT';
0213 C.METHOD_PRMTRS                 = CMd;
0214 
0215 % Need a trick here for setting RJ
0216 Q.IY_UNIT                       = '1';
0217 Q.WSMS_BEFORE_RTE{1}            = 'StringSet(iy_unit,"RJBT")';  
0218 
0219 C.OPT_PROP_PART_AGENDA = { 'ext_matInit', 'abs_vecInit', 'ext_matAddPart', ...
0220                                                          'abs_vecAddPart' };
0221 C.SPT_CALC_AGENDA = { 'opt_prop_sptFromMonoData' };
0222 
0223 
0224 %- Activate cloudbox
0225 %
0226 Q.CLOUDBOX_DO = true;
0227 Q.CLOUDBOX    = C;
0228 
0229 
0230 %= DOIT
0231 %
0232 yd = arts_y( Q, workfolder );
0233 %
0234 if nargout < 3, return, end
0235   
0236 
0237 %= Expand to 3D
0238 %
0239 latlon_grid = [-80:10:80]';
0240 %
0241 Q.ATMOSPHERE_DIM                = 3;
0242 Q.RAW_ATM_EXPAND_1D             = true;
0243 Q.LAT_GRID                      = latlon_grid;
0244 Q.LON_GRID                      = latlon_grid;
0245 %
0246 C.LIMITS                        = [ C.LIMITS [-50 50 -50 50] ];
0247 C.PND_FIELD{1}.grids{2}         = [ -90 -50 -40 40 50 90 ];
0248 C.PND_FIELD{1}.grids{3}         = C.PND_FIELD{1}.grids{2};
0249 p                               = C.PND_FIELD{1}.data;
0250 C.PND_FIELD{1}.data             = zeros( length( C.PND_FIELD{1}.grids{1} ),...
0251                                          length( C.PND_FIELD{1}.grids{2} ),...
0252                                          length( C.PND_FIELD{1}.grids{3} ) );
0253 C.PND_FIELD{1}.data(:,3:4,3:4)  = repmat(p,[1 2 2]);
0254 %
0255 Q.SENSOR_POS(:,2)               = -5;
0256 Q.SENSOR_POS(:,3)               = 1;
0257 Q.SENSOR_LOS(:,2)               = 0;
0258 
0259 
0260 %- Change settings to fit MC
0261 %
0262 C.METHOD                        = 'MC';
0263 C.METHOD_PRMTRS                 = CMm;
0264 %
0265 Q.CLOUDBOX                      = C;
0266 Q.IY_MAIN_AGENDA                = { 'Ignore( rte_pos2 )', 'iyMC', ...
0267                                     'Touch(ppath)' };
0268 Q.IY_AUX_VARS                   = { 'Error (uncorrelated)' };
0269 Q.IY_UNIT                       = 'RJBT';
0270 Q.WSMS_BEFORE_RTE               = {};
0271 
0272   
0273 %= MC
0274 %
0275 [ym,y_aux] = arts_y( Q );
0276 ymerror    = y_aux{1};
0277

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