ARTS  2.3.1285(git:92a29ea9-dirty)
m_montecarlo.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Cory Davis <cory@met.ed.ac.uk>
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 
27 /*===========================================================================
28  === External declarations
29  ===========================================================================*/
30 
31 #include <cmath>
32 #include <ctime>
33 #include <fstream>
34 #include <stdexcept>
35 #include "arts.h"
36 #include "auto_md.h"
37 #include "check_input.h"
38 #include "lin_alg.h"
39 #include "logic.h"
40 #include "math_funcs.h"
41 #include "matpackI.h"
42 #include "mc_interp.h"
43 #include "messages.h"
44 #include "montecarlo.h"
45 #include "physics_funcs.h"
46 #include "ppath.h"
47 #include "refraction.h"
48 #include "rng.h"
49 #include "rte.h"
50 #include "special_interp.h"
51 #include "xml_io.h"
52 
53 extern const Numeric DEG2RAD;
54 extern const Numeric RAD2DEG;
55 extern const Numeric PI;
56 extern const Numeric BOLTZMAN_CONST;
57 extern const Numeric SPEED_OF_LIGHT;
58 
59 /*===========================================================================
60  === The functions (in alphabetical order)
61  ===========================================================================*/
62 
63 /* Workspace method: Doxygen documentation will be auto-generated */
65  //keyword arguments
66  const Numeric& za_sigma,
67  const Numeric& aa_sigma,
68  const Verbosity&) {
69  mc_antenna.set_gaussian(za_sigma, aa_sigma);
70 }
71 
72 /* Workspace method: Doxygen documentation will be auto-generated */
74  //keyword arguments
75  const Numeric& za_fwhm,
76  const Numeric& aa_fwhm,
77  const Verbosity&) {
78  mc_antenna.set_gaussian_fwhm(za_fwhm, aa_fwhm);
79 }
80 
81 /* Workspace method: Doxygen documentation will be auto-generated */
82 void mc_antennaSetPencilBeam(MCAntenna& mc_antenna, const Verbosity&) {
83  mc_antenna.set_pencil_beam();
84 }
85 
86 /* Workspace method: Doxygen documentation will be auto-generated */
88  Vector& y,
89  Index& mc_iteration_count,
90  Vector& mc_error,
91  Tensor3& mc_points,
92  ArrayOfIndex& mc_source_domain,
93  ArrayOfIndex& mc_scat_order,
94  const MCAntenna& mc_antenna,
95  const Vector& f_grid,
96  const Index& f_index,
97  const Matrix& sensor_pos,
98  const Matrix& sensor_los,
99  const Index& stokes_dim,
100  const Index& atmosphere_dim,
101  const Agenda& ppath_step_agenda,
102  const Numeric& ppath_lmax,
103  const Numeric& ppath_lraytrace,
104  const Agenda& iy_space_agenda,
105  const Agenda& surface_rtprop_agenda,
106  const Agenda& propmat_clearsky_agenda,
107  const Vector& p_grid,
108  const Vector& lat_grid,
109  const Vector& lon_grid,
110  const Tensor3& z_field,
111  const Vector& refellipsoid,
112  const Matrix& z_surface,
113  const Tensor3& t_field,
114  const Tensor4& vmr_field,
115  const Index& cloudbox_on,
116  const ArrayOfIndex& cloudbox_limits,
117  const Tensor4& pnd_field,
118  const ArrayOfArrayOfSingleScatteringData& scat_data,
119  const Index& atmfields_checked,
120  const Index& atmgeom_checked,
121  const Index& scat_data_checked,
122  const Index& cloudbox_checked,
123  const String& iy_unit,
124  const Index& mc_seed,
125  const Numeric& std_err,
126  const Index& max_time,
127  const Index& max_iter,
128  const Index& min_iter,
129  const Numeric& taustep_limit,
130  const Index& l_mc_scat_order,
131  const Index& t_interp_order,
132  const Verbosity& verbosity) {
133  // Checks of input
134  //
135  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
136  if (atmfields_checked != 1)
137  throw runtime_error(
138  "The atmospheric fields must be flagged to have "
139  "passed a consistency check (atmfields_checked=1).");
140  if (atmgeom_checked != 1)
141  throw runtime_error(
142  "The atmospheric geometry must be flagged to have "
143  "passed a consistency check (atmgeom_checked=1).");
144 
145  if (!cloudbox_on)
146  throw runtime_error("The cloudbox must be activated (cloudbox_on=1).");
147  if (cloudbox_checked != 1)
148  throw runtime_error(
149  "The cloudbox must be flagged to have "
150  "passed a consistency check (cloudbox_checked=1).");
151  if (scat_data_checked != 1)
152  throw runtime_error(
153  "The scat_data must be flagged to have "
154  "passed a consistency check (scat_data_checked=1).");
155 
156  if (min_iter < 100) throw runtime_error("*mc_min_iter* must be >= 100.");
157 
158  if (max_time < 0 && max_iter < 0 && std_err < 0)
159  throw runtime_error(
160  "At least one of std_err, max_time, and max_iter "
161  "must be positive.");
162 
163  if (l_mc_scat_order <= 0)
164  throw runtime_error("*l_max_scat_order* must be > 0.");
165 
166  if (f_index < 0)
167  throw runtime_error(
168  "The option of f_index < 0 is not handled by this "
169  "method.");
170  if (f_index >= f_grid.nelem())
171  throw runtime_error("*f_index* is outside the range of *f_grid*.");
172 
173  if (atmosphere_dim != 3)
174  throw runtime_error("Only 3D atmospheres are handled. ");
175 
176  if (sensor_pos.ncols() != 3) {
177  ostringstream os;
178  os << "Expected number of columns in sensor_pos: 3.\n";
179  os << "Found: " << sensor_pos.ncols();
180  throw runtime_error(os.str());
181  }
182  if (sensor_pos.nrows() != 1) {
183  ostringstream os;
184  os << "Expected number of rows in sensor_pos: 1.\n";
185  os << "Found: " << sensor_pos.nrows();
186  throw runtime_error(os.str());
187  }
188 
189  if (!(sensor_los.ncols() == 2)) {
190  ostringstream os;
191  os << "Expected number of columns in sensor_los: 2.\n";
192  os << "Found: " << sensor_los.ncols();
193  throw runtime_error(os.str());
194  }
195  if (!(sensor_los.nrows() == 1)) {
196  ostringstream os;
197  os << "Expected number of rows in sensor_los: 1.\n";
198  os << "Found: " << sensor_los.nrows();
199  throw runtime_error(os.str());
200  }
201 
202  Ppath ppath_step;
203  Rng rng; //Random Number generator
204  time_t start_time = time(NULL);
205  Index N_se = pnd_field.nbooks(); //Number of scattering elements
206  Vector pnd_vec(
207  N_se); //Vector of particle number densities used at each point
208  Vector Z11maxvector(
209  N_se); //Vector holding the maximum phase function for each
210 
211  // finding maximum phase function for each scat element
212  Index i_total = -1;
213  Index this_f_index = min(f_index, scat_data[0][0].f_grid.nelem());
214  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
215  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
216  i_total++;
217  Z11maxvector[i_total] = max(scat_data[i_ss][i_se].pha_mat_data(
218  this_f_index, joker, joker, joker, joker, joker, 0));
219  }
220  }
221 
222  rng.seed(mc_seed, verbosity);
223  Numeric g, temperature, albedo, g_los_csc_theta;
224  Matrix A(stokes_dim, stokes_dim), Q(stokes_dim, stokes_dim);
225  Matrix evol_op(stokes_dim, stokes_dim), ext_mat_mono(stokes_dim, stokes_dim);
226  Matrix q(stokes_dim, stokes_dim), newQ(stokes_dim, stokes_dim);
227  Matrix Z(stokes_dim, stokes_dim);
228  Matrix R_ant2enu(3, 3),
229  R_stokes(stokes_dim, stokes_dim); // Needed for antenna rotations
230  q = 0.0;
231  newQ = 0.0;
232  Vector vector1(stokes_dim), abs_vec_mono(stokes_dim), I_i(stokes_dim);
233  Vector Isum(stokes_dim), Isquaredsum(stokes_dim);
234  Index termination_flag = 0;
235  const Numeric f_mono = f_grid[f_index];
236  const Numeric prop_dir =
237  -1.0; // propagation direction opposite of los angles
238 
239  CREATE_OUT0;
240 
241  y.resize(stokes_dim);
242  y = 0;
243 
244  mc_iteration_count = 0;
245  mc_error.resize(stokes_dim);
246  mc_points.resize(p_grid.nelem(), lat_grid.nelem(), lon_grid.nelem());
247  mc_points = 0;
248  mc_scat_order.resize(l_mc_scat_order);
249  mc_scat_order = 0;
250  mc_source_domain.resize(4);
251  mc_source_domain = 0;
252 
253  //local versions of workspace
254  Numeric local_surface_skin_t;
255  Matrix local_iy(1, stokes_dim), local_surface_emission(1, stokes_dim);
256  Matrix local_surface_los;
257  Tensor4 local_surface_rmatrix;
258  Vector local_rte_pos(3); // Fixed this (changed from 2 to 3)
259  Vector local_rte_los(2);
260  Vector new_rte_los(2);
261  Index np;
262  Isum = 0.0;
263  Isquaredsum = 0.0;
264  Numeric std_err_i;
265  bool convert_to_rjbt = false;
266  if (iy_unit == "RJBT") {
267  std_err_i = f_mono * f_mono * 2 * BOLTZMAN_CONST / SPEED_OF_LIGHT /
268  SPEED_OF_LIGHT * std_err;
269  convert_to_rjbt = true;
270  } else if (iy_unit == "1") {
271  std_err_i = std_err;
272  } else {
273  ostringstream os;
274  os << "Invalid value for *iy_unit*:" << iy_unit << ".\n"
275  << "This method allows only the options \"RJBT\" and \"1\".";
276  throw runtime_error(os.str());
277  }
278 
279  // Calculate rotation matrix for boresight
280  rotmat_enu(R_ant2enu, sensor_los(0, joker));
281 
282  //Begin Main Loop
283  //
284  bool keepgoing, oksampling;
285  Index nfails = 0;
286  //
287  while (true) {
288  // Complete content of while inside try/catch to handle occasional
289  // failures in the ppath calculations
290  try {
291  bool inside_cloud;
292 
293  mc_iteration_count += 1;
294  Index scattering_order = 0;
295 
296  keepgoing = true; // indicating whether to continue tracing a photon
297  oksampling = true; // gets false if g becomes zero
298 
299  //Sample a FOV direction
300  Matrix R_prop(3, 3);
301  mc_antenna.draw_los(
302  local_rte_los, R_prop, rng, R_ant2enu, sensor_los(0, joker));
303 
304  // Get stokes rotation matrix for rotating polarization
306  R_stokes, stokes_dim, prop_dir, prop_dir, R_prop, R_ant2enu);
307  id_mat(Q);
308  local_rte_pos = sensor_pos(0, joker);
309  I_i = 0.0;
310 
311  while (keepgoing) {
313  evol_op,
314  abs_vec_mono,
315  temperature,
316  ext_mat_mono,
317  rng,
318  local_rte_pos,
319  local_rte_los,
320  pnd_vec,
321  g,
322  ppath_step,
323  termination_flag,
324  inside_cloud,
325  ppath_step_agenda,
326  ppath_lmax,
327  ppath_lraytrace,
328  taustep_limit,
329  propmat_clearsky_agenda,
330  stokes_dim,
331  f_index,
332  f_grid,
333  p_grid,
334  lat_grid,
335  lon_grid,
336  z_field,
337  refellipsoid,
338  z_surface,
339  t_field,
340  vmr_field,
341  cloudbox_limits,
342  pnd_field,
343  scat_data,
344  verbosity);
345 
346  // GH 2011-09-08: if the lowest layer has large
347  // extent and a thick cloud, g may be 0 due to
348  // underflow, but then I_i should be 0 as well.
349  // Don't turn it into nan for no reason.
350  // If reaching underflow, no point in going on;
351  // hence new photon.
352  // GH 2011-09-14: moved this check to outside the different
353  // scenarios, as this goes wrong regardless of the scenario.
354  if (g == 0) {
355  keepgoing = false;
356  oksampling = false;
357  mc_iteration_count -= 1;
358  out0 << "WARNING: A rejected path sampling (g=0)!\n(if this"
359  << "happens repeatedly, try to decrease *ppath_lmax*)";
360  } else if (termination_flag == 1) {
362  local_iy,
363  Vector(1, f_mono),
364  local_rte_pos,
365  local_rte_los,
366  iy_space_agenda);
367  mult(vector1, evol_op, local_iy(0, joker));
368  mult(I_i, Q, vector1);
369  I_i /= g;
370  keepgoing = false; //stop here. New photon.
371  mc_source_domain[0] += 1;
372  } else if (termination_flag == 2) {
373  //Calculate surface properties
375  local_surface_skin_t,
376  local_surface_emission,
377  local_surface_los,
378  local_surface_rmatrix,
379  Vector(1, f_mono),
380  local_rte_pos,
381  local_rte_los,
382  surface_rtprop_agenda);
383 
384  //if( local_surface_los.nrows() > 1 )
385  // throw runtime_error(
386  // "The method handles only specular reflections." );
387 
388  //deal with blackbody case
389  if (local_surface_los.empty()) {
390  mult(vector1, evol_op, local_surface_emission(0, joker));
391  mult(I_i, Q, vector1);
392  I_i /= g;
393  keepgoing = false;
394  mc_source_domain[1] += 1;
395  } else
396  //decide between reflection and emission
397  {
398  const Numeric rnd = rng.draw();
399 
400  Numeric R11 = 0;
401  for (Index i = 0; i < local_surface_rmatrix.nbooks(); i++) {
402  R11 += local_surface_rmatrix(i, 0, 0, 0);
403  }
404 
405  if (rnd > R11) {
406  //then we have emission
407  mult(vector1, evol_op, local_surface_emission(0, joker));
408  mult(I_i, Q, vector1);
409  I_i /= g * (1 - R11);
410  keepgoing = false;
411  mc_source_domain[1] += 1;
412  } else {
413  //we have reflection
414  // determine which reflection los to use
415  Index i = 0;
416  Numeric rsum = local_surface_rmatrix(i, 0, 0, 0);
417  while (rsum < rnd) {
418  i++;
419  rsum += local_surface_rmatrix(i, 0, 0, 0);
420  }
421 
422  local_rte_los = local_surface_los(i, joker);
423 
424  mult(q, evol_op, local_surface_rmatrix(i, 0, joker, joker));
425  mult(newQ, Q, q);
426  Q = newQ;
427  Q /= g * local_surface_rmatrix(i, 0, 0, 0);
428  }
429  }
430  } else if (inside_cloud) {
431  //we have another scattering/emission point
432  //Estimate single scattering albedo
433  albedo = 1 - abs_vec_mono[0] / ext_mat_mono(0, 0);
434 
435  //determine whether photon is emitted or scattered
436  if (rng.draw() > albedo) {
437  //Calculate emission
438  Numeric planck_value = planck(f_mono, temperature);
439  Vector emission = abs_vec_mono;
440  emission *= planck_value;
441  Vector emissioncontri(stokes_dim);
442  mult(emissioncontri, evol_op, emission);
443  emissioncontri /= (g * (1 - albedo)); //yuck!
444  mult(I_i, Q, emissioncontri);
445  keepgoing = false;
446  mc_source_domain[3] += 1;
447  } else {
448  //we have a scattering event
449  Sample_los(new_rte_los,
450  g_los_csc_theta,
451  Z,
452  rng,
453  local_rte_los,
454  scat_data,
455  f_index,
456  stokes_dim,
457  pnd_vec,
458  Z11maxvector,
459  ext_mat_mono(0, 0) - abs_vec_mono[0],
460  temperature,
461  t_interp_order);
462 
463  Z /= g * g_los_csc_theta * albedo;
464 
465  mult(q, evol_op, Z);
466  mult(newQ, Q, q);
467  Q = newQ;
468  scattering_order += 1;
469  local_rte_los = new_rte_los;
470  }
471  } else {
472  //Must be clear sky emission point
473  //Calculate emission
474  Numeric planck_value = planck(f_mono, temperature);
475  Vector emission = abs_vec_mono;
476  emission *= planck_value;
477  Vector emissioncontri(stokes_dim);
478  mult(emissioncontri, evol_op, emission);
479  emissioncontri /= g;
480  mult(I_i, Q, emissioncontri);
481  keepgoing = false;
482  mc_source_domain[2] += 1;
483  }
484  } // keepgoing
485 
486  if (oksampling) {
487  // Set spome of the bookkeeping variables
488  np = ppath_step.np;
489  mc_points(ppath_step.gp_p[np - 1].idx,
490  ppath_step.gp_lat[np - 1].idx,
491  ppath_step.gp_lon[np - 1].idx) += 1;
492  if (scattering_order < l_mc_scat_order) {
493  mc_scat_order[scattering_order] += 1;
494  }
495 
496  // Rotate into antenna polarization frame
497  Vector I_hold(stokes_dim);
498  mult(I_hold, R_stokes, I_i);
499  Isum += I_i;
500 
501  for (Index j = 0; j < stokes_dim; j++) {
502  assert(!std::isnan(I_i[j]));
503  Isquaredsum[j] += I_i[j] * I_i[j];
504  }
505  y = Isum;
506  y /= (Numeric)mc_iteration_count;
507  for (Index j = 0; j < stokes_dim; j++) {
508  mc_error[j] = sqrt(
509  (Isquaredsum[j] / (Numeric)mc_iteration_count - y[j] * y[j]) /
510  (Numeric)mc_iteration_count);
511  }
512  if (std_err > 0 && mc_iteration_count >= min_iter &&
513  mc_error[0] < std_err_i) {
514  break;
515  }
516  if (max_time > 0 && (Index)(time(NULL) - start_time) >= max_time) {
517  break;
518  }
519  if (max_iter > 0 && mc_iteration_count >= max_iter) {
520  break;
521  }
522  }
523  } // Try
524 
525  catch (const std::runtime_error& e) {
526  mc_iteration_count += 1;
527  nfails += 1;
528  out0 << "WARNING: A MC path sampling failed! Error was:\n";
529  cout << e.what() << endl;
530  if (nfails >= 5) {
531  throw runtime_error(
532  "The MC path sampling has failed five times. A few failures "
533  "should be OK, but this number is suspiciously high and the "
534  "reason to these failures should be tracked down.");
535  }
536  }
537  } // while
538 
539  if (convert_to_rjbt) {
540  for (Index j = 0; j < stokes_dim; j++) {
541  y[j] = invrayjean(y[j], f_mono);
542  mc_error[j] = invrayjean(mc_error[j], f_mono);
543  }
544  }
545 }
546 
547 /* Workspace method: Doxygen documentation will be auto-generated */
548 void MCRadar( // Workspace reference:
549  Workspace& ws,
550  // WS Output:
551  Vector& y,
552  Vector& mc_error,
553  // WS Input:
554  const MCAntenna& mc_antenna,
555  const Vector& f_grid,
556  const Index& f_index,
557  const Matrix& sensor_pos,
558  const Matrix& sensor_los,
559  const Index& stokes_dim,
560  const Index& atmosphere_dim,
561  const Numeric& ppath_lmax,
562  const Agenda& ppath_step_agenda,
563  const Numeric& ppath_lraytrace,
564  const Agenda& propmat_clearsky_agenda,
565  const Vector& p_grid,
566  const Vector& lat_grid,
567  const Vector& lon_grid,
568  const Tensor3& z_field,
569  const Vector& refellipsoid,
570  const Matrix& z_surface,
571  const Tensor3& t_field,
572  const Tensor4& vmr_field,
573  const Index& cloudbox_on,
574  const ArrayOfIndex& cloudbox_limits,
575  const Tensor4& pnd_field,
576  const ArrayOfArrayOfSingleScatteringData& scat_data,
577  const Vector& mc_y_tx,
578  const Vector& range_bins,
579  const Index& atmfields_checked,
580  const Index& atmgeom_checked,
581  const Index& scat_data_checked,
582  const Index& cloudbox_checked,
583  const String& iy_unit,
584  const Index& mc_max_scatorder,
585  const Index& mc_seed,
586  const Index& mc_max_iter,
587  const Numeric& ze_tref,
588  const Numeric& k2,
589  const Index& t_interp_order,
590  // Verbosity object:
591  const Verbosity& verbosity) {
592  CREATE_OUT0;
593 
594  // Important constants
595  const Index nbins = range_bins.nelem() - 1;
596  const Numeric r_min = min(range_bins);
597  const Numeric r_max = max(range_bins);
598 
599  // Basics
600  //
601  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
602  if (stokes_dim < 2)
603  throw runtime_error("This method requires that stokes_dim >= 2");
604  if (atmfields_checked != 1)
605  throw runtime_error(
606  "The atmospheric fields must be flagged to have "
607  "passed a consistency check (atmfields_checked=1).");
608  if (atmgeom_checked != 1)
609  throw runtime_error(
610  "The atmospheric geometry must be flagged to have "
611  "passed a consistency check (atmgeom_checked=1).");
612 
613  if (!cloudbox_on)
614  throw runtime_error("The cloudbox must be activated (cloudbox_on=1).");
615  if (cloudbox_checked != 1)
616  throw runtime_error(
617  "The cloudbox must be flagged to have "
618  "passed a consistency check (cloudbox_checked=1).");
619  if (scat_data_checked != 1)
620  throw runtime_error(
621  "The scat_data must be flagged to have "
622  "passed a consistency check (scat_data_checked=1).");
623 
624  if (f_index < 0)
625  throw runtime_error(
626  "The option of f_index < 0 is not handled by this "
627  "method.");
628  if (f_index >= f_grid.nelem())
629  throw runtime_error("*f_index* is outside the range of *f_grid*.");
630 
631  if (atmosphere_dim != 3)
632  throw runtime_error("Only 3D atmospheres are handled.");
633 
634  if (stokes_dim != mc_y_tx.nelem())
635  throw runtime_error("*mc_y_tx* must have size of *stokes_dim*.");
636 
637  if (sensor_pos.ncols() != 3) {
638  ostringstream os;
639  os << "Expected number of columns in sensor_pos: 3.\n";
640  os << "Found: " << sensor_pos.ncols();
641  throw runtime_error(os.str());
642  }
643 
644  if (sensor_los.ncols() != 2) {
645  ostringstream os;
646  os << "Expected number of columns in sensor_los: 2.\n";
647  os << "Found: " << sensor_los.ncols();
648  throw runtime_error(os.str());
649  }
650 
651  if (mc_max_iter < 0)
652  throw runtime_error(
653  "mc_max_iter must be positive, "
654  "as it is the only limiter.");
655 
656  if (!is_increasing(range_bins))
657  throw runtime_error(
658  "The vector *range_bins* must contain strictly "
659  "increasing values.");
660 
661  if (r_min < 0)
662  throw runtime_error(
663  "The vector *range_bins* is not allowed to contain "
664  "negative distance or round-trip time.");
665 
666  if (mc_antenna.get_type() != ANTENNA_TYPE_GAUSSIAN) {
667  throw runtime_error(
668  "MCRadar only works with "
669  "Gaussian antenna patterns.");
670  }
671 
672  Ppath ppath_step;
673  Rng rng; //Random Number generator
674  Index N_se = pnd_field.nbooks(); //Number of scattering elements
675  Vector pnd_vec(
676  N_se); //Vector of particle number densities used at each point
677  bool anyptype_nonTotRan = is_anyptype_nonTotRan(scat_data);
678  bool is_dist = max(range_bins) > 1; // Is it round trip time or distance
679  rng.seed(mc_seed, verbosity);
680  Numeric ppath_lraytrace_var;
681  //Numeric temperature, albedo;
682  Numeric albedo;
683  Numeric Csca, Cext;
684  Numeric antenna_wgt;
685  Matrix evol_op(stokes_dim, stokes_dim), ext_mat_mono(stokes_dim, stokes_dim);
686  Matrix trans_mat(stokes_dim, stokes_dim);
687  Matrix Z(stokes_dim, stokes_dim);
688  Matrix R_ant2enu(3, 3), R_enu2ant(3, 3), R_stokes(stokes_dim, stokes_dim);
689  Vector abs_vec_mono(stokes_dim), I_i(stokes_dim), I_i_rot(stokes_dim);
690  Vector Isum(nbins * stokes_dim), Isquaredsum(nbins * stokes_dim);
691  Index termination_flag = 0;
692  Vector bin_height(nbins);
693  Vector range_bin_count(nbins);
694  Index mc_iter;
695  Index scat_order;
696 
697  // allocating variables needed for pha_mat extraction (don't want to do this
698  // in every loop step again).
699  ArrayOfArrayOfTensor6 pha_mat_Nse;
700  ArrayOfArrayOfIndex ptypes_Nse;
701  Matrix t_ok;
702  ArrayOfTensor6 pha_mat_ssbulk;
703  ArrayOfIndex ptype_ssbulk;
704  Tensor6 pha_mat_bulk;
705  Index ptype_bulk;
706  Matrix pdir_array(1, 2), idir_array(1, 2);
707  Vector t_array(1);
708  Matrix pnds(N_se, 1);
709 
710  // for pha_mat handling, at the moment we still need scat_data_mono. Hence,
711  // extract that here (but in its local container, not into the WSV
712  // scat_data_mono).
713  //ArrayOfArrayOfSingleScatteringData this_scat_data_mono;
714  //scat_data_monoExtract( this_scat_data_mono, scat_data, f_index, verbosity );
715 
716  const Numeric f_mono = f_grid[f_index];
717  const Numeric tx_dir = 1.0;
718  const Numeric rx_dir = -1.0;
719 
720  for (Index ibin = 0; ibin < nbins; ibin++) {
721  bin_height[ibin] = range_bins[ibin + 1] - range_bins[ibin];
722  }
723  if (!is_dist) {
724  bin_height *= 0.5 * SPEED_OF_LIGHT;
725  }
726 
727  y.resize(nbins * stokes_dim);
728  y = 0;
729 
730  range_bin_count = 0;
731 
732  mc_iter = 0;
733  // this will need to be reshaped differently for range gates
734  mc_error.resize(stokes_dim * nbins);
735 
736  //local versions of workspace
737  Matrix local_iy(1, stokes_dim), local_surface_emission(1, stokes_dim);
738  Matrix local_surface_los;
739  Tensor4 local_surface_rmatrix;
740  Vector local_rte_pos(3);
741  Vector local_rte_los(2);
742  Vector new_rte_los(2);
743  Vector Ipath(stokes_dim), Ihold(stokes_dim), Ipath_norm(stokes_dim);
744  Isum = 0.0;
745  Isquaredsum = 0.0;
746  Numeric s_tot, s_return; // photon distance traveled
747  Numeric t_tot, t_return; // photon time traveled
748  Numeric r_trav, r_bin; // range traveled (1-way distance) or round-trip time
749 
750  Numeric fac;
751  if (iy_unit == "1") {
752  fac = 1.0;
753  }
754  // Conversion from intensity to reflectivity
755  else if (iy_unit == "Ze") {
756  Vector cfac(1);
757  ze_cfac(cfac, Vector(1, f_mono), ze_tref, k2);
758  // Due to different definitions, the factor shall here be scaled with 1/(2pi)
759  fac = cfac[0] / (2 * PI);
760  } else {
761  ostringstream os;
762  os << "Invalid value for *iy_unit*:" << iy_unit << ".\n"
763  << "This method allows only the options \"Ze\" and \"1\".";
764  throw runtime_error(os.str());
765  }
766 
767  // Calculate rotation matrix and polarization bases for boresight
768  rotmat_enu(R_ant2enu, sensor_los(0, joker));
769  R_enu2ant = transpose(R_ant2enu);
770 
771  //Begin Main Loop
772  bool keepgoing, firstpass, integrity;
773  while (mc_iter < mc_max_iter) {
774  bool inside_cloud;
775 
776  mc_iter += 1;
777 
778  integrity = true; // intensity is not nan or below threshold
779  keepgoing = true; // indicating whether to continue tracing a photon
780  firstpass = true; // ensure backscatter is properly calculated
781 
782  //Sample a FOV direction
783  Matrix R_tx(3, 3);
784  mc_antenna.draw_los(
785  local_rte_los, R_tx, rng, R_ant2enu, sensor_los(0, joker));
786  rotmat_stokes(R_stokes, stokes_dim, tx_dir, tx_dir, R_ant2enu, R_tx);
787  mult(Ihold, R_stokes, mc_y_tx);
788 
789  // Initialize other variables
790  local_rte_pos = sensor_pos(0, joker);
791  s_tot = 0.0;
792  t_tot = 0.0;
793  scat_order = 0;
794  while (keepgoing) {
795  Numeric s_path, t_path;
796 
797  mcPathTraceRadar(ws,
798  evol_op,
799  abs_vec_mono,
800  t_array[0],
801  ext_mat_mono,
802  rng,
803  local_rte_pos,
804  local_rte_los,
805  pnd_vec, //pnds(joker,0),
806  s_path,
807  t_path,
808  ppath_step,
809  termination_flag,
810  inside_cloud,
811  ppath_step_agenda,
812  ppath_lmax,
813  ppath_lraytrace,
814  propmat_clearsky_agenda,
815  anyptype_nonTotRan,
816  stokes_dim,
817  f_index,
818  f_grid,
819  Ihold,
820  p_grid,
821  lat_grid,
822  lon_grid,
823  z_field,
824  refellipsoid,
825  z_surface,
826  t_field,
827  vmr_field,
828  cloudbox_limits,
829  pnd_field,
830  scat_data,
831  verbosity);
832  pnds(joker, 0) = pnd_vec;
833  if (!inside_cloud || termination_flag != 0) {
834  keepgoing = false;
835  } else {
836  s_tot += s_path;
837  t_tot += t_path;
838 
839  //
840  Csca = ext_mat_mono(0, 0) - abs_vec_mono[0];
841  Cext = ext_mat_mono(0, 0);
842  if (anyptype_nonTotRan) {
843  const Numeric Irat = Ihold[1] / Ihold[0];
844  Csca += Irat * (ext_mat_mono(1, 0) - abs_vec_mono[1]);
845  Cext += Irat * ext_mat_mono(0, 1);
846  }
847  albedo = Csca / Cext;
848 
849  // Terminate if absorption event, outside cloud, or surface
850  Numeric rn = rng.draw();
851  if (rn > albedo) {
852  keepgoing = false;
853  continue;
854  }
855 
856  Vector rte_los_geom(2);
857 
858  // Compute reflectivity contribution based on local-to-sensor
859  // geometry, path attenuation
860  // Get los angles at atmospheric locale to determine
861  // scattering angle
862  if (firstpass) {
863  // Use this to ensure that the difference in azimuth angle
864  // between incident and scattered lines-of-sight is 180
865  // degrees
866  mirror_los(rte_los_geom, local_rte_los, atmosphere_dim);
867  firstpass = false;
868  } else {
869  // Replace with ppath_agendaExecute??
871  atmosphere_dim,
872  lat_grid,
873  lon_grid,
874  refellipsoid,
875  local_rte_pos,
876  sensor_pos(0, joker),
877  verbosity);
878  }
879 
880  // Get los angles at sensor to determine antenna pattern
881  // weighting of return signal and ppath to determine
882  // propagation path back to sensor
883  // Replace with ppath_agendaExecute??
884  Ppath ppath;
885  Vector rte_los_antenna(2);
886  ppath_lraytrace_var = ppath_lraytrace;
887  Numeric za_accuracy = 2e-5;
888  Numeric pplrt_factor = 5;
889  Numeric pplrt_lowest = 0.5;
890 
891  rte_losGeometricFromRtePosToRtePos2(rte_los_antenna,
892  atmosphere_dim,
893  lat_grid,
894  lon_grid,
895  refellipsoid,
896  sensor_pos(0, joker),
897  local_rte_pos,
898  verbosity);
899 
900  ppathFromRtePos2(ws,
901  ppath,
902  rte_los_antenna,
903  ppath_lraytrace_var,
904  ppath_step_agenda,
905  atmosphere_dim,
906  p_grid,
907  lat_grid,
908  lon_grid,
909  z_field,
910  f_grid,
911  refellipsoid,
912  z_surface,
913  sensor_pos(0, joker),
914  local_rte_pos,
915  ppath_lmax,
916  za_accuracy,
917  pplrt_factor,
918  pplrt_lowest,
919  verbosity);
920 
921  // Return distance
922  const Index np2 = ppath.np;
923  s_return = ppath.end_lstep;
924  t_return = s_return / SPEED_OF_LIGHT;
925  for (Index ip = 1; ip < np2; ip++) {
926  s_return += ppath.lstep[ip - 1];
927  t_return += ppath.lstep[ip - 1] * 0.5 *
928  (ppath.ngroup[ip - 1] + ppath.ngroup[ip]) /
930  }
931 
932  // One-way distance
933  if (is_dist) {
934  r_trav = 0.5 * (s_tot + s_return);
935  }
936 
937  // Round trip travel time
938  else {
939  r_trav = t_tot + t_return;
940  }
941 
942  // Still within max range of radar?
943  if (r_trav <= r_max) {
944  // Compute path extinction as with radio link
946  trans_mat,
947  ppath,
948  propmat_clearsky_agenda,
949  stokes_dim,
950  f_index,
951  f_grid,
952  p_grid,
953  t_field,
954  vmr_field,
955  cloudbox_limits,
956  pnd_field,
957  scat_data,
958  verbosity);
959 
960  // Obtain scattering matrix given incident and scattered angles
961  Matrix P(stokes_dim, stokes_dim);
962 
963  pdir_array(0, joker) = rte_los_geom;
964  idir_array(0, joker) = local_rte_los;
965  pha_mat_NScatElems(pha_mat_Nse,
966  ptypes_Nse,
967  t_ok,
968  scat_data,
969  stokes_dim,
970  t_array,
971  pdir_array,
972  idir_array,
973  f_index,
974  t_interp_order);
975  pha_mat_ScatSpecBulk(pha_mat_ssbulk,
976  ptype_ssbulk,
977  pha_mat_Nse,
978  ptypes_Nse,
979  pnds,
980  t_ok);
981  pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
982  P = pha_mat_bulk(0, 0, 0, 0, joker, joker);
983 
984  P *= 4 * PI;
985  P /= Csca;
986 
987  // Compute reflectivity contribution here
988  mult(Ipath, evol_op, Ihold);
989  Ipath /= Ipath[0];
990  Ipath *= Ihold[0];
991  mult(Ihold, P, Ipath);
992  mult(I_i, trans_mat, Ihold);
993  Ihold = Ipath;
994  if (Ihold[0] < 1e-40 || std::isnan(Ihold[0]) ||
995  std::isnan(Ihold[1]) ||
996  (stokes_dim > 2 && std::isnan(Ihold[2])) ||
997  (stokes_dim > 3 && std::isnan(Ihold[3]))) {
998  integrity = false;
999  }
1000 
1001  if (r_trav > r_min && integrity) {
1002  // Add reflectivity to proper range bin
1003  Index ibin = 0;
1004  r_bin = 0.0;
1005  while (r_bin < r_trav && ibin <= nbins + 1) {
1006  ibin++;
1007  r_bin = range_bins[ibin];
1008  }
1009  ibin -= 1;
1010 
1011  // Calculate rx antenna weight and polarization rotation
1012  Matrix R_rx(3, 3);
1013  rotmat_enu(R_rx, rte_los_antenna);
1014  mc_antenna.return_los(antenna_wgt, R_rx, R_enu2ant);
1015  rotmat_stokes(
1016  R_stokes, stokes_dim, rx_dir, tx_dir, R_rx, R_ant2enu);
1017  mult(I_i_rot, R_stokes, I_i);
1018 
1019  for (Index istokes = 0; istokes < stokes_dim; istokes++) {
1020  Index ibiny = ibin * stokes_dim + istokes;
1021  assert(!std::isnan(I_i_rot[istokes]));
1022  Isum[ibiny] += antenna_wgt * I_i_rot[istokes];
1023  Isquaredsum[ibiny] += antenna_wgt * antenna_wgt *
1024  I_i_rot[istokes] * I_i_rot[istokes];
1025  }
1026  range_bin_count[ibin] += 1;
1027  }
1028 
1029  scat_order++;
1030 
1031  Sample_los_uniform(new_rte_los, rng);
1032  pdir_array(0, joker) = new_rte_los;
1033  // alt:
1034  // Sample_los_uniform( pdir_array(0,joker), rng );
1035  pha_mat_NScatElems(pha_mat_Nse,
1036  ptypes_Nse,
1037  t_ok,
1038  scat_data,
1039  stokes_dim,
1040  t_array,
1041  pdir_array,
1042  idir_array,
1043  f_index,
1044  t_interp_order);
1045  pha_mat_ScatSpecBulk(pha_mat_ssbulk,
1046  ptype_ssbulk,
1047  pha_mat_Nse,
1048  ptypes_Nse,
1049  pnds,
1050  t_ok);
1051  pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
1052  Z = pha_mat_bulk(0, 0, 0, 0, joker, joker);
1053 
1054  Z *= 4 * PI;
1055  Z /= Csca;
1056  mult(Ipath, Z, Ihold);
1057  Ihold = Ipath;
1058  local_rte_los = new_rte_los;
1059  // alt:
1060  //local_rte_los = pdir_array(0,joker);
1061  // or even (but also requires replacements of local_rte_los
1062  // with idir_array throughout the whole loop):
1063  //idir_array = pdir_array;
1064  } else {
1065  // Past farthest range
1066  keepgoing = false;
1067  }
1068  }
1069 
1070  // Some checks
1071  if (scat_order >= mc_max_scatorder) keepgoing = false;
1072  if (!integrity) keepgoing = false;
1073  } // while (inner: keepgoing)
1074 
1075  } // while (outer)
1076 
1077  // Normalize range bins and apply sensor response (polarization)
1078  for (Index ibin = 0; ibin < nbins; ibin++) {
1079  for (Index istokes = 0; istokes < stokes_dim; istokes++) {
1080  Index ibiny = ibin * stokes_dim + istokes;
1081  if (range_bin_count[ibin] > 0) {
1082  y[ibiny] = Isum[ibiny] / ((Numeric)mc_iter) / bin_height[ibin];
1083  mc_error[ibiny] = sqrt((Isquaredsum[ibiny] / (Numeric)mc_iter /
1084  bin_height[ibin] / bin_height[ibin] -
1085  y[ibiny] * y[ibiny]) /
1086  (Numeric)mc_iter);
1087  }
1088  }
1089  }
1090 
1091  y *= fac;
1092  mc_error *= fac;
1093 } // end MCRadar
1094 
1095 /* Workspace method: Doxygen documentation will be auto-generated */
1096 void MCSetSeedFromTime(Index& mc_seed, const Verbosity&) {
1097  mc_seed = (Index)time(NULL);
1098 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void rte_losGeometricFromRtePosToRtePos2(Vector &rte_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Vector &rte_pos, const Vector &rte_pos2, const Verbosity &)
WORKSPACE METHOD: rte_losGeometricFromRtePosToRtePos2.
Definition: m_ppath.cc:1476
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
void surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:25015
const Numeric DEG2RAD
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:195
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:426
Interpolation classes and functions created for use within Monte Carlo scattering simulations...
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.
Declarations having to do with the four output streams.
void set_gaussian(const Numeric &za_sigma, const Numeric &aa_sigma)
set_gaussian.
Definition: mc_antenna.cc:119
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
The Vector class.
Definition: matpackI.h:860
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
The Tensor4 class.
Definition: matpackIV.h:421
void rotmat_stokes(MatrixView R_pra, const Index &stokes_dim, const Numeric &f1_dir, const Numeric &f2_dir, ConstMatrixView R_f1, ConstMatrixView R_f2)
rotmat_stokes.
Definition: mc_antenna.cc:86
Vector lstep
The length between ppath points.
Definition: ppath.h:70
Linear algebra functions.
void ppathFromRtePos2(Workspace &ws, Ppath &ppath, Vector &rte_los, Numeric &ppath_lraytrace, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Vector &rte_pos, const Vector &rte_pos2, const Numeric &ppath_lmax, const Numeric &za_accuracy, const Numeric &pplrt_factor, const Numeric &pplrt_lowest, const Verbosity &verbosity)
WORKSPACE METHOD: ppathFromRtePos2.
Definition: m_ppath.cc:306
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
This file contains basic functions to handle XML data files.
void mc_antennaSetPencilBeam(MCAntenna &mc_antenna, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetPencilBeam.
Definition: m_montecarlo.cc:82
Vector ngroup
The group index of refraction.
Definition: ppath.h:80
#define min(a, b)
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
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
void MCRadar(Workspace &ws, Vector &y, Vector &mc_error, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Numeric &ppath_lmax, const Agenda &ppath_step_agenda, const Numeric &ppath_lraytrace, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Vector &mc_y_tx, const Vector &range_bins, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &cloudbox_checked, const String &iy_unit, const Index &mc_max_scatorder, const Index &mc_seed, const Index &mc_max_iter, const Numeric &ze_tref, const Numeric &k2, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCRadar.
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, ArrayOfIndex &mc_source_domain, ArrayOfIndex &mc_scat_order, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &cloudbox_checked, const String &iy_unit, const Index &mc_seed, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Numeric &taustep_limit, const Index &l_mc_scat_order, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
Definition: m_montecarlo.cc:87
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
Definition: rte.cc:2736
void mc_antennaSetGaussianByFWHM(MCAntenna &mc_antenna, const Numeric &za_fwhm, const Numeric &aa_fwhm, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetGaussianByFWHM.
Definition: m_montecarlo.cc:73
const Numeric SPEED_OF_LIGHT
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
The Tensor3 class.
Definition: matpackIII.h:339
const Numeric PI
The global header file for ARTS.
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.
const Numeric RAD2DEG
_CS_string_type str() const
Definition: sstream.h:491
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:24281
const Joker joker
member functions of the Rng class and gsl_rng code
Numeric vector1(const StokesVector &a, const ConstVectorView &B, const StokesVector &da, const ConstVectorView &dB_dT, const StokesVector &dS, bool dT, Index i) noexcept
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
const Numeric BOLTZMAN_CONST
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
void draw_los(VectorView sampled_rte_los, MatrixView R_los, Rng &rng, ConstMatrixView R_ant2enu, ConstVectorView bore_sight_los) const
draw_los.
Definition: mc_antenna.cc:180
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:51
Implementation of Matrix, Vector, and such stuff.
Header file for special_interp.cc.
Propagation path structure and functions.
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
Header file for logic.cc.
This can be used to make arrays out of anything.
Definition: array.h:40
void mcPathTraceGeneral(Workspace &ws, MatrixView evol_op, Vector &abs_vec_mono, Numeric &temperature, MatrixView ext_mat_mono, Rng &rng, Vector &rte_pos, Vector &rte_los, Vector &pnd_vec, Numeric &g, Ppath &ppath_step, Index &termination_flag, bool &inside_cloud, const Agenda &ppath_step_agenda, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Numeric &f_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Verbosity &verbosity)
mcPathTraceGeneral
Definition: mc_NotUsed.cc:150
Definition: rng.h:554
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
void Sample_los_uniform(VectorView new_rte_los, Rng &rng)
Sample_los_uniform.
Definition: montecarlo.cc:1471
void return_los(Numeric &wgt, ConstMatrixView R_return, ConstMatrixView R_enu2ant) const
return_los
Definition: mc_antenna.cc:143
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
#define max(a, b)
The Tensor6 class.
Definition: matpackVI.h:1088
Index np
Number of points describing the ppath.
Definition: ppath.h:52
AntennaType get_type() const
AntennaType get_type.
Definition: mc_antenna.cc:141
void mc_antennaSetGaussian(MCAntenna &mc_antenna, const Numeric &za_sigma, const Numeric &aa_sigma, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetGaussian.
Definition: m_montecarlo.cc:64
Numeric planck(const Numeric &f, const Numeric &t)
planck
Workspace class.
Definition: workspace_ng.h:40
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
bool is_anyptype_nonTotRan(const ArrayOfArrayOfSingleScatteringData &scat_data)
is_anyptype_nonTotRan.
Definition: montecarlo.cc:694
#define q
Refraction functions.
#define CREATE_OUT0
Definition: messages.h:204
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
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
rotmat_enu.
Definition: mc_antenna.cc:63
void seed(unsigned long int n, const Verbosity &verbosity)
Seeds the Rng with the integer argument.
Definition: rng.cc:54
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
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
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
set_gaussian_fwhm.
Definition: mc_antenna.cc:125
This file contains declerations of functions of physical character.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
Declaration of functions in rte.cc.
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620