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