ARTS  2.4.0(git:4fb77825)
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 
26 /*===========================================================================
27  === External declarations
28  ===========================================================================*/
29 
30 #include <cfloat>
31 #include <sstream>
32 
33 #include "auto_md.h"
34 #include "geodetic.h"
35 #include "mc_interp.h"
36 #include "montecarlo.h"
37 
38 extern const Numeric SPEED_OF_LIGHT;
39 
40 // Some root-finding helper functions (for MCRadar) that don't need
41 // visibility outside this source file
42 //
43 //
44 
59  const Numeric& Q,
60  const Numeric& kI,
61  const Numeric& kQ,
62  const Numeric& s) {
63  Numeric fs;
64 
65  fs = exp(-kI * s) * (I * cosh(kQ * s) + Q * sinh(kQ * s));
66 
67  return fs;
68 }
69 
133  const Numeric& a,
134  const Numeric& b,
135  const Numeric& t,
136  const Numeric& rn,
137  const Numeric& I,
138  const Numeric& Q,
139  const Numeric& KI,
140  const Numeric& KQ) {
141  Numeric c;
142  Numeric d;
143  Numeric e;
144  Numeric fa;
145  Numeric fb;
146  Numeric fc;
147  Numeric m;
148  Numeric macheps;
149  Numeric p;
150  Numeric q;
151  Numeric r;
152  Numeric s;
153  Numeric sa;
154  Numeric tol;
155  //
156  // Make local copies of A and B.
157  //
158  sa = a;
159  sb = b;
160  fa = ext_I(I, Q, KI, KQ, sa); // - rn;
161  fa -= rn;
162  fb = ext_I(I, Q, KI, KQ, sb); // - rn;
163  fb -= rn;
164 
165  c = sa;
166  fc = fa;
167  e = sb - sa;
168  d = e;
169 
170  macheps = DBL_EPSILON;
171 
172  for (;;) {
173  if (abs(fc) < abs(fb)) {
174  sa = sb;
175  sb = c;
176  c = sa;
177  fa = fb;
178  fb = fc;
179  fc = fa;
180  }
181 
182  tol = 2.0 * macheps * abs(sb) + t;
183  m = 0.5 * (c - sb);
184 
185  if (abs(m) <= tol || fb == 0.0) {
186  break;
187  }
188 
189  if (abs(e) < tol || abs(fa) <= abs(fb)) {
190  e = m;
191  d = e;
192  } else {
193  s = fb / fa;
194 
195  if (sa == c) {
196  p = 2.0 * m * s;
197  q = 1.0 - s;
198  } else {
199  q = fa / fc;
200  r = fb / fc;
201  p = s * (2.0 * m * q * (q - r) - (sb - sa) * (r - 1.0));
202  q = (q - 1.0) * (r - 1.0) * (s - 1.0);
203  }
204 
205  if (0.0 < p) {
206  q = -q;
207  } else {
208  p = -p;
209  }
210 
211  s = e;
212  e = d;
213 
214  if (2.0 * p < 3.0 * m * q - abs(tol * q) && p < abs(0.5 * s * q)) {
215  d = p / q;
216  } else {
217  e = m;
218  d = e;
219  }
220  }
221  sa = sb;
222  fa = fb;
223 
224  if (tol < abs(d)) {
225  sb = sb + d;
226  } else if (0.0 < m) {
227  sb = sb + tol;
228  } else {
229  sb = sb - tol;
230  }
231 
232  fb = ext_I(I, Q, KI, KQ, sb);
233  fb -= rn;
234 
235  if ((0.0 < fb && 0.0 < fc) || (fb <= 0.0 && fc <= 0.0)) {
236  c = sa;
237  fc = fa;
238  e = sb - sa;
239  d = e;
240  }
241  }
242 }
243 
245  MatrixView ext_mat_mono,
246  VectorView abs_vec_mono,
247  Numeric& temperature,
249  const Numeric& f_mono,
250  const GridPos& gp_p,
251  const GridPos& gp_lat,
252  const GridPos& gp_lon,
256  const Index ns = vmr_field.nbooks();
257  Vector t_vec(1); //vectors are required by interp_atmfield_gp2itw etc.
258  Vector p_vec(1); //may not be efficient with unecessary vectors
259  Matrix itw_p(1, 2);
260  ArrayOfGridPos ao_gp_p(1), ao_gp_lat(1), ao_gp_lon(1);
261  Matrix vmr_mat(ns, 1), itw_field;
262 
263  //local versions of workspace variables
264  StokesVector local_abs_vec;
265  ArrayOfStokesVector local_nlte_source_dummy;
266  PropagationMatrix local_ext_mat;
267  ArrayOfPropagationMatrix local_propmat_clearsky;
269  local_partial_dummy; // This is right since there should be only clearsky partials
270  ArrayOfStokesVector local_dnlte_dx_source_dummy, local_nlte_dsource_dx_dummy;
271  ao_gp_p[0] = gp_p;
272  ao_gp_lat[0] = gp_lat;
273  ao_gp_lon[0] = gp_lon;
274 
275  // Determine the pressure
276  interpweights(itw_p, ao_gp_p);
277  itw2p(p_vec, p_grid, ao_gp_p, itw_p);
278 
279  // Determine the atmospheric temperature and species VMR
280  //
281  interp_atmfield_gp2itw(itw_field, 3, ao_gp_p, ao_gp_lat, ao_gp_lon);
282  //
284  t_vec, 3, t_field, ao_gp_p, ao_gp_lat, ao_gp_lon, itw_field);
285  //
286  for (Index is = 0; is < ns; is++) {
287  interp_atmfield_by_itw(vmr_mat(is, joker),
288  3,
289  vmr_field(is, joker, joker, joker),
290  ao_gp_p,
291  ao_gp_lat,
292  ao_gp_lon,
293  itw_field);
294  }
295 
296  temperature = t_vec[0];
297 
298  const Vector rtp_mag_dummy(3, 0);
299  const Vector ppath_los_dummy;
300  const EnergyLevelMap nlte_dummy;
301 
302  //calcualte absorption coefficient
304  local_propmat_clearsky,
305  local_nlte_source_dummy,
306  local_partial_dummy,
307  local_dnlte_dx_source_dummy,
308  local_nlte_dsource_dx_dummy,
310  Vector(1, f_mono),
311  rtp_mag_dummy,
312  ppath_los_dummy,
313  p_vec[0],
314  temperature,
315  nlte_dummy,
316  vmr_mat(joker, 0),
318 
320  local_ext_mat, local_abs_vec, local_propmat_clearsky);
321 
322  local_ext_mat.MatrixAtPosition(ext_mat_mono);
323  abs_vec_mono = local_abs_vec.VectorAtPosition();
324 }
325 
327  MatrixView ext_mat_mono,
328  VectorView abs_vec_mono,
329  VectorView pnd_vec,
330  Numeric& temperature,
332  const Index stokes_dim,
333  const Index f_index,
334  const Vector& f_grid,
335  const GridPos& gp_p,
336  const GridPos& gp_lat,
337  const GridPos& gp_lon,
338  ConstVectorView p_grid_cloud,
339  ConstTensor3View t_field_cloud,
340  ConstTensor4View vmr_field_cloud,
341  const Tensor4& pnd_field,
344  const Vector& rte_los)
345 
346 {
347  const Index ns = vmr_field_cloud.nbooks();
348  const Index N_se = pnd_field.nbooks();
349  Matrix pnd_ppath(N_se, 1);
350  Vector t_ppath(1);
351  Vector p_ppath(1); //may not be efficient with unecessary vectors
352  EnergyLevelMap nlte_dummy;
353  Matrix itw_p(1, 2);
354  ArrayOfGridPos ao_gp_p(1), ao_gp_lat(1), ao_gp_lon(1);
355  Matrix vmr_ppath(ns, 1), itw_field;
356 
357  //local versions of workspace variables
359  local_partial_dummy; // This is right since there should be only clearsky partials
361  local_dnlte_dx_source_dummy; // This is right since there should be only clearsky partials
363  local_nlte_dsource_dx_dummy; // This is right since there should be only clearsky partials
364  ArrayOfPropagationMatrix local_propmat_clearsky;
365  ArrayOfStokesVector local_nlte_source_dummy; //FIXME: Do this right?
366  StokesVector local_abs_vec;
367  PropagationMatrix local_ext_mat;
368 
369  ao_gp_p[0] = gp_p;
370  ao_gp_lat[0] = gp_lat;
371  ao_gp_lon[0] = gp_lon;
372 
373  cloud_atm_vars_by_gp(p_ppath,
374  t_ppath,
375  vmr_ppath,
376  pnd_ppath,
377  ao_gp_p,
378  ao_gp_lat,
379  ao_gp_lon,
381  p_grid_cloud,
382  t_field_cloud,
383  vmr_field_cloud,
384  pnd_field);
385  pnd_vec = pnd_ppath(joker, 0);
386  temperature = t_ppath[0];
387 
388  const Vector rtp_mag_dummy(3, 0);
389  const Vector ppath_los_dummy;
390 
391  //rtp_vmr = vmr_ppath(joker,0);
393  local_propmat_clearsky,
394  local_nlte_source_dummy,
395  local_partial_dummy,
396  local_dnlte_dx_source_dummy,
397  local_nlte_dsource_dx_dummy,
399  f_grid[Range(f_index, 1)],
400  rtp_mag_dummy,
401  ppath_los_dummy,
402  p_ppath[0],
403  temperature,
404  nlte_dummy,
405  vmr_ppath(joker, 0),
407 
409  local_ext_mat, local_abs_vec, local_propmat_clearsky);
410 
411  local_ext_mat.MatrixAtPosition(ext_mat_mono);
412  abs_vec_mono = local_abs_vec.VectorAtPosition();
413 
414  ArrayOfArrayOfTensor5 ext_mat_Nse;
415  ArrayOfArrayOfTensor4 abs_vec_Nse;
416  ArrayOfArrayOfIndex ptypes_Nse;
417  Matrix t_ok;
418  ArrayOfTensor5 ext_mat_ssbulk;
419  ArrayOfTensor4 abs_vec_ssbulk;
420  ArrayOfIndex ptype_ssbulk;
421  Tensor5 ext_mat_bulk;
422  Tensor4 abs_vec_bulk;
423  Index ptype_bulk;
424 
425  Vector sca_dir;
426  mirror_los(sca_dir, rte_los, 3);
427  Matrix dir_array(1, 2, 0.);
428  dir_array(0, joker) = sca_dir;
429  //
430  opt_prop_NScatElems(ext_mat_Nse,
431  abs_vec_Nse,
432  ptypes_Nse,
433  t_ok,
434  scat_data,
435  stokes_dim,
436  t_ppath,
437  dir_array,
438  f_index);
439  //
440  opt_prop_ScatSpecBulk(ext_mat_ssbulk,
441  abs_vec_ssbulk,
442  ptype_ssbulk,
443  ext_mat_Nse,
444  abs_vec_Nse,
445  ptypes_Nse,
446  pnd_ppath,
447  t_ok);
448  opt_prop_Bulk(ext_mat_bulk,
449  abs_vec_bulk,
450  ptype_bulk,
451  ext_mat_ssbulk,
452  abs_vec_ssbulk,
453  ptype_ssbulk);
454 
455  ext_mat_mono += ext_mat_bulk(0, 0, 0, joker, joker);
456  abs_vec_mono += abs_vec_bulk(0, 0, 0, joker);
457 }
458 
460  VectorView temperature,
461  MatrixView vmr,
462  MatrixView pnd,
463  const ArrayOfGridPos& gp_p,
464  const ArrayOfGridPos& gp_lat,
465  const ArrayOfGridPos& gp_lon,
467  ConstVectorView p_grid_cloud,
468  ConstTensor3View t_field_cloud,
469  ConstTensor4View vmr_field_cloud,
471  Index np = gp_p.nelem();
472  assert(pressure.nelem() == np);
473  Index ns = vmr_field_cloud.nbooks();
474  Index N_se = pnd_field.nbooks();
475  ArrayOfGridPos gp_p_cloud = gp_p;
476  ArrayOfGridPos gp_lat_cloud = gp_lat;
477  ArrayOfGridPos gp_lon_cloud = gp_lon;
478  Index atmosphere_dim = 3;
479 
480  for (Index i = 0; i < np; i++) {
481  // Calculate grid positions inside the cloud.
482  gp_p_cloud[i].idx -= cloudbox_limits[0];
483  gp_lat_cloud[i].idx -= cloudbox_limits[2];
484  gp_lon_cloud[i].idx -= cloudbox_limits[4];
485  }
486 
487  const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
488  const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
489  const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
490  gridpos_upperend_check(gp_p_cloud[0], n1);
491  gridpos_upperend_check(gp_p_cloud[np - 1], n1);
492  gridpos_upperend_check(gp_lat_cloud[0], n2);
493  gridpos_upperend_check(gp_lat_cloud[np - 1], n2);
494  gridpos_upperend_check(gp_lon_cloud[0], n3);
495  gridpos_upperend_check(gp_lon_cloud[np - 1], n3);
496 
497  // Determine the pressure at each propagation path point
498  Matrix itw_p(np, 2);
499  //
500  //interpweights( itw_p, ppath.gp_p );
501  interpweights(itw_p, gp_p_cloud);
502  itw2p(pressure, p_grid_cloud, gp_p_cloud, itw_p);
503 
504  // Determine the atmospheric temperature and species VMR at
505  // each propagation path point
506  Matrix itw_field;
507  //
509  itw_field, atmosphere_dim, gp_p_cloud, gp_lat_cloud, gp_lon_cloud);
510  //
511  interp_atmfield_by_itw(temperature,
513  t_field_cloud,
514  gp_p_cloud,
515  gp_lat_cloud,
516  gp_lon_cloud,
517  itw_field);
518  //
519  for (Index is = 0; is < ns; is++) {
520  interp_atmfield_by_itw(vmr(is, joker),
522  vmr_field_cloud(is, joker, joker, joker),
523  gp_p_cloud,
524  gp_lat_cloud,
525  gp_lon_cloud,
526  itw_field);
527  }
528 
529  //Determine the particle number density for every scattering element at
530  // each propagation path point
531  for (Index i_se = 0; i_se < N_se; i_se++) {
532  // if grid positions still outside the range the propagation path step
533  // must be outside the cloudbox and pnd is set to zero
534  interp_atmfield_by_itw(pnd(i_se, joker),
536  pnd_field(i_se, joker, joker, joker),
537  gp_p_cloud,
538  gp_lat_cloud,
539  gp_lon_cloud,
540  itw_field);
541  }
542 }
543 
545  MatrixView& trans_mat,
546  const Ppath& ppath,
548  const Index stokes_dim,
549  const Index f_index,
550  const Vector& f_grid,
551  const Vector& p_grid,
552  const Tensor3& t_field,
553  const Tensor4& vmr_field,
555  const Tensor4& pnd_field,
557  const Verbosity& verbosity) {
558  bool inside_cloud;
559  const Index np = ppath.np;
560  ArrayOfMatrix ext_matArray(2);
561  ArrayOfMatrix trans_matArray(2);
562  Index N_se = pnd_field.nbooks(); //Number of scattering elements
563  Vector pnd_vec(N_se);
564  Vector abs_vec_mono(stokes_dim);
566  Matrix ext_mat_mono(stokes_dim, stokes_dim);
567  Matrix incT(stokes_dim, stokes_dim, 0.0);
568  Numeric temperature;
569  Numeric dl = -999;
570 
571  CREATE_OUT0;
572 
573  id_mat(trans_mat);
574 
575  if (np > 1) {
576  // range defining cloudbox
577  Range p_range(cloudbox_limits[0],
578  cloudbox_limits[1] - cloudbox_limits[0] + 1);
579  Range lat_range(cloudbox_limits[2],
580  cloudbox_limits[3] - cloudbox_limits[2] + 1);
581  Range lon_range(cloudbox_limits[4],
582  cloudbox_limits[5] - cloudbox_limits[4] + 1);
583 
584  inside_cloud = is_gp_inside_cloudbox(ppath.gp_p[np - 1],
585  ppath.gp_lat[np - 1],
586  ppath.gp_lon[np - 1],
588  0,
589  3);
590  if (inside_cloud) {
592  ext_mat_mono,
593  abs_vec_mono,
594  pnd_vec,
595  temperature,
597  stokes_dim,
598  f_index,
599  f_grid,
600  ppath.gp_p[np - 1],
601  ppath.gp_lat[np - 1],
602  ppath.gp_lon[np - 1],
603  p_grid[p_range],
604  t_field(p_range, lat_range, lon_range),
605  vmr_field(joker, p_range, lat_range, lon_range),
606  pnd_field,
607  scat_data,
609  ppath.los(np - 1, joker));
610  } else {
612  ext_mat_mono,
613  abs_vec_mono,
614  temperature,
616  f_grid[f_index],
617  ppath.gp_p[np - 1],
618  ppath.gp_lat[np - 1],
619  ppath.gp_lon[np - 1],
620  p_grid,
621  t_field,
622  vmr_field);
623  pnd_vec = 0.0;
624  }
625 
626  trans_matArray[1] = trans_mat;
627  ext_matArray[1] = ext_mat_mono;
628 
629  // Index in ppath of end point considered presently
630  for (Index ip = np - 2; ip >= 0; ip--) {
631  dl = ppath.lstep[ip];
632 
633  ext_matArray[0] = ext_matArray[1];
634  trans_matArray[0] = trans_matArray[1];
635 
636  inside_cloud = is_gp_inside_cloudbox(ppath.gp_p[ip],
637  ppath.gp_lat[ip],
638  ppath.gp_lon[ip],
640  0,
641  3);
642  if (inside_cloud) {
644  ext_mat_mono,
645  abs_vec_mono,
646  pnd_vec,
647  temperature,
649  stokes_dim,
650  f_index,
651  f_grid,
652  ppath.gp_p[ip],
653  ppath.gp_lat[ip],
654  ppath.gp_lon[ip],
655  p_grid[p_range],
656  t_field(p_range, lat_range, lon_range),
657  vmr_field(joker, p_range, lat_range, lon_range),
658  pnd_field,
659  scat_data,
661  ppath.los(ip, joker));
662  } else {
664  ext_mat_mono,
665  abs_vec_mono,
666  temperature,
668  f_grid[f_index],
669  ppath.gp_p[ip],
670  ppath.gp_lat[ip],
671  ppath.gp_lon[ip],
672  p_grid,
673  t_field,
674  vmr_field);
675  pnd_vec = 0.0;
676  }
677 
678  ext_matArray[1] = ext_mat_mono;
679  ext_mat = ext_matArray[0];
680  ext_mat += ext_matArray[1]; // Factor 2 fixed by using dl/2
681  //
682  {
683  Index extmat_case = 0;
684  ext2trans(incT, extmat_case, ext_mat, dl / 2);
685  }
686  //
687  mult(trans_mat, incT, trans_matArray[1]);
688  trans_matArray[1] = trans_mat;
689 
690  } // for( ip... )
691  } // if( np > 1 )
692 }
693 
696  bool is_anyptype_nonTotRan = false;
697  for (Index i_ss = 0;
698  is_anyptype_nonTotRan == false && i_ss < scat_data.nelem();
699  i_ss++) {
700  for (Index i_se = 0;
701  is_anyptype_nonTotRan == false && i_se < scat_data[i_ss].nelem();
702  i_se++) {
703  if (scat_data[i_ss][i_se].ptype > PTYPE_TOTAL_RND) {
704  is_anyptype_nonTotRan = true;
705  }
706  }
707  }
708  return is_anyptype_nonTotRan;
709 }
710 
712  MatrixView evol_op,
713  Vector& abs_vec_mono,
714  Numeric& temperature,
715  MatrixView ext_mat_mono,
716  Rng& rng,
717  Vector& rte_pos,
718  Vector& rte_los,
719  Vector& pnd_vec,
720  Numeric& g,
721  Ppath& ppath_step,
722  Index& termination_flag,
723  bool& inside_cloud,
724  const Agenda& ppath_step_agenda,
725  const Numeric& ppath_lmax,
726  const Numeric& ppath_lraytrace,
727  const Numeric& taustep_limit,
729  const Index stokes_dim,
730  const Index f_index,
731  const Vector& f_grid,
732  const Vector& p_grid,
733  const Vector& lat_grid,
734  const Vector& lon_grid,
735  const Tensor3& z_field,
736  const Vector& refellipsoid,
737  const Matrix& z_surface,
738  const Tensor3& t_field,
739  const Tensor4& vmr_field,
741  const Tensor4& pnd_field,
743  const Verbosity& verbosity) {
744  ArrayOfMatrix evol_opArray(2);
745  ArrayOfMatrix ext_matArray(2);
746  ArrayOfVector abs_vecArray(2);
747  ArrayOfVector pnd_vecArray(2);
749  Matrix incT(stokes_dim, stokes_dim, 0.0);
750  Vector tArray(2);
752  Numeric k;
753  Numeric ds, dl = -999;
754  Index istep = 0; // Counter for number of steps
755  Matrix old_evol_op(stokes_dim, stokes_dim);
756 
757  CREATE_OUT0;
758 
759  //at the start of the path the evolution operator is the identity matrix
760  id_mat(evol_op);
761  evol_opArray[1] = evol_op;
762 
763  // range defining cloudbox
764  Range p_range(cloudbox_limits[0],
765  cloudbox_limits[1] - cloudbox_limits[0] + 1);
766  Range lat_range(cloudbox_limits[2],
767  cloudbox_limits[3] - cloudbox_limits[2] + 1);
768  Range lon_range(cloudbox_limits[4],
769  cloudbox_limits[5] - cloudbox_limits[4] + 1);
770 
771  //initialise Ppath with ppath_start_stepping
773  3,
774  p_grid,
775  lat_grid,
776  lon_grid,
777  z_field,
778  refellipsoid,
779  z_surface,
780  0,
782  false,
783  rte_pos,
784  rte_los,
785  verbosity);
786 
787  // Check if we have already has radiative background
789  termination_flag = ppath_what_background(ppath_step);
790  g = 1;
791  return;
792  }
793 
794  // Index in ppath_step of end point considered presently
795  Index ip = 0;
796 
797  // Is point ip inside the cloudbox?
798  inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
799  ppath_step.gp_lat[ip],
800  ppath_step.gp_lon[ip],
802  0,
803  3);
804 
805  // Determine radiative properties at point
806  if (inside_cloud) {
808  ext_mat_mono,
809  abs_vec_mono,
810  pnd_vec,
811  temperature,
813  stokes_dim,
814  f_index,
815  f_grid,
816  ppath_step.gp_p[0],
817  ppath_step.gp_lat[0],
818  ppath_step.gp_lon[0],
819  p_grid[p_range],
820  t_field(p_range, lat_range, lon_range),
821  vmr_field(joker, p_range, lat_range, lon_range),
822  pnd_field,
823  scat_data,
825  ppath_step.los(0, joker));
826  } else {
828  ext_mat_mono,
829  abs_vec_mono,
830  temperature,
832  f_grid[f_index],
833  ppath_step.gp_p[0],
834  ppath_step.gp_lat[0],
835  ppath_step.gp_lon[0],
836  p_grid,
837  t_field,
838  vmr_field);
839  pnd_vec = 0.0;
840  }
841 
842  // Move the data to end point containers
843  ext_matArray[1] = ext_mat_mono;
844  abs_vecArray[1] = abs_vec_mono;
845  tArray[1] = temperature;
846  pnd_vecArray[1] = pnd_vec;
847 
848  //draw random number to determine end point
849  Numeric r = rng.draw();
850 
851  termination_flag = 0;
852 
853  while ((evol_op(0, 0) > r) && (!termination_flag)) {
854  istep++;
855 
856  if (istep > 100000) {
857  throw runtime_error(
858  "100000 path points have been reached. "
859  "Is this an infinite loop?");
860  }
861 
862  evol_opArray[0] = evol_opArray[1];
863  ext_matArray[0] = ext_matArray[1];
864  abs_vecArray[0] = abs_vecArray[1];
865  tArray[0] = tArray[1];
866  pnd_vecArray[0] = pnd_vecArray[1];
867 
868  // The algorith to meet taustep_lim:
869  // When first calculating a new ppath_step, it assumed that the present
870  // ppath position holds the highest extinction. If the extinction at the
871  // next position is higher, the criterion is checked and a new ppath_step
872  // calculation is triggered if found necessary.
873  // This should work in most cases, but is not 100% safe. Consider a case
874  // with ppath_lmax = -1 and the extinction is zero at all grid box
875  // corners except one. The two "test points" can then both get an
876  // extinction of zero, while in fact is non-zero through the grid box and
877  // the optical depth is underestimated. But this was handled equally bad
878  // before taustep_limit was introduced (2016-10-10, PE)
879  bool oktaustep = false;
880  Index ppath_try = 1;
881  const Index lmax_limit = 10;
882 
883  while (!oktaustep) {
884  // Shall new ppath_step be calculated?
885  if (ip == ppath_step.np - 1) {
886  Numeric lmax = taustep_limit / ext_mat_mono(0, 0);
887  if (ppath_lmax > 0) {
888  lmax = min(ppath_lmax, lmax);
889  }
890  if (lmax < lmax_limit) {
891  lmax = lmax_limit;
892  }
893  //cout << ppath_try << ", lmax = " << lmax << endl;
894  //Print( ppath_step, 0, verbosity );
895 
897  ppath_step,
898  lmax,
900  f_grid[Range(f_index, 1)],
902  ip = 1;
903 
904  inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
905  ppath_step.gp_lat[ip],
906  ppath_step.gp_lon[ip],
908  0,
909  3);
910  } else {
911  ip++;
912  }
913 
914  if (inside_cloud) {
916  ext_mat_mono,
917  abs_vec_mono,
918  pnd_vec,
919  temperature,
921  stokes_dim,
922  f_index,
923  f_grid,
924  ppath_step.gp_p[ip],
925  ppath_step.gp_lat[ip],
926  ppath_step.gp_lon[ip],
927  p_grid[p_range],
928  t_field(p_range, lat_range, lon_range),
929  vmr_field(joker, p_range, lat_range, lon_range),
930  pnd_field,
931  scat_data,
933  ppath_step.los(ip, joker));
934  } else {
936  ext_mat_mono,
937  abs_vec_mono,
938  temperature,
940  f_grid[f_index],
941  ppath_step.gp_p[ip],
942  ppath_step.gp_lat[ip],
943  ppath_step.gp_lon[ip],
944  p_grid,
945  t_field,
946  vmr_field);
947  pnd_vec = 0.0;
948  }
949 
950  dl = ppath_step.lstep[ip - 1];
951 
952  // Check if taustep_limit criterion is met. OK if:
953  // 1. Ppath step already recalculated
954  // 2. New ext_mat <= old one
955  // 3. New ext_mat bigger, but tau of step still below limit
956  if (ppath_try > 1 || ext_mat_mono(0, 0) <= ext_matArray[0](0, 0) ||
957  (ext_mat_mono(0, 0) + ext_matArray[0](0, 0)) * dl / 2 <=
958  taustep_limit) {
959  oktaustep = true;
960  } else {
961  // We trigger a recalculation of ppath_step, from the previous
962  // point
963  ppath_step.np = ip;
964  ip--;
965  ppath_try = 2;
966  // If a background found in first try this has to be reset:
968  }
969  } // while !oktuastep
970 
971  ext_matArray[1] = ext_mat_mono;
972  abs_vecArray[1] = abs_vec_mono;
973  tArray[1] = temperature;
974  pnd_vecArray[1] = pnd_vec;
975  ext_mat = ext_matArray[1];
976  ext_mat += ext_matArray[0]; // Factor 2 fixed by using dl/2
977  //
978  {
979  Index extmat_case = 0;
980  ext2trans(incT, extmat_case, ext_mat, dl / 2);
981  }
982  //
983  mult(evol_op, evol_opArray[0], incT);
984  evol_opArray[1] = evol_op;
985 
986  if (evol_op(0, 0) > r) {
987  // Check whether hit ground or space.
988  // path_step_agenda just detects surface intersections, and
989  // if TOA is reached requires a special check.
990  // But we are already ready if evol_op<=r
991  if (ip == ppath_step.np - 1) {
993  termination_flag = 2;
994  } //we have hit the surface
995  else if (fractional_gp(ppath_step.gp_p[ip]) >=
996  (Numeric)(p_grid.nelem() - 1) - 1e-3) {
997  termination_flag = 1;
998  } //we are at TOA
999  }
1000  }
1001  } // while
1002 
1003  if (termination_flag) { //we must have reached the cloudbox boundary
1004  rte_pos = ppath_step.pos(ip, joker);
1005  rte_los = ppath_step.los(ip, joker);
1006  g = evol_op(0, 0);
1007  } else {
1008  //find position...and evol_op..and everything else required at the new
1009  //scattering/emission point
1010  // GH 2011-09-14:
1011  // log(incT(0,0)) = log(exp(opt_depth_mat(0, 0))) = opt_depth_mat(0, 0)
1012  // Avoid loss of precision, use opt_depth_mat directly
1013  //k=-log(incT(0,0))/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
1014  // PE 2013-05-17, Now the above comes directly from ext_mat:
1015  k = ext_mat(0, 0) / 2; // Factor 2 as sum of end point values
1016  ds = log(evol_opArray[0](0, 0) / r) / k;
1017  g = k * r;
1018  Vector x(2, 0.0);
1019  //interpolate atmospheric variables required later
1020  ArrayOfGridPos gp(1);
1021  x[1] = dl;
1022  Vector itw(2);
1023 
1024  gridpos(gp, x, ds);
1025  assert(gp[0].idx == 0);
1026  interpweights(itw, gp[0]);
1027  interp(ext_mat_mono, itw, ext_matArray, gp[0]);
1028  ext_mat = ext_mat_mono;
1029  ext_mat += ext_matArray[gp[0].idx];
1030  //
1031  {
1032  Index extmat_case = 0;
1033  ext2trans(incT, extmat_case, ext_mat, ds / 2);
1034  }
1035  //
1036  mult(evol_op, evol_opArray[gp[0].idx], incT);
1037  interp(abs_vec_mono, itw, abs_vecArray, gp[0]);
1038  temperature = interp(itw, tArray, gp[0]);
1039  interp(pnd_vec, itw, pnd_vecArray, gp[0]);
1040  for (Index i = 0; i < 2; i++) {
1041  rte_pos[i] = interp(itw, ppath_step.pos(Range(ip - 1, 2), i), gp[0]);
1042  rte_los[i] = interp(itw, ppath_step.los(Range(ip - 1, 2), i), gp[0]);
1043  }
1044  rte_pos[2] = interp(itw, ppath_step.pos(Range(ip - 1, 2), 2), gp[0]);
1045  }
1046 
1047  assert(isfinite(g));
1048 
1049  // A dirty trick to avoid copying ppath_step
1050  const Index np = ip + 1;
1051  ppath_step.np = np;
1052 }
1053 
1055  MatrixView evol_op,
1056  Vector& abs_vec_mono,
1057  Numeric& temperature,
1058  MatrixView ext_mat_mono,
1059  Rng& rng,
1060  Vector& rte_pos,
1061  Vector& rte_los,
1062  Vector& pnd_vec,
1063  Numeric& stot,
1064  Numeric& ttot,
1065  Ppath& ppath_step,
1066  Index& termination_flag,
1067  bool& inside_cloud,
1068  const Agenda& ppath_step_agenda,
1069  const Numeric& ppath_lmax,
1070  const Numeric& ppath_lraytrace,
1072  const bool& anyptype_nonTotRan,
1073  const Index stokes_dim,
1074  const Index f_index,
1075  const Vector& f_grid,
1076  const Vector& Iprop,
1077  const Vector& p_grid,
1078  const Vector& lat_grid,
1079  const Vector& lon_grid,
1080  const Tensor3& z_field,
1081  const Vector& refellipsoid,
1082  const Matrix& z_surface,
1083  const Tensor3& t_field,
1084  const Tensor4& vmr_field,
1086  const Tensor4& pnd_field,
1088  const Verbosity& verbosity) {
1089  ArrayOfMatrix evol_opArray(2);
1090  ArrayOfMatrix ext_matArray(2);
1091  ArrayOfVector abs_vecArray(2);
1092  ArrayOfVector pnd_vecArray(2);
1094  Matrix incT(stokes_dim, stokes_dim, 0.0);
1095  Vector tArray(2);
1097  Numeric kI, kQ;
1098  Numeric ds, dt = -999, dl = -999;
1099  Index istep = 0; // Counter for number of steps
1100  Matrix old_evol_op(stokes_dim, stokes_dim);
1101  Vector local_rte_los(2);
1102 
1103  CREATE_OUT0;
1104 
1105  // Total path length starts at zero
1106  stot = 0.0;
1107  ttot = 0.0;
1108 
1109  //at the start of the path the evolution operator is the identity matrix
1110  id_mat(evol_op);
1111  evol_opArray[1] = evol_op;
1112 
1113  // range defining cloudbox
1114  Range p_range(cloudbox_limits[0],
1115  cloudbox_limits[1] - cloudbox_limits[0] + 1);
1116  Range lat_range(cloudbox_limits[2],
1117  cloudbox_limits[3] - cloudbox_limits[2] + 1);
1118  Range lon_range(cloudbox_limits[4],
1119  cloudbox_limits[5] - cloudbox_limits[4] + 1);
1120 
1121  //initialise Ppath with ppath_start_stepping
1123  3,
1124  p_grid,
1125  lat_grid,
1126  lon_grid,
1127  z_field,
1128  refellipsoid,
1129  z_surface,
1130  0,
1132  false,
1133  rte_pos,
1134  rte_los,
1135  verbosity);
1136 
1137  if (ppath_step.np == 0) {
1138  termination_flag = 1;
1139  return;
1140  }
1141 
1142  // Index in ppath_step of end point considered presently
1143  Index ip = 0;
1144 
1145  // Is point ip inside the cloudbox?
1146  inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
1147  ppath_step.gp_lat[ip],
1148  ppath_step.gp_lon[ip],
1150  0,
1151  3);
1152 
1153  // Determine radiative properties at point
1154  if (inside_cloud) {
1155  local_rte_los[0] = 180 - ppath_step.los(0, 0);
1156  local_rte_los[1] = ppath_step.los(0, 1) - 180;
1158  ext_mat_mono,
1159  abs_vec_mono,
1160  pnd_vec,
1161  temperature,
1163  stokes_dim,
1164  f_index,
1165  f_grid,
1166  ppath_step.gp_p[0],
1167  ppath_step.gp_lat[0],
1168  ppath_step.gp_lon[0],
1169  p_grid[p_range],
1170  t_field(p_range, lat_range, lon_range),
1171  vmr_field(joker, p_range, lat_range, lon_range),
1172  pnd_field,
1173  scat_data,
1175  local_rte_los);
1176  } else {
1178  ext_mat_mono,
1179  abs_vec_mono,
1180  temperature,
1182  f_grid[f_index],
1183  ppath_step.gp_p[0],
1184  ppath_step.gp_lat[0],
1185  ppath_step.gp_lon[0],
1186  p_grid,
1187  t_field,
1188  vmr_field);
1189  pnd_vec = 0.0;
1190  }
1191 
1192  // Move the data to end point containers
1193  ext_matArray[1] = ext_mat_mono;
1194  abs_vecArray[1] = abs_vec_mono;
1195  tArray[1] = temperature;
1196  pnd_vecArray[1] = pnd_vec;
1197 
1198  //draw random number to determine end point
1199  Numeric r = rng.draw();
1200 
1201  termination_flag = 0;
1202 
1203  stot = ppath_step.end_lstep;
1204  ttot = ppath_step.end_lstep / SPEED_OF_LIGHT;
1205 
1206  Numeric evop0 = 1;
1207  while ((evop0 > r) && (!termination_flag)) {
1208  istep++;
1209 
1210  if (istep > 25000) {
1211  throw runtime_error(
1212  "25000 path points have been reached. "
1213  "Is this an infinite loop?");
1214  }
1215 
1216  evol_opArray[0] = evol_opArray[1];
1217  ext_matArray[0] = ext_matArray[1];
1218  abs_vecArray[0] = abs_vecArray[1];
1219  tArray[0] = tArray[1];
1220  pnd_vecArray[0] = pnd_vecArray[1];
1221 
1222  // Shall new ppath_step be calculated?
1223  if (ip >= ppath_step.np - 1) {
1224  ip = 1;
1226  ppath_step,
1227  ppath_lmax,
1229  f_grid[Range(f_index, 1)],
1231 
1232  if (ppath_step.np <= 1) {
1233  termination_flag = 1;
1234  break;
1235  }
1236  inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
1237  ppath_step.gp_lat[ip],
1238  ppath_step.gp_lon[ip],
1240  0,
1241  3);
1242  } else {
1243  ip++;
1244  }
1245 
1246  dl = ppath_step.lstep[ip - 1];
1247  dt = dl * 0.5 * (ppath_step.ngroup[ip - 1] + ppath_step.ngroup[ip]) /
1249  stot += dl;
1250  ttot += dt;
1251  if (inside_cloud) {
1252  local_rte_los[0] = 180 - ppath_step.los(ip, 0);
1253  local_rte_los[1] = ppath_step.los(ip, 1) - 180;
1255  ext_mat_mono,
1256  abs_vec_mono,
1257  pnd_vec,
1258  temperature,
1260  stokes_dim,
1261  f_index,
1262  f_grid,
1263  ppath_step.gp_p[ip],
1264  ppath_step.gp_lat[ip],
1265  ppath_step.gp_lon[ip],
1266  p_grid[p_range],
1267  t_field(p_range, lat_range, lon_range),
1268  vmr_field(joker, p_range, lat_range, lon_range),
1269  pnd_field,
1270  scat_data,
1272  local_rte_los);
1273  } else {
1275  ext_mat_mono,
1276  abs_vec_mono,
1277  temperature,
1279  f_grid[f_index],
1280  ppath_step.gp_p[ip],
1281  ppath_step.gp_lat[ip],
1282  ppath_step.gp_lon[ip],
1283  p_grid,
1284  t_field,
1285  vmr_field);
1286  pnd_vec = 0.0;
1287  }
1288 
1289  ext_matArray[1] = ext_mat_mono;
1290  abs_vecArray[1] = abs_vec_mono;
1291  tArray[1] = temperature;
1292  pnd_vecArray[1] = pnd_vec;
1293  ext_mat = ext_matArray[1];
1294  ext_mat += ext_matArray[0]; // Factor 2 fixed by using dl/2
1295  //
1296  {
1297  Index extmat_case = 0;
1298  ext2trans(incT, extmat_case, ext_mat, dl / 2);
1299  }
1300  //
1301 
1302  mult(evol_op, incT, evol_opArray[0]);
1303  evol_opArray[1] = evol_op;
1304  evop0 = evol_op(0, 0);
1305 
1306  // Handle cross-talk for ptype==30
1307  if (stokes_dim > 1 && anyptype_nonTotRan) {
1308  const Numeric Q1 = evol_op(0, 1) * Iprop[1] / Iprop[0];
1309  evop0 += Q1;
1310  }
1311  if (evop0 > r) {
1312  // Check whether hit ground or space.
1313  // path_step_agenda just detects surface intersections, and
1314  // if TOA is reached requires a special check.
1315  // But we are already ready if evol_op<=r
1316  if (ip >= ppath_step.np - 1) {
1318  termination_flag = 2;
1319  } //we have hit the surface
1320  else if (fractional_gp(ppath_step.gp_p[ip]) >=
1321  (Numeric)(p_grid.nelem() - 1) - 1e-3) {
1322  termination_flag = 1;
1323  } //we are at TOA
1324  }
1325  }
1326  } // while
1327 
1328  if (termination_flag != 0) { //we must have reached the cloudbox boundary
1329  if (ip < ppath_step.np) {
1330  // This is not used if termination flag set,
1331  // so not an issue of ip is too large.
1332  // Need to think about this.
1333  rte_pos = ppath_step.pos(ip, joker);
1334  rte_los = ppath_step.los(ip, joker);
1335  }
1336  } else {
1337  //find position...and evol_op..and everything else required at the new
1338  //scattering/emission point
1339  const Numeric tol = 0.1; // Tolerance of 10 cm
1340  stot -= dl; // Take out last step because we are between stepping points
1341  ttot -= dt; // Take out last step because we are between stepping points
1342  kI = ext_mat(0, 0) / 2; // Factor 2 as sum of end point values
1343  kQ = ext_mat(0, 1) / 2; // Factor 2 as sum of end point values
1344  if (anyptype_nonTotRan) {
1345  const Numeric I1 = evol_opArray[0](0, 0);
1346  const Numeric Q1 = evol_opArray[0](0, 1) * Iprop[1] / Iprop[0];
1347 
1348  // Need to use root finding to solve for ds
1349  brent_zero(
1350  ds, (Numeric)0.0, ppath_step.lstep[ip - 1], tol, r, I1, Q1, kI, kQ);
1351  } else {
1352  // Simple inversion when no cross-talk between I and Q
1353  ds = log(evol_opArray[0](0, 0) / r) / kI;
1354  }
1355  stot += ds;
1356  ttot += ds * dt / dl;
1357  Vector x(2, 0.0);
1358 
1359  //interpolate atmospheric variables required later
1360  ArrayOfGridPos gp(1);
1361  x[1] = dl;
1362  Vector itw(2);
1363  gridpos(gp, x, ds);
1364  assert(gp[0].idx == 0);
1365  interpweights(itw, gp[0]);
1366  interp(ext_mat_mono, itw, ext_matArray, gp[0]);
1367  ext_mat = ext_mat_mono;
1368  ext_mat += ext_matArray[gp[0].idx];
1369  //
1370  {
1371  Index extmat_case = 0;
1372  ext2trans(incT, extmat_case, ext_mat, ds / 2);
1373  }
1374  //
1375  mult(evol_op, incT, evol_opArray[gp[0].idx]);
1376  interp(abs_vec_mono, itw, abs_vecArray, gp[0]);
1377  temperature = interp(itw, tArray, gp[0]);
1378  interp(pnd_vec, itw, pnd_vecArray, gp[0]);
1379  for (Index i = 0; i < 2; i++) {
1380  rte_pos[i] = interp(itw, ppath_step.pos(Range(ip - 1, 2), i), gp[0]);
1381  rte_los[i] = interp(itw, ppath_step.los(Range(ip - 1, 2), i), gp[0]);
1382  }
1383  rte_pos[2] = interp(itw, ppath_step.pos(Range(ip - 1, 2), 2), gp[0]);
1384  }
1385 
1386  // A dirty trick to avoid copying ppath_step
1387  const Index np = ip + 1;
1388  ppath_step.np = np;
1389 }
1390 
1391 void Sample_los(VectorView new_rte_los,
1392  Numeric& g_los_csc_theta,
1393  MatrixView Z,
1394  Rng& rng,
1397  const Index f_index,
1398  const Index stokes_dim,
1399  ConstVectorView pnd_vec,
1400  ConstVectorView Z11maxvector,
1401  const Numeric Csca,
1402  const Numeric rtp_temperature,
1403  const Index t_interp_order) {
1404  Numeric Z11max = 0;
1405  bool tryagain = true;
1406 
1407  Vector sca_dir;
1408  mirror_los(sca_dir, rte_los, 3);
1409 
1410  // Rejection method http://en.wikipedia.org/wiki/Rejection_sampling
1411  Index np = pnd_vec.nelem();
1412  assert(TotalNumberOfElements(scat_data) == np);
1413  for (Index i = 0; i < np; i++) {
1414  Z11max += Z11maxvector[i] * pnd_vec[i];
1415  }
1416 
1418  // allocating variables needed for pha_mat extraction - this seems a little
1419  // disadvantageous. If we move the whole function back into the calling one
1420  // (it's used only at one place anyways), we can do the allocation there once
1421  // and avoid that it is done everytime this function is called within a loop.
1422  ArrayOfArrayOfTensor6 pha_mat_Nse;
1423  ArrayOfArrayOfIndex ptypes_Nse;
1424  Matrix t_ok;
1425  ArrayOfTensor6 pha_mat_ssbulk;
1426  ArrayOfIndex ptype_ssbulk;
1427  Tensor6 pha_mat_bulk;
1428  Index ptype_bulk;
1429  Matrix pdir(1, 2), idir(1, 2);
1430  Vector t(1, rtp_temperature);
1431  Matrix pnds(np, 1);
1432  pnds(joker, 0) = pnd_vec;
1433 
1434  while (tryagain) {
1435  new_rte_los[0] = acos(1 - 2 * rng.draw()) * RAD2DEG;
1436  new_rte_los[1] = rng.draw() * 360 - 180;
1437 
1438  //Calculate Phase matrix////////////////////////////////
1439  Vector inc_dir;
1440  mirror_los(inc_dir, new_rte_los, 3);
1441 
1442  //pha_mat_singleCalc( Z, sca_dir[0], sca_dir[1], inc_dir[0], inc_dir[1],
1443  // scat_data_mono, stokes_dim, pnd_vec, rtp_temperature,
1444  // verbosity );
1445 
1446  pdir(0, joker) = sca_dir;
1447  idir(0, joker) = inc_dir;
1448  pha_mat_NScatElems(pha_mat_Nse,
1449  ptypes_Nse,
1450  t_ok,
1451  scat_data,
1452  stokes_dim,
1453  t,
1454  pdir,
1455  idir,
1456  f_index,
1457  t_interp_order);
1459  pha_mat_ssbulk, ptype_ssbulk, pha_mat_Nse, ptypes_Nse, pnds, t_ok);
1460  pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
1461  Z = pha_mat_bulk(0, 0, 0, 0, joker, joker);
1462 
1463  if (rng.draw() <= Z(0, 0) / Z11max) //then new los is accepted
1464  {
1465  tryagain = false;
1466  }
1467  }
1468  g_los_csc_theta = Z(0, 0) / Csca;
1469 }
1470 
1471 void Sample_los_uniform(VectorView new_rte_los, Rng& rng) {
1472  new_rte_los[1] = rng.draw() * 360 - 180;
1473  new_rte_los[0] = acos(1 - 2 * rng.draw()) * RAD2DEG;
1474 }
PropagationMatrix::MatrixAtPosition
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
Definition: propagationmatrix.cc:1928
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::Ppath::pos
std::size_t pos() const noexcept
Definition: autoarts.h:1300
ARTS::Var::ppath_step
Ppath ppath_step(Workspace &ws) noexcept
Definition: autoarts.h:5213
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
ARTS::Var::ppath_step_agenda
Agenda ppath_step_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5220
auto_md.h
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
mcPathTraceGeneral
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_lmax, const Numeric &ppath_lraytrace, const Numeric &taustep_limit, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Index f_index, const Vector &f_grid, 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 ArrayOfArrayOfSingleScatteringData &scat_data, const Verbosity &verbosity)
mcPathTraceGeneral.
Definition: montecarlo.cc:711
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
mc_interp.h
Interpolation classes and functions created for use within Monte Carlo scattering simulations.
itw2p
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Definition: special_interp.cc:718
RAD2DEG
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
ARTS::Var::pnd_field
Tensor4 pnd_field(Workspace &ws) noexcept
Definition: autoarts.h:5081
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
interp_atmfield_by_itw
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)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Definition: special_interp.cc:101
ARTS::Var::f_index
Index f_index(Workspace &ws) noexcept
Definition: autoarts.h:3464
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
PropagationMatrix
Definition: propagationmatrix.h:87
brent_zero
void brent_zero(Numeric &sb, const Numeric &a, const Numeric &b, const Numeric &t, const Numeric &rn, const Numeric &I, const Numeric &Q, const Numeric &KI, const Numeric &KQ)
brent_zero.
Definition: montecarlo.cc:132
is_gp_inside_cloudbox
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:499
opt_prop_ScatSpecBulk
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
Definition: optproperties.cc:148
fractional_gp
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Definition: interpolation.cc:492
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ArrayOfRetrievalQuantity
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:402
mcPathTraceRadar
void mcPathTraceRadar(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 &stot, Numeric &ttot, Ppath &ppath_step, Index &termination_flag, bool &inside_cloud, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &propmat_clearsky_agenda, const bool &anyptype_nonTotRan, const Index stokes_dim, const Index f_index, const Vector &f_grid, const Vector &Iprop, 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 ArrayOfArrayOfSingleScatteringData &scat_data, const Verbosity &verbosity)
mcPathTraceRadar.
Definition: montecarlo.cc:1054
montecarlo.h
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
Sample_los_uniform
void Sample_los_uniform(VectorView new_rte_los, Rng &rng)
Sample_los_uniform.
Definition: montecarlo.cc:1471
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Agenda
The Agenda class.
Definition: agenda_class.h:44
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
opt_prop_NScatElems
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
Definition: optproperties.cc:248
Array
This can be used to make arrays out of anything.
Definition: array.h:108
Rng::draw
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
is_anyptype_nonTotRan
bool is_anyptype_nonTotRan(const ArrayOfArrayOfSingleScatteringData &scat_data)
is_anyptype_nonTotRan.
Definition: montecarlo.cc:694
min
#define min(a, b)
Definition: legacy_continua.cc:20628
Sample_los
void Sample_los(VectorView new_rte_los, Numeric &g_los_csc_theta, MatrixView Z, Rng &rng, ConstVectorView rte_los, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index f_index, const Index stokes_dim, ConstVectorView pnd_vec, ConstVectorView Z11maxvector, const Numeric Csca, const Numeric rtp_temperature, const Index t_interp_order)
Sample_los.
Definition: montecarlo.cc:1391
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:204
pha_mat_Bulk
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
Definition: optproperties.cc:670
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
gridpos_upperend_check
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
Definition: interpolation.cc:589
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::rte_los
Vector rte_los(Workspace &ws) noexcept
Definition: autoarts.h:5651
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::ppath
Ppath ppath(Workspace &ws) noexcept
Definition: autoarts.h:5139
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
propmat_clearsky_agendaExecute
void propmat_clearsky_agendaExecute(Workspace &ws, ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfStokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_dx_source, ArrayOfStokesVector &nlte_dsource_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:23513
pha_mat_NScatElems
void pha_mat_NScatElems(ArrayOfArrayOfTensor6 &pha_mat, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_index, const Index &t_interp_order)
Phase matrices from all scattering elements.
Definition: optproperties.cc:796
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
geodetic.h
ARTS::Var::refellipsoid
Vector refellipsoid(Workspace &ws) noexcept
Definition: autoarts.h:5519
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
ARTS::Var::ext_mat
PropagationMatrix ext_mat(Workspace &ws) noexcept
Definition: autoarts.h:3390
EnergyLevelMap
Definition: energylevelmap.h:60
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
clear_rt_vars_at_gp
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:244
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
ARTS::Var::ppath_lmax
Numeric ppath_lmax(Workspace &ws) noexcept
Definition: autoarts.h:5183
ARTS::Var::ppath_lraytrace
Numeric ppath_lraytrace(Workspace &ws) noexcept
Definition: autoarts.h:5195
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
Range
The range class.
Definition: matpackI.h:160
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ppath_start_stepping
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)
Initiates a Ppath structure for calculation of a path with ppath_step.
Definition: ppath.cc:4495
interp_atmfield_gp2itw
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
Definition: special_interp.cc:75
ppath_step_agendaExecute
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:24832
opt_prop_sum_propmat_clearsky
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const ArrayOfPropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
Definition: optproperties.cc:2414
q
#define q
Definition: legacy_continua.cc:21712
get_ppath_transmat
void get_ppath_transmat(Workspace &ws, MatrixView &trans_mat, const Ppath &ppath, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Index f_index, const Vector &f_grid, const Vector &p_grid, const Tensor3 &t_field, const Tensor4 &vmr_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Verbosity &verbosity)
get_ppath_transmat.
Definition: montecarlo.cc:544
ARTS::Var::rtp_temperature
Numeric rtp_temperature(Workspace &ws) noexcept
Definition: autoarts.h:5807
cloud_atm_vars_by_gp
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:459
opt_prop_Bulk
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
Definition: optproperties.cc:102
mirror_los
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:2290
ext2trans
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Converts an extinction matrix to a transmission matrix.
Definition: rte.cc:800
Rng
Definition: rng.h:554
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
ext_I
Numeric ext_I(const Numeric &I, const Numeric &Q, const Numeric &kI, const Numeric &kQ, const Numeric &s)
ext_I.
Definition: montecarlo.cc:58
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
ppath_set_background
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
Definition: ppath.cc:1467
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1494
Vector
The Vector class.
Definition: matpackI.h:860
PTYPE_TOTAL_RND
@ PTYPE_TOTAL_RND
Definition: optproperties.h:56
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
ARTS::Var::z_surface
Matrix z_surface(Workspace &ws) noexcept
Definition: autoarts.h:7754
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
StokesVector::VectorAtPosition
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
Definition: propagationmatrix.h:1266
ns
#define ns
Definition: legacy_continua.cc:22183
pha_mat_ScatSpecBulk
void pha_mat_ScatSpecBulk(ArrayOfTensor6 &pha_mat, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor6 &pha_mat_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk phase matrices.
Definition: optproperties.cc:708
cloudy_rt_vars_at_gp
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 Index f_index, const Vector &f_grid, 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 ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfIndex &cloudbox_limits, const Vector &rte_los)
cloudy_rt_vars_at_gp.
Definition: montecarlo.cc:326