ARTS  2.2.66
montecarlo.cc
Go to the documentation of this file.
1 /* copyright (C) 2003-2012 Cory Davis <cory.davis@metservice.com>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
18 
19 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
33 /*===========================================================================
34  === External declarations
35  ===========================================================================*/
36 
37 #include "auto_md.h"
38 #include "geodetic.h"
39 #include "montecarlo.h"
40 #include "mc_interp.h"
41 
42 #ifdef HAVE_SSTREAM
43 #include <sstream>
44 #else
45 #include "sstream.h"
46 #endif
47 
49 
57  MatrixView ext_mat_mono,
58  VectorView abs_vec_mono,
59  Numeric& temperature,
60  const Agenda& propmat_clearsky_agenda,
61  const Numeric& f_mono,
62  const GridPos& gp_p,
63  const GridPos& gp_lat,
64  const GridPos& gp_lon,
65  ConstVectorView p_grid,
66  ConstTensor3View t_field,
67  ConstTensor4View vmr_field)
68 {
69  const Index ns = vmr_field.nbooks();
70  Vector t_vec(1); //vectors are required by interp_atmfield_gp2itw etc.
71  Vector p_vec(1); //may not be efficient with unecessary vectors
72  Matrix itw_p(1,2);
73  ArrayOfGridPos ao_gp_p(1),ao_gp_lat(1),ao_gp_lon(1);
74  Matrix vmr_mat(ns,1), itw_field;
75 
76  //local versions of workspace variables
77  Matrix local_abs_vec;
78  Tensor3 local_ext_mat;
79  Tensor4 local_propmat_clearsky;
80  ao_gp_p[0]=gp_p;
81  ao_gp_lat[0]=gp_lat;
82  ao_gp_lon[0]=gp_lon;
83 
84  // Determine the pressure
85 
86  interpweights( itw_p, ao_gp_p );
87  itw2p( p_vec, p_grid, ao_gp_p, itw_p );
88 
89  // Determine the atmospheric temperature and species VMR
90  //
91  interp_atmfield_gp2itw( itw_field, 3, ao_gp_p, ao_gp_lat, ao_gp_lon );
92  //
93  interp_atmfield_by_itw( t_vec, 3, t_field, ao_gp_p, ao_gp_lat, ao_gp_lon,
94  itw_field );
95  //
96  for( Index is=0; is<ns; is++ )
97  {
98  interp_atmfield_by_itw( vmr_mat(is, joker), 3,
99  vmr_field(is, joker, joker, joker),
100  ao_gp_p, ao_gp_lat,
101  ao_gp_lon, itw_field );
102  }
103 
104 
105  temperature = t_vec[0];
106 
107  const Vector rtp_mag_dummy(3,0);
108  const Vector ppath_los_dummy;
109 
110  //calcualte absorption coefficient
111  propmat_clearsky_agendaExecute(ws, local_propmat_clearsky,
112  Vector(1, f_mono), rtp_mag_dummy,
113  ppath_los_dummy,p_vec[0],
114  temperature, vmr_mat(joker, 0),
115  propmat_clearsky_agenda);
116 
117  opt_prop_sum_propmat_clearsky(local_ext_mat, local_abs_vec, local_propmat_clearsky);
118 
119  ext_mat_mono=local_ext_mat(0, Range(joker), Range(joker));
120  abs_vec_mono=local_abs_vec(0,Range(joker));
121 }
122 
123 
124 
126 
136  MatrixView ext_mat_mono,
137  VectorView abs_vec_mono,
138  VectorView pnd_vec,
139  Numeric& temperature,
140  const Agenda& propmat_clearsky_agenda,
141  const Index stokes_dim,
142  const Numeric& f_mono,
143  const GridPos& gp_p,
144  const GridPos& gp_lat,
145  const GridPos& gp_lon,
146  ConstVectorView p_grid_cloud,
147  ConstTensor3View t_field_cloud,
148  ConstTensor4View vmr_field_cloud,
149  const Tensor4& pnd_field,
150  const ArrayOfSingleScatteringData& scat_data_array_mono,
151  const ArrayOfIndex& cloudbox_limits,
152  const Vector& rte_los,
153  const Verbosity& verbosity
154  )
155 
156 {
157  const Index ns = vmr_field_cloud.nbooks();
158  const Index N_pt = pnd_field.nbooks();
159  Matrix pnd_ppath(N_pt,1);
160  Vector t_ppath(1);
161  Vector p_ppath(1);//may not be efficient with unecessary vectors
162  Matrix itw_p(1,2);
163  ArrayOfGridPos ao_gp_p(1),ao_gp_lat(1),ao_gp_lon(1);
164  Matrix vmr_ppath(ns,1), itw_field;
165  Matrix ext_mat_part(stokes_dim, stokes_dim, 0.0);
166  Vector abs_vec_part(stokes_dim, 0.0);
167  Numeric scat_za,scat_aa;
168 
169  //local versions of workspace variables
170  Tensor4 local_propmat_clearsky;
171  Matrix local_abs_vec;
172  Tensor3 local_ext_mat;
173 
174  ao_gp_p[0]=gp_p;
175  ao_gp_lat[0]=gp_lat;
176  ao_gp_lon[0]=gp_lon;
177 
178 
179  cloud_atm_vars_by_gp(p_ppath,t_ppath,vmr_ppath,pnd_ppath,ao_gp_p,
180  ao_gp_lat,ao_gp_lon,cloudbox_limits,p_grid_cloud,
181  t_field_cloud, vmr_field_cloud,pnd_field);
182 
183  temperature = t_ppath[0];
184 
185  const Vector rtp_mag_dummy(3,0);
186  const Vector ppath_los_dummy;
187 
188  //rtp_vmr = vmr_ppath(joker,0);
189  propmat_clearsky_agendaExecute(ws, local_propmat_clearsky,
190  Vector(1, f_mono), rtp_mag_dummy,
191  ppath_los_dummy,p_ppath[0],
192  temperature,vmr_ppath(joker, 0),
193  propmat_clearsky_agenda);
194 
195  opt_prop_sum_propmat_clearsky(local_ext_mat, local_abs_vec, local_propmat_clearsky);
196 
197 
198  ext_mat_mono=local_ext_mat(0, Range(joker), Range(joker));
199  abs_vec_mono=local_abs_vec(0,Range(joker));
200  ext_mat_part=0.0;
201  abs_vec_part=0.0;
202  scat_za=180-rte_los[0];
203  scat_aa=rte_los[1]+180;
204  //Make sure scat_aa is between -180 and 180
205  if (scat_aa>180){scat_aa-=360;}
206  //
207  //opt_prop_part_agenda.execute( true );
208  //use pnd_ppath and ext_mat_spt to get extmat (and similar for abs_vec
209  pnd_vec=pnd_ppath(joker, 0);
210  opt_propCalc(ext_mat_part,abs_vec_part,scat_za,scat_aa,scat_data_array_mono,
211  stokes_dim, pnd_vec, temperature,verbosity);
212 
213  ext_mat_mono += ext_mat_part;
214  abs_vec_mono += abs_vec_part;
215 
216 }
217 
218 
219 
221 
244  VectorView pressure,
245  VectorView temperature,
246  MatrixView vmr,
247  MatrixView pnd,
248  const ArrayOfGridPos& gp_p,
249  const ArrayOfGridPos& gp_lat,
250  const ArrayOfGridPos& gp_lon,
251  const ArrayOfIndex& cloudbox_limits,
252  ConstVectorView p_grid_cloud,
253  ConstTensor3View t_field_cloud,
254  ConstTensor4View vmr_field_cloud,
255  ConstTensor4View pnd_field
256  )
257 {
258  Index np=gp_p.nelem();
259  assert(pressure.nelem()==np);
260  Index ns=vmr_field_cloud.nbooks();
261  Index N_pt=pnd_field.nbooks();
262  ArrayOfGridPos gp_p_cloud = gp_p;
263  ArrayOfGridPos gp_lat_cloud = gp_lat;
264  ArrayOfGridPos gp_lon_cloud = gp_lon;
265  Index atmosphere_dim=3;
266 
267  for (Index i = 0; i < np; i++ )
268  {
269  // Calculate grid positions inside the cloud.
270  gp_p_cloud[i].idx -= cloudbox_limits[0];
271  gp_lat_cloud[i].idx -= cloudbox_limits[2];
272  gp_lon_cloud[i].idx -= cloudbox_limits[4];
273  }
274 
275  const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
276  const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
277  const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
278  gridpos_upperend_check( gp_p_cloud[0], n1 );
279  gridpos_upperend_check( gp_p_cloud[np-1], n1 );
280  gridpos_upperend_check( gp_lat_cloud[0], n2 );
281  gridpos_upperend_check( gp_lat_cloud[np-1], n2 );
282  gridpos_upperend_check( gp_lon_cloud[0], n3 );
283  gridpos_upperend_check( gp_lon_cloud[np-1], n3 );
284 
285  // Determine the pressure at each propagation path point
286  Matrix itw_p(np,2);
287  //
288  //interpweights( itw_p, ppath.gp_p );
289  interpweights( itw_p, gp_p_cloud );
290  itw2p( pressure, p_grid_cloud, gp_p_cloud, itw_p );
291 
292  // Determine the atmospheric temperature and species VMR at
293  // each propagation path point
294  Matrix itw_field;
295  //
296  interp_atmfield_gp2itw( itw_field, atmosphere_dim,
297  gp_p_cloud, gp_lat_cloud, gp_lon_cloud );
298  //
299  interp_atmfield_by_itw( temperature, atmosphere_dim, t_field_cloud,
300  gp_p_cloud, gp_lat_cloud, gp_lon_cloud,
301  itw_field );
302  //
303  for( Index is=0; is<ns; is++ )
304  {
305  interp_atmfield_by_itw( vmr(is, joker), atmosphere_dim,
306  vmr_field_cloud(is, joker, joker, joker),
307  gp_p_cloud, gp_lat_cloud,
308  gp_lon_cloud, itw_field );
309  }
310 
311  //Determine the particle number density for every particle type at
312  // each propagation path point
313  for( Index ip=0; ip<N_pt; ip++ )
314  {
315  // if grid positions still outside the range the propagation path step
316  // must be outside the cloudbox and pnd is set to zero
317  interp_atmfield_by_itw( pnd(ip, joker), atmosphere_dim,
318  pnd_field(ip, joker, joker, joker),
319  gp_p_cloud, gp_lat_cloud,
320  gp_lon_cloud, itw_field );
321  }
322 }
323 
324 
325 
327 
337 void findZ11max(Vector& Z11maxvector,
338  const ArrayOfSingleScatteringData& scat_data_array_mono)
339 {
340  Index np=scat_data_array_mono.nelem();
341  Z11maxvector.resize(np);
342 
343  for(Index i = 0;i<np;i++)
344  {
345  switch(scat_data_array_mono[i].particle_type){
347  {
348  Z11maxvector[i]=max(scat_data_array_mono[i].pha_mat_data(0,joker,joker,0,0,0,0));
349  }
351  {
352  Z11maxvector[i]=max(scat_data_array_mono[i].pha_mat_data(0,joker,joker,0,joker,0,0));
353  }
354  default:
355  Z11maxvector[i]=max(scat_data_array_mono[i].pha_mat_data(0,joker,joker,joker,joker,joker,0));
356  }
357  }
358 }
359 
360 
361 
362 
364 
373 bool is_anyptype30(const ArrayOfSingleScatteringData& scat_data_array_mono)
374 {
375  Index np=scat_data_array_mono.nelem();
376  bool anyptype30=false;
377  Index i=0;
378  while(i < np && anyptype30==false)
379  {
380  if(scat_data_array_mono[i].particle_type==PARTICLE_TYPE_HORIZ_AL)
381  {
382  anyptype30=true;
383  }
384  i+=1;
385  }
386  return anyptype30;
387 }
388 
389 
390 
392 
412  Workspace& ws,
413  MatrixView evol_op,
414  Vector& abs_vec_mono,
415  Numeric& temperature,
416  MatrixView ext_mat_mono,
417  Rng& rng,
418  Vector& rte_pos,
419  Vector& rte_los,
420  Vector& pnd_vec,
421  Numeric& g,
422  Ppath& ppath_step,
423  Index& termination_flag,
424  bool& inside_cloud,
425  const Agenda& ppath_step_agenda,
426  const Numeric& ppath_lraytrace,
427  const Agenda& propmat_clearsky_agenda,
428  const Index stokes_dim,
429  const Numeric& f_mono,
430  const Vector& p_grid,
431  const Vector& lat_grid,
432  const Vector& lon_grid,
433  const Tensor3& z_field,
434  const Vector& refellipsoid,
435  const Matrix& z_surface,
436  const Tensor3& t_field,
437  const Tensor4& vmr_field,
438  const ArrayOfIndex& cloudbox_limits,
439  const Tensor4& pnd_field,
440  const ArrayOfSingleScatteringData& scat_data_array_mono,
441  const Verbosity& verbosity )
442 {
443  ArrayOfMatrix evol_opArray(2);
444  ArrayOfMatrix ext_matArray(2);
445  ArrayOfVector abs_vecArray(2);
446  ArrayOfVector pnd_vecArray(2);
447  Matrix ext_mat(stokes_dim,stokes_dim);
448  Matrix incT(stokes_dim,stokes_dim,0.0);
449  Vector tArray(2);
450  Vector f_grid(1,f_mono); // Vector version of f_mono
451  Matrix T(stokes_dim,stokes_dim);
452  Numeric k;
453  Numeric ds, dl=-999;
454  Index istep = 0; // Counter for number of steps
455  Matrix old_evol_op(stokes_dim,stokes_dim);
456 
457  CREATE_OUT0;
458 
459  //at the start of the path the evolution operator is the identity matrix
460  id_mat(evol_op);
461  evol_opArray[1]=evol_op;
462 
463  // range defining cloudbox
464  Range p_range( cloudbox_limits[0],
465  cloudbox_limits[1]-cloudbox_limits[0]+1);
466  Range lat_range( cloudbox_limits[2],
467  cloudbox_limits[3]-cloudbox_limits[2]+1 );
468  Range lon_range( cloudbox_limits[4],
469  cloudbox_limits[5]-cloudbox_limits[4]+1 );
470 
471  //initialise Ppath with ppath_start_stepping
472  ppath_start_stepping( ppath_step, 3, p_grid, lat_grid,
473  lon_grid, z_field, refellipsoid, z_surface,
474  0, cloudbox_limits, false,
475  rte_pos, rte_los, verbosity );
476 
477  // Index in ppath_step of end point considered presently
478  Index ip = 0;
479 
480  // Is point ip inside the cloudbox?
481  inside_cloud = is_gp_inside_cloudbox( ppath_step.gp_p[ip],
482  ppath_step.gp_lat[ip],
483  ppath_step.gp_lon[ip],
484  cloudbox_limits, 0, 3 );
485 
486  // Determine radiative properties at point
487  if( inside_cloud )
488  {
489  cloudy_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, pnd_vec,
490  temperature, propmat_clearsky_agenda,
491  stokes_dim, f_mono, ppath_step.gp_p[0],
492  ppath_step.gp_lat[0], ppath_step.gp_lon[0],
493  p_grid[p_range],
494  t_field(p_range,lat_range,lon_range),
495  vmr_field(joker,p_range,lat_range,lon_range),
496  pnd_field, scat_data_array_mono, cloudbox_limits,
497  ppath_step.los(0,joker), verbosity );
498  }
499  else
500  {
501  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
502  propmat_clearsky_agenda, f_mono,
503  ppath_step.gp_p[0], ppath_step.gp_lat[0],
504  ppath_step.gp_lon[0], p_grid, t_field, vmr_field );
505  pnd_vec = 0.0;
506  }
507 
508  // Move the data to end point containers
509  ext_matArray[1] = ext_mat_mono;
510  abs_vecArray[1] = abs_vec_mono;
511  tArray[1] = temperature;
512  pnd_vecArray[1] = pnd_vec;
513 
514  //draw random number to determine end point
515  Numeric r = rng.draw();
516 
517  termination_flag=0;
518 
519  while( (evol_op(0,0)>r) && (!termination_flag) )
520  {
521  istep++;
522 
523  if( istep > 5000 )
524  {
525  throw runtime_error( "5000 path points have been reached. "
526  "Is this an infinite loop?" );
527  }
528 
529  evol_opArray[0] = evol_opArray[1];
530  ext_matArray[0] = ext_matArray[1];
531  abs_vecArray[0] = abs_vecArray[1];
532  tArray[0] = tArray[1];
533  pnd_vecArray[0] = pnd_vecArray[1];
534 
535  // Shall new ppath_step be calculated?
536  if( ip == ppath_step.np-1 )
537  {
538  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field,
539  z_field, vmr_field, f_grid,
540  ppath_step_agenda );
541  //Print( ppath_step, 0, verbosity );
542  ip = 1;
543 
544  inside_cloud = is_gp_inside_cloudbox( ppath_step.gp_p[ip],
545  ppath_step.gp_lat[ip],
546  ppath_step.gp_lon[ip],
547  cloudbox_limits, 0, 3 );
548  }
549  else
550  { ip++; }
551 
552  dl = ppath_step.lstep[ip-1];
553 
554  if( inside_cloud )
555  {
556  cloudy_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, pnd_vec,
557  temperature, propmat_clearsky_agenda,
558  stokes_dim, f_mono, ppath_step.gp_p[ip],
559  ppath_step.gp_lat[ip], ppath_step.gp_lon[ip],
560  p_grid[p_range],
561  t_field(p_range,lat_range,lon_range),
562  vmr_field(joker,p_range,lat_range,lon_range),
563  pnd_field, scat_data_array_mono, cloudbox_limits,
564  ppath_step.los(ip,joker), verbosity );
565  }
566  else
567  {
568  clear_rt_vars_at_gp( ws, ext_mat_mono, abs_vec_mono, temperature,
569  propmat_clearsky_agenda, f_mono,
570  ppath_step.gp_p[ip], ppath_step.gp_lat[ip],
571  ppath_step.gp_lon[ip], p_grid, t_field,
572  vmr_field );
573  pnd_vec = 0.0;
574  }
575 
576  ext_matArray[1] = ext_mat_mono;
577  abs_vecArray[1] = abs_vec_mono;
578  tArray[1] = temperature;
579  pnd_vecArray[1] = pnd_vec;
580  ext_mat = ext_matArray[1];
581  ext_mat += ext_matArray[0]; // Factor 2 fixed by using dl/2
582  //
583  {
584  Index extmat_case = 0;
585  ext2trans( incT, extmat_case, ext_mat, dl/2 );
586  }
587  //
588  mult( evol_op, evol_opArray[0], incT );
589  evol_opArray[1] = evol_op;
590 
591  if( evol_op(0,0)>r )
592  {
593  // Check whether hit ground or space.
594  // path_step_agenda just detects surface intersections, and
595  // if TOA is reached requires a special check.
596  // But we are already ready if evol_op<=r
597  if( ip == ppath_step.np - 1 )
598  {
599  if( ppath_what_background(ppath_step) )
600  { termination_flag = 2; } //we have hit the surface
601  else if( fractional_gp(ppath_step.gp_p[ip]) >=
602  (Numeric)(p_grid.nelem() - 1)-1e-3 )
603  { termination_flag = 1; } //we are at TOA
604  }
605  }
606  } // while
607 
608 
609  if( termination_flag )
610  { //we must have reached the cloudbox boundary
611  rte_pos = ppath_step.pos(ip,joker);
612  rte_los = ppath_step.los(ip,joker);
613  g = evol_op(0,0);
614  }
615  else
616  {
617  //find position...and evol_op..and everything else required at the new
618  //scattering/emission point
619  // GH 2011-09-14:
620  // log(incT(0,0)) = log(exp(opt_depth_mat(0, 0))) = opt_depth_mat(0, 0)
621  // Avoid loss of precision, use opt_depth_mat directly
622  //k=-log(incT(0,0))/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
623  // PE 2013-05-17, Now the above comes directly from ext_mat:
624  k = ext_mat(0,0) / 2; // Factor 2 as sum of end point values
625  ds = log( evol_opArray[0](0,0) / r ) / k;
626  g = k*r;
627  Vector x(2,0.0);
628  //interpolate atmospheric variables required later
629  ArrayOfGridPos gp(1);
630  x[1] = dl;
631  Vector itw(2);
632 
633  gridpos( gp, x, ds );
634  assert( gp[0].idx == 0 );
635  interpweights( itw, gp[0] );
636  interp( ext_mat_mono, itw, ext_matArray, gp[0] );
637  ext_mat = ext_mat_mono;
638  ext_mat += ext_matArray[gp[0].idx];
639  //
640  {
641  Index extmat_case = 0;
642  ext2trans( incT, extmat_case, ext_mat, ds/2 );
643  }
644  //
645  mult( evol_op, evol_opArray[gp[0].idx], incT );
646  interp( abs_vec_mono, itw, abs_vecArray,gp[0] );
647  temperature = interp( itw, tArray, gp[0] );
648  interp( pnd_vec, itw, pnd_vecArray, gp[0] );
649  for( Index i=0; i<2; i++ )
650  {
651  rte_pos[i] = interp( itw, ppath_step.pos(Range(ip-1,2),i), gp[0] );
652  rte_los[i] = interp( itw, ppath_step.los(Range(ip-1,2),i), gp[0] );
653  }
654  rte_pos[2] = interp( itw, ppath_step.pos(Range(ip-1,2),2), gp[0] );
655  }
656 
657  assert(isfinite(g));
658 
659  // A dirty trick to avoid copying ppath_step
660  const Index np = ip+1;
661  ppath_step.np = np;
662 
663 }
664 
665 
666 
668 
685  MatrixView ext_mat_mono,
686  VectorView abs_vec_mono,
687  const Numeric za,
688  const Numeric aa,
689  const ArrayOfSingleScatteringData& scat_data_array_mono,
690  const Index stokes_dim,
691  ConstVectorView pnd_vec,
692  const Numeric rtp_temperature,
693  const Verbosity& verbosity
694  )
695 {
696  assert( stokes_dim>=1 && stokes_dim<=4 );
697  assert( ext_mat_mono.nrows() == stokes_dim );
698  assert( ext_mat_mono.ncols() == stokes_dim );
699  assert( abs_vec_mono.nelem() == stokes_dim );
700 
701  const Index N_pt = scat_data_array_mono.nelem();
702 
703  Matrix ext_mat_mono_spt(stokes_dim,stokes_dim);
704  Vector abs_vec_mono_spt(stokes_dim);
705 
706  ext_mat_mono=0.0;
707  abs_vec_mono=0.0;
708 
709  // Loop over the included particle_types
710  for (Index i_pt = 0; i_pt < N_pt; i_pt++)
711  {
712  if (pnd_vec[i_pt]>0)
713  {
714  Index nT=scat_data_array_mono[i_pt].T_grid.nelem();
715  if( nT > 1 )
716  {
717  //set extrapolfax explicitly. should correspond to the one assumed
718  //by gridpos call in opt_propExtract.
719  Numeric extrapolfac=0.5;
720  Numeric lowlim = scat_data_array_mono[i_pt].T_grid[0]-
721  extrapolfac*(scat_data_array_mono[i_pt].T_grid[1]
722  -scat_data_array_mono[i_pt].T_grid[0]);
723  Numeric uplim = scat_data_array_mono[i_pt].T_grid[nT-1]+
724  extrapolfac*(scat_data_array_mono[i_pt].T_grid[nT-1]
725  -scat_data_array_mono[i_pt].T_grid[nT-2]);
726 
727  if( rtp_temperature<lowlim || rtp_temperature>uplim )
728  {
729  ostringstream os;
730  os << "Atmospheric temperature (" << rtp_temperature
731  << "K) out of valid temperature\n"
732  << "range of particle optical properties ("
733  << lowlim << "-" << uplim
734  << "K) of particle #" << i_pt << ".\n";
735  throw runtime_error( os.str() );
736  }
737  }
738  opt_propExtract( ext_mat_mono_spt, abs_vec_mono_spt,
739  scat_data_array_mono[i_pt], za, aa,
740  rtp_temperature, stokes_dim, verbosity);
741 
742  ext_mat_mono_spt *= pnd_vec[i_pt];
743  abs_vec_mono_spt *= pnd_vec[i_pt];
744  ext_mat_mono += ext_mat_mono_spt;
745  abs_vec_mono += abs_vec_mono_spt;
746  }
747  }
748 }
749 
750 
752  MatrixView ext_mat_mono_spt,
753  VectorView abs_vec_mono_spt,
754  const SingleScatteringData& scat_data,
755  const Numeric za,
756  const Numeric aa _U_, // avoid warning until we use particle_type=10
757  const Numeric rtp_temperature,
758  const Index stokes_dim,
759  const Verbosity& verbosity
760  )
761 {
762  // Temperature grid position
763  GridPos t_gp;
764  if( scat_data.T_grid.nelem() > 1)
765  { gridpos( t_gp, scat_data.T_grid, rtp_temperature ); }
766 
767 
768  switch (scat_data.particle_type){
769 
771  {
772  // This is only included to remove warnings about unused variables
773  // during compilation
774  CREATE_OUT0;
775  out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n";
776  break;
777  }
779  {
780  assert (scat_data.ext_mat_data.ncols() == 1);
781 
782  Vector itw(2);
783  //interpolate over temperature
784  if( scat_data.T_grid.nelem() > 1)
785  {
786  interpweights(itw, t_gp);
787  }
788  // In the case of randomly oriented particles the extinction matrix is
789  // diagonal. The value of each element of the diagonal is the
790  // extinction cross section, which is stored in the database.
791 
792  ext_mat_mono_spt = 0.;
793  abs_vec_mono_spt = 0;
794 
795  if( scat_data.T_grid.nelem() == 1)
796  {
797  ext_mat_mono_spt(0,0) = scat_data.ext_mat_data(0,0,0,0,0);
798  abs_vec_mono_spt[0] = scat_data.abs_vec_data(0,0,0,0,0);
799  }
800  // Temperature interpolation
801  else
802  {
803  ext_mat_mono_spt(0,0) = interp(itw,scat_data.ext_mat_data(0,joker,0,0,0),t_gp);
804  abs_vec_mono_spt[0] = interp(itw,scat_data.abs_vec_data(0,joker,0,0,0),t_gp);
805  }
806 
807  if( stokes_dim == 1 ){
808  break;
809  }
810 
811  ext_mat_mono_spt(1,1) = ext_mat_mono_spt(0,0);
812 
813  if( stokes_dim == 2 ){
814  break;
815  }
816 
817  ext_mat_mono_spt(2,2) = ext_mat_mono_spt(0,0);
818 
819  if( stokes_dim == 3 ){
820  break;
821  }
822 
823  ext_mat_mono_spt(3,3) = ext_mat_mono_spt(0,0);
824  break;
825  }
826 
827  case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis 9/12/03
828  {
829  assert (scat_data.ext_mat_data.ncols() == 3);
830 
831  // In the case of horizontally oriented particles the extinction matrix
832  // has only 3 independent non-zero elements ext_mat_monojj, K12=K21, and
833  // K34=-K43. These values are dependent on the zenith angle of
834  // propagation. The data storage format also makes use of the fact that
835  // in this case K(za_sca)=K(180-za_sca).
836 
837  // 1st interpolate data by za_sca
838  GridPos za_gp;
839  Vector itw(4);
840  Numeric Kjj;
841  Numeric K12;
842  Numeric K34;
843 
844  if( scat_data.T_grid.nelem() == 1)
845  {
846  ostringstream os;
847  os << "Given optical property data are for constant temperature "
848  << "only.\nMC with p30 requires temperature-dependent optical "
849  << "property data\n";
850  throw runtime_error( os.str() );
851  }
852 
853  // This fix for za=90 copied from ext_matTransform in optproperties.cc
854  // (121113, PE).
855  ConstVectorView this_za_datagrid =
856  scat_data.za_grid[ Range( 0, scat_data.ext_mat_data.npages() ) ];
857  if( za > 90 )
858  { gridpos( za_gp, this_za_datagrid, 180-za ); }
859  else
860  { gridpos( za_gp, this_za_datagrid, za ); }
861 
862  interpweights( itw, t_gp, za_gp );
863 
864  ext_mat_mono_spt = 0.0;
865  abs_vec_mono_spt = 0.0;
866  Kjj = interp(itw,scat_data.ext_mat_data(0,joker,joker,0,0),t_gp,za_gp);
867  ext_mat_mono_spt(0,0) = Kjj;
868  abs_vec_mono_spt[0] = interp( itw,
869  scat_data.abs_vec_data(0,joker,joker,0,0), t_gp,za_gp );
870 
871  if( stokes_dim == 1 )
872  { break; }
873 
874  K12=interp(itw,scat_data.ext_mat_data(0,joker,joker,0,1),t_gp,za_gp);
875  ext_mat_mono_spt(1,1)=Kjj;
876  ext_mat_mono_spt(0,1)=K12;
877  ext_mat_mono_spt(1,0)=K12;
878  abs_vec_mono_spt[1] = interp(itw,scat_data.abs_vec_data(0,joker,joker,0,1),
879  t_gp,za_gp);
880 
881  if( stokes_dim == 2 ){
882  break;
883  }
884 
885  ext_mat_mono_spt(2,2)=Kjj;
886 
887  if( stokes_dim == 3 ){
888  break;
889  }
890 
891  K34=interp(itw,scat_data.ext_mat_data(0,joker,joker,0,2),t_gp,za_gp);
892  ext_mat_mono_spt(2,3)=K34;
893  ext_mat_mono_spt(3,2)=-K34;
894  ext_mat_mono_spt(3,3)=Kjj;
895  break;
896 
897  }
898  default:
899  {
900  CREATE_OUT0;
901  out0 << "Not all particle type cases are implemented\n";
902  }
903 
904  }
905 
906 
907 }
908 
909 
911 
929  MatrixView Z,
930  const Numeric za_sca,
931  const Numeric aa_sca,
932  const Numeric za_inc,
933  const Numeric aa_inc,
934  const ArrayOfSingleScatteringData& scat_data_array_mono,
935  const Index stokes_dim,
936  ConstVectorView pnd_vec,
937  const Numeric rtp_temperature,
938  const Verbosity& verbosity
939  )
940 {
941  Index N_pt=pnd_vec.nelem();
942 
943  assert(aa_inc>=-180 && aa_inc<=180);
944  assert(aa_sca>=-180 && aa_sca<=180);
945 
946  Matrix Z_spt(stokes_dim, stokes_dim, 0.);
947  Z=0.0;
948  // this is a loop over the different particle types
949  for (Index i_pt = 0; i_pt < N_pt; i_pt++)
950  {
951  if (pnd_vec[i_pt]>0)
952  {
953  pha_mat_singleExtract(Z_spt,scat_data_array_mono[i_pt],za_sca,aa_sca,za_inc,
954  aa_inc,rtp_temperature,stokes_dim,verbosity);
955  Z_spt*=pnd_vec[i_pt];
956  Z+=Z_spt;
957  }
958  }
959 }
960 
961 
962 
964 
982  MatrixView Z_spt,
983  const SingleScatteringData& scat_data,
984  const Numeric za_sca,
985  const Numeric aa_sca,
986  const Numeric za_inc,
987  const Numeric aa_inc,
988  const Numeric rtp_temperature,
989  const Index stokes_dim,
990  const Verbosity& verbosity
991  )
992 {
993  switch (scat_data.particle_type){
994 
996  {
997  // to remove warnings during compilation.
998  CREATE_OUT0;
999  out0 << "Case PARTICLE_TYPE_GENERAL not yet implemented. \n";
1000  break;
1001  }
1003  {
1004  // Calculate the scattering and interpolate the data on the scattering
1005  // angle:
1006 
1007  Vector pha_mat_int(6);
1008  Numeric theta_rad;
1009 
1010  // Interpolation of the data on the scattering angle:
1011  interp_scat_angle_temperature(pha_mat_int, theta_rad,
1012  scat_data, za_sca, aa_sca,
1013  za_inc, aa_inc, rtp_temperature);
1014 
1015  // Caclulate the phase matrix in the laboratory frame:
1016  pha_mat_labCalc(Z_spt, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1017 
1018  break;
1019  }
1020 
1021  case PARTICLE_TYPE_HORIZ_AL://Added by Cory Davis
1022  //Data is already stored in the laboratory frame, but it is compressed
1023  //a little. Details elsewhere
1024  {
1025  assert (scat_data.pha_mat_data.ncols()==16);
1026  Numeric delta_aa=aa_sca-aa_inc+(aa_sca-aa_inc<-180)*360-
1027  (aa_sca-aa_inc>180)*360;//delta_aa corresponds to the "pages"
1028  //dimension of pha_mat_data
1029  GridPos t_gp;
1030  GridPos za_sca_gp;
1031  GridPos delta_aa_gp;
1032  GridPos za_inc_gp;
1033  Vector itw(16);
1034 
1035  if( scat_data.T_grid.nelem() == 1)
1036  {
1037  ostringstream os;
1038  os << "Given optical property data is for constant temperature only.\n"
1039  "MC with p30 requires temperature-dependent optical property data\n";
1040  throw runtime_error( os.str() );
1041  }
1042 
1043  gridpos(t_gp, scat_data.T_grid, rtp_temperature);
1044  gridpos(delta_aa_gp,scat_data.aa_grid,abs(delta_aa));
1045  if (za_inc>90)
1046  {
1047  gridpos(za_inc_gp,scat_data.za_grid,180-za_inc);
1048  gridpos(za_sca_gp,scat_data.za_grid,180-za_sca);
1049  }
1050  else
1051  {
1052  gridpos(za_inc_gp,scat_data.za_grid,za_inc);
1053  gridpos(za_sca_gp,scat_data.za_grid,za_sca);
1054  }
1055 
1056  interpweights(itw,t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1057 
1058  Z_spt(0,0)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1059  Range(joker),0,0),
1060  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1061  if( stokes_dim == 1 ){
1062  break;
1063  }
1064  Z_spt(0,1)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1065  Range(joker),0,1),
1066  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1067  Z_spt(1,0)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1068  Range(joker),0,4),
1069  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1070  Z_spt(1,1)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1071  Range(joker),0,5),
1072  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1073  if( stokes_dim == 2 ){
1074  break;
1075  }
1076  if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
1077  {
1078  Z_spt(0,2)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1079  Range(joker),0,2),
1080  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1081  Z_spt(1,2)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1082  Range(joker),0,6),
1083  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1084  Z_spt(2,0)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1085  Range(joker),0,8),
1086  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1087  Z_spt(2,1)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1088  Range(joker),0,9),
1089  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1090  }
1091  else
1092  {
1093  Z_spt(0,2)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1094  Range(joker),0,2),
1095  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1096  Z_spt(1,2)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1097  Range(joker),0,6),
1098  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1099  Z_spt(2,0)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1100  Range(joker),0,8),
1101  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1102  Z_spt(2,1)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1103  Range(joker),0,9),
1104  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1105  }
1106  Z_spt(2,2)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1107  Range(joker),0,10),
1108  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1109  if( stokes_dim == 3 ){
1110  break;
1111  }
1112  if ((za_inc<=90 && delta_aa>=0)||(za_inc>90 && delta_aa<0))
1113  {
1114  Z_spt(0,3)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1115  Range(joker),0,3),
1116  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1117  Z_spt(1,3)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1118  Range(joker),0,7),
1119  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1120  Z_spt(3,0)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1121  Range(joker),0,12),
1122  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1123  Z_spt(3,1)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1124  Range(joker),0,13),
1125  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1126  }
1127  else
1128  {
1129  Z_spt(0,3)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1130  Range(joker),0,3),
1131  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1132  Z_spt(1,3)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1133  Range(joker),0,7),
1134  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1135  Z_spt(3,0)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1136  Range(joker),0,12),
1137  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1138  Z_spt(3,1)=-interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1139  Range(joker),0,13),
1140  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1141  }
1142  Z_spt(2,3)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1143  Range(joker),0,11),
1144  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1145  Z_spt(3,2)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1146  Range(joker),0,14),
1147  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1148  Z_spt(3,3)=interp(itw,scat_data.pha_mat_data(0,joker,Range(joker),Range(joker),
1149  Range(joker),0,15),
1150  t_gp,za_sca_gp,delta_aa_gp,za_inc_gp);
1151  break;
1152 
1153  }
1154  default:
1155  CREATE_OUT0;
1156  out0 << "Not all particle type cases are implemented\n";
1157 
1158  }
1159 }
1160 
1161 
1162 
1163 
1165 
1188  VectorView new_rte_los,
1189  Numeric& g_los_csc_theta,
1190  MatrixView Z,
1191  Rng& rng,
1192  ConstVectorView rte_los,
1193  const ArrayOfSingleScatteringData& scat_data_array_mono,
1194  const Index stokes_dim,
1195  ConstVectorView pnd_vec,
1196  const bool anyptype30,
1197  ConstVectorView Z11maxvector,
1198  const Numeric Csca,
1199  const Numeric rtp_temperature,
1200  const Verbosity& verbosity
1201  )
1202 {
1203  Numeric Z11max=0;
1204  bool tryagain=true;
1205  Numeric aa_scat = (rte_los[1]>=0) ?-180+rte_los[1]:180+rte_los[1];
1206 
1207  // Rejection method http://en.wikipedia.org/wiki/Rejection_sampling
1208  if(anyptype30)
1209  {
1210  Index np=pnd_vec.nelem();
1211  assert(scat_data_array_mono.nelem()==np);
1212  for(Index i=0;i<np;i++)
1213  {
1214  Z11max+=Z11maxvector[i]*pnd_vec[i];
1215  }
1216  }
1217  else
1218  {
1219  Matrix dummyZ(stokes_dim,stokes_dim);
1220  //The following is based on the assumption that the maximum value of the
1221  //phase matrix for a given scattered direction is for forward scattering
1222  pha_mat_singleCalc(dummyZ,180-rte_los[0],aa_scat,180-rte_los[0],
1223  aa_scat,scat_data_array_mono,stokes_dim,pnd_vec,rtp_temperature,
1224  verbosity);
1225  Z11max=dummyZ(0,0);
1226  }
1228  while(tryagain)
1229  {
1230  new_rte_los[1] = rng.draw()*360-180;
1231  new_rte_los[0] = acos(1-2*rng.draw())*RAD2DEG;
1232  //Calculate Phase matrix////////////////////////////////
1233  Numeric aa_inc= (new_rte_los[1]>=0) ?
1234  -180+new_rte_los[1]:180+new_rte_los[1];
1235 
1236  pha_mat_singleCalc(Z,180-rte_los[0],aa_scat,180-new_rte_los[0],
1237  aa_inc,scat_data_array_mono,stokes_dim,pnd_vec,rtp_temperature,
1238  verbosity);
1239 
1240  if (rng.draw()<=Z(0,0)/Z11max)//then new los is accepted
1241  {
1242  tryagain=false;
1243  }
1244  }
1245  g_los_csc_theta =Z(0,0)/Csca;
1246 }
Index npages() const
Returns the number of pages.
Definition: matpackV.cc:47
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
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
The VectorView class.
Definition: matpackI.h:372
void pha_mat_labCalc(MatrixView pha_mat_lab, ConstVectorView pha_mat_int, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &theta_rad)
Calculate phase matrix in laboratory coordinate system.
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
Interpolation classes and functions created for use within Monte Carlo scattering simulations...
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
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 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
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 Numeric &ppath_lraytrace, 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 Verbosity &verbosity)
mcPathTraceGeneral
Definition: montecarlo.cc:411
Vector lstep
Definition: ppath.h:70
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Matrix pos
Definition: ppath.h:67
Structure which describes the single scattering properties of a particle or a particle distribution...
Definition: optproperties.h:84
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
void pha_mat_singleExtract(MatrixView Z_spt, const SingleScatteringData &scat_data, const Numeric za_sca, const Numeric aa_sca, const Numeric za_inc, const Numeric aa_inc, const Numeric rtp_temperature, const Index stokes_dim, const Verbosity &verbosity)
Extract the phase matrix from a monochromatic SingleScatteringData object.
Definition: montecarlo.cc:981
Structure to store a grid position.
Definition: interpolation.h:74
A constant view of a Tensor4.
Definition: matpackIV.h:141
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
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
void opt_propExtract(MatrixView ext_mat_mono_spt, VectorView abs_vec_mono_spt, const SingleScatteringData &scat_data, const Numeric za, const Numeric aa, const Numeric rtp_temperature, const Index stokes_dim, const Verbosity &verbosity)
Definition: montecarlo.cc:751
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
#define max(a, b)
Definition: continua.cc:20461
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void propmat_clearsky_agendaExecute(Workspace &ws, Tensor4 &propmat_clearsky, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:13039
#define abs(x)
Definition: continua.cc:20458
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
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 interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
interp_atmfield_gp2itw
double draw()
Definition: rng.cc:120
void interp_scat_angle_temperature(VectorView pha_mat_int, Numeric &theta_rad, const SingleScatteringData &scat_data, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &rtp_temperature)
Definition: mc_interp.cc:183
Definition: rng.h:569
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define ns
Definition: continua.cc:21931
A constant view of a Tensor3.
Definition: matpackIII.h:139
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
A constant view of a Vector.
Definition: matpackI.h:292
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
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
#define _U_
Definition: config.h:167
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:59
bool is_anyptype30(const ArrayOfSingleScatteringData &scat_data_array_mono)
is_anyptype30
Definition: montecarlo.cc:373
#define CREATE_OUT0
Definition: messages.h:211
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:68
ParticleType particle_type
Definition: optproperties.h:85
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.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
ArrayOfGridPos gp_p
Definition: ppath.h:76
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
interp_atmfield_by_itw
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Definition: rte.cc:903
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
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
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
itw2p
void opt_prop_sum_propmat_clearsky(Tensor3 &ext_mat, Matrix &abs_vec, const Tensor4 propmat_clearsky)
Get optical properties from propmat_clearsky.