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