ARTS 2.5.11 (git: 6827797f)
optproperties.cc
Go to the documentation of this file.
1/*===========================================================================
2 === File description
3 ===========================================================================*/
4
16/*===========================================================================
17 === External declarations
18 ===========================================================================*/
19
20#include "optproperties.h"
21#include <cfloat>
22#include <cmath>
23#include <stdexcept>
24#include "array.h"
25#include "arts.h"
26#include "check_input.h"
27#include "arts_conversions.h"
28#include "interpolation.h"
29#include "interp.h"
30#include "logic.h"
31#include "math_funcs.h"
32#include "matpack_data.h"
33#include "messages.h"
34#include "xml_io.h"
35
36inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
37inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
38using Constant::pi;
39inline constexpr Numeric PI=pi;
40
41#define F11 pha_mat_int[0]
42#define F12 pha_mat_int[1]
43#define F22 pha_mat_int[2]
44#define F33 pha_mat_int[3]
45#define F34 pha_mat_int[4]
46#define F44 pha_mat_int[5]
47
49
58/*
59void methodname(//Output
60 type& name,
61 //Input
62 const type& name)
63{
64}
65*/
66
68
88void opt_prop_Bulk( //Output
89 Tensor5& ext_mat, // (nf,nT,ndir,nst,nst)
90 Tensor4& abs_vec, // (nf,nT,ndir,nst)
91 Index& ptype,
92 //Input
93 const ArrayOfTensor5& ext_mat_ss, // [nss](nf,nT,ndir,nst,nst)
94 const ArrayOfTensor4& abs_vec_ss, // [nss](nf,nT,ndir,nst)
95 const ArrayOfIndex& ptypes_ss) {
96 ARTS_ASSERT(ext_mat_ss.nelem() == abs_vec_ss.nelem());
97
98 ext_mat = ext_mat_ss[0];
99 abs_vec = abs_vec_ss[0];
100
101 for (Index i_ss = 1; i_ss < ext_mat_ss.nelem(); i_ss++) {
102 ext_mat += ext_mat_ss[i_ss];
103 abs_vec += abs_vec_ss[i_ss];
104 }
105 ptype = max(ptypes_ss);
106}
107
109
135 ArrayOfTensor5& ext_mat, // [nss](nf,nT,ndir,nst,nst)
136 ArrayOfTensor4& abs_vec, // [nss](nf,nT,ndir,nst)
137 ArrayOfIndex& ptype,
138 //Input
139 const ArrayOfArrayOfTensor5& ext_mat_se, // [nss][nse](nf,nT,ndir,nst,nst)
140 const ArrayOfArrayOfTensor4& abs_vec_se, // [nss][nse](nf,nT,ndir,nst)
141 const ArrayOfArrayOfIndex& ptypes_se,
142 ConstMatrixView pnds,
143 ConstMatrixView t_ok) {
144 ARTS_ASSERT(t_ok.nrows() == pnds.nrows());
145 ARTS_ASSERT(t_ok.ncols() == pnds.ncols());
146 ARTS_ASSERT(TotalNumberOfElements(ext_mat_se) == pnds.nrows());
147 ARTS_ASSERT(TotalNumberOfElements(abs_vec_se) == pnds.nrows());
148 ARTS_ASSERT(ext_mat_se.nelem() == abs_vec_se.nelem());
149
150 const Index nT = pnds.ncols();
151 const Index nf = abs_vec_se[0][0].nbooks();
152 const Index nDir = abs_vec_se[0][0].nrows();
153 const Index stokes_dim = abs_vec_se[0][0].ncols();
154
155 const Index nss = ext_mat_se.nelem();
156 ext_mat.resize(nss);
157 abs_vec.resize(nss);
158 ptype.resize(nss);
159 Tensor4 ext_tmp;
160 Tensor3 abs_tmp;
161
162 Index i_se_flat = 0;
163
164 for (Index i_ss = 0; i_ss < nss; i_ss++) {
165 ARTS_ASSERT(ext_mat_se[i_ss].nelem() == abs_vec_se[i_ss].nelem());
166 ARTS_ASSERT(nT == ext_mat_se[i_ss][0].nbooks());
167 ARTS_ASSERT(nT == abs_vec_se[i_ss][0].npages());
168
169 ext_mat[i_ss].resize(nf, nT, nDir, stokes_dim, stokes_dim);
170 ext_mat[i_ss] = 0.;
171 abs_vec[i_ss].resize(nf, nT, nDir, stokes_dim);
172 abs_vec[i_ss] = 0.;
173
174 for (Index i_se = 0; i_se < ext_mat_se[i_ss].nelem(); i_se++) {
175 ARTS_ASSERT(nT == ext_mat_se[i_ss][i_se].nbooks());
176 ARTS_ASSERT(nT == abs_vec_se[i_ss][i_se].npages());
177
178 for (Index Tind = 0; Tind < nT; Tind++) {
179 if (pnds(i_se_flat, Tind) != 0.) {
180 if (t_ok(i_se_flat, Tind) > 0.) {
181 ext_tmp = ext_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker);
182 ext_tmp *= pnds(i_se_flat, Tind);
183 ext_mat[i_ss](joker, Tind, joker, joker, joker) += ext_tmp;
184
185 abs_tmp = abs_vec_se[i_ss][i_se](joker, Tind, joker, joker);
186 abs_tmp *= pnds(i_se_flat, Tind);
187 abs_vec[i_ss](joker, Tind, joker, joker) += abs_tmp;
188 } else {
190 "Interpolation error for (flat-array) scattering element #",
191 i_se_flat, "\n"
192 "at location/temperature point #", Tind, "\n")
193 }
194 }
195 }
196 i_se_flat++;
197 }
198 ptype[i_ss] = max(ptypes_se[i_ss]);
199 }
200}
201
203
234 ArrayOfArrayOfTensor5& ext_mat, // [nss][nse](nf,nT,ndir,nst,nst)
235 ArrayOfArrayOfTensor4& abs_vec, // [nss][nse](nf,nT,ndir,nst)
236 ArrayOfArrayOfIndex& ptypes,
237 Matrix& t_ok,
238 //Input
239 const ArrayOfArrayOfSingleScatteringData& scat_data,
240 const Index& stokes_dim,
241 const Vector& T_array,
242 const Matrix& dir_array,
243 const Index& f_index,
244 const Index& t_interp_order) {
245 Index f_start, nf;
246 if (f_index < 0) {
247 nf = scat_data[0][0].ext_mat_data.nshelves();
248 f_start = 0;
249 //f_end = f_start+nf;
250 } else {
251 nf = 1;
252 if (scat_data[0][0].ext_mat_data.nshelves() == 1)
253 f_start = 0;
254 else
255 f_start = f_index;
256 //f_end = f_start+nf;
257 }
258
259 const Index nT = T_array.nelem();
260 const Index nDir = dir_array.nrows();
261
262 const Index nss = scat_data.nelem();
263 ext_mat.resize(nss);
264 abs_vec.resize(nss);
265 ptypes.resize(nss);
266
267 const Index Nse_all = TotalNumberOfElements(scat_data);
268 t_ok.resize(Nse_all, nT);
269 Index i_se_flat = 0;
270
271 for (Index i_ss = 0; i_ss < nss; i_ss++) {
272 Index nse = scat_data[i_ss].nelem();
273 ext_mat[i_ss].resize(nse);
274 abs_vec[i_ss].resize(nse);
275 ptypes[i_ss].resize(nse);
276
277 for (Index i_se = 0; i_se < nse; i_se++) {
278 ext_mat[i_ss][i_se].resize(nf, nT, nDir, stokes_dim, stokes_dim);
279 abs_vec[i_ss][i_se].resize(nf, nT, nDir, stokes_dim);
280
281 opt_prop_1ScatElem(ext_mat[i_ss][i_se],
282 abs_vec[i_ss][i_se],
283 ptypes[i_ss][i_se],
284 t_ok(i_se_flat, joker),
285 scat_data[i_ss][i_se],
286 T_array,
287 dir_array,
288 f_start,
289 t_interp_order);
290 i_se_flat++;
291 }
292 }
293 ARTS_ASSERT(i_se_flat == Nse_all);
294}
295
297
314ArrayOfLagrangeInterpolation ssd_tinterp_parameters( //Output
315 VectorView t_ok,
316 Index& this_T_interp_order,
317 //Input
318 ConstVectorView T_grid,
319 const Vector& T_array,
320 const Index& t_interp_order) {
321 const Index nTin = T_grid.nelem();
322 const Index nTout = T_array.nelem();
323
324 this_T_interp_order = -1;
325
326
327 if (nTin > 1) {
328 this_T_interp_order = min(t_interp_order, nTin - 1);
329
330 // we need to check T-grid exceedance. and catch these cases (because T
331 // is assumed to correspond to a location and T-exceedance is ok when pnd=0
332 // there. however, here we do not know pnd.) and report them to have the
333 // calling method dealing with it.
334
335 // we set the extrapolfax explicitly and use it here as well as in
336 // gridpos_poly below.
337 const Numeric extrapolfac = 0.5;
338 const Numeric lowlim = T_grid[0] - extrapolfac * (T_grid[1] - T_grid[0]);
339 const Numeric uplim =
340 T_grid[nTin - 1] + extrapolfac * (T_grid[nTin - 1] - T_grid[nTin - 2]);
341
342 bool any_T_exceed = false;
343 for (Index Tind = 0; Tind < nTout; Tind++) {
344 if (T_array[Tind] < lowlim || T_array[Tind] > uplim) {
345 t_ok[Tind] = -1.;
346 any_T_exceed = true;
347 } else
348 t_ok[Tind] = 1.;
349 }
350
351 if (any_T_exceed) {
352 // Reserve output
353 ArrayOfLagrangeInterpolation T_lag;
354 T_lag.reserve(nTout);
355
356 bool grid_unchecked = true;
357
358 for (Index iT = 0; iT < nTout; iT++) {
359 if (t_ok[iT] < 0) {
360 T_lag.emplace_back(this_T_interp_order);
361 } else {
362 if (grid_unchecked) {
364 "Temperature interpolation in pha_mat_1ScatElem",
365 T_grid,
366 T_array[Range(iT, 1)],
367 this_T_interp_order);
368 grid_unchecked = false;
369 }
370 T_lag.emplace_back(0, T_array[iT], T_grid, this_T_interp_order);
371 }
372 }
373 return T_lag;
374 } else {
375 return my_interp::lagrange_interpolation_list<LagrangeInterpolation>(T_array, T_grid, this_T_interp_order, extrapolfac);
376 }
377 } else {
378 t_ok = 1.;
379 return {};
380 }
381}
382
384
407void opt_prop_1ScatElem( //Output
408 Tensor5View ext_mat, // nf, nT, ndir, nst, nst
409 Tensor4View abs_vec, // nf, nT, ndir, nst
410 Index& ptype,
411 VectorView t_ok,
412 //Input
413 const SingleScatteringData& ssd,
414 const Vector& T_array,
415 const Matrix& dir_array,
416 const Index& f_start,
417 const Index& t_interp_order) {
418 // FIXME: this is prob best done in scat_data_checkedCalc (or
419 // cloudbox_checkedCalc) to have it done once and for all. Here at max ARTS_ASSERT.
420
421 // At very first check validity of the scatt elements ptype (so far we only
422 // handle PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND).
423 /*
424 if( ssd.ptype != PTYPE_TOTAL_RND and ssd.ptype != PTYPE_AZIMUTH_RND )
425 {
426 ostringstream os;
427 os << "Only ptypes " << PTYPE_TOTAL_RND << " and " << PTYPE_AZIMUTH_RND
428 << " can be handled.\n"
429 << "Encountered scattering element with ptype " << ssd.ptype
430 << ", though.";
431 throw runtime_error( os.str() );
432 }
433 */
434
436
437 const Index nf = ext_mat.nshelves();
438 ARTS_ASSERT(abs_vec.nbooks() == nf);
439 if (nf > 1) {
440 ARTS_ASSERT(nf == ssd.f_grid.nelem());
441 }
442
443 const Index nTout = T_array.nelem();
444 ARTS_ASSERT(ext_mat.nbooks() == nTout);
445 ARTS_ASSERT(abs_vec.npages() == nTout);
446 ARTS_ASSERT(t_ok.nelem() == nTout);
447
448 const Index nDir = dir_array.nrows();
449 ARTS_ASSERT(ext_mat.npages() == nDir);
450 ARTS_ASSERT(abs_vec.nrows() == nDir);
451
452 const Index stokes_dim = abs_vec.ncols();
453 ARTS_ASSERT(ext_mat.nrows() == stokes_dim);
454 ARTS_ASSERT(ext_mat.ncols() == stokes_dim);
455
456 ptype = ssd.ptype;
457
458 // Determine T-interpol order as well as interpol positions and weights (they
459 // are the same for all directions (and freqs), ie it is sufficient to
460 // calculate them once).
461 const Index nTin = ssd.T_grid.nelem();
462 Index this_T_interp_order;
463 const auto T_lag = ssd_tinterp_parameters(t_ok,
464 this_T_interp_order,
465 ssd.T_grid,
466 T_array,
467 t_interp_order);
468 const auto T_itw_lag = interpweights(T_lag);
469
470 // Now loop over requested directions (and apply simultaneously for all freqs):
471 // 1) extract/interpolate direction (but not for tot.random)
472 // 2) apply T-interpol
473 if (ptype == PTYPE_TOTAL_RND)
474 if (this_T_interp_order < 0) // just extract (and unpack) ext and abs data
475 // for the fs, the Tin, and the dirin and sort
476 // (copy) into the output fs, Ts, and dirs.
477 {
478 Tensor3 ext_mat_tmp(nf, stokes_dim, stokes_dim);
479 Matrix abs_vec_tmp(nf, stokes_dim);
480 for (Index find = 0; find < nf; find++) {
481 ext_mat_SSD2Stokes(ext_mat_tmp(find, joker, joker),
482 ssd.ext_mat_data(find + f_start, 0, 0, 0, joker),
483 stokes_dim,
484 ptype);
485 abs_vec_SSD2Stokes(abs_vec_tmp(find, joker),
486 ssd.abs_vec_data(find + f_start, 0, 0, 0, joker),
487 stokes_dim,
488 ptype);
489 }
490 for (Index Tind = 0; Tind < nTout; Tind++)
491 for (Index dind = 0; dind < nDir; dind++) {
492 ext_mat(joker, Tind, dind, joker, joker) = ext_mat_tmp;
493 abs_vec(joker, Tind, dind, joker) = abs_vec_tmp;
494 }
495 } else // T-interpolation required (but not dir). To be done on the compact
496 // ssd format.
497 {
498 Tensor4 ext_mat_tmp(nf, nTout, stokes_dim, stokes_dim);
499 Tensor3 abs_vec_tmp(nf, nTout, stokes_dim);
500 Matrix ext_mat_tmp_ssd(nTout, ssd.ext_mat_data.ncols());
501 Matrix abs_vec_tmp_ssd(nTout, ssd.abs_vec_data.ncols());
502 for (Index find = 0; find < nf; find++) {
503 for (Index nst = 0; nst < ext_mat_tmp_ssd.ncols(); nst++) {
504 reinterp(ext_mat_tmp_ssd(joker, nst),
505 ssd.ext_mat_data(find + f_start, joker, 0, 0, nst),
506 T_itw_lag, T_lag);
507 }
508 for (Index Tind = 0; Tind < nTout; Tind++)
509 if (t_ok[Tind] > 0.)
510 ext_mat_SSD2Stokes(ext_mat_tmp(find, Tind, joker, joker),
511 ext_mat_tmp_ssd(Tind, joker),
512 stokes_dim,
513 ptype);
514 else
515 ext_mat_tmp(find, Tind, joker, joker) = 0.;
516
517 for (Index nst = 0; nst < abs_vec_tmp_ssd.ncols(); nst++) {
518 reinterp(abs_vec_tmp_ssd(joker, nst),
519 ssd.abs_vec_data(find + f_start, joker, 0, 0, nst),
520 T_itw_lag,
521 T_lag);
522 }
523 for (Index Tind = 0; Tind < nTout; Tind++)
524 if (t_ok[Tind] > 0.)
525 abs_vec_SSD2Stokes(abs_vec_tmp(find, Tind, joker),
526 abs_vec_tmp_ssd(Tind, joker),
527 stokes_dim,
528 ptype);
529 else
530 abs_vec_tmp(find, Tind, joker) = 0.;
531 }
532
533 for (Index dind = 0; dind < nDir; dind++) {
534 ext_mat(joker, joker, dind, joker, joker) = ext_mat_tmp;
535 abs_vec(joker, joker, dind, joker) = abs_vec_tmp;
536 }
537 }
538
539 else // dir-interpolation for non-tot-random particles
540 {
541 // as we don't allow other than az.random here, ext and abs will only vary
542 // with za, not with aa. Hence, we could be smart here and calc data only
543 // for unique za (and copy the rest). however, this smartness might cost as
544 // well. so for now, we leave that and blindly loop over the given direction
545 // array, and leave smartness in setting up directional array to the calling
546 // methods.
547
548 // derive the direction interpolation weights.
549 ArrayOfGridPos dir_gp(nDir);
550 gridpos(dir_gp, ssd.za_grid, dir_array(joker, 0));
551 Matrix dir_itw(nDir, 2); // only interpolating in za, ie 1D linear interpol
552 interpweights(dir_itw, dir_gp);
553
554 const Index next = ssd.ext_mat_data.ncols();
555 const Index nabs = ssd.abs_vec_data.ncols();
556
557 if (this_T_interp_order < 0) // T only needs to be extracted.
558 {
559 Matrix ext_mat_tmp_ssd(nDir, next);
560 Matrix abs_vec_tmp_ssd(nDir, nabs);
561 Tensor4 ext_mat_tmp(nf, nDir, stokes_dim, stokes_dim);
562 Tensor3 abs_vec_tmp(nf, nDir, stokes_dim);
563 for (Index find = 0; find < nf; find++) {
564 for (Index nst = 0; nst < next; nst++)
565 interp(ext_mat_tmp_ssd(joker, nst),
566 dir_itw,
567 ssd.ext_mat_data(find + f_start, 0, joker, 0, nst),
568 dir_gp);
569 for (Index Dind = 0; Dind < nDir; Dind++)
570 ext_mat_SSD2Stokes(ext_mat_tmp(find, Dind, joker, joker),
571 ext_mat_tmp_ssd(Dind, joker),
572 stokes_dim,
573 ptype);
574
575 for (Index nst = 0; nst < nabs; nst++)
576 interp(abs_vec_tmp_ssd(joker, nst),
577 dir_itw,
578 ssd.abs_vec_data(find + f_start, 0, joker, 0, nst),
579 dir_gp);
580 for (Index Dind = 0; Dind < nDir; Dind++)
581 abs_vec_SSD2Stokes(abs_vec_tmp(find, Dind, joker),
582 abs_vec_tmp_ssd(Dind, joker),
583 stokes_dim,
584 ptype);
585 }
586
587 for (Index Tind = 0; Tind < nTout; Tind++) {
588 ext_mat(joker, Tind, joker, joker, joker) = ext_mat_tmp;
589 abs_vec(joker, Tind, joker, joker) = abs_vec_tmp;
590 }
591 } else // T- and dir-interpolation required. To be done on the compact ssd
592 // format.
593 {
594 Tensor3 ext_mat_tmp_ssd(nTin, nDir, next);
595 Tensor3 abs_vec_tmp_ssd(nTin, nDir, nabs);
596 Matrix ext_mat_tmp(nTout, next);
597 Matrix abs_vec_tmp(nTout, nabs);
598
599 for (Index find = 0; find < nf; find++) {
600 for (Index Tind = 0; Tind < nTin; Tind++) {
601 for (Index nst = 0; nst < next; nst++)
602 interp(ext_mat_tmp_ssd(Tind, joker, nst),
603 dir_itw,
604 ssd.ext_mat_data(find + f_start, Tind, joker, 0, nst),
605 dir_gp);
606 for (Index nst = 0; nst < nabs; nst++)
607 interp(abs_vec_tmp_ssd(Tind, joker, nst),
608 dir_itw,
609 ssd.abs_vec_data(find + f_start, Tind, joker, 0, nst),
610 dir_gp);
611 }
612
613 for (Index Dind = 0; Dind < nDir; Dind++) {
614 for (Index nst = 0; nst < next; nst++) {
615 reinterp(ext_mat_tmp(joker, nst),
616 ext_mat_tmp_ssd(joker, Dind, nst),
617 T_itw_lag,
618 T_lag);
619 }
620 for (Index Tind = 0; Tind < nTout; Tind++)
621 ext_mat_SSD2Stokes(ext_mat(find, Tind, Dind, joker, joker),
622 ext_mat_tmp(Tind, joker),
623 stokes_dim,
624 ptype);
625
626 for (Index nst = 0; nst < nabs; nst++) {
627 reinterp(abs_vec_tmp(joker, nst),
628 abs_vec_tmp_ssd(joker, Dind, nst),
629 T_itw_lag,
630 T_lag);
631 }
632 for (Index Tind = 0; Tind < nTout; Tind++)
633 abs_vec_SSD2Stokes(abs_vec(find, Tind, Dind, joker),
634 abs_vec_tmp(Tind, joker),
635 stokes_dim,
636 ptype);
637 }
638 }
639 }
640 }
641}
642
644
657void ext_mat_SSD2Stokes( //Output
658 MatrixView ext_mat_stokes,
659 //Input
660 ConstVectorView ext_mat_ssd,
661 const Index& stokes_dim,
662 const Index& ptype) {
663 // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
664 // calling methods, though.
666
667 ext_mat_stokes = 0.;
668
669 for (Index ist = 0; ist < stokes_dim; ist++) {
670 ext_mat_stokes(ist, ist) = ext_mat_ssd[0];
671 }
672
673 if (ptype > PTYPE_TOTAL_RND) {
674 switch (stokes_dim) {
675 case 4: {
676 ext_mat_stokes(2, 3) = ext_mat_ssd[2];
677 ext_mat_stokes(3, 2) = -ext_mat_ssd[2];
678 } /* FALLTHROUGH */
679 case 3: {
680 // nothing to be done here. but we need this for executing the below
681 // also in case of stokes_dim==3.
682 }
683 case 2: {
684 ext_mat_stokes(0, 1) = ext_mat_ssd[1];
685 ext_mat_stokes(1, 0) = ext_mat_ssd[1];
686 }
687 }
688 }
689}
690
692
705void abs_vec_SSD2Stokes( //Output
706 VectorView abs_vec_stokes,
707 //Input
708 ConstVectorView abs_vec_ssd,
709 const Index& stokes_dim,
710 const Index& ptype) {
711 // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
712 // calling methods, though.
714
715 abs_vec_stokes = 0.;
716
717 abs_vec_stokes[0] = abs_vec_ssd[0];
718
719 if (ptype > PTYPE_TOTAL_RND and stokes_dim > 1) {
720 abs_vec_stokes[1] = abs_vec_ssd[1];
721 }
722}
723
725
742void pha_mat_Bulk( //Output
743 Tensor6& pha_mat, // (nf,nT,npdir,nidir,nst,nst)
744 Index& ptype,
745 //Input
746 const ArrayOfTensor6& pha_mat_ss, // [nss](nf,nT,npdir,nidir,nst,nst)
747 const ArrayOfIndex& ptypes_ss) {
748 pha_mat = pha_mat_ss[0];
749
750 for (Index i_ss = 1; i_ss < pha_mat_ss.nelem(); i_ss++)
751 pha_mat += pha_mat_ss[i_ss];
752
753 ptype = max(ptypes_ss);
754}
755
757
781 ArrayOfTensor6& pha_mat, // [nss](nf,nT,npdir,nidir,nst,nst)
782 ArrayOfIndex& ptype,
783 //Input
784 const ArrayOfArrayOfTensor6&
785 pha_mat_se, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
786 const ArrayOfArrayOfIndex& ptypes_se,
787 ConstMatrixView pnds,
788 ConstMatrixView t_ok) {
789 ARTS_ASSERT(t_ok.nrows() == pnds.nrows());
790 ARTS_ASSERT(t_ok.ncols() == pnds.ncols());
791 ARTS_ASSERT(TotalNumberOfElements(pha_mat_se) == pnds.nrows());
792
793 const Index nT = pnds.ncols();
794 const Index nf = pha_mat_se[0][0].nvitrines();
795 const Index npDir = pha_mat_se[0][0].nbooks();
796 const Index niDir = pha_mat_se[0][0].npages();
797 const Index stokes_dim = pha_mat_se[0][0].ncols();
798
799 const Index nss = pha_mat_se.nelem();
800 pha_mat.resize(nss);
801 ptype.resize(nss);
802 Tensor5 pha_tmp;
803
804 Index i_se_flat = 0;
805
806 for (Index i_ss = 0; i_ss < nss; i_ss++) {
807 ARTS_ASSERT(nT == pha_mat_se[i_ss][0].nshelves());
808
809 pha_mat[i_ss].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
810 pha_mat[i_ss] = 0.;
811
812 for (Index i_se = 0; i_se < pha_mat_se[i_ss].nelem(); i_se++) {
813 ARTS_ASSERT(nT == pha_mat_se[i_ss][i_se].nshelves());
814
815 for (Index Tind = 0; Tind < nT; Tind++) {
816 if (pnds(i_se_flat, Tind) != 0.) {
817 if (t_ok(i_se_flat, Tind) > 0.) {
818 pha_tmp =
819 pha_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker, joker);
820 pha_tmp *= pnds(i_se_flat, Tind);
821 pha_mat[i_ss](joker, Tind, joker, joker, joker, joker) += pha_tmp;
822 } else {
824 "Interpolation error for (flat-array) scattering element #",
825 i_se_flat, "\n"
826 "at location/temperature point #", Tind, "\n")
827 }
828 }
829 }
830 i_se_flat++;
831 }
832 ptype[i_ss] = max(ptypes_se[i_ss]);
833 }
834}
835
837
867void pha_mat_NScatElems( //Output
868 ArrayOfArrayOfTensor6& pha_mat, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
869 ArrayOfArrayOfIndex& ptypes,
870 Matrix& t_ok,
871 //Input
872 const ArrayOfArrayOfSingleScatteringData& scat_data,
873 const Index& stokes_dim,
874 const Vector& T_array,
875 const Matrix& pdir_array,
876 const Matrix& idir_array,
877 const Index& f_index,
878 const Index& t_interp_order) {
879 Index f_start, nf;
880 if (f_index < 0) {
881 nf = scat_data[0][0].pha_mat_data.nlibraries();
882 f_start = 0;
883 //f_end = f_start+nf;
884 } else {
885 nf = 1;
886 if (scat_data[0][0].pha_mat_data.nlibraries() == 1)
887 f_start = 0;
888 else
889 f_start = f_index;
890 //f_end = f_start+nf;
891 }
892
893 const Index nT = T_array.nelem();
894 const Index npDir = pdir_array.nrows();
895 const Index niDir = idir_array.nrows();
896
897 const Index nss = scat_data.nelem();
898 pha_mat.resize(nss);
899 ptypes.resize(nss);
900
901 const Index Nse_all = TotalNumberOfElements(scat_data);
902 t_ok.resize(Nse_all, nT);
903 Index i_se_flat = 0;
904
905 for (Index i_ss = 0; i_ss < nss; i_ss++) {
906 Index nse = scat_data[i_ss].nelem();
907 pha_mat[i_ss].resize(nse);
908 ptypes[i_ss].resize(nse);
909
910 for (Index i_se = 0; i_se < nse; i_se++) {
911 pha_mat[i_ss][i_se].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
912
913 pha_mat_1ScatElem(pha_mat[i_ss][i_se],
914 ptypes[i_ss][i_se],
915 t_ok(i_se_flat, joker),
916 scat_data[i_ss][i_se],
917 T_array,
918 pdir_array,
919 idir_array,
920 f_start,
921 t_interp_order);
922 i_se_flat++;
923 }
924 }
925 ARTS_ASSERT(i_se_flat == Nse_all);
926}
927
929
952void pha_mat_1ScatElem( //Output
953 Tensor6View pha_mat, // nf, nT, npdir, nidir, nst, nst
954 Index& ptype,
955 VectorView t_ok,
956 //Input
957 const SingleScatteringData& ssd,
958 const Vector& T_array,
959 const Matrix& pdir_array,
960 const Matrix& idir_array,
961 const Index& f_start,
962 const Index& t_interp_order) {
964
965 const Index nf = pha_mat.nvitrines();
966 if (nf > 1) {
967 ARTS_ASSERT(nf == ssd.f_grid.nelem());
968 }
969
970 const Index nTout = T_array.nelem();
971 ARTS_ASSERT(pha_mat.nshelves() == nTout);
972 ARTS_ASSERT(t_ok.nelem() == nTout);
973
974 const Index npDir = pdir_array.nrows();
975 ARTS_ASSERT(pha_mat.nbooks() == npDir);
976
977 const Index niDir = idir_array.nrows();
978 ARTS_ASSERT(pha_mat.npages() == niDir);
979
980 const Index stokes_dim = pha_mat.ncols();
981 ARTS_ASSERT(pha_mat.nrows() == stokes_dim);
982
983 ptype = ssd.ptype;
984
985 // Determine T-interpol order as well as interpol positions and weights (they
986 // are the same for all directions (and freqs), ie it is sufficient to
987 // calculate them once).
988 const Index nTin = ssd.T_grid.nelem();
989 Index this_T_interp_order;
990 const auto T_lag = ssd_tinterp_parameters(t_ok,
991 this_T_interp_order,
992 ssd.T_grid,
993 T_array,
994 t_interp_order);
995 const auto T_itw_lag = interpweights(T_lag);
996
997 // Now loop over requested directions (and apply simultaneously for all freqs):
998 // 1) interpolate direction
999 // 2) apply T-interpol
1000 if (ptype == PTYPE_TOTAL_RND) {
1001 // determine how many of the compact stokes elements we will need.
1002 // restrict interpolations to those.
1003 Index npha;
1004 if (stokes_dim == 1)
1005 npha = 1;
1006 else if (stokes_dim < 4) // stokes_dim==2 || stokes_dim==3
1007 npha = 4;
1008 else
1009 npha = 6;
1010 if (this_T_interp_order < 0) // just extract (and unpack) pha data for the
1011 // fs and Tin, and interpolate in sca-angs, and
1012 // sort (copy) into the output fs, Ts, and dirs.
1013 {
1014 for (Index pdir = 0; pdir < npDir; pdir++)
1015 for (Index idir = 0; idir < niDir; idir++) {
1016 // calc scat ang theta from incident and prop dirs
1017 Numeric theta = scat_angle(pdir_array(pdir, 0),
1018 pdir_array(pdir, 1),
1019 idir_array(idir, 0),
1020 idir_array(idir, 1));
1021
1022 // get scat angle interpolation weights
1023 GridPos dir_gp;
1024 gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
1025 Vector dir_itw(2);
1026 interpweights(dir_itw, dir_gp);
1027
1028 Vector pha_mat_int(npha, 0.);
1029 Matrix pha_mat_tmp(stokes_dim, stokes_dim);
1030 for (Index find = 0; find < nf; find++) {
1031 // perform the scat angle interpolation
1032 for (Index nst = 0; nst < npha; nst++)
1033 pha_mat_int[nst] = interp(
1034 dir_itw,
1035 ssd.pha_mat_data(find + f_start, 0, joker, 0, 0, 0, nst),
1036 dir_gp);
1037
1038 // convert from scat to lab frame
1039 pha_mat_labCalc(pha_mat_tmp(joker, joker),
1040 pha_mat_int,
1041 pdir_array(pdir, 0),
1042 pdir_array(pdir, 1),
1043 idir_array(idir, 0),
1044 idir_array(idir, 1),
1045 theta);
1046
1047 for (Index Tind = 0; Tind < nTout; Tind++)
1048 pha_mat(find, Tind, pdir, idir, joker, joker) = pha_mat_tmp;
1049 }
1050 }
1051 } else // T-interpolation required. To be done on the compact ssd format.
1052 {
1053 for (Index pdir = 0; pdir < npDir; pdir++)
1054 for (Index idir = 0; idir < niDir; idir++) {
1055 // calc scat ang theta from incident and prop dirs
1056 Numeric theta = scat_angle(pdir_array(pdir, 0),
1057 pdir_array(pdir, 1),
1058 idir_array(idir, 0),
1059 idir_array(idir, 1));
1060
1061 // get scat angle interpolation weights
1062 GridPos dir_gp;
1063 gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
1064 Vector dir_itw(2);
1065 interpweights(dir_itw, dir_gp);
1066
1067 Matrix pha_mat_int(nTin, npha, 0.);
1068 Matrix pha_mat_tmp(nTout, npha, 0.);
1069 for (Index find = 0; find < nf; find++) {
1070 for (Index Tind = 0; Tind < nTin; Tind++)
1071 // perform the scat angle interpolation
1072 for (Index nst = 0; nst < npha; nst++) {
1073 pha_mat_int(Tind, nst) = interp(
1074 dir_itw,
1075 ssd.pha_mat_data(find + f_start, Tind, joker, 0, 0, 0, nst),
1076 dir_gp);
1077 }
1078 // perform the T-interpolation
1079 for (Index nst = 0; nst < npha; nst++) {
1080 reinterp(pha_mat_tmp(joker, nst),
1081 pha_mat_int(joker, nst),
1082 T_itw_lag,
1083 T_lag);
1084 }
1085 // FIXME: it's probably better to do the frame conversion first,
1086 // then the T-interpolation (which is better depends on how many T-
1087 // aka vertical grid points we have. for single point, T-interpol
1088 // first is better, for a full vertical grid, frame conversion first
1089 // should be faster.)
1090 for (Index Tind = 0; Tind < nTout; Tind++) {
1091 // convert from scat to lab frame
1092 pha_mat_labCalc(pha_mat(find, Tind, pdir, idir, joker, joker),
1093 pha_mat_tmp(Tind, joker),
1094 pdir_array(pdir, 0),
1095 pdir_array(pdir, 1),
1096 idir_array(idir, 0),
1097 idir_array(idir, 1),
1098 theta);
1099 }
1100 }
1101 }
1102 }
1103 } else // dir-interpolation for non-tot-random particles
1104 // Data is already stored in the laboratory frame,
1105 // but it is compressed a little. Details elsewhere.
1106 {
1107 Index nDir = npDir * niDir;
1108 Vector adelta_aa(nDir);
1109 Matrix delta_aa(npDir, niDir);
1110 ArrayOfGridPos daa_gp(nDir), pza_gp(nDir), iza_gp(nDir);
1111 ArrayOfGridPos pza_gp_tmp(npDir), iza_gp_tmp(niDir);
1112
1113 gridpos(pza_gp_tmp, ssd.za_grid, pdir_array(joker, 0));
1114 gridpos(iza_gp_tmp, ssd.za_grid, idir_array(joker, 0));
1115
1116 Index j = 0;
1117 for (Index pdir = 0; pdir < npDir; pdir++) {
1118 for (Index idir = 0; idir < niDir; idir++) {
1119 delta_aa(pdir, idir) =
1120 pdir_array(pdir, 1) - idir_array(idir, 1) +
1121 (pdir_array(pdir, 1) - idir_array(idir, 1) < -180) * 360 -
1122 (pdir_array(pdir, 1) - idir_array(idir, 1) > 180) * 360;
1123 adelta_aa[j] = abs(delta_aa(pdir, idir));
1124 pza_gp[j] = pza_gp_tmp[pdir];
1125 iza_gp[j] = iza_gp_tmp[idir];
1126 j++;
1127 }
1128 }
1129
1130 gridpos(daa_gp, ssd.aa_grid, adelta_aa);
1131
1132 Matrix dir_itw(nDir, 8);
1133 interpweights(dir_itw, pza_gp, daa_gp, iza_gp);
1134
1135 if (this_T_interp_order < 0) // T only needs to be extracted.
1136 {
1137 Tensor3 pha_mat_int(nDir, stokes_dim, stokes_dim, 0.);
1138 Tensor4 pha_mat_tmp(npDir, niDir, stokes_dim, stokes_dim);
1139
1140 for (Index find = 0; find < nf; find++) {
1141 // perform the (tri-linear) angle interpolation. but only for the
1142 // pha_mat elements that we actually need.
1143
1144 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1145 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1146 interp(
1147 pha_mat_int(joker, ist1, ist2),
1148 dir_itw,
1149 ssd.pha_mat_data(
1150 find + f_start, 0, joker, joker, joker, 0, ist1 * 4 + ist2),
1151 pza_gp,
1152 daa_gp,
1153 iza_gp);
1154
1155 // sort direction-combined 1D-array back into prop and incident
1156 // direction matrix
1157 Index i = 0;
1158 for (Index pdir = 0; pdir < npDir; pdir++)
1159 for (Index idir = 0; idir < niDir; idir++) {
1160 pha_mat_tmp(pdir, idir, joker, joker) =
1161 pha_mat_int(i, joker, joker);
1162 i++;
1163 }
1164
1165 if (stokes_dim > 2) {
1166 for (Index pdir = 0; pdir < npDir; pdir++)
1167 for (Index idir = 0; idir < niDir; idir++)
1168 if (delta_aa(pdir, idir) < 0.) {
1169 pha_mat_tmp(pdir, idir, 0, 2) *= -1;
1170 pha_mat_tmp(pdir, idir, 1, 2) *= -1;
1171 pha_mat_tmp(pdir, idir, 2, 0) *= -1;
1172 pha_mat_tmp(pdir, idir, 2, 1) *= -1;
1173 }
1174 }
1175
1176 if (stokes_dim > 3) {
1177 for (Index pdir = 0; pdir < npDir; pdir++)
1178 for (Index idir = 0; idir < niDir; idir++)
1179 if (delta_aa(pdir, idir) < 0.) {
1180 pha_mat_tmp(pdir, idir, 0, 3) *= -1;
1181 pha_mat_tmp(pdir, idir, 1, 3) *= -1;
1182 pha_mat_tmp(pdir, idir, 3, 0) *= -1;
1183 pha_mat_tmp(pdir, idir, 3, 1) *= -1;
1184 }
1185 }
1186
1187 for (Index Tind = 0; Tind < nTout; Tind++)
1188 pha_mat(find, Tind, joker, joker, joker, joker) = pha_mat_tmp;
1189 }
1190 }
1191
1192 else // T- and dir-interpolation required. To be done on the compact ssd
1193 // format.
1194 {
1195 Tensor4 pha_mat_int(nTin, nDir, stokes_dim, stokes_dim, 0.);
1196
1197 for (Index find = 0; find < nf; find++) {
1198 // perform the (tri-linear) angle interpolation. but only for the
1199 // pha_mat elements that we actually need.
1200 for (Index Tind = 0; Tind < nTin; Tind++) {
1201 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1202 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1203 interp(pha_mat_int(Tind, joker, ist1, ist2),
1204 dir_itw,
1205 ssd.pha_mat_data(find + f_start,
1206 Tind,
1207 joker,
1208 joker,
1209 joker,
1210 0,
1211 ist1 * 4 + ist2),
1212 pza_gp,
1213 daa_gp,
1214 iza_gp);
1215 }
1216
1217 // perform the T-interpolation and simultaneously sort
1218 // direction-combined 1D-array back into prop and incident direction
1219 // matrix.
1220 Index i = 0;
1221 for (Index pdir = 0; pdir < npDir; pdir++)
1222 for (Index idir = 0; idir < niDir; idir++) {
1223 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1224 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1225 reinterp(pha_mat(find, joker, pdir, idir, ist1, ist2),
1226 pha_mat_int(joker, i, ist1, ist2),
1227 T_itw_lag,
1228 T_lag);
1229 i++;
1230 }
1231
1232 if (stokes_dim > 2) {
1233 for (Index pdir = 0; pdir < npDir; pdir++)
1234 for (Index idir = 0; idir < niDir; idir++)
1235 if (delta_aa(pdir, idir) < 0.) {
1236 pha_mat(find, joker, pdir, idir, 0, 2) *= -1;
1237 pha_mat(find, joker, pdir, idir, 1, 2) *= -1;
1238 pha_mat(find, joker, pdir, idir, 2, 0) *= -1;
1239 pha_mat(find, joker, pdir, idir, 2, 1) *= -1;
1240 }
1241 }
1242
1243 if (stokes_dim > 3) {
1244 for (Index pdir = 0; pdir < npDir; pdir++)
1245 for (Index idir = 0; idir < niDir; idir++)
1246 if (delta_aa(pdir, idir) < 0.) {
1247 pha_mat(find, joker, pdir, idir, 0, 3) *= -1;
1248 pha_mat(find, joker, pdir, idir, 1, 3) *= -1;
1249 pha_mat(find, joker, pdir, idir, 3, 0) *= -1;
1250 pha_mat(find, joker, pdir, idir, 3, 1) *= -1;
1251 }
1252 }
1253 }
1254 }
1255 }
1256}
1257
1259
1280void abs_vecTransform( //Output and Input
1281 StokesVector& abs_vec_lab,
1282 //Input
1283 ConstTensor3View abs_vec_data,
1284 ConstVectorView za_datagrid,
1285 ConstVectorView aa_datagrid _U_,
1286 const PType& ptype,
1287 const Numeric& za_sca _U_,
1288 const Numeric& aa_sca _U_,
1289 const Verbosity& verbosity) {
1290 const Index stokes_dim = abs_vec_lab.StokesDimensions();
1291 ARTS_ASSERT(abs_vec_lab.NumberOfFrequencies() == 1);
1292
1293 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1294 "The dimension of the stokes vector \n"
1295 "must be 1,2,3 or 4");
1296
1297 switch (ptype) {
1298 case PTYPE_GENERAL: {
1299 /*
1300 TO ANY DEVELOPER:
1301 current usage of coordinate systems in scattering solvers (RT and SSD
1302 extraction) and general radiative transfer is not consistent. Not an
1303 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1304 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1305 adding PTYPE_GENERAL support (see AUG appendix for more info).
1306 */
1307
1309 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1310 break;
1311 }
1312 case PTYPE_TOTAL_RND: {
1313 // The first element of the vector corresponds to the absorption
1314 // coefficient which is stored in the database, the others are 0.
1315
1316 abs_vec_lab.SetZero();
1317
1318 abs_vec_lab.Kjj()[0] = abs_vec_data(0, 0, 0);
1319 break;
1320 }
1321
1322 case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1323 {
1324 ARTS_ASSERT(abs_vec_data.ncols() == 2);
1325
1326 // In the case of azimuthally randomly oriented particles, only the first
1327 // two elements of the absorption coefficient vector are non-zero.
1328 // These values are dependent on the zenith angle of propagation.
1329
1330 // 1st interpolate data by za_sca
1331 GridPos gp;
1332 Vector itw(2);
1333
1334 gridpos(gp, za_datagrid, za_sca);
1335 interpweights(itw, gp);
1336
1337 abs_vec_lab.SetZero();
1338
1339 abs_vec_lab.Kjj()[0] = interp(itw, abs_vec_data(Range(joker), 0, 0), gp);
1340
1341 if (stokes_dim == 1) {
1342 break;
1343 }
1344 abs_vec_lab.K12()[0] = interp(itw, abs_vec_data(Range(joker), 0, 1), gp);
1345 break;
1346 }
1347 default: {
1349 out0 << "Not all ptype cases are implemented\n";
1350 }
1351 }
1352}
1353
1355
1376void ext_matTransform( //Output and Input
1377 PropagationMatrix& ext_mat_lab,
1378 //Input
1379 ConstTensor3View ext_mat_data,
1380 ConstVectorView za_datagrid,
1381 ConstVectorView aa_datagrid _U_,
1382 const PType& ptype,
1383 const Numeric& za_sca,
1384 const Numeric& aa_sca _U_,
1385 const Verbosity& verbosity) {
1386 const Index stokes_dim = ext_mat_lab.StokesDimensions();
1387 ARTS_ASSERT(ext_mat_lab.NumberOfFrequencies() == 1);
1388
1389 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1390 "The dimension of the stokes vector \n"
1391 "must be 1,2,3 or 4");
1392
1393 switch (ptype) {
1394 case PTYPE_GENERAL: {
1395 /*
1396 TO ANY DEVELOPER:
1397 current usage of coordinate systems in scattering solvers (RT and SSD
1398 extraction) and general radiative transfer is not consistent. Not an
1399 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1400 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1401 adding PTYPE_GENERAL support (see AUG appendix for more info).
1402 */
1403
1405 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1406 break;
1407 }
1408 case PTYPE_TOTAL_RND: {
1409 ARTS_ASSERT(ext_mat_data.ncols() == 1);
1410
1411 // In the case of randomly oriented particles the extinction matrix is
1412 // diagonal. The value of each element of the diagonal is the
1413 // extinction cross section, which is stored in the database.
1414
1415 ext_mat_lab.SetZero();
1416
1417 ext_mat_lab.Kjj() = ext_mat_data(0, 0, 0);
1418 break;
1419 }
1420
1421 case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1422 {
1423 ARTS_ASSERT(ext_mat_data.ncols() == 3);
1424
1425 // In the case of azimuthally randomly oriented particles, the extinction
1426 // matrix has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
1427 // These values are dependent on the zenith angle of propagation.
1428
1429 // 1st interpolate data by za_sca
1430 GridPos gp;
1431 Vector itw(2);
1432 Numeric Kjj;
1433 Numeric K12;
1434 Numeric K34;
1435
1436 gridpos(gp, za_datagrid, za_sca);
1437 interpweights(itw, gp);
1438
1439 ext_mat_lab.SetZero();
1440
1441 Kjj = interp(itw, ext_mat_data(Range(joker), 0, 0), gp);
1442 ext_mat_lab.Kjj() = Kjj;
1443
1444 if (stokes_dim < 2) {
1445 break;
1446 }
1447
1448 K12 = interp(itw, ext_mat_data(Range(joker), 0, 1), gp);
1449 ext_mat_lab.K12()[0] = K12;
1450
1451 if (stokes_dim < 4) {
1452 break;
1453 }
1454
1455 K34 = interp(itw, ext_mat_data(Range(joker), 0, 2), gp);
1456 ext_mat_lab.K34()[0] = K34;
1457 break;
1458 }
1459 default: {
1461 out0 << "Not all ptype cases are implemented\n";
1462 }
1463 }
1464}
1465
1467
1494void pha_matTransform( //Output
1495 MatrixView pha_mat_lab,
1496 //Input
1497 ConstTensor5View pha_mat_data,
1498 ConstVectorView za_datagrid,
1499 ConstVectorView aa_datagrid,
1500 const PType& ptype,
1501 const Index& za_sca_idx,
1502 const Index& aa_sca_idx,
1503 const Index& za_inc_idx,
1504 const Index& aa_inc_idx,
1505 ConstVectorView za_grid,
1506 ConstVectorView aa_grid,
1507 const Verbosity& verbosity) {
1508 const Index stokes_dim = pha_mat_lab.ncols();
1509
1510 Numeric za_sca = za_grid[za_sca_idx];
1511 Numeric aa_sca = aa_grid[aa_sca_idx];
1512 Numeric za_inc = za_grid[za_inc_idx];
1513 Numeric aa_inc = aa_grid[aa_inc_idx];
1514
1515 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1516 "The dimension of the stokes vector \n"
1517 "must be 1,2,3 or 4");
1518
1519 switch (ptype) {
1520 case PTYPE_GENERAL: {
1521 /*
1522 TO ANY DEVELOPER:
1523 current usage of coordinate systems in scattering solvers (RT and SSD
1524 extraction) and general radiative transfer is not consistent. Not an
1525 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1526 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1527 adding PTYPE_GENERAL support (see AUG appendix for more info).
1528 */
1529
1531 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1532 break;
1533 }
1534 case PTYPE_TOTAL_RND: {
1535 // Calculate the scattering angle and interpolate the data in it:
1536 Numeric theta_rad = scat_angle(za_sca, aa_sca, za_inc, aa_inc);
1537 const Numeric theta = RAD2DEG * theta_rad;
1538
1539 // Interpolation of the data on the scattering angle:
1540 Vector pha_mat_int(6);
1541 interpolate_scat_angle(pha_mat_int, pha_mat_data, za_datagrid, theta);
1542
1543 // Calculate the phase matrix in the laboratory frame:
1545 pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1546
1547 break;
1548 }
1549
1550 case PTYPE_AZIMUTH_RND: //Added by Cory Davis
1551 //Data is already stored in the laboratory frame,
1552 //but it is compressed a little. Details elsewhere.
1553 {
1554 ARTS_ASSERT(pha_mat_data.ncols() == 16);
1555 ARTS_ASSERT(pha_mat_data.npages() == za_datagrid.nelem());
1556 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
1557 (aa_sca - aa_inc > 180) *
1558 360; //delta_aa corresponds to the "books"
1559 //dimension of pha_mat_data
1560 GridPos za_sca_gp;
1561 GridPos delta_aa_gp;
1562 GridPos za_inc_gp;
1563 Vector itw(8);
1564
1565 gridpos(delta_aa_gp, aa_datagrid, abs(delta_aa));
1566 gridpos(za_inc_gp, za_datagrid, za_inc);
1567 gridpos(za_sca_gp, za_datagrid, za_sca);
1568
1569 interpweights(itw, za_sca_gp, delta_aa_gp, za_inc_gp);
1570
1571 pha_mat_lab(0, 0) =
1572 interp(itw,
1573 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 0),
1574 za_sca_gp,
1575 delta_aa_gp,
1576 za_inc_gp);
1577 if (stokes_dim == 1) {
1578 break;
1579 }
1580 pha_mat_lab(0, 1) =
1581 interp(itw,
1582 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 1),
1583 za_sca_gp,
1584 delta_aa_gp,
1585 za_inc_gp);
1586 pha_mat_lab(1, 0) =
1587 interp(itw,
1588 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 4),
1589 za_sca_gp,
1590 delta_aa_gp,
1591 za_inc_gp);
1592 pha_mat_lab(1, 1) =
1593 interp(itw,
1594 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 5),
1595 za_sca_gp,
1596 delta_aa_gp,
1597 za_inc_gp);
1598 if (stokes_dim == 2) {
1599 break;
1600 }
1601 if (delta_aa >= 0) {
1602 pha_mat_lab(0, 2) = interp(
1603 itw,
1604 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1605 za_sca_gp,
1606 delta_aa_gp,
1607 za_inc_gp);
1608 pha_mat_lab(1, 2) = interp(
1609 itw,
1610 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1611 za_sca_gp,
1612 delta_aa_gp,
1613 za_inc_gp);
1614 pha_mat_lab(2, 0) = interp(
1615 itw,
1616 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1617 za_sca_gp,
1618 delta_aa_gp,
1619 za_inc_gp);
1620 pha_mat_lab(2, 1) = interp(
1621 itw,
1622 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1623 za_sca_gp,
1624 delta_aa_gp,
1625 za_inc_gp);
1626 } else {
1627 pha_mat_lab(0, 2) = -interp(
1628 itw,
1629 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1630 za_sca_gp,
1631 delta_aa_gp,
1632 za_inc_gp);
1633 pha_mat_lab(1, 2) = -interp(
1634 itw,
1635 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1636 za_sca_gp,
1637 delta_aa_gp,
1638 za_inc_gp);
1639 pha_mat_lab(2, 0) = -interp(
1640 itw,
1641 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1642 za_sca_gp,
1643 delta_aa_gp,
1644 za_inc_gp);
1645 pha_mat_lab(2, 1) = -interp(
1646 itw,
1647 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1648 za_sca_gp,
1649 delta_aa_gp,
1650 za_inc_gp);
1651 }
1652 pha_mat_lab(2, 2) = interp(
1653 itw,
1654 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 10),
1655 za_sca_gp,
1656 delta_aa_gp,
1657 za_inc_gp);
1658 if (stokes_dim == 3) {
1659 break;
1660 }
1661 if (delta_aa >= 0) {
1662 pha_mat_lab(0, 3) = interp(
1663 itw,
1664 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1665 za_sca_gp,
1666 delta_aa_gp,
1667 za_inc_gp);
1668 pha_mat_lab(1, 3) = interp(
1669 itw,
1670 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1671 za_sca_gp,
1672 delta_aa_gp,
1673 za_inc_gp);
1674 pha_mat_lab(3, 0) = interp(
1675 itw,
1676 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1677 za_sca_gp,
1678 delta_aa_gp,
1679 za_inc_gp);
1680 pha_mat_lab(3, 1) = interp(
1681 itw,
1682 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1683 za_sca_gp,
1684 delta_aa_gp,
1685 za_inc_gp);
1686 } else {
1687 pha_mat_lab(0, 3) = -interp(
1688 itw,
1689 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1690 za_sca_gp,
1691 delta_aa_gp,
1692 za_inc_gp);
1693 pha_mat_lab(1, 3) = -interp(
1694 itw,
1695 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1696 za_sca_gp,
1697 delta_aa_gp,
1698 za_inc_gp);
1699 pha_mat_lab(3, 0) = -interp(
1700 itw,
1701 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1702 za_sca_gp,
1703 delta_aa_gp,
1704 za_inc_gp);
1705 pha_mat_lab(3, 1) = -interp(
1706 itw,
1707 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1708 za_sca_gp,
1709 delta_aa_gp,
1710 za_inc_gp);
1711 }
1712 pha_mat_lab(2, 3) = interp(
1713 itw,
1714 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 11),
1715 za_sca_gp,
1716 delta_aa_gp,
1717 za_inc_gp);
1718 pha_mat_lab(3, 2) = interp(
1719 itw,
1720 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 14),
1721 za_sca_gp,
1722 delta_aa_gp,
1723 za_inc_gp);
1724 pha_mat_lab(3, 3) = interp(
1725 itw,
1726 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 15),
1727 za_sca_gp,
1728 delta_aa_gp,
1729 za_inc_gp);
1730 break;
1731 }
1732
1733 default: {
1735 out0 << "Not all ptype cases are implemented\n";
1736 }
1737 }
1738}
1739
1741
1770 MatrixView ext_mat,
1771 //Input
1772 ConstVectorView abs_vec,
1773 const Index& stokes_dim) {
1774 ARTS_ASSERT(stokes_dim >= 1 && stokes_dim <= 4);
1775 ARTS_ASSERT(ext_mat.nrows() == stokes_dim);
1776 ARTS_ASSERT(ext_mat.ncols() == stokes_dim);
1777 ARTS_ASSERT(abs_vec.nelem() == stokes_dim);
1778
1779 // first: diagonal elements
1780 for (Index is = 0; is < stokes_dim; is++) {
1781 ext_mat(is, is) += abs_vec[0];
1782 }
1783 // second: off-diagonal elements, namely first row and column
1784 for (Index is = 1; is < stokes_dim; is++) {
1785 ext_mat(0, is) += abs_vec[is];
1786 ext_mat(is, 0) += abs_vec[is];
1787 }
1788}
1789
1791
1804Numeric scat_angle(const Numeric& za_sca,
1805 const Numeric& aa_sca,
1806 const Numeric& za_inc,
1807 const Numeric& aa_inc) {
1808 Numeric theta_rad;
1809 Numeric ANG_TOL = 1e-7;
1810
1811 // CPD 5/10/03.
1812 // Two special cases are implemented here to avoid NaNs that can sometimes
1813 // occur in in the acos() formula in forward and backscattering cases.
1814 //
1815 // GH 2011-05-31
1816 // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
1817
1818 if ((abs(aa_sca - aa_inc) < ANG_TOL) ||
1819 (abs(abs(aa_sca - aa_inc) - 360) < ANG_TOL)) {
1820 theta_rad = Conversion::deg2rad(abs(za_sca - za_inc));
1821 } else if (abs(abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
1822 theta_rad = Conversion::deg2rad(za_sca + za_inc);
1823 if (theta_rad > PI) {
1824 theta_rad = 2 * PI - theta_rad;
1825 }
1826 } else {theta_rad =
1827 acos(Conversion::cosd(za_sca) * Conversion::cosd(za_inc) +
1828 Conversion::sind(za_sca) * Conversion::sind(za_inc) *
1829 Conversion::cosd(aa_sca - aa_inc));
1830 }
1831 return theta_rad;
1832}
1833
1835
1859 VectorView pha_mat_int,
1860 //Input:
1861 ConstTensor5View pha_mat_data,
1862 ConstVectorView za_datagrid,
1863 const Numeric theta) {
1864 GridPos thet_gp;
1865 gridpos(thet_gp, za_datagrid, theta);
1866
1867 Vector itw(2);
1868 interpweights(itw, thet_gp);
1869
1870 ARTS_ASSERT(pha_mat_data.ncols() == 6);
1871 for (Index i = 0; i < 6; i++) {
1872 pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), thet_gp);
1873 }
1874}
1875
1877
1902void pha_mat_labCalc( //Output:
1903 MatrixView pha_mat_lab,
1904 //Input:
1905 ConstVectorView pha_mat_int,
1906 const Numeric& za_sca,
1907 const Numeric& aa_sca,
1908 const Numeric& za_inc,
1909 const Numeric& aa_inc,
1910 const Numeric& theta_rad) {
1911 const Index stokes_dim = pha_mat_lab.ncols();
1912
1913 ARTS_USER_ERROR_IF (std::isnan(F11),
1914 "NaN value(s) detected in *pha_mat_labCalc* (0,0). Could the "
1915 "input data contain NaNs? Please check with *scat_dataCheck*. If "
1916 "input data are OK and you critically need the ongoing calculations, "
1917 "try to change the observation LOS slightly. If you can reproduce "
1918 "this error, please contact Patrick in order to help tracking down "
1919 "the reason to this problem. If you see this message occasionally "
1920 "when doing MC calculations, it should not be critical. This path "
1921 "sampling will be rejected and replaced with a new one.");
1922
1923 // For stokes_dim = 1, we only need Z11=F11:
1924 pha_mat_lab(0, 0) = F11;
1925
1926 if (stokes_dim > 1) {
1927 Numeric za_sca_rad = Conversion::deg2rad(za_sca);
1928 Numeric za_inc_rad = Conversion::deg2rad(za_inc);
1929 Numeric aa_sca_rad = Conversion::deg2rad(aa_sca);
1930 Numeric aa_inc_rad = Conversion::deg2rad(aa_inc);
1931
1932 const Numeric ANGTOL_RAD = 1e-6; //CPD: this constant is used to adjust
1933 //zenith angles close to 0 and PI. This is
1934 //also used to avoid float == float statements.
1935
1936 //
1937 // Several cases have to be considered:
1938 //
1939
1940 if ((abs(theta_rad) < ANGTOL_RAD) // forward scattering
1941 || (abs(theta_rad - Constant::pi) < ANGTOL_RAD) // backward scattering
1942 ||
1943 (abs(aa_inc_rad - aa_sca_rad) < ANGTOL_RAD) // inc and sca on meridian
1944 || (abs(abs(aa_inc_rad - aa_sca_rad) - Constant::two_pi) < ANGTOL_RAD) // "
1945 || (abs(abs(aa_inc_rad - aa_sca_rad) - Constant::pi) < ANGTOL_RAD) // "
1946 ) {
1947 pha_mat_lab(0, 1) = F12;
1948 pha_mat_lab(1, 0) = F12;
1949 pha_mat_lab(1, 1) = F22;
1950
1951 if (stokes_dim > 2) {
1952 pha_mat_lab(0, 2) = 0;
1953 pha_mat_lab(1, 2) = 0;
1954 pha_mat_lab(2, 0) = 0;
1955 pha_mat_lab(2, 1) = 0;
1956 pha_mat_lab(2, 2) = F33;
1957
1958 if (stokes_dim > 3) {
1959 pha_mat_lab(0, 3) = 0;
1960 pha_mat_lab(1, 3) = 0;
1961 pha_mat_lab(2, 3) = F34;
1962 pha_mat_lab(3, 0) = 0;
1963 pha_mat_lab(3, 1) = 0;
1964 pha_mat_lab(3, 2) = -F34;
1965 pha_mat_lab(3, 3) = F44;
1966 }
1967 }
1968 }
1969
1970 else {
1971 Numeric sigma1;
1972 Numeric sigma2;
1973
1974 Numeric s1, s2;
1975
1976 // In these cases we have to take limiting values.
1977
1978 if (za_inc_rad < ANGTOL_RAD) {
1979 sigma1 = pi + aa_sca_rad - aa_inc_rad;
1980 sigma2 = 0;
1981 } else if (za_inc_rad > pi - ANGTOL_RAD) {
1982 sigma1 = aa_sca_rad - aa_inc_rad;
1983 sigma2 = pi;
1984 } else if (za_sca_rad < ANGTOL_RAD) {
1985 sigma1 = 0;
1986 sigma2 = pi + aa_sca_rad - aa_inc_rad;
1987 } else if (za_sca_rad > pi - ANGTOL_RAD) {
1988 sigma1 = pi;
1989 sigma2 = aa_sca_rad - aa_inc_rad;
1990 } else {
1991 s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad)) /
1992 (sin(za_inc_rad) * sin(theta_rad));
1993 s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos(theta_rad)) /
1994 (sin(za_sca_rad) * sin(theta_rad));
1995
1996 sigma1 = acos(s1);
1997 sigma2 = acos(s2);
1998
1999 // Arccos is only defined in the range from -1 ... 1
2000 // Numerical problems can appear for values close to 1 or -1
2001 // this (also) catches the case when inc and sca are on one meridian
2002 if (std::isnan(sigma1) || std::isnan(sigma2)) {
2003 if (abs(s1 - 1) < ANGTOL_RAD) sigma1 = 0;
2004 if (abs(s1 + 1) < ANGTOL_RAD) sigma1 = pi;
2005 if (abs(s2 - 1) < ANGTOL_RAD) sigma2 = 0;
2006 if (abs(s2 + 1) < ANGTOL_RAD) sigma2 = pi;
2007 }
2008 }
2009
2010 const Numeric C1 = cos(2 * sigma1);
2011 const Numeric C2 = cos(2 * sigma2);
2012
2013 const Numeric S1 = sin(2 * sigma1);
2014 const Numeric S2 = sin(2 * sigma2);
2015
2016 pha_mat_lab(0, 1) = C1 * F12;
2017 pha_mat_lab(1, 0) = C2 * F12;
2018 pha_mat_lab(1, 1) = C1 * C2 * F22 - S1 * S2 * F33;
2019
2020 //ARTS_ASSERT(!std::isnan(pha_mat_lab(0,1)));
2021 //ARTS_ASSERT(!std::isnan(pha_mat_lab(1,0)));
2022 //ARTS_ASSERT(!std::isnan(pha_mat_lab(1,1)));
2023 ARTS_USER_ERROR_IF (std::isnan(pha_mat_lab(0, 1)) || std::isnan(pha_mat_lab(1, 0)) ||
2024 std::isnan(pha_mat_lab(1, 1)),
2025 "NaN value(s) detected in *pha_mat_labCalc* (0/1,1). Could the "
2026 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2027 "input data are OK and you critically need the ongoing calculations, "
2028 "try to change the observation LOS slightly. If you can reproduce "
2029 "this error, please contact Patrick in order to help tracking down "
2030 "the reason to this problem. If you see this message occasionally "
2031 "when doing MC calculations, it should not be critical. This path "
2032 "sampling will be rejected and replaced with a new one.");
2033
2034 if (stokes_dim > 2) {
2035 /*CPD: For skokes_dim > 2 some of the transformation formula
2036 for each element have a different sign depending on whether or
2037 not 0<aa_scat-aa_inc<180. For details see pages 94 and 95 of
2038 Mishchenkos chapter in :
2039 Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic
2040 scattering by nonspherical particles. In Exploring the Atmosphere
2041 by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag,
2042 Berlin, pp. 77-127.
2043 This is available at http://www.giss.nasa.gov/~crmim/publications/ */
2044 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
2045 (aa_sca - aa_inc > 180) * 360;
2046 if (delta_aa >= 0) {
2047 pha_mat_lab(0, 2) = S1 * F12;
2048 pha_mat_lab(1, 2) = S1 * C2 * F22 + C1 * S2 * F33;
2049 pha_mat_lab(2, 0) = -S2 * F12;
2050 pha_mat_lab(2, 1) = -C1 * S2 * F22 - S1 * C2 * F33;
2051 } else {
2052 pha_mat_lab(0, 2) = -S1 * F12;
2053 pha_mat_lab(1, 2) = -S1 * C2 * F22 - C1 * S2 * F33;
2054 pha_mat_lab(2, 0) = S2 * F12;
2055 pha_mat_lab(2, 1) = C1 * S2 * F22 + S1 * C2 * F33;
2056 }
2057 pha_mat_lab(2, 2) = -S1 * S2 * F22 + C1 * C2 * F33;
2058
2059 if (stokes_dim > 3) {
2060 if (delta_aa >= 0) {
2061 pha_mat_lab(1, 3) = S2 * F34;
2062 pha_mat_lab(3, 1) = S1 * F34;
2063 } else {
2064 pha_mat_lab(1, 3) = -S2 * F34;
2065 pha_mat_lab(3, 1) = -S1 * F34;
2066 }
2067 pha_mat_lab(0, 3) = 0;
2068 pha_mat_lab(2, 3) = C2 * F34;
2069 pha_mat_lab(3, 0) = 0;
2070 pha_mat_lab(3, 2) = -C1 * F34;
2071 pha_mat_lab(3, 3) = F44;
2072 }
2073 }
2074 }
2075 }
2076}
2077
2078ostream& operator<<(ostream& os, const SingleScatteringData& /*ssd*/) {
2079 os << "SingleScatteringData: Output operator not implemented";
2080 return os;
2081}
2082
2083ostream& operator<<(ostream& os, const ScatteringMetaData& /*ssd*/) {
2084 os << "ScatteringMetaData: Output operator not implemented";
2085 return os;
2086}
2087
2089
2105 PropagationMatrix& ext_mat,
2106 StokesVector& abs_vec,
2107 //Input:
2108 const PropagationMatrix& propmat_clearsky) {
2109 const Index stokes_dim = propmat_clearsky.StokesDimensions();
2110 const Index freq_dim = propmat_clearsky.NumberOfFrequencies();
2111
2112 // old abs_vecInit
2113 abs_vec = StokesVector(freq_dim, stokes_dim);
2114 abs_vec.SetZero();
2115
2116 ext_mat = PropagationMatrix(freq_dim, stokes_dim);
2117 ext_mat.SetZero();
2118
2119 // old ext_matAddGas and abs_vecAddGas for 0 vector and matrix
2120 abs_vec += propmat_clearsky;
2121 ext_mat += propmat_clearsky;
2122}
2123
2125
2135PType PTypeFromString(const String& ptype_string) {
2136 PType ptype;
2137 if (ptype_string == "general")
2138 ptype = PTYPE_GENERAL;
2139 else if (ptype_string == "totally_random")
2140 ptype = PTYPE_TOTAL_RND;
2141 else if (ptype_string == "azimuthally_random")
2142 ptype = PTYPE_AZIMUTH_RND;
2143 else {
2145 "Unknown ptype: ", ptype_string, "\n"
2146 "Valid types are: general, totally_random and "
2147 "azimuthally_random.")
2148 }
2149
2150 return ptype;
2151}
2152
2154
2164PType PType2FromString(const String& ptype_string) {
2165 PType ptype;
2166 if (ptype_string == "general")
2167 ptype = PTYPE_GENERAL;
2168 else if (ptype_string == "macroscopically_isotropic")
2169 ptype = PTYPE_TOTAL_RND;
2170 else if (ptype_string == "horizontally_aligned")
2171 ptype = PTYPE_AZIMUTH_RND;
2172 else {
2174 "Unknown ptype: ", ptype_string, "\n"
2175 "Valid types are: general, macroscopically_isotropic and "
2176 "horizontally_aligned.")
2177 }
2178
2179 return ptype;
2180}
2181
2183
2192 String ptype_string;
2193
2194 switch (ptype) {
2195 case PTYPE_GENERAL:
2196 ptype_string = "general";
2197 break;
2198 case PTYPE_TOTAL_RND:
2199 ptype_string = "totally_random";
2200 break;
2201 case PTYPE_AZIMUTH_RND:
2202 ptype_string = "azimuthally_random";
2203 break;
2204 default:
2206 "Internal error: Cannot map PType enum value ",
2207 ptype, " to String.")
2208 break;
2209 }
2210
2211 return ptype_string;
2212}
2213
2215
2223 // First check that input fulfills requirements on older data formats:
2224 // 1) Is za_grid symmetric and includes 90deg?
2225 Index nza = ssd.za_grid.nelem();
2226 for (Index i = 0; i < nza / 2; i++) {
2227 ARTS_USER_ERROR_IF (!is_same_within_epsilon(
2228 180. - ssd.za_grid[nza - 1 - i], ssd.za_grid[i], 2 * DBL_EPSILON),
2229 "Zenith grid of azimuthally_random single scattering data\n"
2230 "is not symmetric with respect to 90degree.")
2231 }
2232 ARTS_USER_ERROR_IF (!is_same_within_epsilon(ssd.za_grid[nza / 2], 90., 2 * DBL_EPSILON),
2233 "Zenith grid of azimuthally_random single scattering data\n"
2234 "does not contain 90 degree grid point.")
2235
2236 // 2) Are data sizes correct?
2237 ostringstream os_pha_mat;
2238 os_pha_mat << "pha_mat ";
2239 ostringstream os_ext_mat;
2240 os_ext_mat << "ext_mat ";
2241 ostringstream os_abs_vec;
2242 os_abs_vec << "abs_vec ";
2243 chk_size(os_pha_mat.str(),
2244 ssd.pha_mat_data,
2245 ssd.f_grid.nelem(),
2246 ssd.T_grid.nelem(),
2247 ssd.za_grid.nelem(),
2248 ssd.aa_grid.nelem(),
2249 ssd.za_grid.nelem() / 2 + 1,
2250 1,
2251 16);
2252
2253 chk_size(os_ext_mat.str(),
2254 ssd.ext_mat_data,
2255 ssd.f_grid.nelem(),
2256 ssd.T_grid.nelem(),
2257 ssd.za_grid.nelem() / 2 + 1,
2258 1,
2259 3);
2260
2261 chk_size(os_abs_vec.str(),
2262 ssd.abs_vec_data,
2263 ssd.f_grid.nelem(),
2264 ssd.T_grid.nelem(),
2265 ssd.za_grid.nelem() / 2 + 1,
2266 1,
2267 2);
2268
2269 // Now that we are sure that za_grid is properly symmetric, we just need to
2270 // copy over the data (ie no interpolation).
2271 Tensor5 tmpT5 = ssd.abs_vec_data;
2272 ssd.abs_vec_data.resize(tmpT5.nshelves(),
2273 tmpT5.nbooks(),
2274 ssd.za_grid.nelem(),
2275 tmpT5.nrows(),
2276 tmpT5.ncols());
2277 ssd.abs_vec_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2278 for (Index i = 0; i < nza / 2; i++) {
2279 ssd.abs_vec_data(joker, joker, nza - 1 - i, joker, joker) =
2280 tmpT5(joker, joker, i, joker, joker);
2281 }
2282
2283 tmpT5 = ssd.ext_mat_data;
2284 ssd.ext_mat_data.resize(tmpT5.nshelves(),
2285 tmpT5.nbooks(),
2286 ssd.za_grid.nelem(),
2287 tmpT5.nrows(),
2288 tmpT5.ncols());
2289 ssd.ext_mat_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2290 for (Index i = 0; i < nza / 2; i++) {
2291 ssd.ext_mat_data(joker, joker, nza - 1 - i, joker, joker) =
2292 tmpT5(joker, joker, i, joker, joker);
2293 }
2294
2295 Tensor7 tmpT7 = ssd.pha_mat_data;
2296 ssd.pha_mat_data.resize(tmpT7.nlibraries(),
2297 tmpT7.nvitrines(),
2298 tmpT7.nshelves(),
2299 tmpT7.nbooks(),
2300 ssd.za_grid.nelem(),
2301 tmpT7.nrows(),
2302 tmpT7.ncols());
2303 ssd.pha_mat_data(
2304 joker, joker, joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT7;
2305
2306 // scatt. matrix elements 13,23,31,32 and 14,24,41,42 (=elements 2,6,8,9 and
2307 // 3,7,12,13 in ARTS' flattened format, respectively) change sign.
2308 tmpT7(joker, joker, joker, joker, joker, joker, Range(2, 2)) *= -1.;
2309 tmpT7(joker, joker, joker, joker, joker, joker, Range(6, 4)) *= -1.;
2310 tmpT7(joker, joker, joker, joker, joker, joker, Range(12, 2)) *= -1.;
2311
2312 // For second half of incident polar angles (>90deg), we need to mirror the
2313 // original data in both incident and scattered polar angle around 90deg "planes".
2314 for (Index i = 0; i < nza / 2; i++)
2315 for (Index j = 0; j < nza; j++)
2316 ssd.pha_mat_data(
2317 joker, joker, nza - 1 - j, joker, nza - 1 - i, joker, joker) =
2318 tmpT7(joker, joker, j, joker, i, joker, joker);
2319}
2320
2322
2331 const String& particle_ssdmethod_string) {
2332 ParticleSSDMethod particle_ssdmethod;
2333 if (particle_ssdmethod_string == "tmatrix")
2334 particle_ssdmethod = PARTICLE_SSDMETHOD_TMATRIX;
2335 else {
2337 "Unknown particle SSD method: ",
2338 particle_ssdmethod_string, "\n"
2339 "Valid methods: tmatrix")
2340 }
2341
2342 return particle_ssdmethod;
2343}
2344
2346
2354String PTypeToString(const ParticleSSDMethod& particle_ssdmethod) {
2355 String particle_ssdmethod_string;
2356
2357 switch (particle_ssdmethod) {
2359 particle_ssdmethod_string = "tmatrix";
2360 break;
2361 default:
2363 "Internal error: Cannot map ParticleSSDMethod enum value ",
2364 particle_ssdmethod, " to String.")
2365 break;
2366 }
2367
2368 return particle_ssdmethod_string;
2369}
2370
2371
2373
2400void ext_abs_pfun_from_tro(MatrixView ext_data,
2401 MatrixView abs_data,
2402 Tensor3View pfun_data,
2403 const ArrayOfSingleScatteringData& scat_data,
2404 const Index& iss,
2405 ConstMatrixView pnd_data,
2406 ArrayOfIndex& cloudbox_limits,
2407 ConstVectorView T_grid,
2408 ConstVectorView sa_grid)
2409{
2410 // Sizes
2411 const Index nse = scat_data.nelem();
2412 const Index nf = scat_data[0].f_grid.nelem();
2413 [[maybe_unused]] const Index nt = T_grid.nelem();
2414 const Index nsa = sa_grid.nelem();
2415 const Index ncl = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2416
2417 ARTS_ASSERT(ext_data.nrows() == nf);
2418 ARTS_ASSERT(ext_data.ncols() == nt);
2419 ARTS_ASSERT(abs_data.nrows() == nf);
2420 ARTS_ASSERT(abs_data.ncols() == nt);
2421 ARTS_ASSERT(pfun_data.npages() == nf);
2422 ARTS_ASSERT(pfun_data.nrows() == nt);
2423 ARTS_ASSERT(pfun_data.ncols() == nsa);
2424 ARTS_ASSERT(cloudbox_limits.nelem() == 2);
2425 ARTS_ASSERT(pnd_data.nrows() == nse);
2426 ARTS_ASSERT(pnd_data.ncols() == ncl);
2427
2428 // Check that all data is TRO
2429 {
2430 bool all_totrand = true;
2431 for (Index ie = 0; ie < nse; ie++) {
2432 if (scat_data[ie].ptype != PTYPE_TOTAL_RND)
2433 all_totrand = false;
2434 }
2435 ARTS_USER_ERROR_IF (!all_totrand,
2436 "This method demands that all scat_data are TRO");
2437 }
2438
2439 // Help variables to hold non-zero data inside the cloudbox
2440 const Index cl_start = cloudbox_limits[0];
2441 Vector T_values(ncl);
2442 ArrayOfIndex cboxlayer(ncl);
2443
2444 // Loop scattering elements
2445 for (Index ie = 0; ie < nse; ie++) {
2446 // Allowed temperature range
2447 const Index last = scat_data[ie].T_grid.nelem() - 1;
2448 const Numeric tmin = 1.5*scat_data[ie].T_grid[0] -
2449 0.5*scat_data[ie].T_grid[1];
2450 const Numeric tmax = 1.5*scat_data[ie].T_grid[last] -
2451 0.5*scat_data[ie].T_grid[last-1];
2452
2453 Index nvals = 0;
2454 for(Index icl=0; icl<ncl; ++icl){
2455 // Nothing to do if PND is zero
2456 if (abs(pnd_data(ie,icl)) > 1e-3) {
2457 const Numeric Tthis = T_grid[cl_start + icl];
2458 ARTS_USER_ERROR_IF(Tthis < tmin || Tthis > tmax,
2459 "Temperature interpolation error for scattering element ", ie,
2460 " of species ", iss, "\nYour temperature: ", Tthis, " K\n"
2461 "Allowed range of temperatures: ", tmin, " - ", tmax, " K");
2462 T_values[nvals] = Tthis;
2463 cboxlayer[nvals] = icl;
2464 nvals++;
2465 }
2466 }
2467
2468 if (nvals > 0) {
2469 // Temperature-only interpolation weights
2470 ArrayOfGridPos gp_t(nvals);
2471 gridpos(gp_t, scat_data[ie].T_grid, T_values[Range(0,nvals)]);
2472 Matrix itw1(nvals, 2);
2473 interpweights(itw1, gp_t);
2474
2475 // Temperature + scattering angle interpolation weights
2476 ArrayOfGridPos gp_sa(nsa);
2477 gridpos(gp_sa, scat_data[ie].za_grid, sa_grid);
2478 Tensor3 itw2(nvals, nsa, 4);
2479 interpweights(itw2, gp_t, gp_sa);
2480
2481 // Loop frequencies
2482 for (Index iv = 0; iv < nf; iv++) {
2483 // Interpolate
2484 Vector ext1(nvals), abs1(nvals);
2485 Matrix pfu1(nvals,nsa);
2486 interp(ext1, itw1, scat_data[ie].ext_mat_data(iv,joker,0,0,0), gp_t);
2487 interp(abs1, itw1, scat_data[ie].abs_vec_data(iv,joker,0,0,0), gp_t);
2488 interp(pfu1, itw2, scat_data[ie].pha_mat_data(iv,joker,joker,0,0,0,0),
2489 gp_t, gp_sa);
2490
2491 // Add to container variables
2492 for (Index i = 0; i < nvals; i++) {
2493 const Index ic = cboxlayer[i];
2494 const Index it = cl_start + ic;
2495 ext_data(iv,it) += pnd_data(ie,ic) * ext1[i];
2496 abs_data(iv,it) += pnd_data(ie,ic) * abs1[i];
2497 for (Index ia = 0; ia < nsa; ia++) {
2498 pfun_data(iv,it,ia) += pnd_data(ie,ic) * pfu1(i,ia);
2499 }
2500 }
2501 }
2502 }
2503 }
2504}
This file contains the definition of Array.
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:224
base max(const Array< base > &x)
Max function.
Definition: array.h:128
base min(const Array< base > &x)
Min function.
Definition: array.h:144
The global header file for ARTS.
Common ARTS conversions.
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:887
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:402
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
#define _U_
Definition: config.h:177
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
#define ARTS_USER_ERROR(...)
Definition: debug.h:151
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:142
Declarations having to do with the four output streams.
#define CREATE_OUT0
Definition: messages.h:186
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric two_pi
Two times pi.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto sind(auto x) noexcept
Returns the sine of the deg2rad of the input.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
PType PType2FromString(const String &ptype_string)
Convert ptype name to enum value.
#define F22
#define F44
void ext_mat_SSD2Stokes(MatrixView ext_mat_stokes, ConstVectorView ext_mat_ssd, const Index &stokes_dim, const Index &ptype)
Extinction matrix scat_data to stokes format conversion.
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.
constexpr Numeric DEG2RAD
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
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_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
Numeric scat_angle(const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc)
Calculates the scattering angle.
ostream & operator<<(ostream &os, const SingleScatteringData &)
void pha_mat_labCalc(MatrixView pha_mat_lab, ConstVectorView pha_mat_int, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &theta_rad)
Calculate phase matrix in laboratory coordinate system.
void abs_vecTransform(StokesVector &abs_vec_lab, ConstTensor3View abs_vec_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of absorption vector.
ArrayOfLagrangeInterpolation ssd_tinterp_parameters(VectorView t_ok, Index &this_T_interp_order, ConstVectorView T_grid, const Vector &T_array, const Index &t_interp_order)
Determine T-interpol parameters for a specific scattering element.
constexpr Numeric RAD2DEG
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
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 ext_matTransform(PropagationMatrix &ext_mat_lab, ConstTensor3View ext_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of extinction matrix.
#define F34
#define F11
void ConvertAzimuthallyRandomSingleScatteringData(SingleScatteringData &ssd)
Convert azimuthally-random oriented SingleScatteringData to latest version.
void abs_vec_SSD2Stokes(VectorView abs_vec_stokes, ConstVectorView abs_vec_ssd, const Index &stokes_dim, const Index &ptype)
Absorption vector scat_data to stokes format conversion.
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
ParticleSSDMethod ParticleSSDMethodFromString(const String &particle_ssdmethod_string)
Convert particle ssd method name to enum value.
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.
#define F33
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.
void interpolate_scat_angle(VectorView pha_mat_int, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, const Numeric theta)
Interpolate data on the scattering angle.
void opt_prop_1ScatElem(Tensor5View ext_mat, Tensor4View abs_vec, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &dir_array, const Index &f_start, const Index &t_interp_order)
Preparing extinction and absorption from one scattering element.
#define F12
constexpr Numeric PI
void pha_matTransform(MatrixView pha_mat_lab, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Index &za_sca_idx, const Index &aa_sca_idx, const Index &za_inc_idx, const Index &aa_inc_idx, ConstVectorView za_grid, ConstVectorView aa_grid, const Verbosity &verbosity)
Transformation of phase matrix.
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
void ext_abs_pfun_from_tro(MatrixView ext_data, MatrixView abs_data, Tensor3View pfun_data, const ArrayOfSingleScatteringData &scat_data, const Index &iss, ConstMatrixView pnd_data, ArrayOfIndex &cloudbox_limits, ConstVectorView T_grid, ConstVectorView sa_grid)
Extinction, absorption and phase function for one scattering species, 1D and TRO.
Scattering database structure and functions.
PType
An attribute to classify the particle type (ptype) of a SingleScatteringData.
Definition: optproperties.h:35
@ PTYPE_GENERAL
Definition: optproperties.h:36
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:37
@ PTYPE_TOTAL_RND
Definition: optproperties.h:38
ParticleSSDMethod
An attribute to classify the method to be used for SingleScatteringData.
Definition: optproperties.h:46
@ PARTICLE_SSDMETHOD_TMATRIX
Definition: optproperties.h:48
Structure to store a grid position.
Definition: interpolation.h:56
This file contains basic functions to handle XML data files.