ARTS  2.2.66
mc_NotUsed.cc
Go to the documentation of this file.
1  wsv_data.push_back
2  (WsvRecord
3  ( NAME( "mc_z_field_is_1D" ),
5  (
6  "Flag for MCGeneral and associated methods.\n"
7  "\n"
8  "If fields outside the cloud box are 1D, this flag can be set to 1\n"
9  "and the calculations will be more rapid.\n"
10  "\n"
11  "Usage: Set by the user.\n"
12  ),
13  GROUP( "Index" )));
14 
15 
16 
17  md_data_raw.push_back
18  ( MdRecord
19  ( NAME( "mc_IWP_cloud_opt_pathCalc" ),
21  (
22  "Calculates the FOV averaged ice water path and cloud optical path\n"
23  "for a given viewing direction\n"
24  ),
25  AUTHORS( "Cory Davis" ),
26  OUT( "mc_IWP", "mc_cloud_opt_path", "mc_IWP_error",
27  "mc_cloud_opt_path_error", "mc_iteration_count" ),
28  GOUT(),
29  GOUT_TYPE(),
30  GOUT_DESC(),
31  IN( "mc_antenna", "sensor_pos", "sensor_los", "ppath_step_agenda",
32  "p_grid", "lat_grid", "lon_grid", "refellipsoid", "z_surface",
33  "z_field", "t_field", "vmr_field", "edensity_field",
34  "f_grid", "f_index", "cloudbox_limits", "pnd_field",
35  "scat_data_array_mono", "particle_masses", "mc_seed" ),
36  GIN( "max_iter" ),
37  GIN_TYPE( "Index" ),
38  GIN_DEFAULT( NODEF ),
39  GIN_DESC( "Maximum number of iterations." )
40  ));
41 
42 /* Workspace method: Doxygen documentation will be auto-generated */
44  //output
45  Numeric& mc_IWP,
46  Numeric& mc_cloud_opt_path,
47  Numeric& mc_IWP_error,
48  Numeric& mc_cloud_opt_path_error,
49  Index& mc_iteration_count,
50  //input
51  const MCAntenna& mc_antenna,
52  const Matrix& sensor_pos,
53  const Matrix& sensor_los,
54  const Agenda& ppath_step_agenda,
55  const Vector& p_grid,
56  const Vector& lat_grid,
57  const Vector& lon_grid,
58  const Vector& refellipsoid,
59  const Matrix& z_surface,
60  const Tensor3& z_field,
61  const Tensor3& t_field,
62  const Tensor4& vmr_field,
63  const Tensor3& edensity_field,
64  const Vector& f_grid,
65  const Index& f_index,
66  const ArrayOfIndex& cloudbox_limits,
67  const Tensor4& pnd_field,
68  const ArrayOfSingleScatteringData& scat_data_array_mono,
69  const Matrix& particle_masses,
70  const Index& mc_seed,
71  //Keyword params
72  const Index& max_iter,
73  const Verbosity& verbosity)
74 {
75  const Numeric f_mono = f_grid[f_index];
76 
77  if( particle_masses.ncols() != 1 )
78  throw runtime_error(
79  "*particle_masses* is only allowed to have a single column" );
80 
81  //if antenna is pencil beam just need a single path integral, otherwise
82  //for now do monte carlo integration
83  if (mc_antenna.get_type()==ANTENNA_TYPE_PENCIL_BEAM)
84  {
85  iwp_cloud_opt_pathCalc(ws, mc_IWP,mc_cloud_opt_path,sensor_pos(0,joker),
86  sensor_los(0,joker),
87  ppath_step_agenda,p_grid,lat_grid,lon_grid,
88  refellipsoid,z_surface,z_field,t_field,vmr_field,
89  edensity_field, f_mono, cloudbox_limits,
90  pnd_field,scat_data_array_mono,particle_masses,
91  verbosity);
92  }
93  else
94  {
95  Numeric iwp,cloud_opt_path;
96  Numeric iwp_squared=0;
97  Numeric tau_squared=0;
98  mc_iteration_count=0;
99  mc_IWP=0;
100  mc_cloud_opt_path=0;
101  Vector local_rte_los(2);
102  Rng rng;
103  rng.seed(mc_seed, verbosity);
104  //Begin Main Loop
105  while (mc_iteration_count<max_iter)
106  {
107  mc_antenna.draw_los(local_rte_los,rng,sensor_los(0,joker));
108  mc_iteration_count+=1;
109  iwp_cloud_opt_pathCalc(ws, iwp,cloud_opt_path,sensor_pos(0,joker),
110  local_rte_los,
111  ppath_step_agenda,p_grid,lat_grid,lon_grid,
112  refellipsoid,z_surface,z_field,t_field,
113  vmr_field,edensity_field, f_mono,
114  cloudbox_limits,
115  pnd_field,scat_data_array_mono,particle_masses,
116  verbosity);
117  mc_IWP+=iwp;
118  iwp_squared+=iwp*iwp;
119  mc_cloud_opt_path+=cloud_opt_path;
120  tau_squared+=cloud_opt_path*cloud_opt_path;
121  }
122  mc_IWP/=(Numeric)mc_iteration_count;
123  mc_cloud_opt_path/=(Numeric)mc_iteration_count;
124  mc_IWP_error=sqrt((iwp_squared/(Numeric)mc_iteration_count-mc_IWP*mc_IWP)
125  /(Numeric)mc_iteration_count);
126  mc_cloud_opt_path_error=sqrt((tau_squared/(Numeric)mc_iteration_count-
127  mc_cloud_opt_path*mc_cloud_opt_path)
128  /(Numeric)mc_iteration_count);
129 
130  }
131 
132 }
133 
134 
135 
136 
137 
138 // The version that was used up to 2012-11-20:
140 
151  MatrixView evol_op,
152  Vector& abs_vec_mono,
153  Numeric& temperature,
154  MatrixView ext_mat_mono,
155  Rng& rng,
156  Vector& rte_pos,
157  Vector& rte_los,
158  Vector& pnd_vec,
159  Numeric& g,
160  Ppath& ppath_step,
161  Index& termination_flag,
162  bool& inside_cloud,
163  const Agenda& ppath_step_agenda,
164  const Agenda& propmat_clearsky_agenda,
165  const Index stokes_dim,
166  const Numeric& f_mono,
167  const Vector& p_grid,
168  const Vector& lat_grid,
169  const Vector& lon_grid,
170  const Tensor3& z_field,
171  const Vector& refellipsoid,
172  const Matrix& z_surface,
173  const Tensor3& t_field,
174  const Tensor4& vmr_field,
175  const Tensor3& edensity_field,
176  const ArrayOfIndex& cloudbox_limits,
177  const Tensor4& pnd_field,
178  const ArrayOfSingleScatteringData& scat_data_array_mono,
179  const Verbosity& verbosity)
180  // 2011-06-17 GH commented out, unused?
181  // const Index z_field_is_1D)
182 
183 {
184  ArrayOfMatrix evol_opArray(2);
185  ArrayOfMatrix ext_matArray(2);
186  ArrayOfVector abs_vecArray(2),pnd_vecArray(2);
187  Vector tArray(2);
188  Vector cum_l_step;
189  Matrix T(stokes_dim,stokes_dim);
190  Numeric k;
191  Numeric ds;
192  Index np=0;
193  Index istep = 0; // Counter for number of steps
194  Matrix opt_depth_mat(stokes_dim,stokes_dim),incT(stokes_dim,stokes_dim,0.0);
195  Matrix old_evol_op(stokes_dim,stokes_dim);
196  //g=0;k=0;ds=0;
197  //at the start of the path the evolution operator is the identity matrix
198  id_mat(evol_op);
199  evol_opArray[1]=evol_op;
200  //initialise Ppath with ppath_start_stepping
201  //Index rubbish=z_field_is_1D;rubbish+=1; // 2011-06-17 GH commented out, unused?
202  ppath_start_stepping( ppath_step, 3, p_grid, lat_grid,
203  lon_grid, z_field, refellipsoid, z_surface,
204  0, cloudbox_limits, false,
205  rte_pos, rte_los, verbosity );
206  //Use cloudbox_ppath_start_stepping to avoid unnecessary z_at_latlon calls.
207  //cloudbox_ppath_start_stepping( ppath_step, 3, p_grid, lat_grid,
208  // lon_grid, z_field, refellipsoid, z_surface, rte_pos,
209  // rte_los ,z_field_is_1D);
210  Range p_range(cloudbox_limits[0],
211  cloudbox_limits[1]-cloudbox_limits[0]+1);
212  Range lat_range(cloudbox_limits[2],
213  cloudbox_limits[3]-cloudbox_limits[2]+1);
214  Range lon_range(cloudbox_limits[4],
215  cloudbox_limits[5]-cloudbox_limits[4]+1);
216  termination_flag=0;
217 
218  inside_cloud=is_inside_cloudbox( ppath_step, cloudbox_limits,true );
219 
220  if (inside_cloud)
221  {
222  cloudy_rt_vars_at_gp(ws,ext_mat_mono,abs_vec_mono,pnd_vec,temperature,
223  propmat_clearsky_agenda,
224  stokes_dim, f_mono, ppath_step.gp_p[0], ppath_step.gp_lat[0],
225  ppath_step.gp_lon[0],p_grid[p_range],
226  t_field(p_range,lat_range,lon_range),
227  vmr_field(joker,p_range,lat_range,lon_range),pnd_field,
228  scat_data_array_mono, cloudbox_limits,ppath_step.los(0,joker),
229  verbosity);
230  }
231  else
232  {
233  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
234  propmat_clearsky_agenda, f_mono,
235  ppath_step.gp_p[0], ppath_step.gp_lat[0], ppath_step.gp_lon[0],
236  p_grid, t_field, vmr_field );
237  pnd_vec=0.0;
238  }
239  ext_matArray[1]=ext_mat_mono;
240  abs_vecArray[1]=abs_vec_mono;
241  tArray[1]=temperature;
242  pnd_vecArray[1]=pnd_vec;
243  //draw random number to determine end point
244  Numeric r = rng.draw();
245  while ((evol_op(0,0)>r)&(!termination_flag))
246  {
247  istep++;
248  //we consider more than 5000
249  // path points to be an indication on that the calcululations have
250  // got stuck in an infinite loop.
251  if( istep > 5000 )
252  {
253  throw runtime_error(
254  "5000 path points have been reached. Is this an infinite loop?" );
255  }
256  evol_opArray[0]=evol_opArray[1];
257  ext_matArray[0]=ext_matArray[1];
258  abs_vecArray[0]=abs_vecArray[1];
259  tArray[0]=tArray[1];
260  pnd_vecArray[0]=pnd_vecArray[1];
261  //perform single path step using ppath_step_geom_3d
262  ppath_step_geom_3d(ppath_step, lat_grid, lon_grid, z_field,
263  refellipsoid, z_surface, -1);
264 
265  // For debugging:
266  // Print( ppath_step, 0, verbosity );
267 
268  np=ppath_step.np;//I think this should always be 2.
269  cum_l_stepCalc(cum_l_step, ppath_step);
270  //path_step should now have two elements.
271  //calculate evol_op
272  inside_cloud=is_inside_cloudbox( ppath_step, cloudbox_limits, true );
273  if (inside_cloud)
274  {
275  cloudy_rt_vars_at_gp(ws,ext_mat_mono,abs_vec_mono,pnd_vec,temperature,
276  propmat_clearsky_agenda, stokes_dim, f_mono,
277  ppath_step.gp_p[np-1],ppath_step.gp_lat[np-1],
278  ppath_step.gp_lon[np-1],p_grid[p_range],
279  t_field(p_range,lat_range,lon_range),
280  vmr_field(joker,p_range,lat_range,lon_range),pnd_field,
281  scat_data_array_mono, cloudbox_limits,ppath_step.los(np-1,joker),
282  verbosity);
283  }
284  else
285  {
286  clear_rt_vars_at_gp(ws, ext_mat_mono,abs_vec_mono,temperature,
287  propmat_clearsky_agenda, f_mono,
288  ppath_step.gp_p[np-1], ppath_step.gp_lat[np-1],
289  ppath_step.gp_lon[np-1], p_grid, t_field, vmr_field);
290  pnd_vec=0.0;
291  }
292  ext_matArray[1]=ext_mat_mono;
293  abs_vecArray[1]=abs_vec_mono;
294  tArray[1]=temperature;
295  pnd_vecArray[1]=pnd_vec;
296  opt_depth_mat=ext_matArray[1];
297  opt_depth_mat+=ext_matArray[0];
298  opt_depth_mat*=-cum_l_step[np-1]/2;
299  incT=0;
300  if ( stokes_dim == 1 )
301  {
302  incT(0,0)=exp(opt_depth_mat(0,0));
303  }
304  else if ( is_diagonal( opt_depth_mat))
305  {
306  for ( Index j=0;j<stokes_dim;j++)
307  {
308  incT(j,j)=exp(opt_depth_mat(j,j));
309  }
310  }
311  else
312  {
313  //matrix_exp(incT,opt_depth_mat,10);
314  matrix_exp_p30(incT,opt_depth_mat);
315  }
316  mult(evol_op,evol_opArray[0],incT);
317  evol_opArray[1]=evol_op;
318 
319  //check whether hit ground or space
320  Index bg = ppath_what_background(ppath_step);
321  if ( bg==2 )
322  {
323  //we have hit the surface
324  termination_flag=2;
325 
326  }
327  //else if ( is_gridpos_at_index_i( ppath_step.gp_p[np-1], p_grid.nelem() - 1 ) )
328  else if (fractional_gp(ppath_step.gp_p[np-1])>=(Numeric)(p_grid.nelem() - 1)-1e-3)
329  {
330  //we have left the top of the atmosphere
331  termination_flag=1;
332 
333  }
334 
335 
336  }
337  if (termination_flag)
338  {//we must have reached the cloudbox boundary
339  //
340  rte_pos = ppath_step.pos(np-1,Range(0,3));
341  rte_los = ppath_step.los(np-1,joker);
342  g=evol_op(0,0);
343  }
344  else
345  {
346  //find position...and evol_op..and everything else required at the new
347  //scattering/emission point
348  // GH 2011-09-14:
349  // log(incT(0,0)) = log(exp(opt_depth_mat(0, 0))) = opt_depth_mat(0, 0)
350  // Avoid loss of precision, use opt_depth_mat directly
351  //k=-log(incT(0,0))/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
352  k=-opt_depth_mat(0,0)/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
353  ds=log(evol_opArray[0](0,0)/r)/k;
354  g=k*r;
355  Vector x(2,0.0);
356  //interpolate atmospheric variables required later
357  //length 2
358  ArrayOfGridPos gp(1);
359  x[1]=cum_l_step[np-1];
360  Vector itw(2);
361 
362  gridpos(gp, x, ds);
363  assert(gp[0].idx==0);
364  interpweights(itw,gp[0]);
365  interp(ext_mat_mono, itw, ext_matArray, gp[0]);
366  opt_depth_mat = ext_mat_mono;
367  opt_depth_mat+=ext_matArray[gp[0].idx];
368  opt_depth_mat*=-ds/2;
369  if ( stokes_dim == 1 )
370  {
371  incT(0,0)=exp(opt_depth_mat(0,0));
372  }
373  else if ( is_diagonal( opt_depth_mat))
374  {
375  for ( Index i=0;i<stokes_dim;i++)
376  {
377  incT(i,i)=exp(opt_depth_mat(i,i));
378  }
379  }
380  else
381  {
382  //matrix_exp(incT,opt_depth_mat,10);
383  matrix_exp_p30(incT,opt_depth_mat);
384  }
385  mult(evol_op,evol_opArray[gp[0].idx],incT);
386  interp(abs_vec_mono, itw, abs_vecArray,gp[0]);
387  temperature=interp(itw,tArray,gp[0]);
388  interp(pnd_vec, itw,pnd_vecArray,gp[0]);
389  if (np > 2)
390  {
391  gridpos(gp,cum_l_step,ds);
392  interpweights(itw,gp[0]);
393  }
394  for (Index i=0; i<2; i++)
395  {
396  rte_pos[i] = interp(itw,ppath_step.pos(Range(joker),i),gp[0]);
397  rte_los[i] = interp(itw,ppath_step.los(Range(joker),i),gp[0]);
398  }
399  rte_pos[2] = interp(itw,ppath_step.pos(Range(joker),2),gp[0]);
400  }
401  assert(isfinite(g));
402 }
403 
404 
405 
406  md_data_raw.push_back
407  ( MdRecord
408  ( NAME( "MCIPA" ),
410  ( "A specialised 3D reversed Monte Carlo radiative algorithm, that\n"
411  "mimics independent pixel appoximation simulations.\n"
412  ),
413  AUTHORS( "Cory Davis" ),
414  OUT( "y", "mc_iteration_count", "mc_error", "mc_points" ),
415  GOUT(),
416  GOUT_TYPE(),
417  GOUT_DESC(),
418  IN( "mc_antenna", "f_grid", "f_index", "sensor_pos", "sensor_los",
419  "stokes_dim", "atmosphere_dim", "iy_space_agenda",
420  "surface_rtprop_agenda",
421  "propmat_clearsky_agenda", "ppath_step_agenda", "p_grid", "lat_grid",
422  "lon_grid", "z_field", "refellipsoid", "z_surface", "t_field",
423  "vmr_field", "edensity_field", "cloudbox_limits", "pnd_field",
424  "scat_data_array_mono", "mc_seed", "iy_unit",
425  "mc_std_err", "mc_max_time", "mc_max_iter", "mc_min_iter",
426  "mc_z_field_is_1D" ),
427  GIN(),
428  GIN_TYPE(),
429  GIN_DEFAULT(),
430  GIN_DESC()
431  ));
432 /* Workspace method: Doxygen documentation will be auto-generated */
433 void MCIPA(Workspace& ws,
434  Vector& y,
435  Index& mc_iteration_count,
436  Vector& mc_error,
437  Tensor3& mc_points,
438  const MCAntenna& mc_antenna,
439  const Vector& f_grid,
440  const Index& f_index,
441  const Matrix& sensor_pos,
442  const Matrix& sensor_los,
443  const Index& stokes_dim,
444  const Index& atmosphere_dim,
445  const Agenda& iy_space_agenda,
446  const Agenda& surface_rtprop_agenda,
447  const Agenda& propmat_clearsky_agenda,
448  const Agenda& ppath_step_agenda,
449  const Vector& p_grid,
450  const Vector& lat_grid,
451  const Vector& lon_grid,
452  const Tensor3& z_field,
453  const Vector& refellipsoid,
454  const Matrix& z_surface,
455  const Tensor3& t_field,
456  const Tensor4& vmr_field,
457  const Tensor3& edensity_field,
458  const ArrayOfIndex& cloudbox_limits,
459  const Tensor4& pnd_field,
460  const ArrayOfSingleScatteringData& scat_data_array_mono,
461  const Index& mc_seed,
462  const String& y_unit,
463  //Keyword params
464  const Numeric& std_err,
465  const Index& max_time,
466  const Index& max_iter,
467  const Index& min_iter,
468  const Index& z_field_is_1D,
469  const Verbosity& verbosity)
470 {
471  if( min_iter < 100 )
472  { throw runtime_error( "*mc_min_iter* must be >= 100." ); }
473 
474  //Check keyword input
475  if (max_time<0 && max_iter<0 && std_err<0){
476  throw runtime_error( "At least one of std_err, max_time, and max_iter must be positive" );
477  }
478 
479  if (! (atmosphere_dim == 3))
480  {
481  throw runtime_error(
482  "For montecarlo, atmosphere_dim must be 3.");
483  }
484 
485  Ppath ppath_step;
486  Rng rng; //Random Nuimber generator
487  time_t start_time=time(NULL);
488  Index N_pt=pnd_field.nbooks();//Number of particle types
489  Vector pnd_vec(N_pt); //Vector of particle number densities used at each point
490  Vector Z11maxvector;//Vector holding the maximum phase function for each
491  bool anyptype30=is_anyptype30(scat_data_array_mono);
492  if (anyptype30)
493  {
494  findZ11max(Z11maxvector,scat_data_array_mono);
495  }
496  rng.seed(mc_seed, verbosity);
497  bool keepgoing,inside_cloud; // flag indicating whether to stop tracing a photons path
498  Numeric g,temperature,albedo,g_los_csc_theta;
499  Matrix A(stokes_dim,stokes_dim),Q(stokes_dim,stokes_dim),evol_op(stokes_dim,stokes_dim),
500  ext_mat_mono(stokes_dim,stokes_dim),q(stokes_dim,stokes_dim),newQ(stokes_dim,stokes_dim),
501  Z(stokes_dim,stokes_dim);
502  q=0.0;newQ=0.0;
503  mc_iteration_count=0;
504  Vector vector1(stokes_dim),abs_vec_mono(stokes_dim),I_i(stokes_dim),Isum(stokes_dim),
505  Isquaredsum(stokes_dim);
506  y.resize(stokes_dim);
507  y=0;
508  Index termination_flag=0;
509  mc_error.resize(stokes_dim);
510  //local versions of workspace
511  Matrix local_iy(1,stokes_dim),local_surface_emission(1,stokes_dim),local_surface_los;
512  Tensor4 local_surface_rmatrix;
513  Vector local_rte_pos(2);
514  Vector local_rte_los(2);
515  Vector new_rte_los(2);
516  Index np;
517  mc_points.resize(p_grid.nelem(),lat_grid.nelem(),lon_grid.nelem());
518  mc_points=0;
519  Isum=0.0;Isquaredsum=0.0;
520  const Numeric f_mono = f_grid[f_index];
521  Numeric std_err_i;
522  bool convert_to_rjbt=false;
523  if ( y_unit == "RJBT" )
524  {
525  std_err_i=f_grid[0]*f_grid[0]*2*BOLTZMAN_CONST/SPEED_OF_LIGHT/SPEED_OF_LIGHT*
526  std_err;
527  convert_to_rjbt=true;
528  }
529  else if ( y_unit == "1" )
530  {
531  std_err_i=std_err;
532  }
533  else
534  {
535  ostringstream os;
536  os << "Invalid value for *y_unit*:" << y_unit <<".\n"
537  << "This method allows only the options \"RJBT\" and \"1\".";
538  throw runtime_error( os.str() );
539  }
540 
541  //Begin Main Loop
542  while (true)
543  {
544  mc_iteration_count+=1;
545  keepgoing=true; //flag indicating whether to continue tracing a photon path
546  //Sample a FOV direction
547  mc_antenna.draw_los(local_rte_los,rng,sensor_los(0,joker));
548  id_mat(Q);
549  local_rte_pos=sensor_pos(0,joker);
550  I_i=0.0;
551 
552  //Obtain a reference propagation path to use for obtaining optical properties
553  //for the IPA method.
554  Ppath ppath;
555  ppath_calc( ws, ppath, ppath_step_agenda, 3,
556  p_grid, lat_grid, lon_grid, t_field, z_field, vmr_field,
557  edensity_field, Vector(1, f_mono), refellipsoid,
558  z_surface, 0, cloudbox_limits, local_rte_pos, local_rte_los,
559  0, verbosity );
560 
561  while (keepgoing)
562  {
563  //modified path tracing routine for independent pixel approximation
564  mcPathTraceIPA( ws,
565  evol_op, abs_vec_mono, temperature, ext_mat_mono, rng,
566  local_rte_pos, local_rte_los, pnd_vec, g,
567  termination_flag, inside_cloud,
568  propmat_clearsky_agenda,
569  stokes_dim, f_mono, p_grid,
570  lat_grid, lon_grid, z_field, refellipsoid, z_surface,
571  t_field, vmr_field, cloudbox_limits, pnd_field,
572  scat_data_array_mono, z_field_is_1D, ppath,
573  verbosity );
574 
575  np=ppath.np;
576  mc_points(ppath.gp_p[np-1].idx,ppath.gp_lat[np-1].idx,
577  ppath.gp_lon[np-1].idx)+=1;
578  if (termination_flag==1)
579  {
580  iy_space_agendaExecute( ws, local_iy, Vector(1,f_grid[f_index]),
581  local_rte_pos, local_rte_los,
582  iy_space_agenda );
583  mult(vector1,evol_op,local_iy(0,joker));
584  mult(I_i,Q,vector1);
585  I_i/=g;
586  keepgoing=false; //stop here. New photon.
587  }
588  else if (termination_flag==2)
589  {
590  //Choose appropriate lat and lon grid points for IPA
591  //if ppath (the reference ppath) hits the srface just take the last point.
592  GridPos latgp;
593  GridPos longp;
594  if (ppath.background=="surface")
595  {
596  latgp=ppath.gp_lat[ppath.np-1];
597  longp=ppath.gp_lon[ppath.np-1];
598  }
599  else
600  {
601  //Use lat and lon at the lowest z (ie. the tangent point)
602  Numeric latt, lont, zmin=9e99;
603  for( Index i=0; i<ppath.np; i++ )
604  {
605  if( ppath.pos(i,0) < zmin )
606  {
607  zmin = ppath.pos(i,0);
608  latt = ppath.pos(i,1);
609  lont = ppath.pos(i,2);
610  }
611  }
612  gridpos(latgp,lat_grid,latt);
613  gridpos(longp,lon_grid,lont);
614  }
615  //decide whether we have reflection or emission
617  local_surface_emission, local_surface_los,
618  local_surface_rmatrix, Vector(1,f_grid[f_index]),
619  local_rte_pos, local_rte_los, surface_rtprop_agenda );
620  //deal with blackbody case
621  if (local_surface_los.nrows()==0)
622  {
623  mult(vector1,evol_op,local_surface_emission(0,joker));
624  mult(I_i,Q,vector1);
625  I_i/=g;
626  keepgoing=false;
627  }
628  else
629  //decide between reflection and emission
630  {
631  Numeric R11=local_surface_rmatrix(0,0,0,0);
632  if (rng.draw()>R11)
633  {
634  //then we have emission
635  //Matrix oneminusR(stokes_dim,stokes_dim);
636  //id_mat(oneminusR);
637  //oneminusR-=local_surface_rmatrix(0,0,joker,joker);
638  //oneminusR/=1-R11;
639  //mult(vector1,oneminusR,local_surface_emission(0,joker));
640  mult(vector1,evol_op,local_surface_emission(0,joker));
641  mult(I_i,Q,vector1);
642  I_i/=g*(1-R11);
643  keepgoing=false;
644  }
645  else
646  {
647  //we have reflection
648  local_rte_los=local_surface_los(0,joker);
649 
650  mult(q,evol_op,local_surface_rmatrix(0,0,joker,joker));
651  mult(newQ,Q,q);
652  Q=newQ;
653  Q/=g*R11;
654  }
655  }
656  }
657  else if (inside_cloud)
658  {
659  //we have another scattering/emission point
660  //Estimate single scattering albedo
661  albedo=1-abs_vec_mono[0]/ext_mat_mono(0,0);
662  //cout<<"albedo = "<<albedo<<" ext_mat_mono(0,0) = "<<ext_mat_mono(0,0)<<" abs_vec_mono[0] = "<<abs_vec_mono[0]<<"\n";
663  //determine whether photon is emitted or scattered
664  if (rng.draw()>albedo)
665  {
666  //Calculate emission
667  Numeric planck_value = planck( f_grid[0], temperature );
668  Vector emission=abs_vec_mono;
669  emission*=planck_value;
670  Vector emissioncontri(stokes_dim);
671  mult(emissioncontri,evol_op,emission);
672  emissioncontri/=(g*(1-albedo));//yuck!
673  mult(I_i,Q,emissioncontri);
674  keepgoing=false;
675  //cout << "emission contri" << I_i[0] << "\n";
676  }
677  else
678  {
679  //we have a scattering event
680  //Sample new line of sight.
681 
682  Sample_los (new_rte_los,g_los_csc_theta,Z,rng,local_rte_los,
683  scat_data_array_mono,stokes_dim,
684  pnd_vec,anyptype30,Z11maxvector,ext_mat_mono(0,0)-abs_vec_mono[0],temperature,
685  verbosity);
686 
687  Z/=g*g_los_csc_theta*albedo;
688 
689  mult(q,evol_op,Z);
690  mult(newQ,Q,q);
691  Q=newQ;
692  //scattering_order+=1;
693  local_rte_los=new_rte_los;
694  //if (silent==0){cout <<"mc_iteration_count = "<<mc_iteration_count <<
695  // ", scattering_order = " <<scattering_order <<"\n";}
696  }
697  }
698  else
699  {
700  //Must be clear sky emission point
701  //Calculate emission
702  Numeric planck_value = planck( f_grid[0], temperature );
703  Vector emission=abs_vec_mono;
704  emission*=planck_value;
705  Vector emissioncontri(stokes_dim);
706  mult(emissioncontri,evol_op,emission);
707  emissioncontri/=g;
708  mult(I_i,Q,emissioncontri);
709  keepgoing=false;
710  //cout << "emission contri" << I_i[0] << "\n";
711  }
712  }
713  Isum += I_i;
714  for(Index j=0; j<stokes_dim; j++)
715  {
716  assert(!isnan(I_i[j]));
717  Isquaredsum[j] += I_i[j]*I_i[j];
718  }
719  y=Isum;
720  y/=(Numeric)mc_iteration_count;
721  for(Index j=0; j<stokes_dim; j++)
722  {
723  mc_error[j]=sqrt((Isquaredsum[j]/(Numeric)mc_iteration_count-y[j]*y[j])/(Numeric)mc_iteration_count);
724  }
725  if (std_err>0 && mc_iteration_count>=min_iter && mc_error[0]<std_err_i){break;}
726  if (max_time>0 && (Index)(time(NULL)-start_time)>=max_time){break;}
727  if (max_iter>0 && mc_iteration_count>=max_iter){break;}
728  }
729  if ( convert_to_rjbt )
730  {
731  for(Index j=0; j<stokes_dim; j++)
732  {
733  y[j]=invrayjean(y[j],f_grid[0]);
734  mc_error[j]=invrayjean(mc_error[j],f_grid[0]);
735  }
736  }
737 }
738 
739 
740  md_data_raw.push_back
741  ( MdRecord
742  ( NAME( "pha_matExtractManually" ),
744  (
745  "A simple function for manually extract a single phase matrix.\n"
746  "\n"
747  "The function returns the phase matrix for a single particle, for\n"
748  "scattering from (za_in,aa_in) to (za_out,aa_out).\n"
749  "\n"
750  "Only a single particle type is handled and *scat_data_array* must\n"
751  "have length 1. The frequency is selected by *f_grid* and *f_index*.\n"
752  "The temperature is set by *rtp_temperature*.\n"
753  ),
754  AUTHORS( "Patrick Eriksson" ),
755  OUT( ),
756  GOUT( "pha_mat_single" ),
757  GOUT_TYPE( "Matrix" ),
758  GOUT_DESC(
759  "Phase matrix for a single frequency and combination of angles" ),
760  IN( "f_grid", "f_index", "stokes_dim", "scat_data_array",
761  "rtp_temperature" ),
762  GIN( "za_out", "aa_out", "za_in", "aa_in" ),
763  GIN_TYPE( "Numeric", "Numeric", "Numeric", "Numeric" ),
765  GIN_DESC( "Outgoing zenith angle", "Outgoing azimuth angle",
766  "Incoming zenith angle", "Incoming azimuth angle" )
767  ));
768 
769 /* Workspace method: Doxygen documentation will be auto-generated */
771  const Vector& f_grid,
772  const Index& f_index,
773  const Index& stokes_dim,
774  const ArrayOfSingleScatteringData& scat_data_array,
775  const Numeric& rtp_temperature,
776  const Numeric& za_out,
777  const Numeric& aa_out,
778  const Numeric& za_in,
779  const Numeric& aa_in,
780  const Verbosity& verbosity)
781 {
782  if( scat_data_array.nelem() > 1 )
783  throw runtime_error( "Only one element in *scat_data_array* is allowed." );
784 
785  ArrayOfSingleScatteringData scat_data;
786  scat_data_array_monoCalc( scat_data, scat_data_array, f_grid, f_index, verbosity);
787 
788  Vector pnd( 1, 1.0 );
789 
790  pha_mat.resize( stokes_dim, stokes_dim );
791 
792  pha_mat_singleCalc( pha_mat, za_out, aa_out, za_in, aa_in,
793  scat_data, stokes_dim, pnd, rtp_temperature,
794  verbosity );
795 }
796 
797 
799 
811  Vector& cum_l_step,
812  const Ppath& ppath )
813 {
814  cum_l_step.resize(ppath.np);
815  Numeric cumsum = 0.0;
816  cum_l_step[0] = 0.0;
817 
818  for (Index i=0; i<ppath.np-1; i++)
819  {
820  cumsum += ppath.lstep[i];
821  cum_l_step[i+1] = cumsum;
822  }
823 }
824 
825 
827 
840  MatrixView evol_op,
841  Vector& abs_vec_mono,
842  Numeric& temperature,
843  MatrixView ext_mat_mono,
844  Rng& rng,
845  Vector& rte_pos,
846  Vector& rte_los,
847  Vector& pnd_vec,
848  Numeric& g,
849  Index& termination_flag,
850  bool& inside_cloud,
851  const Agenda& propmat_clearsky_agenda,
852  const Index stokes_dim,
853  const Numeric& f_mono,
854  const Vector& p_grid,
855  const Vector& lat_grid,
856  const Vector& lon_grid,
857  const Tensor3& z_field,
858  const Vector& refellipsoid,
859  const Matrix& z_surface,
860  const Tensor3& t_field,
861  const Tensor4& vmr_field,
862  const ArrayOfIndex& cloudbox_limits,
863  const Tensor4& pnd_field,
864  const ArrayOfSingleScatteringData& scat_data_array_mono,
865  const Index z_field_is_1D,
866  const Ppath& ppath,
867  const Verbosity& verbosity)
868 {
869 
870  //Internal declarations
871  ArrayOfMatrix evol_opArray(2);
872  ArrayOfMatrix ext_matArray(2);
873  ArrayOfVector abs_vecArray(2),pnd_vecArray(2);
874  GridPos gp_p,gp_lat,gp_lon;
875  Vector itw(4);
876  Vector tArray(2);
877  Vector cum_l_step;
878  Vector z_grid(p_grid.nelem());
879  Matrix T(stokes_dim,stokes_dim);
880  Numeric k;
881  Numeric x,y,z;
882  Numeric ds,lstep=0.,dx,dy,dz;
883  Index istep = 0; // Counter for number of steps
884  Matrix opt_depth_mat(stokes_dim,stokes_dim),incT(stokes_dim,stokes_dim,0.0);
885  Matrix old_evol_op(stokes_dim,stokes_dim);
886  Numeric rv_ellips=0.;
887  Numeric rv_surface=0.;
888  Numeric alt;
889  Numeric gpnum;
890  const Index np_ipa=ppath.np;
891  Index i_closest=0;
892  Numeric gp_diff;
893  Vector lon_iw(2),lat_iw(2);
894  Vector l_vec(2);
895  GridPos gp;
896  Vector itw1d(2);
897 
898  //Initialisations
899 
900  if (z_field_is_1D)
901  {
902  z_grid=z_field(joker,0,0);
903  rv_ellips=refellipsoid[0];
904  rv_surface=rv_ellips+z_surface(0,0);
905  }
906 
907  id_mat(evol_op);
908  evol_opArray[1]=evol_op;
909  //initialise Ppath with ppath_start_stepping
910  //
911  Ppath ppath_step;
912  //
913  ppath_start_stepping( ppath_step, 3, p_grid, lat_grid,
914  lon_grid, z_field, refellipsoid, z_surface,
915  0, cloudbox_limits, false,
916  rte_pos, rte_los, verbosity );
917 
918  gp_p=ppath_step.gp_p[0];
919  gp_lat=ppath_step.gp_lat[0];
920  gp_lon=ppath_step.gp_lon[0];
921  rte_pos=ppath_step.pos(0,joker);
922  rte_los=ppath_step.los(0,joker);
923 
924  // For simplicity, rte_pos holds the radius (not the altitude) until end of
925  // function
926  const Numeric rre = refell2d( refellipsoid, lat_grid, gp_lat );
927  rte_pos[0] += rre;
928 
929  Range p_range(cloudbox_limits[0],
930  cloudbox_limits[1]-cloudbox_limits[0]+1);
931  Range lat_range(cloudbox_limits[2],
932  cloudbox_limits[3]-cloudbox_limits[2]+1);
933  Range lon_range(cloudbox_limits[4],
934  cloudbox_limits[5]-cloudbox_limits[4]+1);
935  termination_flag=0;
936 
937  inside_cloud=is_inside_cloudbox( ppath_step, cloudbox_limits,false );
938 
939  if (inside_cloud)
940  {
941  cloudy_rt_vars_at_gp(ws, ext_mat_mono,abs_vec_mono,pnd_vec,temperature,
942  propmat_clearsky_agenda, stokes_dim, f_mono,
943  gp_p, gp_lat, gp_lon,p_grid[p_range],
944  t_field(p_range,lat_range,lon_range),
945  vmr_field(joker,p_range,lat_range,lon_range),
946  pnd_field,scat_data_array_mono, cloudbox_limits,rte_los,
947  verbosity);
948  }
949  else
950  {
951  clear_rt_vars_at_gp( ws, ext_mat_mono,abs_vec_mono,temperature,
952  propmat_clearsky_agenda, f_mono,
953  gp_p, gp_lat, gp_lon, p_grid, t_field, vmr_field);
954  pnd_vec=0.0;
955  }
956  ext_matArray[1]=ext_mat_mono;
957  abs_vecArray[1]=abs_vec_mono;
958  tArray[1]=temperature;
959  pnd_vecArray[1]=pnd_vec;
960  //draw random number to determine end point
961  Numeric r = rng.draw();
962  //
963  Numeric ppc, lat0, lon0, za0, aa0;
964  //
965  while ((evol_op(0,0)>r)&(!termination_flag))
966  {
967  //check we are not in an infinite loop (assert for ease of debugging
968  istep++;
969  assert(istep<=5000);
970 
971  //these array variables hold values from the two ends of a path step.
972  //Here we make the values for the last point of the previous step, the values
973  //for the first poinst of the current step.
974  evol_opArray[0]=evol_opArray[1];
975  ext_matArray[0]=ext_matArray[1];
976  abs_vecArray[0]=abs_vecArray[1];
977  tArray[0]=tArray[1];
978  pnd_vecArray[0]=pnd_vecArray[1];
979 
980  //Choose sensible path step length. This is done by considering the
981  //dimensions of the current grid cell and the LOS direction. The aim is
982  //to move only one grid cell.
983  Numeric Dy=(lat_grid[gp_lat.idx+1]-lat_grid[gp_lat.idx])*DEG2RAD*
984  refell2r(refellipsoid,lat_grid[gp_lat.idx]);
985  Numeric Dx=(lon_grid[gp_lon.idx+1]-lon_grid[gp_lon.idx])*
986  DEG2RAD*refell2r(refellipsoid,lat_grid[gp_lat.idx])*
987  cos(DEG2RAD*rte_pos[1]);
988  Numeric Dz=z_field(gp_p.idx+1,gp_lat.idx,gp_lon.idx)-
989  z_field(gp_p.idx,gp_lat.idx,gp_lon.idx);
990  Numeric Dxy=sqrt(Dx*Dx+Dy*Dy);
991  if (Dxy/abs(tan(rte_los[0]*DEG2RAD))>Dz){dz=Dz;}
992  else {dz=Dxy/abs(tan(rte_los[0]*DEG2RAD));}
993  Numeric dxy;
994  if(dz*abs(tan(rte_los[0]*DEG2RAD))>Dxy){dxy=Dxy;}
995  else{dxy=dz*abs(tan(rte_los[0]*DEG2RAD));}
996  lstep=sqrt(dxy*dxy+dz*dz);
997  //calculate new position
998  //current position and direction vector
999  poslos2cart( x, y, z, dx, dy, dz, rte_pos[0], rte_pos[1], rte_pos[2],
1000  rte_los[0], rte_los[1] );
1001  lat0 = rte_pos[1];
1002  lon0 = rte_pos[2];
1003  za0 = rte_los[0];
1004  aa0 = rte_los[1];
1005  ppc = rte_pos[0]*sin(DEG2RAD*rte_los[0]);
1006  //new_position
1007  x+=lstep*dx;
1008  y+=lstep*dy;
1009  z+=lstep*dz;
1010  //back to spherical coords
1011  cart2poslos(rte_pos[0],rte_pos[1],rte_pos[2],rte_los[0],rte_los[1],
1012  x,y,z,dx,dy,dz,ppc,lat0,lon0,za0,aa0);
1013  //get new grid_positions
1014  gridpos( gp_lat, lat_grid, rte_pos[1] );
1015  gridpos( gp_lon, lon_grid, rte_pos[2] );
1016  if (!z_field_is_1D)
1017  {
1018  z_at_latlon( z_grid, p_grid, lat_grid, lon_grid, z_field, gp_lat, gp_lon );
1019  interpweights( itw, gp_lat, gp_lon );
1020  rv_ellips = refell2d(refellipsoid,lat_grid,gp_lat);
1021  rv_surface = rv_ellips + interp( itw, z_surface, gp_lat, gp_lon );
1022  }
1023  alt = rte_pos[0] - rv_ellips;
1024  //Have we left the top of the atmosphere?
1025  if ((alt>=z_grid[p_grid.nelem()-1])&(rte_los[0]<=90))
1026  {
1027  termination_flag=1;
1028  //shift relavent position variables to the appropriate point
1029  //z is linearly interpolated with respect to lat and lon
1030  //I don't think this is overly important. Just change gp_p and rte_pos
1031  alt=z_grid[p_grid.nelem()-1];
1032  rte_pos[0]=alt+rv_ellips;
1033  }
1034  //Have we hit the surface?
1035  if( (rte_pos[0] <= rv_surface)&(rte_los[0]>=90) )
1036  {
1037  termination_flag=2;
1038  rte_pos[0]=rv_surface;
1039  alt = rte_pos[0] - rv_ellips;
1040  }
1041  gridpos( gp_p, z_grid, alt );
1042 
1043  //now project the new point back to the original IPA path
1044  //find lat and lon gridpoints on reference ppath_ipa
1045  gpnum=fractional_gp(gp_p);
1046  //search through ppath_ipa to find the closest grid point
1047  gp_diff=abs(fractional_gp(ppath.gp_p[0])-gpnum);
1048  for (Index i=1;i<np_ipa;i++)
1049  {
1050  Numeric new_gp_diff=abs(fractional_gp(ppath.gp_p[i])-gpnum);
1051  if (new_gp_diff<=gp_diff)
1052  {
1053  gp_diff=new_gp_diff;
1054  i_closest=i;
1055  }
1056  else
1057  {
1058  //getting further away - can stop now
1059  break;
1060  }
1061  }
1062  gp_lat=ppath.gp_lat[i_closest];
1063  gp_lon=ppath.gp_lon[i_closest];
1064  interpweights(lon_iw,gp_lon);
1065  interpweights(lat_iw,gp_lat);
1066  //
1067  if (!z_field_is_1D)
1068  {
1069  interpweights( itw, gp_lat, gp_lon );
1070  rv_ellips = refell2d(refellipsoid,lat_grid,gp_lat);
1071  }
1072  rte_pos[0]=rv_ellips+interp_atmfield_by_gp(3,z_field,
1073  gp_p,gp_lat,gp_lon);
1074  rte_pos[1]=interp(lat_iw, lat_grid,gp_lat);
1075  rte_pos[2]=interp(lon_iw, lon_grid,gp_lon);
1076 
1077  //calculate evolution operator
1078  inside_cloud=is_gp_inside_cloudbox(gp_p,gp_lat,gp_lon,
1079  cloudbox_limits, false );
1080  //calculate RT variables at new point
1081  if (inside_cloud)
1082  {
1084  ext_mat_mono, abs_vec_mono, pnd_vec, temperature,
1085  propmat_clearsky_agenda, stokes_dim,
1086  f_mono, gp_p,gp_lat, gp_lon, p_grid[p_range],
1087  t_field(p_range,lat_range,lon_range),
1088  vmr_field(joker,p_range,lat_range,lon_range),
1089  pnd_field, scat_data_array_mono, cloudbox_limits,rte_los,
1090  verbosity);
1091  }
1092  else
1093  {
1094  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
1095  propmat_clearsky_agenda, f_mono,
1096  gp_p, gp_lat, gp_lon, p_grid, t_field, vmr_field);
1097  pnd_vec=0.0;
1098  }
1099  //put these variables in the last Array slot
1100  ext_matArray[1]=ext_mat_mono;
1101  abs_vecArray[1]=abs_vec_mono;
1102  tArray[1]=temperature;
1103  pnd_vecArray[1]=pnd_vec;
1104  opt_depth_mat=ext_matArray[1];
1105  //now calculate evolution operator
1106  opt_depth_mat+=ext_matArray[0];
1107  opt_depth_mat*=-lstep/2;
1108  incT=0;
1109  if ( stokes_dim == 1 )
1110  {
1111  incT(0,0)=exp(opt_depth_mat(0,0));
1112  }
1113  else if ( is_diagonal( opt_depth_mat))
1114  {
1115  for ( Index j=0;j<stokes_dim;j++)
1116  {
1117  incT(j,j)=exp(opt_depth_mat(j,j));
1118  }
1119  }
1120  else
1121  {
1122  matrix_exp_p30(incT,opt_depth_mat);
1123  }
1124  mult(evol_op,evol_opArray[0],incT);
1125  assert(incT(0,0)<1);
1126  assert(evol_op(0,0)<1);
1127  evol_opArray[1]=evol_op;
1128 
1129  }
1130  if (termination_flag)
1131  {
1132  //we must have reached the surface or the TOA
1133  g=evol_op(0,0);
1134  }
1135  else
1136  {
1137  //then we have met the monte carlo pathlength criteria somewhere inside
1138  //the atmosphere and between the two points corresponding to the *Array
1139  //variables.
1140 
1141  // If istep is 0 means we never entered the above while loop and
1142  // thus x, y, z, ds, dx, dy, dz are uninitialized
1143  assert(istep);
1144 
1145  //find position corresponding to the pathlength criteria
1146  k=-log(incT(0,0))/lstep;
1147  //length of path segment required by critera
1148  ds=log(evol_opArray[0](0,0)/r)/k;
1149  g=k*r;
1150  l_vec=0.0;
1151  l_vec[1]=lstep;
1152  //gridpos and interpolations for required step length
1153  gridpos(gp, l_vec, ds);
1154  interpweights(itw1d,gp);
1155  //interpolate optical properties at required point
1156  interp(ext_mat_mono, itw1d, ext_matArray, gp);
1157  opt_depth_mat = ext_mat_mono;
1158  opt_depth_mat+=ext_matArray[gp.idx];
1159  opt_depth_mat*=-ds/2;
1160  if ( stokes_dim == 1 )
1161  {
1162  incT(0,0)=exp(opt_depth_mat(0,0));
1163  }
1164  else if ( is_diagonal( opt_depth_mat))
1165  {
1166  for ( Index i=0;i<stokes_dim;i++)
1167  {
1168  incT(i,i)=exp(opt_depth_mat(i,i));
1169  }
1170  }
1171  else
1172  {
1173  //matrix_exp(incT,opt_depth_mat,10);
1174  matrix_exp_p30(incT,opt_depth_mat);
1175  }
1176  mult(evol_op,evol_opArray[gp.idx],incT);
1177  interp(abs_vec_mono, itw1d, abs_vecArray,gp);
1178  temperature=interp(itw1d,tArray, gp);
1179  interp(pnd_vec, itw1d, pnd_vecArray, gp);
1180  //move cartesion coordinates back to required point.
1181  x+=(ds-lstep)*dx;
1182  y+=(ds-lstep)*dy;
1183  z+=(ds-lstep)*dz;
1184  //and convert back to spherical.
1185  cart2poslos(rte_pos[0],rte_pos[1],rte_pos[2],rte_los[0],rte_los[1],
1186  x,y,z,dx,dy,dz,ppc,lat0,lon0,za0,aa0);
1187  }
1188 
1189  // Go back to altitude in rte_pos
1190  rte_pos[0] -= rre;
1191 
1192 }
1193 
1194 
1195 
1197 
1208  Numeric& iwp,
1209  Numeric& cloud_opt_path,
1210  //input
1211  const Vector& rte_pos,
1212  const Vector& rte_los,
1213  const Agenda& ppath_step_agenda,
1214  const Vector& p_grid,
1215  const Vector& lat_grid,
1216  const Vector& lon_grid,
1217  const Vector& refellipsoid,
1218  const Matrix& z_surface,
1219  const Tensor3& z_field,
1220  const Tensor3& t_field,
1221  const Tensor4& vmr_field,
1222  const Tensor3& edensity_field,
1223  const Numeric& f_mono,
1224  const ArrayOfIndex& cloudbox_limits,
1225  const Tensor4& pnd_field,
1226  const ArrayOfSingleScatteringData& scat_data_array_mono,
1227  const Matrix& particle_masses,
1228  const Verbosity& verbosity)
1229 {
1230  //internal declarations
1231  Ppath ppath;
1232  Vector local_rte_pos=rte_pos;
1233  Vector local_rte_los=rte_los;
1234  iwp=0;
1235  cloud_opt_path=0;
1236  //calculate ppath to cloudbox boundary
1237  ppath_calc( ws, ppath, ppath_step_agenda, 3,
1238  p_grid, lat_grid, lon_grid, t_field, z_field, vmr_field,
1239  edensity_field, Vector(1, f_mono), refellipsoid,
1240  z_surface, 1, cloudbox_limits, local_rte_pos, local_rte_los, 0,
1241  verbosity );
1242  //if this ppath hit a cloud, now take ppath inside cloud
1243  if (ppath_what_background(ppath)>2)
1244  {
1245  //move to the intersection of the cloudbox boundary and the
1246  //propagation path
1247  local_rte_pos=ppath.pos(ppath.np-1,joker);
1248  local_rte_los=ppath.los(ppath.np-1,joker);
1249 
1250  Range p_range(cloudbox_limits[0],
1251  cloudbox_limits[1]-cloudbox_limits[0]+1);
1252  Range lat_range(cloudbox_limits[2],
1253  cloudbox_limits[3]-cloudbox_limits[2]+1);
1254  Range lon_range(cloudbox_limits[4],
1255  cloudbox_limits[5]-cloudbox_limits[4]+1);
1256 
1257  ppath_calc( ws, ppath, ppath_step_agenda, 3,
1258  p_grid, lat_grid, lon_grid, t_field, z_field, vmr_field,
1259  edensity_field, Vector(1, f_mono), refellipsoid,
1260  z_surface, 1, cloudbox_limits, local_rte_pos, local_rte_los,
1261  1, verbosity );
1262 
1263  Matrix pnd_ppath(particle_masses.nrows(),ppath.np);
1264  Vector t_ppath(ppath.np);
1265  Vector p_ppath(ppath.np);
1266  Matrix vmr_ppath(vmr_field.nbooks(),ppath.np);
1267  Matrix ext_mat_part(1, 1, 0.0);
1268  Vector abs_vec_part(1, 0.0);
1269 
1270  cloud_atm_vars_by_gp(p_ppath,t_ppath,vmr_ppath,pnd_ppath,ppath.gp_p,
1271  ppath.gp_lat,ppath.gp_lon,cloudbox_limits,p_grid[p_range],
1272  t_field(p_range,lat_range,lon_range),
1273  vmr_field(joker,p_range,lat_range,lon_range),
1274  pnd_field);
1275 
1276  Vector k_vec(ppath.np);
1277  Vector iwc_vec(ppath.np);
1278  //calculate integrands
1279  for (Index i = 0; i < ppath.np ; ++i)
1280  {
1281  opt_propCalc(ext_mat_part,abs_vec_part,ppath.los(i,0),ppath.los(i,1),scat_data_array_mono,
1282  1, pnd_ppath(joker, i), t_ppath[i], verbosity);
1283  k_vec[i]=ext_mat_part(0,0);
1284  Vector pnd_vec=pnd_ppath(joker, i);
1285  assert(particle_masses.ncols()==1);
1286  assert(pnd_vec.nelem()==particle_masses.nrows());
1287  iwc_vec[i]=pnd_vec*particle_masses(joker,0);//hopefully this is the dot product
1288  }
1289  //integrate IWP and optical properties
1290  for (Index i = 0; i < ppath.np-1 ; ++i)
1291  {
1292  //Add to path integrals
1293  iwp+=0.5*(iwc_vec[i+1]+iwc_vec[i])*ppath.lstep[i];
1294  cloud_opt_path+=0.5*(k_vec[i+1]+k_vec[i])*ppath.lstep[i];
1295  }
1296  }
1297 }
1298 
1299 
1300 
1302 
1330 //interpolates TArray and PPath to give T and rte_pos(los) at a given pathlength
1332  Vector& K_abs,
1333  Numeric& temperature,
1334  MatrixView& K,
1335  Vector& rte_pos,
1336  Vector& rte_los,
1337  VectorView& pnd_vec,
1338  const ArrayOfMatrix& TArray,
1339  const ArrayOfMatrix& ext_matArray,
1340  const ArrayOfVector& abs_vecArray,
1341  const Vector& t_ppath,
1342  const Matrix& pnd_ppath,
1343  const Vector& cum_l_step,
1344  const Numeric& pathlength,
1345  const Index& stokes_dim,
1346  const Ppath& ppath
1347  )
1348 {
1349  //Internal Declarations
1350  Matrix incT(stokes_dim,stokes_dim,0.0);
1351  Matrix opt_depth_mat(stokes_dim,stokes_dim);
1352  Vector itw(2);
1353  Numeric delta_s;
1354  Index N_pt=pnd_vec.nelem();
1355  ArrayOfGridPos gp(1);
1356 
1357  //interpolate transmittance matrix
1358  gridpos(gp, cum_l_step, pathlength);
1359  interpweights(itw,gp[0]);
1360  interp(K, itw,ext_matArray,gp[0]);
1361  delta_s = pathlength - cum_l_step[gp[0].idx];
1362  opt_depth_mat = K;
1363  opt_depth_mat+=ext_matArray[gp[0].idx];
1364  opt_depth_mat*=-delta_s/2;
1365  if ( stokes_dim == 1 )
1366  {
1367  incT(0,0)=exp(opt_depth_mat(0,0));
1368  }
1369  else if ( is_diagonal( opt_depth_mat))
1370  {
1371  for ( Index i=0;i<stokes_dim;i++)
1372  {
1373  incT(i,i)=exp(opt_depth_mat(i,i));
1374  }
1375  }
1376  else
1377  {
1378  //matrix_exp(incT,opt_depth_mat,10);
1379  matrix_exp_p30(incT,opt_depth_mat);
1380  }
1381  mult(T,TArray[gp[0].idx],incT);
1382 
1383  interp(K_abs, itw, abs_vecArray,gp[0]);
1384 
1385  temperature=interp(itw,t_ppath,gp[0]);
1386 
1387  for (Index i=0;i<N_pt;i++)
1388  {
1389  pnd_vec[i]=interp(itw,pnd_ppath(i,Range(joker)),gp[0]);
1390  }
1391 
1392  for (Index i=0; i<2; i++)
1393  {
1394  rte_pos[i] = interp(itw,ppath.pos(Range(joker),i),gp[0]);
1395  rte_los[i] = interp(itw,ppath.los(Range(joker),i),gp[0]);
1396  }
1397  rte_pos[2] = interp(itw,ppath.pos(Range(joker),2),gp[0]);
1398 }
1399 
1400 
1401 
1403  Workspace& ws,
1404  MatrixView evol_op,
1405  Vector& abs_vec_mono,
1406  Numeric& temperature,
1407  MatrixView ext_mat_mono,
1408  Rng& rng,
1409  Vector& rte_pos,
1410  Vector& rte_los,
1411  Vector& pnd_vec,
1412  Numeric& g,
1413  Ppath& ppath_step,
1414  Index& termination_flag,
1415  bool& inside_cloud,
1416  const Agenda& ppath_step_agenda,
1417  const Numeric& ppath_lraytrace,
1418  const Agenda& propmat_clearsky_agenda,
1419  const Index stokes_dim,
1420  const Numeric& f_mono,
1421  const Vector& p_grid,
1422  const Vector& lat_grid,
1423  const Vector& lon_grid,
1424  const Tensor3& z_field,
1425  const Vector& refellipsoid,
1426  const Matrix& z_surface,
1427  const Tensor3& t_field,
1428  const Tensor4& vmr_field,
1429  const ArrayOfIndex& cloudbox_limits,
1430  const Tensor4& pnd_field,
1431  const ArrayOfSingleScatteringData& scat_data_array_mono,
1432  const Verbosity& verbosity )
1433 {
1434  ArrayOfMatrix evol_opArray(2);
1435  ArrayOfMatrix ext_matArray(2);
1436  ArrayOfVector abs_vecArray(2);
1437  ArrayOfVector pnd_vecArray(2);
1438  Matrix opt_depth_mat(stokes_dim,stokes_dim);
1439  Matrix incT(stokes_dim,stokes_dim,0.0);
1440  Vector tArray(2);
1441  Vector f_grid(1,f_mono); // Vector version of f_mono
1442  Matrix T(stokes_dim,stokes_dim);
1443  Numeric k;
1444  Numeric ds, dl=-999;
1445  Index istep = 0; // Counter for number of steps
1446  Matrix old_evol_op(stokes_dim,stokes_dim);
1447 
1448  CREATE_OUT0;
1449 
1450  //at the start of the path the evolution operator is the identity matrix
1451  id_mat(evol_op);
1452  evol_opArray[1]=evol_op;
1453 
1454  // range defining cloudbox
1455  Range p_range( cloudbox_limits[0],
1456  cloudbox_limits[1]-cloudbox_limits[0]+1);
1457  Range lat_range( cloudbox_limits[2],
1458  cloudbox_limits[3]-cloudbox_limits[2]+1 );
1459  Range lon_range( cloudbox_limits[4],
1460  cloudbox_limits[5]-cloudbox_limits[4]+1 );
1461 
1462  //initialise Ppath with ppath_start_stepping
1463  ppath_start_stepping( ppath_step, 3, p_grid, lat_grid,
1464  lon_grid, z_field, refellipsoid, z_surface,
1465  0, cloudbox_limits, false,
1466  rte_pos, rte_los, verbosity );
1467 
1468  // Index in ppath_step of end point considered presently
1469  Index ip = 0;
1470 
1471  // Is point ip inside the cloudbox?
1472  inside_cloud = is_gp_inside_cloudbox( ppath_step.gp_p[ip],
1473  ppath_step.gp_lat[ip],
1474  ppath_step.gp_lon[ip],
1475  cloudbox_limits, 0, 3 );
1476 
1477  // Determine radiative properties at point
1478  if( inside_cloud )
1479  {
1480  cloudy_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, pnd_vec,
1481  temperature, propmat_clearsky_agenda,
1482  stokes_dim, f_mono, ppath_step.gp_p[0],
1483  ppath_step.gp_lat[0], ppath_step.gp_lon[0],
1484  p_grid[p_range],
1485  t_field(p_range,lat_range,lon_range),
1486  vmr_field(joker,p_range,lat_range,lon_range),
1487  pnd_field, scat_data_array_mono, cloudbox_limits,
1488  ppath_step.los(0,joker), verbosity );
1489  }
1490  else
1491  {
1492  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
1493  propmat_clearsky_agenda, f_mono,
1494  ppath_step.gp_p[0], ppath_step.gp_lat[0],
1495  ppath_step.gp_lon[0], p_grid, t_field, vmr_field );
1496  pnd_vec = 0.0;
1497  }
1498 
1499  // Move the data to end point containers
1500  ext_matArray[1] = ext_mat_mono;
1501  abs_vecArray[1] = abs_vec_mono;
1502  tArray[1] = temperature;
1503  pnd_vecArray[1] = pnd_vec;
1504 
1505  //draw random number to determine end point
1506  Numeric r = rng.draw();
1507 
1508  termination_flag=0;
1509 
1510  while( (evol_op(0,0)>r) && (!termination_flag) )
1511  {
1512  istep++;
1513 
1514  if( istep > 5000 )
1515  {
1516  throw runtime_error( "5000 path points have been reached. "
1517  "Is this an infinite loop?" );
1518  }
1519 
1520  evol_opArray[0] = evol_opArray[1];
1521  ext_matArray[0] = ext_matArray[1];
1522  abs_vecArray[0] = abs_vecArray[1];
1523  tArray[0] = tArray[1];
1524  pnd_vecArray[0] = pnd_vecArray[1];
1525 
1526  // Shall new ppath_step be calculated?
1527  if( ip == ppath_step.np-1 )
1528  {
1529  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field,
1530  z_field, vmr_field, f_grid,
1531  ppath_step_agenda );
1532  //Print( ppath_step, 0, verbosity );
1533  ip = 1;
1534 
1535  inside_cloud = is_gp_inside_cloudbox( ppath_step.gp_p[ip],
1536  ppath_step.gp_lat[ip],
1537  ppath_step.gp_lon[ip],
1538  cloudbox_limits, 0, 3 );
1539  }
1540  else
1541  { ip++; }
1542 
1543  dl = ppath_step.lstep[ip-1];
1544 
1545  if( inside_cloud )
1546  {
1547  cloudy_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, pnd_vec,
1548  temperature, propmat_clearsky_agenda,
1549  stokes_dim, f_mono, ppath_step.gp_p[ip],
1550  ppath_step.gp_lat[ip], ppath_step.gp_lon[ip],
1551  p_grid[p_range],
1552  t_field(p_range,lat_range,lon_range),
1553  vmr_field(joker,p_range,lat_range,lon_range),
1554  pnd_field, scat_data_array_mono, cloudbox_limits,
1555  ppath_step.los(ip,joker), verbosity );
1556  }
1557  else
1558  {
1559  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
1560  propmat_clearsky_agenda, f_mono,
1561  ppath_step.gp_p[ip], ppath_step.gp_lat[ip],
1562  ppath_step.gp_lon[ip], p_grid, t_field,
1563  vmr_field );
1564  pnd_vec = 0.0;
1565  }
1566 
1567  ext_matArray[1] = ext_mat_mono;
1568  abs_vecArray[1] = abs_vec_mono;
1569  tArray[1] = temperature;
1570  pnd_vecArray[1] = pnd_vec;
1571  opt_depth_mat = ext_matArray[1];
1572  opt_depth_mat += ext_matArray[0];
1573  opt_depth_mat *= -dl/2;
1574  incT = 0;
1575 
1576  if( opt_depth_mat(0,0) < -4 )
1577  {
1578  out0 << "WARNING: A MC path step of high optical depth ("
1579  << abs(opt_depth_mat(0,0)) << ")!\n";
1580  }
1581 
1582  if( stokes_dim == 1 )
1583  { incT(0,0) = exp( opt_depth_mat(0,0) ); }
1584  else if( is_diagonal( opt_depth_mat ) )
1585  {
1586  for ( Index j=0;j<stokes_dim;j++)
1587  { incT(j,j) = exp( opt_depth_mat(j,j) ); }
1588  }
1589  else
1590  { matrix_exp_p30( incT, opt_depth_mat ); }
1591  mult( evol_op, evol_opArray[0], incT );
1592  evol_opArray[1] = evol_op;
1593 
1594  if( evol_op(0,0)>r )
1595  {
1596  // Check whether hit ground or space.
1597  // path_step_agenda just detects surface intersections, and
1598  // if TOA is reached requires a special check.
1599  // But we are already ready if evol_op<=r
1600  if( ip == ppath_step.np - 1 )
1601  {
1602  if( ppath_what_background(ppath_step) )
1603  { termination_flag = 2; } //we have hit the surface
1604  else if( fractional_gp(ppath_step.gp_p[ip]) >=
1605  (Numeric)(p_grid.nelem() - 1)-1e-3 )
1606  { termination_flag = 1; } //we are at TOA
1607  }
1608  }
1609  } // while
1610 
1611 
1612  if( termination_flag )
1613  { //we must have reached the cloudbox boundary
1614  rte_pos = ppath_step.pos(ip,joker);
1615  rte_los = ppath_step.los(ip,joker);
1616  g = evol_op(0,0);
1617  }
1618  else
1619  {
1620  //find position...and evol_op..and everything else required at the new
1621  //scattering/emission point
1622  // GH 2011-09-14:
1623  // log(incT(0,0)) = log(exp(opt_depth_mat(0, 0))) = opt_depth_mat(0, 0)
1624  // Avoid loss of precision, use opt_depth_mat directly
1625  //k=-log(incT(0,0))/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
1626  k = -opt_depth_mat(0,0) / dl;
1627  ds = log( evol_opArray[0](0,0) / r ) / k;
1628  g = k*r;
1629  Vector x(2,0.0);
1630  //interpolate atmospheric variables required later
1631  ArrayOfGridPos gp(1);
1632  x[1] = dl;
1633  Vector itw(2);
1634 
1635  gridpos( gp, x, ds );
1636  assert( gp[0].idx == 0 );
1637  interpweights( itw, gp[0] );
1638  interp( ext_mat_mono, itw, ext_matArray, gp[0] );
1639  opt_depth_mat = ext_mat_mono;
1640  opt_depth_mat += ext_matArray[gp[0].idx];
1641  opt_depth_mat *= -ds/2;
1642  if( stokes_dim == 1 )
1643  { incT(0,0) = exp( opt_depth_mat(0,0) ); }
1644  else if( is_diagonal( opt_depth_mat) )
1645  {
1646  for ( Index i=0;i<stokes_dim;i++)
1647  { incT(i,i) = exp( opt_depth_mat(i,i) ); }
1648  }
1649  else
1650  { matrix_exp_p30( incT, opt_depth_mat ); }
1651  mult( evol_op, evol_opArray[gp[0].idx], incT );
1652  interp( abs_vec_mono, itw, abs_vecArray,gp[0] );
1653  temperature = interp( itw, tArray, gp[0] );
1654  interp( pnd_vec, itw, pnd_vecArray, gp[0] );
1655  for( Index i=0; i<2; i++ )
1656  {
1657  rte_pos[i] = interp( itw, ppath_step.pos(Range(ip-1,2),i), gp[0] );
1658  rte_los[i] = interp( itw, ppath_step.los(Range(ip-1,2),i), gp[0] );
1659  }
1660  rte_pos[2] = interp( itw, ppath_step.pos(Range(ip-1,2),2), gp[0] );
1661  }
1662 
1663  assert(isfinite(g));
1664 
1665  // A dirty trick to avoid copying ppath_step
1666  const Index np = ip+1;
1667  ppath_step.np = np;
1668 
1669 }
1670 
1671 
1672 
1674 
1686  ConstMatrixView A)
1687 {
1688  Index m=A.nrows();
1689  assert( A.ncols()==m );
1690  M=0;
1691  Numeric a=A(0,0);
1692  Numeric b=A(0,1);
1693  M(0,0)=cosh(b);
1694  //M(1,1)=cosh(b);
1695  M(1,1)=M(0,0);
1696  M(0,1)=sinh(b);
1697  //M(1,0)=sinh(b);
1698  M(1,0)=M(0,1);
1699  if ( m>2 )
1700  {
1701  Numeric c=A(2,3);
1702  M(2,2)=cos(c);
1703  if ( m > 3 )
1704  {
1705  M(2,3)=sin(c);
1706  //M(3,2)=-sin(c);
1707  M(3,2)=-M(2,3);
1708  M(3,3)=M(2,2);
1709  //M(3,3)=cos(c); // Added by GH 2011-06-15 as per e-mail 2011-06-13
1710  }
1711  }
1712  M*=exp(a);
1713 }
1714 
1715 
1716 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
#define GIN_TYPE
Definition: methods.cc:53
ArrayOfGridPos gp_lat
Definition: ppath.h:77
void clear_rt_vars_at_gp(Workspace &ws, MatrixView ext_mat_mono, VectorView abs_vec_mono, Numeric &temperature, const Agenda &propmat_clearsky_agenda, const Numeric &f_mono, const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, ConstVectorView p_grid, ConstTensor3View t_field, ConstTensor4View vmr_field)
clear_rt_vars_at_gp
Definition: montecarlo.cc:56
#define GIN
Definition: methods.cc:52
The VectorView class.
Definition: matpackI.h:372
The Agenda class.
Definition: agenda_class.h:44
#define M
Definition: rng.cc:196
Index nelem() const
Number of elements.
Definition: array.h:176
void interp_atmfield_by_gp(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
interp_atmfield_by_gp
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
#define q
Definition: continua.cc:21469
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
void findZ11max(Vector &Z11maxvector, const ArrayOfSingleScatteringData &scat_data_array_mono)
findZ11max
Definition: montecarlo.cc:337
Matrix los
Definition: ppath.h:68
void ppath_start_stepping(Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const bool &ppath_inside_cloudbox_do, ConstVectorView rte_pos, ConstVectorView rte_los, const Verbosity &verbosity)
ppath_start_stepping
Definition: ppath.cc:4606
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void cum_l_stepCalc(Vector &cum_l_step, const Ppath &ppath)
cum_l_stepCalc
Definition: mc_NotUsed.cc:810
void draw_los(VectorView &sampled_rte_los, Rng &rng, ConstVectorView bore_sight_los) const
draws a line of sight by sampling the antenna response function
Definition: mc_antenna.cc:167
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lraytrace, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:14720
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
Vector lstep
Definition: ppath.h:70
void ppath_step_geom_3d(Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax)
ppath_step_geom_3d
Definition: ppath.cc:3553
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Matrix pos
Definition: ppath.h:67
void pha_matExtractManually(Matrix &pha_mat, const Vector &f_grid, const Index &f_index, const Index &stokes_dim, const ArrayOfSingleScatteringData &scat_data_array, const Numeric &rtp_temperature, const Numeric &za_out, const Numeric &aa_out, const Numeric &za_in, const Numeric &aa_in, const Verbosity &verbosity)
Definition: mc_NotUsed.cc:770
#define GIN_DEFAULT
Definition: methods.cc:54
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1244
All information for one workspace method.
Definition: methods.h:42
#define GROUP(x)
Definition: workspace.cc:47
void mcPathTraceIPA(Workspace &ws, MatrixView evol_op, Vector &abs_vec_mono, Numeric &temperature, MatrixView ext_mat_mono, Rng &rng, Vector &rte_pos, Vector &rte_los, Vector &pnd_vec, Numeric &g, Index &termination_flag, bool &inside_cloud, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Numeric &f_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index z_field_is_1D, const Ppath &ppath, const Verbosity &verbosity)
mcPathTraceIPA
Definition: mc_NotUsed.cc:839
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Structure to store a grid position.
Definition: interpolation.h:74
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
String background
Definition: ppath.h:63
#define NODEF
Definition: methods.h:36
void Sample_los(VectorView new_rte_los, Numeric &g_los_csc_theta, MatrixView Z, Rng &rng, ConstVectorView rte_los, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const bool anyptype30, ConstVectorView Z11maxvector, const Numeric Csca, const Numeric rtp_temperature, const Verbosity &verbosity)
Sample_los.
Definition: montecarlo.cc:1187
void scat_data_array_monoCalc(ArrayOfSingleScatteringData &scat_data_array_mono, const ArrayOfSingleScatteringData &scat_data_array, const Vector &f_grid, const Index &f_index, const Verbosity &)
WORKSPACE METHOD: scat_data_array_monoCalc.
void cart2poslos(Numeric &r, Numeric &lat, Numeric &za, const Numeric &x, const Numeric &z, const Numeric &dx, const Numeric &dz, const Numeric &ppc, const Numeric &lat0, const Numeric &za0)
cart2poslos
Definition: geodetic.cc:125
void matrix_exp_p30(MatrixView M, ConstMatrixView A)
matrix_exp_p30
Definition: mc_NotUsed.cc:1685
#define OUT
Definition: methods.cc:47
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
bool is_inside_cloudbox(const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, const bool include_boundaries)
Definition: cloudbox.cc:892
#define NAME(x)
Definition: agendas.cc:34
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:13886
#define abs(x)
Definition: continua.cc:20458
bool is_diagonal(ConstMatrixView A)
Checks if a square matrix is diagonal.
Definition: logic.cc:411
const Joker joker
Index idx
Definition: interpolation.h:75
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define dx
Definition: continua.cc:21928
void opt_propCalc(MatrixView ext_mat_mono, VectorView abs_vec_mono, const Numeric za, const Numeric aa, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const Numeric rtp_temperature, const Verbosity &verbosity)
opt_propCalc
Definition: montecarlo.cc:684
The Matrix class.
Definition: matpackI.h:788
void interpTArray(Matrix &T, Vector &K_abs, Numeric &temperature, MatrixView &K, Vector &rte_pos, Vector &rte_los, VectorView &pnd_vec, const ArrayOfMatrix &TArray, const ArrayOfMatrix &ext_matArray, const ArrayOfVector &abs_vecArray, const Vector &t_ppath, const Matrix &pnd_ppath, const Vector &cum_l_step, const Numeric &pathlength, const Index &stokes_dim, const Ppath &ppath)
interpTarray
Definition: mc_NotUsed.cc:1331
void z_at_latlon(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, const GridPos &gp_lat, const GridPos &gp_lon)
z_at_latlon
#define DESCRIPTION(x)
Definition: agendas.cc:35
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:56
#define GOUT_DESC
Definition: methods.cc:50
void iwp_cloud_opt_pathCalc(Workspace &ws, Numeric &iwp, Numeric &cloud_opt_path, const Vector &rte_pos, const Vector &rte_los, const Agenda &ppath_step_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const Numeric &f_mono, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Matrix &particle_masses, const Verbosity &verbosity)
iwp_cloud_opt_pathCalc
Definition: mc_NotUsed.cc:1207
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
#define GIN_DESC
Definition: methods.cc:55
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
double draw()
Definition: rng.cc:120
This can be used to make arrays out of anything.
Definition: array.h:40
void MCIPA(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &ppath_step_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index &mc_seed, const String &y_unit, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Index &z_field_is_1D, const Verbosity &verbosity)
Definition: mc_NotUsed.cc:433
void mcPathTraceGeneral(Workspace &ws, MatrixView evol_op, Vector &abs_vec_mono, Numeric &temperature, MatrixView ext_mat_mono, Rng &rng, Vector &rte_pos, Vector &rte_los, Vector &pnd_vec, Numeric &g, Ppath &ppath_step, Index &termination_flag, bool &inside_cloud, const Agenda &ppath_step_agenda, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Numeric &f_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Verbosity &verbosity)
mcPathTraceGeneral
Definition: mc_NotUsed.cc:150
Definition: rng.h:569
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
void cloud_atm_vars_by_gp(VectorView pressure, VectorView temperature, MatrixView vmr, MatrixView pnd, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, ConstVectorView p_grid_cloud, ConstTensor3View t_field_cloud, ConstTensor4View vmr_field_cloud, ConstTensor4View pnd_field)
cloud_atm_vars_by_gp
Definition: montecarlo.cc:243
Index np
Definition: ppath.h:61
ArrayOfGridPos gp_lon
Definition: ppath.h:78
void cloudy_rt_vars_at_gp(Workspace &ws, MatrixView ext_mat_mono, VectorView abs_vec_mono, VectorView pnd_vec, Numeric &temperature, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Numeric &f_mono, const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, ConstVectorView p_grid_cloud, ConstTensor3View t_field_cloud, ConstTensor4View vmr_field_cloud, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const ArrayOfIndex &cloudbox_limits, const Vector &rte_los, const Verbosity &verbosity)
cloudy_rt_vars_at_gp
Definition: montecarlo.cc:135
A constant view of a Matrix.
Definition: matpackI.h:596
Numeric planck(const Numeric &f, const Numeric &t)
planck
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
bool is_anyptype30(const ArrayOfSingleScatteringData &scat_data_array_mono)
is_anyptype30
Definition: montecarlo.cc:373
#define CREATE_OUT0
Definition: messages.h:211
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
ppath_calc
Definition: ppath.cc:5344
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
#define AUTHORS
Definition: methods.cc:46
#define GOUT
Definition: methods.cc:48
void seed(unsigned long int n, const Verbosity &verbosity)
Definition: rng.cc:72
void mc_IWP_cloud_opt_pathCalc(Workspace &ws, Numeric &mc_IWP, Numeric &mc_cloud_opt_path, Numeric &mc_IWP_error, Numeric &mc_cloud_opt_path_error, Index &mc_iteration_count, const MCAntenna &mc_antenna, const Matrix &sensor_pos, const Matrix &sensor_los, const Agenda &ppath_step_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const Vector &f_grid, const Index &f_index, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Matrix &particle_masses, const Index &mc_seed, const Index &max_iter, const Verbosity &verbosity)
Definition: mc_NotUsed.cc:43
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
This class contains all static information for one workspace variable.
Definition: wsv_aux.h:54
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
Definition: cloudbox.cc:808
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void surface_rtprop_agendaExecute(Workspace &ws, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:15046
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:387
const Numeric SPEED_OF_LIGHT
#define GOUT_TYPE
Definition: methods.cc:49
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
ArrayOfGridPos gp_p
Definition: ppath.h:76
#define IN
Definition: methods.cc:51
void pha_mat_singleCalc(MatrixView Z, const Numeric za_sca, const Numeric aa_sca, const Numeric za_inc, const Numeric aa_inc, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const Numeric rtp_temperature, const Verbosity &verbosity)
pha_mat_singleCalc
Definition: montecarlo.cc:928
AntennaType get_type(void) const
returns the antenna type
Definition: mc_antenna.cc:150
const Array< MdRecord > md_data_raw
Lookup information for workspace methods.
Definition: methods.cc:39
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580