ARTS 2.5.4 (git: 31ce4f0e)
montecarlo.cc
Go to the documentation of this file.
1/* copyright (C) 2003-2012 Cory Davis <cory.davis@metservice.com>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
26/*===========================================================================
27 === External declarations
28 ===========================================================================*/
29
30#include <cfloat>
31#include <sstream>
32
33#include "arts_constants.h"
34#include "arts_conversions.h"
35#include "auto_md.h"
36#include "geodetic.h"
37#include "mc_interp.h"
38#include "montecarlo.h"
39
41
42// Some root-finding helper functions (for MCRadar) that don't need
43// visibility outside this source file
44//
45//
46
61 const Numeric& Q,
62 const Numeric& kI,
63 const Numeric& kQ,
64 const Numeric& s) {
65 Numeric fs;
66
67 fs = exp(-kI * s) * (I * cosh(kQ * s) + Q * sinh(kQ * s));
68
69 return fs;
70}
71
135 const Numeric& a,
136 const Numeric& b,
137 const Numeric& t,
138 const Numeric& rn,
139 const Numeric& I,
140 const Numeric& Q,
141 const Numeric& KI,
142 const Numeric& KQ) {
143 Numeric c;
144 Numeric d;
145 Numeric e;
146 Numeric fa;
147 Numeric fb;
148 Numeric fc;
149 Numeric m;
150 Numeric macheps;
151 Numeric p;
152 Numeric q;
153 Numeric r;
154 Numeric s;
155 Numeric sa;
156 Numeric tol;
157 //
158 // Make local copies of A and B.
159 //
160 sa = a;
161 sb = b;
162 fa = ext_I(I, Q, KI, KQ, sa); // - rn;
163 fa -= rn;
164 fb = ext_I(I, Q, KI, KQ, sb); // - rn;
165 fb -= rn;
166
167 c = sa;
168 fc = fa;
169 e = sb - sa;
170 d = e;
171
172 macheps = DBL_EPSILON;
173
174 for (;;) {
175 if (abs(fc) < abs(fb)) {
176 sa = sb;
177 sb = c;
178 c = sa;
179 fa = fb;
180 fb = fc;
181 fc = fa;
182 }
183
184 tol = 2.0 * macheps * abs(sb) + t;
185 m = 0.5 * (c - sb);
186
187 if (abs(m) <= tol || fb == 0.0) {
188 break;
189 }
190
191 if (abs(e) < tol || abs(fa) <= abs(fb)) {
192 e = m;
193 d = e;
194 } else {
195 s = fb / fa;
196
197 if (sa == c) {
198 p = 2.0 * m * s;
199 q = 1.0 - s;
200 } else {
201 q = fa / fc;
202 r = fb / fc;
203 p = s * (2.0 * m * q * (q - r) - (sb - sa) * (r - 1.0));
204 q = (q - 1.0) * (r - 1.0) * (s - 1.0);
205 }
206
207 if (0.0 < p) {
208 q = -q;
209 } else {
210 p = -p;
211 }
212
213 s = e;
214 e = d;
215
216 if (2.0 * p < 3.0 * m * q - abs(tol * q) && p < abs(0.5 * s * q)) {
217 d = p / q;
218 } else {
219 e = m;
220 d = e;
221 }
222 }
223 sa = sb;
224 fa = fb;
225
226 if (tol < abs(d)) {
227 sb = sb + d;
228 } else if (0.0 < m) {
229 sb = sb + tol;
230 } else {
231 sb = sb - tol;
232 }
233
234 fb = ext_I(I, Q, KI, KQ, sb);
235 fb -= rn;
236
237 if ((0.0 < fb && 0.0 < fc) || (fb <= 0.0 && fc <= 0.0)) {
238 c = sa;
239 fc = fa;
240 e = sb - sa;
241 d = e;
242 }
243 }
244}
245
247 MatrixView ext_mat_mono,
248 VectorView abs_vec_mono,
249 Numeric& temperature,
250 const Agenda& propmat_clearsky_agenda,
251 const Numeric& f_mono,
252 const GridPos& gp_p,
253 const GridPos& gp_lat,
254 const GridPos& gp_lon,
255 ConstVectorView p_grid,
256 ConstTensor3View t_field,
257 ConstTensor4View vmr_field) {
258 const Index ns = vmr_field.nbooks();
259 Vector t_vec(1); //vectors are required by interp_atmfield_gp2itw etc.
260 Vector p_vec(1); //may not be efficient with unecessary vectors
261 Matrix itw_p(1, 2);
262 ArrayOfGridPos ao_gp_p(1), ao_gp_lat(1), ao_gp_lon(1);
263 Matrix vmr_mat(ns, 1), itw_field;
264
265 //local versions of workspace variables
266 StokesVector local_abs_vec;
267 StokesVector local_nlte_source_dummy;
268 PropagationMatrix local_ext_mat;
269 PropagationMatrix local_propmat_clearsky;
271 local_partial_dummy; // This is right since there should be only clearsky partials
272 ArrayOfStokesVector local_dnlte_source_dx_dummy;
273 ao_gp_p[0] = gp_p;
274 ao_gp_lat[0] = gp_lat;
275 ao_gp_lon[0] = gp_lon;
276
277 // Determine the pressure
278 interpweights(itw_p, ao_gp_p);
279 itw2p(p_vec, p_grid, ao_gp_p, itw_p);
280
281 // Determine the atmospheric temperature and species VMR
282 //
283 interp_atmfield_gp2itw(itw_field, 3, ao_gp_p, ao_gp_lat, ao_gp_lon);
284 //
286 t_vec, 3, t_field, ao_gp_p, ao_gp_lat, ao_gp_lon, itw_field);
287 //
288 for (Index is = 0; is < ns; is++) {
289 interp_atmfield_by_itw(vmr_mat(is, joker),
290 3,
291 vmr_field(is, joker, joker, joker),
292 ao_gp_p,
293 ao_gp_lat,
294 ao_gp_lon,
295 itw_field);
296 }
297
298 temperature = t_vec[0];
299
300 const Vector rtp_mag_dummy(3, 0);
301 const Vector ppath_los_dummy;
302 const EnergyLevelMap nlte_dummy;
303
304 //calcualte absorption coefficient
306 local_propmat_clearsky,
307 local_nlte_source_dummy,
308 local_partial_dummy,
309 local_dnlte_source_dx_dummy,
311 {},
312 Vector(1, f_mono),
313 rtp_mag_dummy,
314 ppath_los_dummy,
315 p_vec[0],
316 temperature,
317 nlte_dummy,
318 vmr_mat(joker, 0),
319 propmat_clearsky_agenda);
320
322 local_ext_mat, local_abs_vec, local_propmat_clearsky);
323
324 local_ext_mat.MatrixAtPosition(ext_mat_mono);
325 abs_vec_mono = local_abs_vec.VectorAtPosition();
326}
327
329 MatrixView ext_mat_mono,
330 VectorView abs_vec_mono,
331 VectorView pnd_vec,
332 Numeric& temperature,
333 const Agenda& propmat_clearsky_agenda,
334 const Index stokes_dim,
335 const Index f_index,
336 const Vector& f_grid,
337 const GridPos& gp_p,
338 const GridPos& gp_lat,
339 const GridPos& gp_lon,
340 ConstVectorView p_grid_cloud,
341 ConstTensor3View t_field_cloud,
342 ConstTensor4View vmr_field_cloud,
343 const Tensor4& pnd_field,
344 const ArrayOfArrayOfSingleScatteringData& scat_data,
345 const ArrayOfIndex& cloudbox_limits,
346 const Vector& rte_los)
347
348{
349 const Index ns = vmr_field_cloud.nbooks();
350 const Index N_se = pnd_field.nbooks();
351 Matrix pnd_ppath(N_se, 1);
352 Vector t_ppath(1);
353 Vector p_ppath(1); //may not be efficient with unecessary vectors
354 EnergyLevelMap nlte_dummy;
355 Matrix itw_p(1, 2);
356 ArrayOfGridPos ao_gp_p(1), ao_gp_lat(1), ao_gp_lon(1);
357 Matrix vmr_ppath(ns, 1), itw_field;
358
359 //local versions of workspace variables
361 local_partial_dummy; // This is right since there should be only clearsky partials
362 ArrayOfStokesVector local_dnlte_source_dx_dummy; // This is right since there should be only clearsky partials
363 PropagationMatrix local_propmat_clearsky;
364 StokesVector local_nlte_source_dummy; //FIXME: Do this right?
365 StokesVector local_abs_vec;
366 PropagationMatrix local_ext_mat;
367
368 ao_gp_p[0] = gp_p;
369 ao_gp_lat[0] = gp_lat;
370 ao_gp_lon[0] = gp_lon;
371
372 cloud_atm_vars_by_gp(p_ppath,
373 t_ppath,
374 vmr_ppath,
375 pnd_ppath,
376 ao_gp_p,
377 ao_gp_lat,
378 ao_gp_lon,
379 cloudbox_limits,
380 p_grid_cloud,
381 t_field_cloud,
382 vmr_field_cloud,
383 pnd_field);
384 pnd_vec = pnd_ppath(joker, 0);
385 temperature = t_ppath[0];
386
387 const Vector rtp_mag_dummy(3, 0);
388 const Vector ppath_los_dummy;
389
390 //rtp_vmr = vmr_ppath(joker,0);
392 local_propmat_clearsky,
393 local_nlte_source_dummy,
394 local_partial_dummy,
395 local_dnlte_source_dx_dummy,
397 {},
398 f_grid[Range(f_index, 1)],
399 rtp_mag_dummy,
400 ppath_los_dummy,
401 p_ppath[0],
402 temperature,
403 nlte_dummy,
404 vmr_ppath(joker, 0),
405 propmat_clearsky_agenda);
406
408 local_ext_mat, local_abs_vec, local_propmat_clearsky);
409
410 local_ext_mat.MatrixAtPosition(ext_mat_mono);
411 abs_vec_mono = local_abs_vec.VectorAtPosition();
412
413 ArrayOfArrayOfTensor5 ext_mat_Nse;
414 ArrayOfArrayOfTensor4 abs_vec_Nse;
415 ArrayOfArrayOfIndex ptypes_Nse;
416 Matrix t_ok;
417 ArrayOfTensor5 ext_mat_ssbulk;
418 ArrayOfTensor4 abs_vec_ssbulk;
419 ArrayOfIndex ptype_ssbulk;
420 Tensor5 ext_mat_bulk;
421 Tensor4 abs_vec_bulk;
422 Index ptype_bulk;
423
424 Vector sca_dir;
425 mirror_los(sca_dir, rte_los, 3);
426 Matrix dir_array(1, 2, 0.);
427 dir_array(0, joker) = sca_dir;
428 //
429 opt_prop_NScatElems(ext_mat_Nse,
430 abs_vec_Nse,
431 ptypes_Nse,
432 t_ok,
433 scat_data,
434 stokes_dim,
435 t_ppath,
436 dir_array,
437 f_index);
438 //
439 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
440 abs_vec_ssbulk,
441 ptype_ssbulk,
442 ext_mat_Nse,
443 abs_vec_Nse,
444 ptypes_Nse,
445 pnd_ppath,
446 t_ok);
447 opt_prop_Bulk(ext_mat_bulk,
448 abs_vec_bulk,
449 ptype_bulk,
450 ext_mat_ssbulk,
451 abs_vec_ssbulk,
452 ptype_ssbulk);
453
454 ext_mat_mono += ext_mat_bulk(0, 0, 0, joker, joker);
455 abs_vec_mono += abs_vec_bulk(0, 0, 0, joker);
456}
457
459 VectorView temperature,
460 MatrixView vmr,
461 MatrixView pnd,
462 const ArrayOfGridPos& gp_p,
463 const ArrayOfGridPos& gp_lat,
464 const ArrayOfGridPos& gp_lon,
465 const ArrayOfIndex& cloudbox_limits,
466 ConstVectorView p_grid_cloud,
467 ConstTensor3View t_field_cloud,
468 ConstTensor4View vmr_field_cloud,
469 ConstTensor4View pnd_field) {
470 Index np = gp_p.nelem();
471 ARTS_ASSERT(pressure.nelem() == np);
472 Index ns = vmr_field_cloud.nbooks();
473 Index N_se = pnd_field.nbooks();
474 ArrayOfGridPos gp_p_cloud = gp_p;
475 ArrayOfGridPos gp_lat_cloud = gp_lat;
476 ArrayOfGridPos gp_lon_cloud = gp_lon;
477 Index atmosphere_dim = 3;
478
479 for (Index i = 0; i < np; i++) {
480 // Calculate grid positions inside the cloud.
481 gp_p_cloud[i].idx -= cloudbox_limits[0];
482 gp_lat_cloud[i].idx -= cloudbox_limits[2];
483 gp_lon_cloud[i].idx -= cloudbox_limits[4];
484 }
485
486 const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
487 const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
488 const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
489 gridpos_upperend_check(gp_p_cloud[0], n1);
490 gridpos_upperend_check(gp_p_cloud[np - 1], n1);
491 gridpos_upperend_check(gp_lat_cloud[0], n2);
492 gridpos_upperend_check(gp_lat_cloud[np - 1], n2);
493 gridpos_upperend_check(gp_lon_cloud[0], n3);
494 gridpos_upperend_check(gp_lon_cloud[np - 1], n3);
495
496 // Determine the pressure at each propagation path point
497 Matrix itw_p(np, 2);
498 //
499 //interpweights( itw_p, ppath.gp_p );
500 interpweights(itw_p, gp_p_cloud);
501 itw2p(pressure, p_grid_cloud, gp_p_cloud, itw_p);
502
503 // Determine the atmospheric temperature and species VMR at
504 // each propagation path point
505 Matrix itw_field;
506 //
508 itw_field, atmosphere_dim, gp_p_cloud, gp_lat_cloud, gp_lon_cloud);
509 //
510 interp_atmfield_by_itw(temperature,
511 atmosphere_dim,
512 t_field_cloud,
513 gp_p_cloud,
514 gp_lat_cloud,
515 gp_lon_cloud,
516 itw_field);
517 //
518 for (Index is = 0; is < ns; is++) {
520 atmosphere_dim,
521 vmr_field_cloud(is, joker, joker, joker),
522 gp_p_cloud,
523 gp_lat_cloud,
524 gp_lon_cloud,
525 itw_field);
526 }
527
528 //Determine the particle number density for every scattering element at
529 // each propagation path point
530 for (Index i_se = 0; i_se < N_se; i_se++) {
531 // if grid positions still outside the range the propagation path step
532 // must be outside the cloudbox and pnd is set to zero
533 interp_atmfield_by_itw(pnd(i_se, joker),
534 atmosphere_dim,
535 pnd_field(i_se, joker, joker, joker),
536 gp_p_cloud,
537 gp_lat_cloud,
538 gp_lon_cloud,
539 itw_field);
540 }
541}
542
543void ext_mat_case(Index& icase,
544 ConstMatrixView ext_mat,
545 const Index stokes_dim) {
546 if (icase == 0) {
547 icase = 1; // Start guess is diagonal
548
549 //--- Scalar case ----------------------------------------------------------
550 if (stokes_dim == 1) {
551 }
552
553 //--- Vector RT ------------------------------------------------------------
554 else {
555 // Check symmetries and analyse structure of exp_mat:
556 ARTS_ASSERT(ext_mat(1, 1) == ext_mat(0, 0));
557 ARTS_ASSERT(ext_mat(1, 0) == ext_mat(0, 1));
558
559 if (ext_mat(1, 0) != 0) {
560 icase = 2;
561 }
562
563 if (stokes_dim >= 3) {
564 ARTS_ASSERT(ext_mat(2, 2) == ext_mat(0, 0));
565 ARTS_ASSERT(ext_mat(2, 1) == -ext_mat(1, 2));
566 ARTS_ASSERT(ext_mat(2, 0) == ext_mat(0, 2));
567
568 if (ext_mat(2, 0) != 0 || ext_mat(2, 1) != 0) {
569 icase = 3;
570 }
571
572 if (stokes_dim > 3) {
573 ARTS_ASSERT(ext_mat(3, 3) == ext_mat(0, 0));
574 ARTS_ASSERT(ext_mat(3, 2) == -ext_mat(2, 3));
575 ARTS_ASSERT(ext_mat(3, 1) == -ext_mat(1, 3));
576 ARTS_ASSERT(ext_mat(3, 0) == ext_mat(0, 3));
577
578 if (icase < 3) // if icase==3, already at most complex case
579 {
580 if (ext_mat(3, 0) != 0 || ext_mat(3, 1) != 0) {
581 icase = 3;
582 } else if (ext_mat(3, 2) != 0) {
583 icase = 2;
584 }
585 }
586 }
587 }
588 }
589 }
590}
591
617void ext2trans(MatrixView trans_mat,
618 Index& icase,
619 ConstMatrixView ext_mat,
620 const Numeric& lstep) {
621 const Index stokes_dim = ext_mat.ncols();
622
623 ARTS_ASSERT(ext_mat.nrows() == stokes_dim);
624 ARTS_ASSERT(trans_mat.nrows() == stokes_dim && trans_mat.ncols() == stokes_dim);
625
626 // Theoretically ext_mat(0,0) >= 0, but to demand this can cause problems for
627 // iterative retrievals, and the assert is skipped. Negative should be a
628 // result of negative vmr, and this issue is checked in basics_checkedCalc.
629 //ARTS_ASSERT( ext_mat(0,0) >= 0 );
630
631 ARTS_ASSERT(icase >= 0 && icase <= 3);
632 ARTS_ASSERT(!is_singular(ext_mat));
633 ARTS_ASSERT(lstep >= 0);
634
635 // Analyse ext_mat?
636 ext_mat_case(icase, ext_mat, stokes_dim);
637
638 // Calculation options:
639 if (icase == 1) {
640 trans_mat = 0;
641 trans_mat(0, 0) = exp(-ext_mat(0, 0) * lstep);
642 for (Index i = 1; i < stokes_dim; i++) {
643 trans_mat(i, i) = trans_mat(0, 0);
644 }
645 }
646
647 else if (icase == 2 && stokes_dim < 3) {
648 // Expressions below are found in "Polarization in Spectral Lines" by
649 // Landi Degl'Innocenti and Landolfi (2004).
650 const Numeric tI = exp(-ext_mat(0, 0) * lstep);
651 const Numeric HQ = ext_mat(0, 1) * lstep;
652 trans_mat(0, 0) = tI * cosh(HQ);
653 trans_mat(1, 1) = trans_mat(0, 0);
654 trans_mat(1, 0) = -tI * sinh(HQ);
655 trans_mat(0, 1) = trans_mat(1, 0);
656 /* Does not work for stokes_dim==3, and commnted out 180502 by PE:
657 if( stokes_dim >= 3 )
658 {
659 trans_mat(2,0) = 0;
660 trans_mat(2,1) = 0;
661 trans_mat(0,2) = 0;
662 trans_mat(1,2) = 0;
663 const Numeric RQ = ext_mat(2,3) * lstep;
664 trans_mat(2,2) = tI * cos( RQ );
665 if( stokes_dim > 3 )
666 {
667 trans_mat(3,0) = 0;
668 trans_mat(3,1) = 0;
669 trans_mat(0,3) = 0;
670 trans_mat(1,3) = 0;
671 trans_mat(3,3) = trans_mat(2,2);
672 trans_mat(3,2) = tI * sin( RQ );
673 trans_mat(2,3) = -trans_mat(3,2);
674 }
675 }
676 */
677 } else {
678 Matrix ext_mat_ds = ext_mat;
679 ext_mat_ds *= -lstep;
680 //
681 constexpr Index q = 10; // index for the precision of the matrix exp function
682 matrix_exp(trans_mat, ext_mat_ds, q);
683 }
684}
685
687 MatrixView& trans_mat,
688 const Ppath& ppath,
689 const Agenda& propmat_clearsky_agenda,
690 const Index stokes_dim,
691 const Index f_index,
692 const Vector& f_grid,
693 const Vector& p_grid,
694 const Tensor3& t_field,
695 const Tensor4& vmr_field,
696 const ArrayOfIndex& cloudbox_limits,
697 const Tensor4& pnd_field,
698 const ArrayOfArrayOfSingleScatteringData& scat_data,
699 const Verbosity& verbosity) {
700 bool inside_cloud;
701 const Index np = ppath.np;
702 ArrayOfMatrix ext_matArray(2);
703 ArrayOfMatrix trans_matArray(2);
704 Index N_se = pnd_field.nbooks(); //Number of scattering elements
705 Vector pnd_vec(N_se);
706 Vector abs_vec_mono(stokes_dim);
707 Matrix ext_mat(stokes_dim, stokes_dim);
708 Matrix ext_mat_mono(stokes_dim, stokes_dim);
709 Matrix incT(stokes_dim, stokes_dim, 0.0);
710 Numeric temperature;
711 Numeric dl = -999;
712
714
715 id_mat(trans_mat);
716
717 if (np > 1) {
718 // range defining cloudbox
719 Range p_range(cloudbox_limits[0],
720 cloudbox_limits[1] - cloudbox_limits[0] + 1);
721 Range lat_range(cloudbox_limits[2],
722 cloudbox_limits[3] - cloudbox_limits[2] + 1);
723 Range lon_range(cloudbox_limits[4],
724 cloudbox_limits[5] - cloudbox_limits[4] + 1);
725
726 inside_cloud = is_gp_inside_cloudbox(ppath.gp_p[np - 1],
727 ppath.gp_lat[np - 1],
728 ppath.gp_lon[np - 1],
729 cloudbox_limits,
730 0,
731 3);
732 if (inside_cloud) {
734 ext_mat_mono,
735 abs_vec_mono,
736 pnd_vec,
737 temperature,
738 propmat_clearsky_agenda,
739 stokes_dim,
740 f_index,
741 f_grid,
742 ppath.gp_p[np - 1],
743 ppath.gp_lat[np - 1],
744 ppath.gp_lon[np - 1],
745 p_grid[p_range],
746 t_field(p_range, lat_range, lon_range),
747 vmr_field(joker, p_range, lat_range, lon_range),
748 pnd_field,
749 scat_data,
750 cloudbox_limits,
751 ppath.los(np - 1, joker));
752 } else {
754 ext_mat_mono,
755 abs_vec_mono,
756 temperature,
757 propmat_clearsky_agenda,
758 f_grid[f_index],
759 ppath.gp_p[np - 1],
760 ppath.gp_lat[np - 1],
761 ppath.gp_lon[np - 1],
762 p_grid,
763 t_field,
764 vmr_field);
765 pnd_vec = 0.0;
766 }
767
768 trans_matArray[1] = trans_mat;
769 ext_matArray[1] = ext_mat_mono;
770
771 // Index in ppath of end point considered presently
772 for (Index ip = np - 2; ip >= 0; ip--) {
773 dl = ppath.lstep[ip];
774
775 ext_matArray[0] = ext_matArray[1];
776 trans_matArray[0] = trans_matArray[1];
777
778 inside_cloud = is_gp_inside_cloudbox(ppath.gp_p[ip],
779 ppath.gp_lat[ip],
780 ppath.gp_lon[ip],
781 cloudbox_limits,
782 0,
783 3);
784 if (inside_cloud) {
786 ext_mat_mono,
787 abs_vec_mono,
788 pnd_vec,
789 temperature,
790 propmat_clearsky_agenda,
791 stokes_dim,
792 f_index,
793 f_grid,
794 ppath.gp_p[ip],
795 ppath.gp_lat[ip],
796 ppath.gp_lon[ip],
797 p_grid[p_range],
798 t_field(p_range, lat_range, lon_range),
799 vmr_field(joker, p_range, lat_range, lon_range),
800 pnd_field,
801 scat_data,
802 cloudbox_limits,
803 ppath.los(ip, joker));
804 } else {
806 ext_mat_mono,
807 abs_vec_mono,
808 temperature,
809 propmat_clearsky_agenda,
810 f_grid[f_index],
811 ppath.gp_p[ip],
812 ppath.gp_lat[ip],
813 ppath.gp_lon[ip],
814 p_grid,
815 t_field,
816 vmr_field);
817 pnd_vec = 0.0;
818 }
819
820 ext_matArray[1] = ext_mat_mono;
821 ext_mat = ext_matArray[0];
822 ext_mat += ext_matArray[1]; // Factor 2 fixed by using dl/2
823 //
824 {
825 Index extmat_case = 0;
826 ext2trans(incT, extmat_case, ext_mat, dl / 2);
827 }
828 //
829 mult(trans_mat, incT, trans_matArray[1]);
830 trans_matArray[1] = trans_mat;
831
832 } // for( ip... )
833 } // if( np > 1 )
834}
835
837 const ArrayOfArrayOfSingleScatteringData& scat_data) {
838 bool is_anyptype_nonTotRan = false;
839 for (Index i_ss = 0;
840 is_anyptype_nonTotRan == false && i_ss < scat_data.nelem();
841 i_ss++) {
842 for (Index i_se = 0;
843 is_anyptype_nonTotRan == false && i_se < scat_data[i_ss].nelem();
844 i_se++) {
845 if (scat_data[i_ss][i_se].ptype > PTYPE_TOTAL_RND) {
847 }
848 }
849 }
851}
852
854 MatrixView evol_op,
855 Vector& abs_vec_mono,
856 Numeric& temperature,
857 MatrixView ext_mat_mono,
858 Rng& rng,
859 Vector& rte_pos,
860 Vector& rte_los,
861 Vector& pnd_vec,
862 Numeric& g,
863 Ppath& ppath_step,
864 Index& termination_flag,
865 bool& inside_cloud,
866 const Agenda& ppath_step_agenda,
867 const Numeric& ppath_lmax,
868 const Numeric& ppath_lraytrace,
869 const Numeric& taustep_limit,
870 const Agenda& propmat_clearsky_agenda,
871 const Index stokes_dim,
872 const Index f_index,
873 const Vector& f_grid,
874 const Vector& p_grid,
875 const Vector& lat_grid,
876 const Vector& lon_grid,
877 const Tensor3& z_field,
878 const Vector& refellipsoid,
879 const Matrix& z_surface,
880 const Tensor3& t_field,
881 const Tensor4& vmr_field,
882 const ArrayOfIndex& cloudbox_limits,
883 const Tensor4& pnd_field,
884 const ArrayOfArrayOfSingleScatteringData& scat_data,
885 const Verbosity& verbosity) {
886 ArrayOfMatrix evol_opArray(2);
887 ArrayOfMatrix ext_matArray(2);
888 ArrayOfVector abs_vecArray(2);
889 ArrayOfVector pnd_vecArray(2);
890 Matrix ext_mat(stokes_dim, stokes_dim);
891 Matrix incT(stokes_dim, stokes_dim, 0.0);
892 Vector tArray(2);
893 Matrix T(stokes_dim, stokes_dim);
894 Numeric k;
895 Numeric ds, dl = -999;
896 Index istep = 0; // Counter for number of steps
897 Matrix old_evol_op(stokes_dim, stokes_dim);
898
900
901 //at the start of the path the evolution operator is the identity matrix
902 id_mat(evol_op);
903 evol_opArray[1] = evol_op;
904
905 // range defining cloudbox
906 Range p_range(cloudbox_limits[0],
907 cloudbox_limits[1] - cloudbox_limits[0] + 1);
908 Range lat_range(cloudbox_limits[2],
909 cloudbox_limits[3] - cloudbox_limits[2] + 1);
910 Range lon_range(cloudbox_limits[4],
911 cloudbox_limits[5] - cloudbox_limits[4] + 1);
912
913 //initialise Ppath with ppath_start_stepping
914 ppath_start_stepping(ppath_step,
915 3,
916 p_grid,
917 lat_grid,
918 lon_grid,
919 z_field,
920 refellipsoid,
921 z_surface,
922 0,
923 cloudbox_limits,
924 false,
925 rte_pos,
926 rte_los,
927 verbosity);
928
929 // Check if we have already has radiative background
930 if (ppath_what_background(ppath_step)) {
931 termination_flag = ppath_what_background(ppath_step);
932 g = 1;
933 return;
934 }
935
936 // Index in ppath_step of end point considered presently
937 Index ip = 0;
938
939 // Is point ip inside the cloudbox?
940 inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
941 ppath_step.gp_lat[ip],
942 ppath_step.gp_lon[ip],
943 cloudbox_limits,
944 0,
945 3);
946
947 // Determine radiative properties at point
948 if (inside_cloud) {
950 ext_mat_mono,
951 abs_vec_mono,
952 pnd_vec,
953 temperature,
954 propmat_clearsky_agenda,
955 stokes_dim,
956 f_index,
957 f_grid,
958 ppath_step.gp_p[0],
959 ppath_step.gp_lat[0],
960 ppath_step.gp_lon[0],
961 p_grid[p_range],
962 t_field(p_range, lat_range, lon_range),
963 vmr_field(joker, p_range, lat_range, lon_range),
964 pnd_field,
965 scat_data,
966 cloudbox_limits,
967 ppath_step.los(0, joker));
968 } else {
970 ext_mat_mono,
971 abs_vec_mono,
972 temperature,
973 propmat_clearsky_agenda,
974 f_grid[f_index],
975 ppath_step.gp_p[0],
976 ppath_step.gp_lat[0],
977 ppath_step.gp_lon[0],
978 p_grid,
979 t_field,
980 vmr_field);
981 pnd_vec = 0.0;
982 }
983
984 // Move the data to end point containers
985 ext_matArray[1] = ext_mat_mono;
986 abs_vecArray[1] = abs_vec_mono;
987 tArray[1] = temperature;
988 pnd_vecArray[1] = pnd_vec;
989
990 //draw random number to determine end point
991 Numeric r = rng.draw();
992
993 termination_flag = 0;
994
995 while ((evol_op(0, 0) > r) && (!termination_flag)) {
996 istep++;
997
998 ARTS_USER_ERROR_IF (istep > 100000,
999 "100000 path points have been reached. "
1000 "Is this an infinite loop?");
1001
1002 evol_opArray[0] = evol_opArray[1];
1003 ext_matArray[0] = ext_matArray[1];
1004 abs_vecArray[0] = abs_vecArray[1];
1005 tArray[0] = tArray[1];
1006 pnd_vecArray[0] = pnd_vecArray[1];
1007
1008 // The algorith to meet taustep_lim:
1009 // When first calculating a new ppath_step, it assumed that the present
1010 // ppath position holds the highest extinction. If the extinction at the
1011 // next position is higher, the criterion is checked and a new ppath_step
1012 // calculation is triggered if found necessary.
1013 // This should work in most cases, but is not 100% safe. Consider a case
1014 // with ppath_lmax = -1 and the extinction is zero at all grid box
1015 // corners except one. The two "test points" can then both get an
1016 // extinction of zero, while in fact is non-zero through the grid box and
1017 // the optical depth is underestimated. But this was handled equally bad
1018 // before taustep_limit was introduced (2016-10-10, PE)
1019 bool oktaustep = false;
1020 Index ppath_try = 1;
1021 const Index lmax_limit = 10;
1022
1023 while (!oktaustep) {
1024 // Shall new ppath_step be calculated?
1025 if (ip == ppath_step.np - 1) {
1026 Numeric lmax = taustep_limit / ext_mat_mono(0, 0);
1027 if (ppath_lmax > 0) {
1028 lmax = min(ppath_lmax, lmax);
1029 }
1030 if (lmax < lmax_limit) {
1031 lmax = lmax_limit;
1032 }
1033 //cout << ppath_try << ", lmax = " << lmax << endl;
1034 //Print( ppath_step, 0, verbosity );
1035
1037 ppath_step,
1038 lmax,
1039 ppath_lraytrace,
1040 f_grid[Range(f_index, 1)],
1041 ppath_step_agenda);
1042 ip = 1;
1043
1044 inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
1045 ppath_step.gp_lat[ip],
1046 ppath_step.gp_lon[ip],
1047 cloudbox_limits,
1048 0,
1049 3);
1050 } else {
1051 ip++;
1052 }
1053
1054 if (inside_cloud) {
1056 ext_mat_mono,
1057 abs_vec_mono,
1058 pnd_vec,
1059 temperature,
1060 propmat_clearsky_agenda,
1061 stokes_dim,
1062 f_index,
1063 f_grid,
1064 ppath_step.gp_p[ip],
1065 ppath_step.gp_lat[ip],
1066 ppath_step.gp_lon[ip],
1067 p_grid[p_range],
1068 t_field(p_range, lat_range, lon_range),
1069 vmr_field(joker, p_range, lat_range, lon_range),
1070 pnd_field,
1071 scat_data,
1072 cloudbox_limits,
1073 ppath_step.los(ip, joker));
1074 } else {
1076 ext_mat_mono,
1077 abs_vec_mono,
1078 temperature,
1079 propmat_clearsky_agenda,
1080 f_grid[f_index],
1081 ppath_step.gp_p[ip],
1082 ppath_step.gp_lat[ip],
1083 ppath_step.gp_lon[ip],
1084 p_grid,
1085 t_field,
1086 vmr_field);
1087 pnd_vec = 0.0;
1088 }
1089
1090 dl = ppath_step.lstep[ip - 1];
1091
1092 // Check if taustep_limit criterion is met. OK if:
1093 // 1. Ppath step already recalculated
1094 // 2. New ext_mat <= old one
1095 // 3. New ext_mat bigger, but tau of step still below limit
1096 if (ppath_try > 1 || ext_mat_mono(0, 0) <= ext_matArray[0](0, 0) ||
1097 (ext_mat_mono(0, 0) + ext_matArray[0](0, 0)) * dl / 2 <=
1098 taustep_limit) {
1099 oktaustep = true;
1100 } else {
1101 // We trigger a recalculation of ppath_step, from the previous
1102 // point
1103 ppath_step.np = ip;
1104 ip--;
1105 ppath_try = 2;
1106 // If a background found in first try this has to be reset:
1107 ppath_set_background(ppath_step, 0);
1108 }
1109 } // while !oktuastep
1110
1111 ext_matArray[1] = ext_mat_mono;
1112 abs_vecArray[1] = abs_vec_mono;
1113 tArray[1] = temperature;
1114 pnd_vecArray[1] = pnd_vec;
1115 ext_mat = ext_matArray[1];
1116 ext_mat += ext_matArray[0]; // Factor 2 fixed by using dl/2
1117 //
1118 {
1119 Index extmat_case = 0;
1120 ext2trans(incT, extmat_case, ext_mat, dl / 2);
1121 }
1122 //
1123 mult(evol_op, evol_opArray[0], incT);
1124 evol_opArray[1] = evol_op;
1125
1126 if (evol_op(0, 0) > r) {
1127 // Check whether hit ground or space.
1128 // path_step_agenda just detects surface intersections, and
1129 // if TOA is reached requires a special check.
1130 // But we are already ready if evol_op<=r
1131 if (ip == ppath_step.np - 1) {
1132 if (ppath_what_background(ppath_step)) {
1133 termination_flag = 2;
1134 } //we have hit the surface
1135 else if (fractional_gp(ppath_step.gp_p[ip]) >=
1136 (Numeric)(p_grid.nelem() - 1) - 1e-3) {
1137 termination_flag = 1;
1138 } //we are at TOA
1139 }
1140 }
1141 } // while
1142
1143 if (termination_flag) { //we must have reached the cloudbox boundary
1144 rte_pos = ppath_step.pos(ip, joker);
1145 rte_los = ppath_step.los(ip, joker);
1146 g = evol_op(0, 0);
1147 } else {
1148 //find position...and evol_op..and everything else required at the new
1149 //scattering/emission point
1150 // GH 2011-09-14:
1151 // log(incT(0,0)) = log(exp(opt_depth_mat(0, 0))) = opt_depth_mat(0, 0)
1152 // Avoid loss of precision, use opt_depth_mat directly
1153 //k=-log(incT(0,0))/cum_l_step[np-1];//K=K11 only for diagonal ext_mat
1154 // PE 2013-05-17, Now the above comes directly from ext_mat:
1155 k = ext_mat(0, 0) / 2; // Factor 2 as sum of end point values
1156 ds = log(evol_opArray[0](0, 0) / r) / k;
1157 g = k * r;
1158 Vector x(2, 0.0);
1159 //interpolate atmospheric variables required later
1160 ArrayOfGridPos gp(1);
1161 x[1] = dl;
1162 Vector itw(2);
1163
1164 gridpos(gp, x, ds);
1165 ARTS_ASSERT(gp[0].idx == 0);
1166 interpweights(itw, gp[0]);
1167 interp(ext_mat_mono, itw, ext_matArray, gp[0]);
1168 ext_mat = ext_mat_mono;
1169 ext_mat += ext_matArray[gp[0].idx];
1170 //
1171 {
1172 Index extmat_case = 0;
1173 ext2trans(incT, extmat_case, ext_mat, ds / 2);
1174 }
1175 //
1176 mult(evol_op, evol_opArray[gp[0].idx], incT);
1177 interp(abs_vec_mono, itw, abs_vecArray, gp[0]);
1178 temperature = interp(itw, tArray, gp[0]);
1179 interp(pnd_vec, itw, pnd_vecArray, gp[0]);
1180 for (Index i = 0; i < 2; i++) {
1181 rte_pos[i] = interp(itw, ppath_step.pos(Range(ip - 1, 2), i), gp[0]);
1182 rte_los[i] = interp(itw, ppath_step.los(Range(ip - 1, 2), i), gp[0]);
1183 }
1184 rte_pos[2] = interp(itw, ppath_step.pos(Range(ip - 1, 2), 2), gp[0]);
1185 }
1186
1187 ARTS_ASSERT(isfinite(g));
1188
1189 // A dirty trick to avoid copying ppath_step
1190 const Index np = ip + 1;
1191 ppath_step.np = np;
1192}
1193
1195 MatrixView evol_op,
1196 Vector& abs_vec_mono,
1197 Numeric& temperature,
1198 MatrixView ext_mat_mono,
1199 Rng& rng,
1200 Vector& rte_pos,
1201 Vector& rte_los,
1202 Vector& pnd_vec,
1203 Numeric& stot,
1204 Numeric& ttot,
1205 Ppath& ppath_step,
1206 Index& termination_flag,
1207 bool& inside_cloud,
1208 const Agenda& ppath_step_agenda,
1209 const Numeric& ppath_lmax,
1210 const Numeric& ppath_lraytrace,
1211 const Agenda& propmat_clearsky_agenda,
1212 const bool& anyptype_nonTotRan,
1213 const Index stokes_dim,
1214 const Index f_index,
1215 const Vector& f_grid,
1216 const Vector& Iprop,
1217 const Vector& p_grid,
1218 const Vector& lat_grid,
1219 const Vector& lon_grid,
1220 const Tensor3& z_field,
1221 const Vector& refellipsoid,
1222 const Matrix& z_surface,
1223 const Tensor3& t_field,
1224 const Tensor4& vmr_field,
1225 const ArrayOfIndex& cloudbox_limits,
1226 const Tensor4& pnd_field,
1227 const ArrayOfArrayOfSingleScatteringData& scat_data,
1228 const Verbosity& verbosity) {
1229 ArrayOfMatrix evol_opArray(2);
1230 ArrayOfMatrix ext_matArray(2);
1231 ArrayOfVector abs_vecArray(2);
1232 ArrayOfVector pnd_vecArray(2);
1233 Matrix ext_mat(stokes_dim, stokes_dim);
1234 Matrix incT(stokes_dim, stokes_dim, 0.0);
1235 Vector tArray(2);
1236 Matrix T(stokes_dim, stokes_dim);
1237 Numeric kI, kQ;
1238 Numeric ds, dt = -999, dl = -999;
1239 Index istep = 0; // Counter for number of steps
1240 Matrix old_evol_op(stokes_dim, stokes_dim);
1241 Vector local_rte_los(2);
1242
1244
1245 // Total path length starts at zero
1246 stot = 0.0;
1247 ttot = 0.0;
1248
1249 //at the start of the path the evolution operator is the identity matrix
1250 id_mat(evol_op);
1251 evol_opArray[1] = evol_op;
1252
1253 // range defining cloudbox
1254 Range p_range(cloudbox_limits[0],
1255 cloudbox_limits[1] - cloudbox_limits[0] + 1);
1256 Range lat_range(cloudbox_limits[2],
1257 cloudbox_limits[3] - cloudbox_limits[2] + 1);
1258 Range lon_range(cloudbox_limits[4],
1259 cloudbox_limits[5] - cloudbox_limits[4] + 1);
1260
1261 //initialise Ppath with ppath_start_stepping
1262 ppath_start_stepping(ppath_step,
1263 3,
1264 p_grid,
1265 lat_grid,
1266 lon_grid,
1267 z_field,
1268 refellipsoid,
1269 z_surface,
1270 0,
1271 cloudbox_limits,
1272 false,
1273 rte_pos,
1274 rte_los,
1275 verbosity);
1276
1277 if (ppath_step.np == 0) {
1278 termination_flag = 1;
1279 return;
1280 }
1281
1282 // Index in ppath_step of end point considered presently
1283 Index ip = 0;
1284
1285 // Is point ip inside the cloudbox?
1286 inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
1287 ppath_step.gp_lat[ip],
1288 ppath_step.gp_lon[ip],
1289 cloudbox_limits,
1290 0,
1291 3);
1292
1293 // Determine radiative properties at point
1294 if (inside_cloud) {
1295 local_rte_los[0] = 180 - ppath_step.los(0, 0);
1296 local_rte_los[1] = ppath_step.los(0, 1) - 180;
1298 ext_mat_mono,
1299 abs_vec_mono,
1300 pnd_vec,
1301 temperature,
1302 propmat_clearsky_agenda,
1303 stokes_dim,
1304 f_index,
1305 f_grid,
1306 ppath_step.gp_p[0],
1307 ppath_step.gp_lat[0],
1308 ppath_step.gp_lon[0],
1309 p_grid[p_range],
1310 t_field(p_range, lat_range, lon_range),
1311 vmr_field(joker, p_range, lat_range, lon_range),
1312 pnd_field,
1313 scat_data,
1314 cloudbox_limits,
1315 local_rte_los);
1316 } else {
1318 ext_mat_mono,
1319 abs_vec_mono,
1320 temperature,
1321 propmat_clearsky_agenda,
1322 f_grid[f_index],
1323 ppath_step.gp_p[0],
1324 ppath_step.gp_lat[0],
1325 ppath_step.gp_lon[0],
1326 p_grid,
1327 t_field,
1328 vmr_field);
1329 pnd_vec = 0.0;
1330 }
1331
1332 // Move the data to end point containers
1333 ext_matArray[1] = ext_mat_mono;
1334 abs_vecArray[1] = abs_vec_mono;
1335 tArray[1] = temperature;
1336 pnd_vecArray[1] = pnd_vec;
1337
1338 //draw random number to determine end point
1339 Numeric r = rng.draw();
1340
1341 termination_flag = 0;
1342
1343 stot = ppath_step.end_lstep;
1344 ttot = ppath_step.end_lstep / SPEED_OF_LIGHT;
1345
1346 Numeric evop0 = 1;
1347 while ((evop0 > r) && (!termination_flag)) {
1348 istep++;
1349
1350 ARTS_USER_ERROR_IF (istep > 25000,
1351 "25000 path points have been reached. "
1352 "Is this an infinite loop?");
1353
1354 evol_opArray[0] = evol_opArray[1];
1355 ext_matArray[0] = ext_matArray[1];
1356 abs_vecArray[0] = abs_vecArray[1];
1357 tArray[0] = tArray[1];
1358 pnd_vecArray[0] = pnd_vecArray[1];
1359
1360 // Shall new ppath_step be calculated?
1361 if (ip >= ppath_step.np - 1) {
1362 ip = 1;
1364 ppath_step,
1365 ppath_lmax,
1366 ppath_lraytrace,
1367 f_grid[Range(f_index, 1)],
1368 ppath_step_agenda);
1369
1370 if (ppath_step.np <= 1) {
1371 termination_flag = 1;
1372 break;
1373 }
1374 inside_cloud = is_gp_inside_cloudbox(ppath_step.gp_p[ip],
1375 ppath_step.gp_lat[ip],
1376 ppath_step.gp_lon[ip],
1377 cloudbox_limits,
1378 0,
1379 3);
1380 } else {
1381 ip++;
1382 }
1383
1384 dl = ppath_step.lstep[ip - 1];
1385 dt = dl * 0.5 * (ppath_step.ngroup[ip - 1] + ppath_step.ngroup[ip]) /
1387 stot += dl;
1388 ttot += dt;
1389 if (inside_cloud) {
1390 local_rte_los[0] = 180 - ppath_step.los(ip, 0);
1391 local_rte_los[1] = ppath_step.los(ip, 1) - 180;
1393 ext_mat_mono,
1394 abs_vec_mono,
1395 pnd_vec,
1396 temperature,
1397 propmat_clearsky_agenda,
1398 stokes_dim,
1399 f_index,
1400 f_grid,
1401 ppath_step.gp_p[ip],
1402 ppath_step.gp_lat[ip],
1403 ppath_step.gp_lon[ip],
1404 p_grid[p_range],
1405 t_field(p_range, lat_range, lon_range),
1406 vmr_field(joker, p_range, lat_range, lon_range),
1407 pnd_field,
1408 scat_data,
1409 cloudbox_limits,
1410 local_rte_los);
1411 } else {
1413 ext_mat_mono,
1414 abs_vec_mono,
1415 temperature,
1416 propmat_clearsky_agenda,
1417 f_grid[f_index],
1418 ppath_step.gp_p[ip],
1419 ppath_step.gp_lat[ip],
1420 ppath_step.gp_lon[ip],
1421 p_grid,
1422 t_field,
1423 vmr_field);
1424 pnd_vec = 0.0;
1425 }
1426
1427 ext_matArray[1] = ext_mat_mono;
1428 abs_vecArray[1] = abs_vec_mono;
1429 tArray[1] = temperature;
1430 pnd_vecArray[1] = pnd_vec;
1431 ext_mat = ext_matArray[1];
1432 ext_mat += ext_matArray[0]; // Factor 2 fixed by using dl/2
1433 //
1434 {
1435 Index extmat_case = 0;
1436 ext2trans(incT, extmat_case, ext_mat, dl / 2);
1437 }
1438 //
1439
1440 mult(evol_op, incT, evol_opArray[0]);
1441 evol_opArray[1] = evol_op;
1442 evop0 = evol_op(0, 0);
1443
1444 // Handle cross-talk for ptype==30
1445 if (stokes_dim > 1 && anyptype_nonTotRan) {
1446 const Numeric Q1 = evol_op(0, 1) * Iprop[1] / Iprop[0];
1447 evop0 += Q1;
1448 }
1449 if (evop0 > r) {
1450 // Check whether hit ground or space.
1451 // path_step_agenda just detects surface intersections, and
1452 // if TOA is reached requires a special check.
1453 // But we are already ready if evol_op<=r
1454 if (ip >= ppath_step.np - 1) {
1455 if (ppath_what_background(ppath_step)) {
1456 termination_flag = 2;
1457 } //we have hit the surface
1458 else if (fractional_gp(ppath_step.gp_p[ip]) >=
1459 (Numeric)(p_grid.nelem() - 1) - 1e-3) {
1460 termination_flag = 1;
1461 } //we are at TOA
1462 }
1463 }
1464 } // while
1465
1466 if (termination_flag != 0) { //we must have reached the cloudbox boundary
1467 if (ip < ppath_step.np) {
1468 // This is not used if termination flag set,
1469 // so not an issue of ip is too large.
1470 // Need to think about this.
1471 rte_pos = ppath_step.pos(ip, joker);
1472 rte_los = ppath_step.los(ip, joker);
1473 }
1474 } else {
1475 //find position...and evol_op..and everything else required at the new
1476 //scattering/emission point
1477 const Numeric tol = 0.1; // Tolerance of 10 cm
1478 stot -= dl; // Take out last step because we are between stepping points
1479 ttot -= dt; // Take out last step because we are between stepping points
1480 kI = ext_mat(0, 0) / 2; // Factor 2 as sum of end point values
1481 kQ = ext_mat(0, 1) / 2; // Factor 2 as sum of end point values
1482 if (anyptype_nonTotRan) {
1483 const Numeric I1 = evol_opArray[0](0, 0);
1484 const Numeric Q1 = evol_opArray[0](0, 1) * Iprop[1] / Iprop[0];
1485
1486 // Need to use root finding to solve for ds
1487 brent_zero(
1488 ds, (Numeric)0.0, ppath_step.lstep[ip - 1], tol, r, I1, Q1, kI, kQ);
1489 } else {
1490 // Simple inversion when no cross-talk between I and Q
1491 ds = log(evol_opArray[0](0, 0) / r) / kI;
1492 }
1493 stot += ds;
1494 ttot += ds * dt / dl;
1495 Vector x(2, 0.0);
1496
1497 //interpolate atmospheric variables required later
1498 ArrayOfGridPos gp(1);
1499 x[1] = dl;
1500 Vector itw(2);
1501 gridpos(gp, x, ds);
1502 ARTS_ASSERT(gp[0].idx == 0);
1503 interpweights(itw, gp[0]);
1504 interp(ext_mat_mono, itw, ext_matArray, gp[0]);
1505 ext_mat = ext_mat_mono;
1506 ext_mat += ext_matArray[gp[0].idx];
1507 //
1508 {
1509 Index extmat_case = 0;
1510 ext2trans(incT, extmat_case, ext_mat, ds / 2);
1511 }
1512 //
1513 mult(evol_op, incT, evol_opArray[gp[0].idx]);
1514 interp(abs_vec_mono, itw, abs_vecArray, gp[0]);
1515 temperature = interp(itw, tArray, gp[0]);
1516 interp(pnd_vec, itw, pnd_vecArray, gp[0]);
1517 for (Index i = 0; i < 2; i++) {
1518 rte_pos[i] = interp(itw, ppath_step.pos(Range(ip - 1, 2), i), gp[0]);
1519 rte_los[i] = interp(itw, ppath_step.los(Range(ip - 1, 2), i), gp[0]);
1520 }
1521 rte_pos[2] = interp(itw, ppath_step.pos(Range(ip - 1, 2), 2), gp[0]);
1522 }
1523
1524 // A dirty trick to avoid copying ppath_step
1525 const Index np = ip + 1;
1526 ppath_step.np = np;
1527}
1528
1529void Sample_los(VectorView new_rte_los,
1530 Numeric& g_los_csc_theta,
1531 MatrixView Z,
1532 Rng& rng,
1533 ConstVectorView rte_los,
1534 const ArrayOfArrayOfSingleScatteringData& scat_data,
1535 const Index f_index,
1536 const Index stokes_dim,
1537 ConstVectorView pnd_vec,
1538 ConstVectorView Z11maxvector,
1539 const Numeric Csca,
1540 const Numeric rtp_temperature,
1541 const Index t_interp_order) {
1542 Numeric Z11max = 0;
1543 bool tryagain = true;
1544
1545 Vector sca_dir;
1546 mirror_los(sca_dir, rte_los, 3);
1547
1548 // Rejection method http://en.wikipedia.org/wiki/Rejection_sampling
1549 Index np = pnd_vec.nelem();
1550 ARTS_ASSERT(TotalNumberOfElements(scat_data) == np);
1551 for (Index i = 0; i < np; i++) {
1552 Z11max += Z11maxvector[i] * pnd_vec[i];
1553 }
1554
1556 // allocating variables needed for pha_mat extraction - this seems a little
1557 // disadvantageous. If we move the whole function back into the calling one
1558 // (it's used only at one place anyways), we can do the allocation there once
1559 // and avoid that it is done everytime this function is called within a loop.
1560 ArrayOfArrayOfTensor6 pha_mat_Nse;
1561 ArrayOfArrayOfIndex ptypes_Nse;
1562 Matrix t_ok;
1563 ArrayOfTensor6 pha_mat_ssbulk;
1564 ArrayOfIndex ptype_ssbulk;
1565 Tensor6 pha_mat_bulk;
1566 Index ptype_bulk;
1567 Matrix pdir(1, 2), idir(1, 2);
1568 Vector t(1, rtp_temperature);
1569 Matrix pnds(np, 1);
1570 pnds(joker, 0) = pnd_vec;
1571
1572 while (tryagain) {
1573 new_rte_los[0] = Conversion::acosd(1 - 2 * rng.draw());
1574 new_rte_los[1] = rng.draw() * 360 - 180;
1575
1576 //Calculate Phase matrix////////////////////////////////
1577 Vector inc_dir;
1578 mirror_los(inc_dir, new_rte_los, 3);
1579
1580 //pha_mat_singleCalc( Z, sca_dir[0], sca_dir[1], inc_dir[0], inc_dir[1],
1581 // scat_data_mono, stokes_dim, pnd_vec, rtp_temperature,
1582 // verbosity );
1583
1584 pdir(0, joker) = sca_dir;
1585 idir(0, joker) = inc_dir;
1586 pha_mat_NScatElems(pha_mat_Nse,
1587 ptypes_Nse,
1588 t_ok,
1589 scat_data,
1590 stokes_dim,
1591 t,
1592 pdir,
1593 idir,
1594 f_index,
1595 t_interp_order);
1597 pha_mat_ssbulk, ptype_ssbulk, pha_mat_Nse, ptypes_Nse, pnds, t_ok);
1598 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
1599 Z = pha_mat_bulk(0, 0, 0, 0, joker, joker);
1600
1601 if (rng.draw() <= Z(0, 0) / Z11max) //then new los is accepted
1602 {
1603 tryagain = false;
1604 }
1605 }
1606 g_los_csc_theta = Z(0, 0) / Csca;
1607}
1608
1609void Sample_los_uniform(VectorView new_rte_los, Rng& rng) {
1610 new_rte_los[1] = rng.draw() * 360 - 180;
1611 new_rte_los[0] = Conversion::acosd(1 - 2 * rng.draw());
1612}
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:241
Constants of physical expressions as constexpr.
Common ARTS conversions.
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:25619
void propmat_clearsky_agendaExecute(Workspace &ws, PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfSpeciesTag &select_abs_species, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:25655
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Tensor3.
Definition: matpackIII.h:130
A constant view of a Tensor4.
Definition: matpackIV.h:131
Index nbooks() const noexcept
Definition: matpackIV.h:139
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
The range class.
Definition: matpackI.h:159
Definition: rng.h:554
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
Stokes vector is as Propagation matrix but only has 4 possible values.
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
The Tensor5 class.
Definition: matpackV.h:516
The Tensor6 class.
Definition: matpackVI.h:1099
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
Workspace class.
Definition: workspace_ng.h:53
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
Definition: cloudbox.cc:461
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:539
#define abs(x)
#define q
#define ns
#define min(a, b)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:396
bool is_singular(ConstMatrixView A)
Checks if a square matrix is singular.
Definition: logic.cc:298
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
const Joker joker
Interpolation classes and functions created for use within Monte Carlo scattering simulations.
#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 ext_mat_case(Index &icase, ConstMatrixView ext_mat, const Index stokes_dim)
Definition: montecarlo.cc:543
void Sample_los_uniform(VectorView new_rte_los, Rng &rng)
Sample_los_uniform.
Definition: montecarlo.cc:1609
void clear_rt_vars_at_gp(Workspace &ws, MatrixView ext_mat_mono, VectorView abs_vec_mono, Numeric &temperature, const Agenda &propmat_clearsky_agenda, const Numeric &f_mono, const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, ConstVectorView p_grid, ConstTensor3View t_field, ConstTensor4View vmr_field)
clear_rt_vars_at_gp.
Definition: montecarlo.cc:246
constexpr Numeric SPEED_OF_LIGHT
Definition: montecarlo.cc:40
void cloud_atm_vars_by_gp(VectorView pressure, VectorView temperature, MatrixView vmr, MatrixView pnd, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, ConstVectorView p_grid_cloud, ConstTensor3View t_field_cloud, ConstTensor4View vmr_field_cloud, ConstTensor4View pnd_field)
cloud_atm_vars_by_gp.
Definition: montecarlo.cc:458
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Converts an extinction matrix to a transmission matrix.
Definition: montecarlo.cc:617
void brent_zero(Numeric &sb, const Numeric &a, const Numeric &b, const Numeric &t, const Numeric &rn, const Numeric &I, const Numeric &Q, const Numeric &KI, const Numeric &KQ)
brent_zero.
Definition: montecarlo.cc:134
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
Numeric ext_I(const Numeric &I, const Numeric &Q, const Numeric &kI, const Numeric &kQ, const Numeric &s)
ext_I.
Definition: montecarlo.cc:60
void cloudy_rt_vars_at_gp(Workspace &ws, MatrixView ext_mat_mono, VectorView abs_vec_mono, VectorView pnd_vec, Numeric &temperature, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Index f_index, const Vector &f_grid, const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, ConstVectorView p_grid_cloud, ConstTensor3View t_field_cloud, ConstTensor4View vmr_field_cloud, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfIndex &cloudbox_limits, const Vector &rte_los)
cloudy_rt_vars_at_gp.
Definition: montecarlo.cc:328
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 k
Boltzmann constant convenience name [J/K].
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 Numeric e
Elementary charge convenience name [C].
auto acosd(auto x) noexcept
Returns rad2deg of the arc-cosine of the input.
constexpr Numeric dl(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
Numeric Q(Numeric T, const Species::IsotopeRecord &ir)
Definition: partfun.cc:14
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const PropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
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 opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
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.
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1477
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
Definition: ppath.cc:1451
void ppath_start_stepping(Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const bool &ppath_inside_cloudbox_do, ConstVectorView rte_pos, ConstVectorView rte_los, const Verbosity &verbosity)
Initiates a Ppath structure for calculation of a path with ppath_step.
Definition: ppath.cc:4477
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:1812
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Structure to store a grid position.
Definition: interpolation.h:73
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Matrix los
Line-of-sight at each ppath point.
Definition: ppath_struct.h:35
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
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
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
#define d
#define a
#define c
#define b