ARTS 2.5.0 (git: 9ee3ac6c)
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
53extern const Numeric DEG2RAD;
54extern const Numeric RAD2DEG;
55extern const Numeric PI;
56extern const Numeric BOLTZMAN_CONST;
57extern 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 */
82void 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
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 ARTS_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
548
549/* Workspace method: Doxygen documentation will be auto-generated */
550void MCRadar( // Workspace reference:
551 Workspace& ws,
552 // WS Output:
553 Vector& y,
554 Vector& mc_error,
555 // WS Input:
556 const MCAntenna& mc_antenna,
557 const Vector& f_grid,
558 const Index& f_index,
559 const Matrix& sensor_pos,
560 const Matrix& sensor_los,
561 const Index& stokes_dim,
562 const Index& atmosphere_dim,
563 const Numeric& ppath_lmax,
564 const Agenda& ppath_step_agenda,
565 const Numeric& ppath_lraytrace,
566 const Agenda& propmat_clearsky_agenda,
567 const Vector& p_grid,
568 const Vector& lat_grid,
569 const Vector& lon_grid,
570 const Tensor3& z_field,
571 const Vector& refellipsoid,
572 const Matrix& z_surface,
573 const Tensor3& t_field,
574 const Tensor4& vmr_field,
575 const Index& cloudbox_on,
576 const ArrayOfIndex& cloudbox_limits,
577 const Tensor4& pnd_field,
578 const ArrayOfArrayOfSingleScatteringData& scat_data,
579 const Vector& mc_y_tx,
580 const Vector& range_bins,
581 const Index& atmfields_checked,
582 const Index& atmgeom_checked,
583 const Index& scat_data_checked,
584 const Index& cloudbox_checked,
585 const String& iy_unit_radar,
586 const Index& mc_max_scatorder,
587 const Index& mc_seed,
588 const Index& mc_max_iter,
589 const Numeric& ze_tref,
590 const Numeric& k2,
591 const Index& t_interp_order,
592 // Verbosity object:
593 const Verbosity& verbosity) {
595
596 // Important constants
597 const Index nbins = range_bins.nelem() - 1;
598 const Numeric r_min = min(range_bins);
599 const Numeric r_max = max(range_bins);
600
601 // Basics
602 //
603 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
604 if (stokes_dim < 2)
605 throw runtime_error("This method requires that stokes_dim >= 2");
606 if (atmfields_checked != 1)
607 throw runtime_error(
608 "The atmospheric fields must be flagged to have "
609 "passed a consistency check (atmfields_checked=1).");
610 if (atmgeom_checked != 1)
611 throw runtime_error(
612 "The atmospheric geometry must be flagged to have "
613 "passed a consistency check (atmgeom_checked=1).");
614
615 if (!cloudbox_on)
616 throw runtime_error("The cloudbox must be activated (cloudbox_on=1).");
617 if (cloudbox_checked != 1)
618 throw runtime_error(
619 "The cloudbox must be flagged to have "
620 "passed a consistency check (cloudbox_checked=1).");
621 if (scat_data_checked != 1)
622 throw runtime_error(
623 "The scat_data must be flagged to have "
624 "passed a consistency check (scat_data_checked=1).");
625
626 if (f_index < 0)
627 throw runtime_error(
628 "The option of f_index < 0 is not handled by this "
629 "method.");
630 if (f_index >= f_grid.nelem())
631 throw runtime_error("*f_index* is outside the range of *f_grid*.");
632
633 if (atmosphere_dim != 3)
634 throw runtime_error("Only 3D atmospheres are handled.");
635
636 if (stokes_dim != mc_y_tx.nelem())
637 throw runtime_error("*mc_y_tx* must have size of *stokes_dim*.");
638
639 if (sensor_pos.ncols() != 3) {
640 ostringstream os;
641 os << "Expected number of columns in sensor_pos: 3.\n";
642 os << "Found: " << sensor_pos.ncols();
643 throw runtime_error(os.str());
644 }
645
646 if (sensor_los.ncols() != 2) {
647 ostringstream os;
648 os << "Expected number of columns in sensor_los: 2.\n";
649 os << "Found: " << sensor_los.ncols();
650 throw runtime_error(os.str());
651 }
652
653 if (mc_max_iter < 0)
654 throw runtime_error(
655 "mc_max_iter must be positive, "
656 "as it is the only limiter.");
657
658 if (!is_increasing(range_bins))
659 throw runtime_error(
660 "The vector *range_bins* must contain strictly "
661 "increasing values.");
662
663 if (r_min < 0)
664 throw runtime_error(
665 "The vector *range_bins* is not allowed to contain "
666 "negative distance or round-trip time.");
667
668 if (mc_antenna.get_type() != ANTENNA_TYPE_GAUSSIAN) {
669 throw runtime_error(
670 "MCRadar only works with "
671 "Gaussian antenna patterns.");
672 }
673
674 Ppath ppath_step;
675 Rng rng; //Random Number generator
676 Index N_se = pnd_field.nbooks(); //Number of scattering elements
677 Vector pnd_vec(
678 N_se); //Vector of particle number densities used at each point
679 bool anyptype_nonTotRan = is_anyptype_nonTotRan(scat_data);
680 bool is_dist = max(range_bins) > 1; // Is it round trip time or distance
681 rng.seed(mc_seed, verbosity);
682 Numeric ppath_lraytrace_var;
683 //Numeric temperature, albedo;
684 Numeric albedo;
685 Numeric Csca, Cext;
686 Numeric antenna_wgt;
687 Matrix evol_op(stokes_dim, stokes_dim), ext_mat_mono(stokes_dim, stokes_dim);
688 Matrix trans_mat(stokes_dim, stokes_dim);
689 Matrix Z(stokes_dim, stokes_dim);
690 Matrix R_ant2enu(3, 3), R_enu2ant(3, 3), R_stokes(stokes_dim, stokes_dim);
691 Vector abs_vec_mono(stokes_dim), I_i(stokes_dim), I_i_rot(stokes_dim);
692 Vector Isum(nbins * stokes_dim), Isquaredsum(nbins * stokes_dim);
693 Index termination_flag = 0;
694 Vector bin_height(nbins);
695 Vector range_bin_count(nbins);
696 Index mc_iter;
697 Index scat_order;
698
699 // allocating variables needed for pha_mat extraction (don't want to do this
700 // in every loop step again).
701 ArrayOfArrayOfTensor6 pha_mat_Nse;
702 ArrayOfArrayOfIndex ptypes_Nse;
703 Matrix t_ok;
704 ArrayOfTensor6 pha_mat_ssbulk;
705 ArrayOfIndex ptype_ssbulk;
706 Tensor6 pha_mat_bulk;
707 Index ptype_bulk;
708 Matrix pdir_array(1, 2), idir_array(1, 2);
709 Vector t_array(1);
710 Matrix pnds(N_se, 1);
711
712 // for pha_mat handling, at the moment we still need scat_data_mono. Hence,
713 // extract that here (but in its local container, not into the WSV
714 // scat_data_mono).
715 //ArrayOfArrayOfSingleScatteringData this_scat_data_mono;
716 //scat_data_monoExtract( this_scat_data_mono, scat_data, f_index, verbosity );
717
718 const Numeric f_mono = f_grid[f_index];
719 const Numeric tx_dir = 1.0;
720 const Numeric rx_dir = -1.0;
721
722 for (Index ibin = 0; ibin < nbins; ibin++) {
723 bin_height[ibin] = range_bins[ibin + 1] - range_bins[ibin];
724 }
725 if (!is_dist) {
726 bin_height *= 0.5 * SPEED_OF_LIGHT;
727 }
728
729 y.resize(nbins * stokes_dim);
730 y = 0;
731
732 range_bin_count = 0;
733
734 mc_iter = 0;
735 // this will need to be reshaped differently for range gates
736 mc_error.resize(stokes_dim * nbins);
737
738 //local versions of workspace
739 Matrix local_iy(1, stokes_dim), local_surface_emission(1, stokes_dim);
740 Matrix local_surface_los;
741 Tensor4 local_surface_rmatrix;
742 Vector local_rte_pos(3);
743 Vector local_rte_los(2);
744 Vector new_rte_los(2);
745 Vector Ipath(stokes_dim), Ihold(stokes_dim), Ipath_norm(stokes_dim);
746 Isum = 0.0;
747 Isquaredsum = 0.0;
748 Numeric s_tot, s_return; // photon distance traveled
749 Numeric t_tot, t_return; // photon time traveled
750 Numeric r_trav, r_bin; // range traveled (1-way distance) or round-trip time
751
752 Numeric fac;
753 if (iy_unit_radar == "1") {
754 fac = 1.0;
755 }
756 // Conversion from intensity to reflectivity
757 else if (iy_unit_radar == "Ze") {
758 Vector cfac(1);
759 ze_cfac(cfac, Vector(1, f_mono), ze_tref, k2);
760 // Due to different definitions, the factor shall here be scaled with 1/(2pi)
761 fac = cfac[0] / (2 * PI);
762 } else {
763 ostringstream os;
764 os << "Invalid value for *iy_unit_radar*:" << iy_unit_radar << ".\n"
765 << "This method allows only the options \"Ze\" and \"1\".";
766 throw runtime_error(os.str());
767 }
768
769 // Calculate rotation matrix and polarization bases for boresight
770 rotmat_enu(R_ant2enu, sensor_los(0, joker));
771 R_enu2ant = transpose(R_ant2enu);
772
773 //Begin Main Loop
774 bool keepgoing, firstpass, integrity;
775 while (mc_iter < mc_max_iter) {
776 bool inside_cloud;
777
778 mc_iter += 1;
779
780 integrity = true; // intensity is not nan or below threshold
781 keepgoing = true; // indicating whether to continue tracing a photon
782 firstpass = true; // ensure backscatter is properly calculated
783
784 //Sample a FOV direction
785 Matrix R_tx(3, 3);
786 mc_antenna.draw_los(
787 local_rte_los, R_tx, rng, R_ant2enu, sensor_los(0, joker));
788 rotmat_stokes(R_stokes, stokes_dim, tx_dir, tx_dir, R_ant2enu, R_tx);
789 mult(Ihold, R_stokes, mc_y_tx);
790
791 // Initialize other variables
792 local_rte_pos = sensor_pos(0, joker);
793 s_tot = 0.0;
794 t_tot = 0.0;
795 scat_order = 0;
796 while (keepgoing) {
797 Numeric s_path, t_path;
798
800 evol_op,
801 abs_vec_mono,
802 t_array[0],
803 ext_mat_mono,
804 rng,
805 local_rte_pos,
806 local_rte_los,
807 pnd_vec, //pnds(joker,0),
808 s_path,
809 t_path,
810 ppath_step,
811 termination_flag,
812 inside_cloud,
813 ppath_step_agenda,
814 ppath_lmax,
815 ppath_lraytrace,
816 propmat_clearsky_agenda,
817 anyptype_nonTotRan,
818 stokes_dim,
819 f_index,
820 f_grid,
821 Ihold,
822 p_grid,
823 lat_grid,
824 lon_grid,
825 z_field,
826 refellipsoid,
827 z_surface,
828 t_field,
829 vmr_field,
830 cloudbox_limits,
831 pnd_field,
832 scat_data,
833 verbosity);
834 pnds(joker, 0) = pnd_vec;
835 if (!inside_cloud || termination_flag != 0) {
836 keepgoing = false;
837 } else {
838 s_tot += s_path;
839 t_tot += t_path;
840
841 //
842 Csca = ext_mat_mono(0, 0) - abs_vec_mono[0];
843 Cext = ext_mat_mono(0, 0);
844 if (anyptype_nonTotRan) {
845 const Numeric Irat = Ihold[1] / Ihold[0];
846 Csca += Irat * (ext_mat_mono(1, 0) - abs_vec_mono[1]);
847 Cext += Irat * ext_mat_mono(0, 1);
848 }
849 albedo = Csca / Cext;
850
851 // Terminate if absorption event, outside cloud, or surface
852 Numeric rn = rng.draw();
853 if (rn > albedo) {
854 keepgoing = false;
855 continue;
856 }
857
858 Vector rte_los_geom(2);
859
860 // Compute reflectivity contribution based on local-to-sensor
861 // geometry, path attenuation
862 // Get los angles at atmospheric locale to determine
863 // scattering angle
864 if (firstpass) {
865 // Use this to ensure that the difference in azimuth angle
866 // between incident and scattered lines-of-sight is 180
867 // degrees
868 mirror_los(rte_los_geom, local_rte_los, atmosphere_dim);
869 firstpass = false;
870 } else {
871 // Replace with ppath_agendaExecute??
873 atmosphere_dim,
874 lat_grid,
875 lon_grid,
876 refellipsoid,
877 local_rte_pos,
878 sensor_pos(0, joker),
879 verbosity);
880 }
881
882 // Get los angles at sensor to determine antenna pattern
883 // weighting of return signal and ppath to determine
884 // propagation path back to sensor
885 // Replace with ppath_agendaExecute??
886 Ppath ppath;
887 Vector rte_los_antenna(2);
888 ppath_lraytrace_var = ppath_lraytrace;
889 Numeric za_accuracy = 2e-5;
890 Numeric pplrt_factor = 5;
891 Numeric pplrt_lowest = 0.5;
892
894 atmosphere_dim,
895 lat_grid,
896 lon_grid,
897 refellipsoid,
898 sensor_pos(0, joker),
899 local_rte_pos,
900 verbosity);
901
903 ppath,
904 rte_los_antenna,
905 ppath_lraytrace_var,
906 ppath_step_agenda,
907 atmosphere_dim,
908 p_grid,
909 lat_grid,
910 lon_grid,
911 z_field,
912 f_grid,
913 refellipsoid,
914 z_surface,
915 sensor_pos(0, joker),
916 local_rte_pos,
917 ppath_lmax,
918 za_accuracy,
919 pplrt_factor,
920 pplrt_lowest,
921 verbosity);
922
923 // Return distance
924 const Index np2 = ppath.np;
925 s_return = ppath.end_lstep;
926 t_return = s_return / SPEED_OF_LIGHT;
927 for (Index ip = 1; ip < np2; ip++) {
928 s_return += ppath.lstep[ip - 1];
929 t_return += ppath.lstep[ip - 1] * 0.5 *
930 (ppath.ngroup[ip - 1] + ppath.ngroup[ip]) /
932 }
933
934 // One-way distance
935 if (is_dist) {
936 r_trav = 0.5 * (s_tot + s_return);
937 }
938
939 // Round trip travel time
940 else {
941 r_trav = t_tot + t_return;
942 }
943
944 // Still within max range of radar?
945 if (r_trav <= r_max) {
946 // Compute path extinction as with radio link
948 trans_mat,
949 ppath,
950 propmat_clearsky_agenda,
951 stokes_dim,
952 f_index,
953 f_grid,
954 p_grid,
955 t_field,
956 vmr_field,
957 cloudbox_limits,
958 pnd_field,
959 scat_data,
960 verbosity);
961
962 // Obtain scattering matrix given incident and scattered angles
963 Matrix P(stokes_dim, stokes_dim);
964
965 pdir_array(0, joker) = rte_los_geom;
966 idir_array(0, joker) = local_rte_los;
967 pha_mat_NScatElems(pha_mat_Nse,
968 ptypes_Nse,
969 t_ok,
970 scat_data,
971 stokes_dim,
972 t_array,
973 pdir_array,
974 idir_array,
975 f_index,
976 t_interp_order);
977 pha_mat_ScatSpecBulk(pha_mat_ssbulk,
978 ptype_ssbulk,
979 pha_mat_Nse,
980 ptypes_Nse,
981 pnds,
982 t_ok);
983 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
984 P = pha_mat_bulk(0, 0, 0, 0, joker, joker);
985
986 P *= 4 * PI;
987 P /= Csca;
988
989 // Compute reflectivity contribution here
990 mult(Ipath, evol_op, Ihold);
991 Ipath /= Ipath[0];
992 Ipath *= Ihold[0];
993 mult(Ihold, P, Ipath);
994 mult(I_i, trans_mat, Ihold);
995 Ihold = Ipath;
996 if (Ihold[0] < 1e-40 || std::isnan(Ihold[0]) ||
997 std::isnan(Ihold[1]) ||
998 (stokes_dim > 2 && std::isnan(Ihold[2])) ||
999 (stokes_dim > 3 && std::isnan(Ihold[3]))) {
1000 integrity = false;
1001 }
1002
1003 if (r_trav > r_min && integrity) {
1004 // Add reflectivity to proper range bin
1005 Index ibin = 0;
1006 r_bin = 0.0;
1007 while (r_bin < r_trav && ibin <= nbins + 1) {
1008 ibin++;
1009 r_bin = range_bins[ibin];
1010 }
1011 ibin -= 1;
1012
1013 // Calculate rx antenna weight and polarization rotation
1014 Matrix R_rx(3, 3);
1015 rotmat_enu(R_rx, rte_los_antenna);
1016 mc_antenna.return_los(antenna_wgt, R_rx, R_enu2ant);
1018 R_stokes, stokes_dim, rx_dir, tx_dir, R_rx, R_ant2enu);
1019 mult(I_i_rot, R_stokes, I_i);
1020
1021 for (Index istokes = 0; istokes < stokes_dim; istokes++) {
1022 Index ibiny = ibin * stokes_dim + istokes;
1023 ARTS_ASSERT(!std::isnan(I_i_rot[istokes]));
1024 Isum[ibiny] += antenna_wgt * I_i_rot[istokes];
1025 Isquaredsum[ibiny] += antenna_wgt * antenna_wgt *
1026 I_i_rot[istokes] * I_i_rot[istokes];
1027 }
1028 range_bin_count[ibin] += 1;
1029 }
1030
1031 scat_order++;
1032
1033 Sample_los_uniform(new_rte_los, rng);
1034 pdir_array(0, joker) = new_rte_los;
1035 // alt:
1036 // Sample_los_uniform( pdir_array(0,joker), rng );
1037 pha_mat_NScatElems(pha_mat_Nse,
1038 ptypes_Nse,
1039 t_ok,
1040 scat_data,
1041 stokes_dim,
1042 t_array,
1043 pdir_array,
1044 idir_array,
1045 f_index,
1046 t_interp_order);
1047 pha_mat_ScatSpecBulk(pha_mat_ssbulk,
1048 ptype_ssbulk,
1049 pha_mat_Nse,
1050 ptypes_Nse,
1051 pnds,
1052 t_ok);
1053 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
1054 Z = pha_mat_bulk(0, 0, 0, 0, joker, joker);
1055
1056 Z *= 4 * PI;
1057 Z /= Csca;
1058 mult(Ipath, Z, Ihold);
1059 Ihold = Ipath;
1060 local_rte_los = new_rte_los;
1061 // alt:
1062 //local_rte_los = pdir_array(0,joker);
1063 // or even (but also requires replacements of local_rte_los
1064 // with idir_array throughout the whole loop):
1065 //idir_array = pdir_array;
1066 } else {
1067 // Past farthest range
1068 keepgoing = false;
1069 }
1070 }
1071
1072 // Some checks
1073 if (scat_order >= mc_max_scatorder) keepgoing = false;
1074 if (!integrity) keepgoing = false;
1075 } // while (inner: keepgoing)
1076
1077 } // while (outer)
1078
1079 // Normalize range bins and apply sensor response (polarization)
1080 for (Index ibin = 0; ibin < nbins; ibin++) {
1081 for (Index istokes = 0; istokes < stokes_dim; istokes++) {
1082 Index ibiny = ibin * stokes_dim + istokes;
1083 if (range_bin_count[ibin] > 0) {
1084 y[ibiny] = Isum[ibiny] / ((Numeric)mc_iter) / bin_height[ibin];
1085 mc_error[ibiny] = sqrt((Isquaredsum[ibiny] / (Numeric)mc_iter /
1086 bin_height[ibin] / bin_height[ibin] -
1087 y[ibiny] * y[ibiny]) /
1088 (Numeric)mc_iter);
1089 }
1090 }
1091 }
1092
1093 y *= fac;
1094 mc_error *= fac;
1095} // end MCRadar
1096
1097/* Workspace method: Doxygen documentation will be auto-generated */
1098void MCSetSeedFromTime(Index& mc_seed, const Verbosity&) {
1099 mc_seed = (Index)time(NULL);
1100}
The global header file for ARTS.
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:23658
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:24392
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 &verbosity)
WORKSPACE METHOD: rte_losGeometricFromRtePosToRtePos2.
Definition: m_ppath.cc:1898
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:737
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:44
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:430
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
An Antenna object used by MCGeneral.
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
AntennaType get_type() const
AntennaType get_type.
Definition: mc_antenna.cc:141
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:119
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
set_gaussian_fwhm.
Definition: mc_antenna.cc:125
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
The Matrix class.
Definition: matpackI.h:1225
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:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
The Tensor4 class.
Definition: matpackIV.h:421
The Tensor6 class.
Definition: matpackVI.h:1088
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define q
#define min(a, b)
#define max(a, b)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2235
Linear algebra functions.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
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.
void mc_antennaSetPencilBeam(MCAntenna &mc_antenna, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetPencilBeam.
Definition: m_montecarlo.cc:82
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
const Numeric RAD2DEG
const Numeric PI
const Numeric SPEED_OF_LIGHT
const Numeric BOLTZMAN_CONST
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
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
const Numeric DEG2RAD
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
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
Implementation of Matrix, Vector, and such stuff.
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
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:86
void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
rotmat_enu.
Definition: mc_antenna.cc:63
@ ANTENNA_TYPE_GAUSSIAN
Definition: mc_antenna.h:43
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:204
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:1198
bool is_anyptype_nonTotRan(const ArrayOfArrayOfSingleScatteringData &scat_data)
is_anyptype_nonTotRan.
Definition: montecarlo.cc:840
void Sample_los_uniform(VectorView new_rte_los, Rng &rng)
Sample_los_uniform.
Definition: montecarlo.cc:1613
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:1533
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:857
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:690
constexpr Numeric Q(Numeric T, const Species::IsotopeRecord &ir)
Definition: partfun.h:139
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
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:739
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:2389
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:2010
Declaration of functions in rte.cc.
Header file for special_interp.cc.
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath.h:86
Index np
Number of points describing the ppath.
Definition: ppath.h:52
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath.h:84
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
Vector lstep
The length between ppath points.
Definition: ppath.h:70
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
Vector ngroup
The group index of refraction.
Definition: ppath.h:80
This file contains basic functions to handle XML data files.