ARTS 2.5.4 (git: 31ce4f0e)
optproperties.cc
Go to the documentation of this file.
1/* Copyright (C) 2003-2012 Claudia Emde <claudia.emde@dlr.de>
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
18/*===========================================================================
19 === File description
20 ===========================================================================*/
21
33/*===========================================================================
34 === External declarations
35 ===========================================================================*/
36
37#include "optproperties.h"
38#include <cfloat>
39#include <cmath>
40#include <stdexcept>
41#include "array.h"
42#include "arts.h"
43#include "check_input.h"
44#include "arts_conversions.h"
45#include "interpolation.h"
47#include "logic.h"
48#include "math_funcs.h"
49#include "matpackVII.h"
50#include "messages.h"
51#include "xml_io.h"
52
55using Constant::pi;
56inline constexpr Numeric PI=pi;
57
58#define F11 pha_mat_int[0]
59#define F12 pha_mat_int[1]
60#define F22 pha_mat_int[2]
61#define F33 pha_mat_int[3]
62#define F34 pha_mat_int[4]
63#define F44 pha_mat_int[5]
64
66
75/*
76void methodname(//Output
77 type& name,
78 //Input
79 const type& name)
80{
81}
82*/
83
85
105void opt_prop_Bulk( //Output
106 Tensor5& ext_mat, // (nf,nT,ndir,nst,nst)
107 Tensor4& abs_vec, // (nf,nT,ndir,nst)
108 Index& ptype,
109 //Input
110 const ArrayOfTensor5& ext_mat_ss, // [nss](nf,nT,ndir,nst,nst)
111 const ArrayOfTensor4& abs_vec_ss, // [nss](nf,nT,ndir,nst)
112 const ArrayOfIndex& ptypes_ss) {
113 ARTS_ASSERT(ext_mat_ss.nelem() == abs_vec_ss.nelem());
114
115 ext_mat = ext_mat_ss[0];
116 abs_vec = abs_vec_ss[0];
117
118 for (Index i_ss = 1; i_ss < ext_mat_ss.nelem(); i_ss++) {
119 ext_mat += ext_mat_ss[i_ss];
120 abs_vec += abs_vec_ss[i_ss];
121 }
122 ptype = max(ptypes_ss);
123}
124
126
152 ArrayOfTensor5& ext_mat, // [nss](nf,nT,ndir,nst,nst)
153 ArrayOfTensor4& abs_vec, // [nss](nf,nT,ndir,nst)
154 ArrayOfIndex& ptype,
155 //Input
156 const ArrayOfArrayOfTensor5& ext_mat_se, // [nss][nse](nf,nT,ndir,nst,nst)
157 const ArrayOfArrayOfTensor4& abs_vec_se, // [nss][nse](nf,nT,ndir,nst)
158 const ArrayOfArrayOfIndex& ptypes_se,
159 ConstMatrixView pnds,
160 ConstMatrixView t_ok) {
161 ARTS_ASSERT(t_ok.nrows() == pnds.nrows());
162 ARTS_ASSERT(t_ok.ncols() == pnds.ncols());
163 ARTS_ASSERT(TotalNumberOfElements(ext_mat_se) == pnds.nrows());
164 ARTS_ASSERT(TotalNumberOfElements(abs_vec_se) == pnds.nrows());
165 ARTS_ASSERT(ext_mat_se.nelem() == abs_vec_se.nelem());
166
167 const Index nT = pnds.ncols();
168 const Index nf = abs_vec_se[0][0].nbooks();
169 const Index nDir = abs_vec_se[0][0].nrows();
170 const Index stokes_dim = abs_vec_se[0][0].ncols();
171
172 const Index nss = ext_mat_se.nelem();
173 ext_mat.resize(nss);
174 abs_vec.resize(nss);
175 ptype.resize(nss);
176 Tensor4 ext_tmp;
177 Tensor3 abs_tmp;
178
179 Index i_se_flat = 0;
180
181 for (Index i_ss = 0; i_ss < nss; i_ss++) {
182 ARTS_ASSERT(ext_mat_se[i_ss].nelem() == abs_vec_se[i_ss].nelem());
183 ARTS_ASSERT(nT == ext_mat_se[i_ss][0].nbooks());
184 ARTS_ASSERT(nT == abs_vec_se[i_ss][0].npages());
185
186 ext_mat[i_ss].resize(nf, nT, nDir, stokes_dim, stokes_dim);
187 ext_mat[i_ss] = 0.;
188 abs_vec[i_ss].resize(nf, nT, nDir, stokes_dim);
189 abs_vec[i_ss] = 0.;
190
191 for (Index i_se = 0; i_se < ext_mat_se[i_ss].nelem(); i_se++) {
192 ARTS_ASSERT(nT == ext_mat_se[i_ss][i_se].nbooks());
193 ARTS_ASSERT(nT == abs_vec_se[i_ss][i_se].npages());
194
195 for (Index Tind = 0; Tind < nT; Tind++) {
196 if (pnds(i_se_flat, Tind) != 0.) {
197 if (t_ok(i_se_flat, Tind) > 0.) {
198 ext_tmp = ext_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker);
199 ext_tmp *= pnds(i_se_flat, Tind);
200 ext_mat[i_ss](joker, Tind, joker, joker, joker) += ext_tmp;
201
202 abs_tmp = abs_vec_se[i_ss][i_se](joker, Tind, joker, joker);
203 abs_tmp *= pnds(i_se_flat, Tind);
204 abs_vec[i_ss](joker, Tind, joker, joker) += abs_tmp;
205 } else {
207 "Interpolation error for (flat-array) scattering element #",
208 i_se_flat, "\n"
209 "at location/temperature point #", Tind, "\n")
210 }
211 }
212 }
213 i_se_flat++;
214 }
215 ptype[i_ss] = max(ptypes_se[i_ss]);
216 }
217}
218
220
251 ArrayOfArrayOfTensor5& ext_mat, // [nss][nse](nf,nT,ndir,nst,nst)
252 ArrayOfArrayOfTensor4& abs_vec, // [nss][nse](nf,nT,ndir,nst)
253 ArrayOfArrayOfIndex& ptypes,
254 Matrix& t_ok,
255 //Input
256 const ArrayOfArrayOfSingleScatteringData& scat_data,
257 const Index& stokes_dim,
258 const Vector& T_array,
259 const Matrix& dir_array,
260 const Index& f_index,
261 const Index& t_interp_order) {
262 Index f_start, nf;
263 if (f_index < 0) {
264 nf = scat_data[0][0].ext_mat_data.nshelves();
265 f_start = 0;
266 //f_end = f_start+nf;
267 } else {
268 nf = 1;
269 if (scat_data[0][0].ext_mat_data.nshelves() == 1)
270 f_start = 0;
271 else
272 f_start = f_index;
273 //f_end = f_start+nf;
274 }
275
276 const Index nT = T_array.nelem();
277 const Index nDir = dir_array.nrows();
278
279 const Index nss = scat_data.nelem();
280 ext_mat.resize(nss);
281 abs_vec.resize(nss);
282 ptypes.resize(nss);
283
284 const Index Nse_all = TotalNumberOfElements(scat_data);
285 t_ok.resize(Nse_all, nT);
286 Index i_se_flat = 0;
287
288 for (Index i_ss = 0; i_ss < nss; i_ss++) {
289 Index nse = scat_data[i_ss].nelem();
290 ext_mat[i_ss].resize(nse);
291 abs_vec[i_ss].resize(nse);
292 ptypes[i_ss].resize(nse);
293
294 for (Index i_se = 0; i_se < nse; i_se++) {
295 ext_mat[i_ss][i_se].resize(nf, nT, nDir, stokes_dim, stokes_dim);
296 abs_vec[i_ss][i_se].resize(nf, nT, nDir, stokes_dim);
297
298 opt_prop_1ScatElem(ext_mat[i_ss][i_se],
299 abs_vec[i_ss][i_se],
300 ptypes[i_ss][i_se],
301 t_ok(i_se_flat, joker),
302 scat_data[i_ss][i_se],
303 T_array,
304 dir_array,
305 f_start,
306 t_interp_order);
307 i_se_flat++;
308 }
309 }
310 ARTS_ASSERT(i_se_flat == Nse_all);
311}
312
314
332 VectorView t_ok,
333 Index& this_T_interp_order,
334 //Input
335 ConstVectorView T_grid,
336 const Vector& T_array,
337 const Index& t_interp_order) {
338 const Index nTin = T_grid.nelem();
339 const Index nTout = T_array.nelem();
340
341 this_T_interp_order = -1;
342
343
344 if (nTin > 1) {
345 this_T_interp_order = min(t_interp_order, nTin - 1);
346
347 // we need to check T-grid exceedance. and catch these cases (because T
348 // is assumed to correspond to a location and T-exceedance is ok when pnd=0
349 // there. however, here we do not know pnd.) and report them to have the
350 // calling method dealing with it.
351
352 // we set the extrapolfax explicitly and use it here as well as in
353 // gridpos_poly below.
354 const Numeric extrapolfac = 0.5;
355 const Numeric lowlim = T_grid[0] - extrapolfac * (T_grid[1] - T_grid[0]);
356 const Numeric uplim =
357 T_grid[nTin - 1] + extrapolfac * (T_grid[nTin - 1] - T_grid[nTin - 2]);
358
359 bool any_T_exceed = false;
360 for (Index Tind = 0; Tind < nTout; Tind++) {
361 if (T_array[Tind] < lowlim || T_array[Tind] > uplim) {
362 t_ok[Tind] = -1.;
363 any_T_exceed = true;
364 } else
365 t_ok[Tind] = 1.;
366 }
367
368 if (any_T_exceed) {
369 // Reserve output
371 T_lag.reserve(nTout);
372
373 bool grid_unchecked = true;
374
375 for (Index iT = 0; iT < nTout; iT++) {
376 if (t_ok[iT] < 0) {
377 T_lag.emplace_back();
378 } else {
379 if (grid_unchecked) {
381 "Temperature interpolation in pha_mat_1ScatElem",
382 T_grid,
383 T_array[Range(iT, 1)],
384 this_T_interp_order);
385 grid_unchecked = false;
386 }
387 T_lag.emplace_back(0, T_array[iT], T_grid, this_T_interp_order, false);
388 }
389 }
390 return T_lag;
391 } else {
392 return Interpolation::LagrangeVector(T_array, T_grid, this_T_interp_order, extrapolfac);
393 }
394 } else {
395 t_ok = 1.;
396 return {};
397 }
398}
399
401
424void opt_prop_1ScatElem( //Output
425 Tensor5View ext_mat, // nf, nT, ndir, nst, nst
426 Tensor4View abs_vec, // nf, nT, ndir, nst
427 Index& ptype,
428 VectorView t_ok,
429 //Input
430 const SingleScatteringData& ssd,
431 const Vector& T_array,
432 const Matrix& dir_array,
433 const Index& f_start,
434 const Index& t_interp_order) {
435 // FIXME: this is prob best done in scat_data_checkedCalc (or
436 // cloudbox_checkedCalc) to have it done once and for all. Here at max ARTS_ASSERT.
437
438 // At very first check validity of the scatt elements ptype (so far we only
439 // handle PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND).
440 /*
441 if( ssd.ptype != PTYPE_TOTAL_RND and ssd.ptype != PTYPE_AZIMUTH_RND )
442 {
443 ostringstream os;
444 os << "Only ptypes " << PTYPE_TOTAL_RND << " and " << PTYPE_AZIMUTH_RND
445 << " can be handled.\n"
446 << "Encountered scattering element with ptype " << ssd.ptype
447 << ", though.";
448 throw runtime_error( os.str() );
449 }
450 */
451
453
454 const Index nf = ext_mat.nshelves();
455 ARTS_ASSERT(abs_vec.nbooks() == nf);
456 if (nf > 1) {
457 ARTS_ASSERT(nf == ssd.f_grid.nelem());
458 }
459
460 const Index nTout = T_array.nelem();
461 ARTS_ASSERT(ext_mat.nbooks() == nTout);
462 ARTS_ASSERT(abs_vec.npages() == nTout);
463 ARTS_ASSERT(t_ok.nelem() == nTout);
464
465 const Index nDir = dir_array.nrows();
466 ARTS_ASSERT(ext_mat.npages() == nDir);
467 ARTS_ASSERT(abs_vec.nrows() == nDir);
468
469 const Index stokes_dim = abs_vec.ncols();
470 ARTS_ASSERT(ext_mat.nrows() == stokes_dim);
471 ARTS_ASSERT(ext_mat.ncols() == stokes_dim);
472
473 ptype = ssd.ptype;
474
475 // Determine T-interpol order as well as interpol positions and weights (they
476 // are the same for all directions (and freqs), ie it is sufficient to
477 // calculate them once).
478 const Index nTin = ssd.T_grid.nelem();
479 Index this_T_interp_order;
480 const auto T_lag = ssd_tinterp_parameters(t_ok,
481 this_T_interp_order,
482 ssd.T_grid,
483 T_array,
484 t_interp_order);
485 const auto T_itw_lag = interpweights(T_lag);
486
487 // Now loop over requested directions (and apply simultaneously for all freqs):
488 // 1) extract/interpolate direction (but not for tot.random)
489 // 2) apply T-interpol
490 if (ptype == PTYPE_TOTAL_RND)
491 if (this_T_interp_order < 0) // just extract (and unpack) ext and abs data
492 // for the fs, the Tin, and the dirin and sort
493 // (copy) into the output fs, Ts, and dirs.
494 {
495 Tensor3 ext_mat_tmp(nf, stokes_dim, stokes_dim);
496 Matrix abs_vec_tmp(nf, stokes_dim);
497 for (Index find = 0; find < nf; find++) {
498 ext_mat_SSD2Stokes(ext_mat_tmp(find, joker, joker),
499 ssd.ext_mat_data(find + f_start, 0, 0, 0, joker),
500 stokes_dim,
501 ptype);
502 abs_vec_SSD2Stokes(abs_vec_tmp(find, joker),
503 ssd.abs_vec_data(find + f_start, 0, 0, 0, joker),
504 stokes_dim,
505 ptype);
506 }
507 for (Index Tind = 0; Tind < nTout; Tind++)
508 for (Index dind = 0; dind < nDir; dind++) {
509 ext_mat(joker, Tind, dind, joker, joker) = ext_mat_tmp;
510 abs_vec(joker, Tind, dind, joker) = abs_vec_tmp;
511 }
512 } else // T-interpolation required (but not dir). To be done on the compact
513 // ssd format.
514 {
515 Tensor4 ext_mat_tmp(nf, nTout, stokes_dim, stokes_dim);
516 Tensor3 abs_vec_tmp(nf, nTout, stokes_dim);
517 Matrix ext_mat_tmp_ssd(nTout, ssd.ext_mat_data.ncols());
518 Matrix abs_vec_tmp_ssd(nTout, ssd.abs_vec_data.ncols());
519 for (Index find = 0; find < nf; find++) {
520 for (Index nst = 0; nst < ext_mat_tmp_ssd.ncols(); nst++) {
521 reinterp(ext_mat_tmp_ssd(joker, nst),
522 ssd.ext_mat_data(find + f_start, joker, 0, 0, nst),
523 T_itw_lag, T_lag);
524 }
525 for (Index Tind = 0; Tind < nTout; Tind++)
526 if (t_ok[Tind] > 0.)
527 ext_mat_SSD2Stokes(ext_mat_tmp(find, Tind, joker, joker),
528 ext_mat_tmp_ssd(Tind, joker),
529 stokes_dim,
530 ptype);
531 else
532 ext_mat_tmp(find, Tind, joker, joker) = 0.;
533
534 for (Index nst = 0; nst < abs_vec_tmp_ssd.ncols(); nst++) {
535 reinterp(abs_vec_tmp_ssd(joker, nst),
536 ssd.abs_vec_data(find + f_start, joker, 0, 0, nst),
537 T_itw_lag,
538 T_lag);
539 }
540 for (Index Tind = 0; Tind < nTout; Tind++)
541 if (t_ok[Tind] > 0.)
542 abs_vec_SSD2Stokes(abs_vec_tmp(find, Tind, joker),
543 abs_vec_tmp_ssd(Tind, joker),
544 stokes_dim,
545 ptype);
546 else
547 abs_vec_tmp(find, Tind, joker) = 0.;
548 }
549
550 for (Index dind = 0; dind < nDir; dind++) {
551 ext_mat(joker, joker, dind, joker, joker) = ext_mat_tmp;
552 abs_vec(joker, joker, dind, joker) = abs_vec_tmp;
553 }
554 }
555
556 else // dir-interpolation for non-tot-random particles
557 {
558 // as we don't allow other than az.random here, ext and abs will only vary
559 // with za, not with aa. Hence, we could be smart here and calc data only
560 // for unique za (and copy the rest). however, this smartness might cost as
561 // well. so for now, we leave that and blindly loop over the given direction
562 // array, and leave smartness in setting up directional array to the calling
563 // methods.
564
565 // derive the direction interpolation weights.
566 ArrayOfGridPos dir_gp(nDir);
567 gridpos(dir_gp, ssd.za_grid, dir_array(joker, 0));
568 Matrix dir_itw(nDir, 2); // only interpolating in za, ie 1D linear interpol
569 interpweights(dir_itw, dir_gp);
570
571 const Index next = ssd.ext_mat_data.ncols();
572 const Index nabs = ssd.abs_vec_data.ncols();
573
574 if (this_T_interp_order < 0) // T only needs to be extracted.
575 {
576 Matrix ext_mat_tmp_ssd(nDir, next);
577 Matrix abs_vec_tmp_ssd(nDir, nabs);
578 Tensor4 ext_mat_tmp(nf, nDir, stokes_dim, stokes_dim);
579 Tensor3 abs_vec_tmp(nf, nDir, stokes_dim);
580 for (Index find = 0; find < nf; find++) {
581 for (Index nst = 0; nst < next; nst++)
582 interp(ext_mat_tmp_ssd(joker, nst),
583 dir_itw,
584 ssd.ext_mat_data(find + f_start, 0, joker, 0, nst),
585 dir_gp);
586 for (Index Dind = 0; Dind < nDir; Dind++)
587 ext_mat_SSD2Stokes(ext_mat_tmp(find, Dind, joker, joker),
588 ext_mat_tmp_ssd(Dind, joker),
589 stokes_dim,
590 ptype);
591
592 for (Index nst = 0; nst < nabs; nst++)
593 interp(abs_vec_tmp_ssd(joker, nst),
594 dir_itw,
595 ssd.abs_vec_data(find + f_start, 0, joker, 0, nst),
596 dir_gp);
597 for (Index Dind = 0; Dind < nDir; Dind++)
598 abs_vec_SSD2Stokes(abs_vec_tmp(find, Dind, joker),
599 abs_vec_tmp_ssd(Dind, joker),
600 stokes_dim,
601 ptype);
602 }
603
604 for (Index Tind = 0; Tind < nTout; Tind++) {
605 ext_mat(joker, Tind, joker, joker, joker) = ext_mat_tmp;
606 abs_vec(joker, Tind, joker, joker) = abs_vec_tmp;
607 }
608 } else // T- and dir-interpolation required. To be done on the compact ssd
609 // format.
610 {
611 Tensor3 ext_mat_tmp_ssd(nTin, nDir, next);
612 Tensor3 abs_vec_tmp_ssd(nTin, nDir, nabs);
613 Matrix ext_mat_tmp(nTout, next);
614 Matrix abs_vec_tmp(nTout, nabs);
615
616 for (Index find = 0; find < nf; find++) {
617 for (Index Tind = 0; Tind < nTin; Tind++) {
618 for (Index nst = 0; nst < next; nst++)
619 interp(ext_mat_tmp_ssd(Tind, joker, nst),
620 dir_itw,
621 ssd.ext_mat_data(find + f_start, Tind, joker, 0, nst),
622 dir_gp);
623 for (Index nst = 0; nst < nabs; nst++)
624 interp(abs_vec_tmp_ssd(Tind, joker, nst),
625 dir_itw,
626 ssd.abs_vec_data(find + f_start, Tind, joker, 0, nst),
627 dir_gp);
628 }
629
630 for (Index Dind = 0; Dind < nDir; Dind++) {
631 for (Index nst = 0; nst < next; nst++) {
632 reinterp(ext_mat_tmp(joker, nst),
633 ext_mat_tmp_ssd(joker, Dind, nst),
634 T_itw_lag,
635 T_lag);
636 }
637 for (Index Tind = 0; Tind < nTout; Tind++)
638 ext_mat_SSD2Stokes(ext_mat(find, Tind, Dind, joker, joker),
639 ext_mat_tmp(Tind, joker),
640 stokes_dim,
641 ptype);
642
643 for (Index nst = 0; nst < nabs; nst++) {
644 reinterp(abs_vec_tmp(joker, nst),
645 abs_vec_tmp_ssd(joker, Dind, nst),
646 T_itw_lag,
647 T_lag);
648 }
649 for (Index Tind = 0; Tind < nTout; Tind++)
650 abs_vec_SSD2Stokes(abs_vec(find, Tind, Dind, joker),
651 abs_vec_tmp(Tind, joker),
652 stokes_dim,
653 ptype);
654 }
655 }
656 }
657 }
658}
659
661
674void ext_mat_SSD2Stokes( //Output
675 MatrixView ext_mat_stokes,
676 //Input
677 ConstVectorView ext_mat_ssd,
678 const Index& stokes_dim,
679 const Index& ptype) {
680 // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
681 // calling methods, though.
683
684 ext_mat_stokes = 0.;
685
686 for (Index ist = 0; ist < stokes_dim; ist++) {
687 ext_mat_stokes(ist, ist) = ext_mat_ssd[0];
688 }
689
690 if (ptype > PTYPE_TOTAL_RND) {
691 switch (stokes_dim) {
692 case 4: {
693 ext_mat_stokes(2, 3) = ext_mat_ssd[2];
694 ext_mat_stokes(3, 2) = -ext_mat_ssd[2];
695 } /* FALLTHROUGH */
696 case 3: {
697 // nothing to be done here. but we need this for executing the below
698 // also in case of stokes_dim==3.
699 }
700 case 2: {
701 ext_mat_stokes(0, 1) = ext_mat_ssd[1];
702 ext_mat_stokes(1, 0) = ext_mat_ssd[1];
703 }
704 }
705 }
706}
707
709
722void abs_vec_SSD2Stokes( //Output
723 VectorView abs_vec_stokes,
724 //Input
725 ConstVectorView abs_vec_ssd,
726 const Index& stokes_dim,
727 const Index& ptype) {
728 // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
729 // calling methods, though.
731
732 abs_vec_stokes = 0.;
733
734 abs_vec_stokes[0] = abs_vec_ssd[0];
735
736 if (ptype > PTYPE_TOTAL_RND and stokes_dim > 1) {
737 abs_vec_stokes[1] = abs_vec_ssd[1];
738 }
739}
740
742
759void pha_mat_Bulk( //Output
760 Tensor6& pha_mat, // (nf,nT,npdir,nidir,nst,nst)
761 Index& ptype,
762 //Input
763 const ArrayOfTensor6& pha_mat_ss, // [nss](nf,nT,npdir,nidir,nst,nst)
764 const ArrayOfIndex& ptypes_ss) {
765 pha_mat = pha_mat_ss[0];
766
767 for (Index i_ss = 1; i_ss < pha_mat_ss.nelem(); i_ss++)
768 pha_mat += pha_mat_ss[i_ss];
769
770 ptype = max(ptypes_ss);
771}
772
774
798 ArrayOfTensor6& pha_mat, // [nss](nf,nT,npdir,nidir,nst,nst)
799 ArrayOfIndex& ptype,
800 //Input
802 pha_mat_se, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
803 const ArrayOfArrayOfIndex& ptypes_se,
804 ConstMatrixView pnds,
805 ConstMatrixView t_ok) {
806 ARTS_ASSERT(t_ok.nrows() == pnds.nrows());
807 ARTS_ASSERT(t_ok.ncols() == pnds.ncols());
808 ARTS_ASSERT(TotalNumberOfElements(pha_mat_se) == pnds.nrows());
809
810 const Index nT = pnds.ncols();
811 const Index nf = pha_mat_se[0][0].nvitrines();
812 const Index npDir = pha_mat_se[0][0].nbooks();
813 const Index niDir = pha_mat_se[0][0].npages();
814 const Index stokes_dim = pha_mat_se[0][0].ncols();
815
816 const Index nss = pha_mat_se.nelem();
817 pha_mat.resize(nss);
818 ptype.resize(nss);
819 Tensor5 pha_tmp;
820
821 Index i_se_flat = 0;
822
823 for (Index i_ss = 0; i_ss < nss; i_ss++) {
824 ARTS_ASSERT(nT == pha_mat_se[i_ss][0].nshelves());
825
826 pha_mat[i_ss].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
827 pha_mat[i_ss] = 0.;
828
829 for (Index i_se = 0; i_se < pha_mat_se[i_ss].nelem(); i_se++) {
830 ARTS_ASSERT(nT == pha_mat_se[i_ss][i_se].nshelves());
831
832 for (Index Tind = 0; Tind < nT; Tind++) {
833 if (pnds(i_se_flat, Tind) != 0.) {
834 if (t_ok(i_se_flat, Tind) > 0.) {
835 pha_tmp =
836 pha_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker, joker);
837 pha_tmp *= pnds(i_se_flat, Tind);
838 pha_mat[i_ss](joker, Tind, joker, joker, joker, joker) += pha_tmp;
839 } else {
841 "Interpolation error for (flat-array) scattering element #",
842 i_se_flat, "\n"
843 "at location/temperature point #", Tind, "\n")
844 }
845 }
846 }
847 i_se_flat++;
848 }
849 ptype[i_ss] = max(ptypes_se[i_ss]);
850 }
851}
852
854
884void pha_mat_NScatElems( //Output
885 ArrayOfArrayOfTensor6& pha_mat, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
886 ArrayOfArrayOfIndex& ptypes,
887 Matrix& t_ok,
888 //Input
889 const ArrayOfArrayOfSingleScatteringData& scat_data,
890 const Index& stokes_dim,
891 const Vector& T_array,
892 const Matrix& pdir_array,
893 const Matrix& idir_array,
894 const Index& f_index,
895 const Index& t_interp_order) {
896 Index f_start, nf;
897 if (f_index < 0) {
898 nf = scat_data[0][0].pha_mat_data.nlibraries();
899 f_start = 0;
900 //f_end = f_start+nf;
901 } else {
902 nf = 1;
903 if (scat_data[0][0].pha_mat_data.nlibraries() == 1)
904 f_start = 0;
905 else
906 f_start = f_index;
907 //f_end = f_start+nf;
908 }
909
910 const Index nT = T_array.nelem();
911 const Index npDir = pdir_array.nrows();
912 const Index niDir = idir_array.nrows();
913
914 const Index nss = scat_data.nelem();
915 pha_mat.resize(nss);
916 ptypes.resize(nss);
917
918 const Index Nse_all = TotalNumberOfElements(scat_data);
919 t_ok.resize(Nse_all, nT);
920 Index i_se_flat = 0;
921
922 for (Index i_ss = 0; i_ss < nss; i_ss++) {
923 Index nse = scat_data[i_ss].nelem();
924 pha_mat[i_ss].resize(nse);
925 ptypes[i_ss].resize(nse);
926
927 for (Index i_se = 0; i_se < nse; i_se++) {
928 pha_mat[i_ss][i_se].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
929
930 pha_mat_1ScatElem(pha_mat[i_ss][i_se],
931 ptypes[i_ss][i_se],
932 t_ok(i_se_flat, joker),
933 scat_data[i_ss][i_se],
934 T_array,
935 pdir_array,
936 idir_array,
937 f_start,
938 t_interp_order);
939 i_se_flat++;
940 }
941 }
942 ARTS_ASSERT(i_se_flat == Nse_all);
943}
944
946
969void pha_mat_1ScatElem( //Output
970 Tensor6View pha_mat, // nf, nT, npdir, nidir, nst, nst
971 Index& ptype,
972 VectorView t_ok,
973 //Input
974 const SingleScatteringData& ssd,
975 const Vector& T_array,
976 const Matrix& pdir_array,
977 const Matrix& idir_array,
978 const Index& f_start,
979 const Index& t_interp_order) {
981
982 const Index nf = pha_mat.nvitrines();
983 if (nf > 1) {
984 ARTS_ASSERT(nf == ssd.f_grid.nelem());
985 }
986
987 const Index nTout = T_array.nelem();
988 ARTS_ASSERT(pha_mat.nshelves() == nTout);
989 ARTS_ASSERT(t_ok.nelem() == nTout);
990
991 const Index npDir = pdir_array.nrows();
992 ARTS_ASSERT(pha_mat.nbooks() == npDir);
993
994 const Index niDir = idir_array.nrows();
995 ARTS_ASSERT(pha_mat.npages() == niDir);
996
997 const Index stokes_dim = pha_mat.ncols();
998 ARTS_ASSERT(pha_mat.nrows() == stokes_dim);
999
1000 ptype = ssd.ptype;
1001
1002 // Determine T-interpol order as well as interpol positions and weights (they
1003 // are the same for all directions (and freqs), ie it is sufficient to
1004 // calculate them once).
1005 const Index nTin = ssd.T_grid.nelem();
1006 Index this_T_interp_order;
1007 const auto T_lag = ssd_tinterp_parameters(t_ok,
1008 this_T_interp_order,
1009 ssd.T_grid,
1010 T_array,
1011 t_interp_order);
1012 const auto T_itw_lag = interpweights(T_lag);
1013
1014 // Now loop over requested directions (and apply simultaneously for all freqs):
1015 // 1) interpolate direction
1016 // 2) apply T-interpol
1017 if (ptype == PTYPE_TOTAL_RND) {
1018 // determine how many of the compact stokes elements we will need.
1019 // restrict interpolations to those.
1020 Index npha;
1021 if (stokes_dim == 1)
1022 npha = 1;
1023 else if (stokes_dim < 4) // stokes_dim==2 || stokes_dim==3
1024 npha = 4;
1025 else
1026 npha = 6;
1027 if (this_T_interp_order < 0) // just extract (and unpack) pha data for the
1028 // fs and Tin, and interpolate in sca-angs, and
1029 // sort (copy) into the output fs, Ts, and dirs.
1030 {
1031 for (Index pdir = 0; pdir < npDir; pdir++)
1032 for (Index idir = 0; idir < niDir; idir++) {
1033 // calc scat ang theta from incident and prop dirs
1034 Numeric theta = scat_angle(pdir_array(pdir, 0),
1035 pdir_array(pdir, 1),
1036 idir_array(idir, 0),
1037 idir_array(idir, 1));
1038
1039 // get scat angle interpolation weights
1040 GridPos dir_gp;
1041 gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
1042 Vector dir_itw(2);
1043 interpweights(dir_itw, dir_gp);
1044
1045 Vector pha_mat_int(npha, 0.);
1046 Matrix pha_mat_tmp(stokes_dim, stokes_dim);
1047 for (Index find = 0; find < nf; find++) {
1048 // perform the scat angle interpolation
1049 for (Index nst = 0; nst < npha; nst++)
1050 pha_mat_int[nst] = interp(
1051 dir_itw,
1052 ssd.pha_mat_data(find + f_start, 0, joker, 0, 0, 0, nst),
1053 dir_gp);
1054
1055 // convert from scat to lab frame
1056 pha_mat_labCalc(pha_mat_tmp(joker, joker),
1057 pha_mat_int,
1058 pdir_array(pdir, 0),
1059 pdir_array(pdir, 1),
1060 idir_array(idir, 0),
1061 idir_array(idir, 1),
1062 theta);
1063
1064 for (Index Tind = 0; Tind < nTout; Tind++)
1065 pha_mat(find, Tind, pdir, idir, joker, joker) = pha_mat_tmp;
1066 }
1067 }
1068 } else // T-interpolation required. To be done on the compact ssd format.
1069 {
1070 for (Index pdir = 0; pdir < npDir; pdir++)
1071 for (Index idir = 0; idir < niDir; idir++) {
1072 // calc scat ang theta from incident and prop dirs
1073 Numeric theta = scat_angle(pdir_array(pdir, 0),
1074 pdir_array(pdir, 1),
1075 idir_array(idir, 0),
1076 idir_array(idir, 1));
1077
1078 // get scat angle interpolation weights
1079 GridPos dir_gp;
1080 gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
1081 Vector dir_itw(2);
1082 interpweights(dir_itw, dir_gp);
1083
1084 Matrix pha_mat_int(nTin, npha, 0.);
1085 Matrix pha_mat_tmp(nTout, npha, 0.);
1086 for (Index find = 0; find < nf; find++) {
1087 for (Index Tind = 0; Tind < nTin; Tind++)
1088 // perform the scat angle interpolation
1089 for (Index nst = 0; nst < npha; nst++) {
1090 pha_mat_int(Tind, nst) = interp(
1091 dir_itw,
1092 ssd.pha_mat_data(find + f_start, Tind, joker, 0, 0, 0, nst),
1093 dir_gp);
1094 }
1095 // perform the T-interpolation
1096 for (Index nst = 0; nst < npha; nst++) {
1097 reinterp(pha_mat_tmp(joker, nst),
1098 pha_mat_int(joker, nst),
1099 T_itw_lag,
1100 T_lag);
1101 }
1102 // FIXME: it's probably better to do the frame conversion first,
1103 // then the T-interpolation (which is better depends on how many T-
1104 // aka vertical grid points we have. for single point, T-interpol
1105 // first is better, for a full vertical grid, frame conversion first
1106 // should be faster.)
1107 for (Index Tind = 0; Tind < nTout; Tind++) {
1108 // convert from scat to lab frame
1109 pha_mat_labCalc(pha_mat(find, Tind, pdir, idir, joker, joker),
1110 pha_mat_tmp(Tind, joker),
1111 pdir_array(pdir, 0),
1112 pdir_array(pdir, 1),
1113 idir_array(idir, 0),
1114 idir_array(idir, 1),
1115 theta);
1116 }
1117 }
1118 }
1119 }
1120 } else // dir-interpolation for non-tot-random particles
1121 // Data is already stored in the laboratory frame,
1122 // but it is compressed a little. Details elsewhere.
1123 {
1124 Index nDir = npDir * niDir;
1125 Vector adelta_aa(nDir);
1126 Matrix delta_aa(npDir, niDir);
1127 ArrayOfGridPos daa_gp(nDir), pza_gp(nDir), iza_gp(nDir);
1128 ArrayOfGridPos pza_gp_tmp(npDir), iza_gp_tmp(niDir);
1129
1130 gridpos(pza_gp_tmp, ssd.za_grid, pdir_array(joker, 0));
1131 gridpos(iza_gp_tmp, ssd.za_grid, idir_array(joker, 0));
1132
1133 Index j = 0;
1134 for (Index pdir = 0; pdir < npDir; pdir++) {
1135 for (Index idir = 0; idir < niDir; idir++) {
1136 delta_aa(pdir, idir) =
1137 pdir_array(pdir, 1) - idir_array(idir, 1) +
1138 (pdir_array(pdir, 1) - idir_array(idir, 1) < -180) * 360 -
1139 (pdir_array(pdir, 1) - idir_array(idir, 1) > 180) * 360;
1140 adelta_aa[j] = abs(delta_aa(pdir, idir));
1141 pza_gp[j] = pza_gp_tmp[pdir];
1142 iza_gp[j] = iza_gp_tmp[idir];
1143 j++;
1144 }
1145 }
1146
1147 gridpos(daa_gp, ssd.aa_grid, adelta_aa);
1148
1149 Matrix dir_itw(nDir, 8);
1150 interpweights(dir_itw, pza_gp, daa_gp, iza_gp);
1151
1152 if (this_T_interp_order < 0) // T only needs to be extracted.
1153 {
1154 Tensor3 pha_mat_int(nDir, stokes_dim, stokes_dim, 0.);
1155 Tensor4 pha_mat_tmp(npDir, niDir, stokes_dim, stokes_dim);
1156
1157 for (Index find = 0; find < nf; find++) {
1158 // perform the (tri-linear) angle interpolation. but only for the
1159 // pha_mat elements that we actually need.
1160
1161 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1162 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1163 interp(
1164 pha_mat_int(joker, ist1, ist2),
1165 dir_itw,
1166 ssd.pha_mat_data(
1167 find + f_start, 0, joker, joker, joker, 0, ist1 * 4 + ist2),
1168 pza_gp,
1169 daa_gp,
1170 iza_gp);
1171
1172 // sort direction-combined 1D-array back into prop and incident
1173 // direction matrix
1174 Index i = 0;
1175 for (Index pdir = 0; pdir < npDir; pdir++)
1176 for (Index idir = 0; idir < niDir; idir++) {
1177 pha_mat_tmp(pdir, idir, joker, joker) =
1178 pha_mat_int(i, joker, joker);
1179 i++;
1180 }
1181
1182 if (stokes_dim > 2) {
1183 for (Index pdir = 0; pdir < npDir; pdir++)
1184 for (Index idir = 0; idir < niDir; idir++)
1185 if (delta_aa(pdir, idir) < 0.) {
1186 pha_mat_tmp(pdir, idir, 0, 2) *= -1;
1187 pha_mat_tmp(pdir, idir, 1, 2) *= -1;
1188 pha_mat_tmp(pdir, idir, 2, 0) *= -1;
1189 pha_mat_tmp(pdir, idir, 2, 1) *= -1;
1190 }
1191 }
1192
1193 if (stokes_dim > 3) {
1194 for (Index pdir = 0; pdir < npDir; pdir++)
1195 for (Index idir = 0; idir < niDir; idir++)
1196 if (delta_aa(pdir, idir) < 0.) {
1197 pha_mat_tmp(pdir, idir, 0, 3) *= -1;
1198 pha_mat_tmp(pdir, idir, 1, 3) *= -1;
1199 pha_mat_tmp(pdir, idir, 3, 0) *= -1;
1200 pha_mat_tmp(pdir, idir, 3, 1) *= -1;
1201 }
1202 }
1203
1204 for (Index Tind = 0; Tind < nTout; Tind++)
1205 pha_mat(find, Tind, joker, joker, joker, joker) = pha_mat_tmp;
1206 }
1207 }
1208
1209 else // T- and dir-interpolation required. To be done on the compact ssd
1210 // format.
1211 {
1212 Tensor4 pha_mat_int(nTin, nDir, stokes_dim, stokes_dim, 0.);
1213
1214 for (Index find = 0; find < nf; find++) {
1215 // perform the (tri-linear) angle interpolation. but only for the
1216 // pha_mat elements that we actually need.
1217 for (Index Tind = 0; Tind < nTin; Tind++) {
1218 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1219 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1220 interp(pha_mat_int(Tind, joker, ist1, ist2),
1221 dir_itw,
1222 ssd.pha_mat_data(find + f_start,
1223 Tind,
1224 joker,
1225 joker,
1226 joker,
1227 0,
1228 ist1 * 4 + ist2),
1229 pza_gp,
1230 daa_gp,
1231 iza_gp);
1232 }
1233
1234 // perform the T-interpolation and simultaneously sort
1235 // direction-combined 1D-array back into prop and incident direction
1236 // matrix.
1237 Index i = 0;
1238 for (Index pdir = 0; pdir < npDir; pdir++)
1239 for (Index idir = 0; idir < niDir; idir++) {
1240 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1241 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1242 reinterp(pha_mat(find, joker, pdir, idir, ist1, ist2),
1243 pha_mat_int(joker, i, ist1, ist2),
1244 T_itw_lag,
1245 T_lag);
1246 i++;
1247 }
1248
1249 if (stokes_dim > 2) {
1250 for (Index pdir = 0; pdir < npDir; pdir++)
1251 for (Index idir = 0; idir < niDir; idir++)
1252 if (delta_aa(pdir, idir) < 0.) {
1253 pha_mat(find, joker, pdir, idir, 0, 2) *= -1;
1254 pha_mat(find, joker, pdir, idir, 1, 2) *= -1;
1255 pha_mat(find, joker, pdir, idir, 2, 0) *= -1;
1256 pha_mat(find, joker, pdir, idir, 2, 1) *= -1;
1257 }
1258 }
1259
1260 if (stokes_dim > 3) {
1261 for (Index pdir = 0; pdir < npDir; pdir++)
1262 for (Index idir = 0; idir < niDir; idir++)
1263 if (delta_aa(pdir, idir) < 0.) {
1264 pha_mat(find, joker, pdir, idir, 0, 3) *= -1;
1265 pha_mat(find, joker, pdir, idir, 1, 3) *= -1;
1266 pha_mat(find, joker, pdir, idir, 3, 0) *= -1;
1267 pha_mat(find, joker, pdir, idir, 3, 1) *= -1;
1268 }
1269 }
1270 }
1271 }
1272 }
1273}
1274
1276
1303void FouComp_1ScatElem( //Output
1304 Tensor7View pha_mat_fou, // nf, nT, npdir, nidir, nst, nst, m
1305 Index& ptype,
1306 VectorView t_ok,
1307 //Input
1308 const SingleScatteringData& ssd,
1309 const Vector& T_array,
1310 const Vector& pdir_array,
1311 const Vector& idir_array,
1312 const Index& f_start,
1313 const Index& t_interp_order,
1314 const Index& naa_totran) {
1316
1317 const Index nf = pha_mat_fou.nlibraries();
1318 if (nf > 1) {
1319 ARTS_ASSERT(nf == ssd.f_grid.nelem());
1320 }
1321
1322 const Index nTout = T_array.nelem();
1323 ARTS_ASSERT(pha_mat_fou.nvitrines() == nTout);
1324 ARTS_ASSERT(t_ok.nelem() == nTout);
1325
1326 const Index npDir = pdir_array.nelem();
1327 ARTS_ASSERT(pha_mat_fou.nshelves() == npDir);
1328 const Index niDir = idir_array.nelem();
1329 ARTS_ASSERT(pha_mat_fou.nbooks() == niDir);
1330
1331 const Index stokes_dim = pha_mat_fou.nrows();
1332 ARTS_ASSERT(pha_mat_fou.npages() == stokes_dim);
1333 // currently code is only prepared for stokes_dim up to 2 (nothing else needed in
1334 // RT4 and generally in azimuth-symmetrical system)
1335 ARTS_ASSERT(stokes_dim < 3);
1336
1337 const Index nmodes = pha_mat_fou.ncols();
1338 // currently code is only prepared for fourier mode 0 (nothing else needed in
1339 // RT4 and generally in azimuth-symmetrical system)
1340 ARTS_ASSERT(nmodes == 0);
1341
1342 ptype = ssd.ptype;
1343
1344 // Determine T-interpol order as well as interpol positions and weights (they
1345 // are the same for all directions (and freqs), ie it is sufficient to
1346 // calculate them once).
1347 const Index nTin = ssd.T_grid.nelem();
1348 Index this_T_interp_order;
1349 const auto T_lag = ssd_tinterp_parameters(t_ok,
1350 this_T_interp_order,
1351 ssd.T_grid,
1352 T_array,
1353 t_interp_order);
1354 const auto T_itw_lag = interpweights(T_lag);
1355
1356 // 1) derive Fourier component(s)
1357 // 2) apply T-interpol
1358 if (ptype == PTYPE_TOTAL_RND) {
1359 // DCalculate azimuth angles and their integration weights for Fourier
1360 // component derivation (they are only determined by naa_totran).
1361 ARTS_USER_ERROR_IF (naa_totran < 3,
1362 "Azimuth grid size for scatt matrix extraction"
1363 " (*naa_totran*) must be >3.\n"
1364 "Yours is ", naa_totran, ".\n")
1365 Vector aa_grid;
1366 nlinspace(aa_grid, 0, 180, naa_totran);
1367 Numeric daa_totran =
1368 1. / float(naa_totran - 1); // 2*180./360./(naa_totran-1)
1369 Vector theta(naa_totran);
1370 ArrayOfGridPos theta_gp;
1371 Matrix theta_itw(naa_totran, 2);
1372
1373 Index npha;
1374 if (stokes_dim == 1)
1375 npha = 1;
1376 else if (stokes_dim < 4) // stokes_dim==2 || stokes_dim==3
1377 npha = 4;
1378 else
1379 npha = 6;
1380
1381 Matrix pha_mat(stokes_dim, stokes_dim);
1382 Matrix pha_mat_angint(naa_totran, npha);
1383 Tensor3 Fou_int(nTin, stokes_dim, stokes_dim);
1384
1385 for (Index idir = 0; idir < niDir; idir++)
1386 for (Index pdir = 0; pdir < npDir; pdir++) {
1387 for (Index iaa = 0; iaa < naa_totran; iaa++)
1388 // calc scat ang theta from incident and prop dirs
1389 theta[iaa] =
1390 scat_angle(pdir_array[pdir], aa_grid[iaa], idir_array[idir], 0.);
1391
1392 // get scat angle interpolation weights
1393 gridpos(theta_gp, ssd.za_grid, theta * RAD2DEG);
1394 interpweights(theta_itw, theta_gp);
1395
1396 for (Index find = 0; find < nf; find++) {
1397 Fou_int = 0.;
1398 for (Index Tind = 0; Tind < nTin; Tind++) {
1399 // perform the scat angle interpolation
1400 for (Index nst = 0; nst < npha; nst++)
1401 interp(
1402 pha_mat_angint(joker, nst),
1403 theta_itw,
1404 ssd.pha_mat_data(find + f_start, Tind, joker, 0, 0, 0, nst),
1405 theta_gp);
1406 for (Index iaa = 0; iaa < naa_totran; iaa++) {
1407 // convert from scat to lab frame
1408 pha_mat_labCalc(pha_mat,
1409 pha_mat_angint(iaa, joker),
1410 pdir_array[pdir],
1411 aa_grid[iaa],
1412 idir_array[idir],
1413 0.,
1414 theta[iaa]);
1415 // and sum up/integrate
1416 // FIXME: can the integration probably be done in lab frame?
1417 // test!
1418 if (iaa == 0 || iaa == naa_totran - 1)
1419 pha_mat *= (daa_totran / 2.);
1420 else
1421 pha_mat *= daa_totran;
1422 Fou_int(Tind, joker, joker) += pha_mat;
1423 }
1424 }
1425
1426 if (this_T_interp_order <
1427 0) // T only needs to be sorted into pha_mat_fou.
1428 {
1429 for (Index Tind = 0; Tind < nTout; Tind++)
1430 pha_mat_fou(find, Tind, pdir, idir, joker, joker, 0) =
1431 Fou_int(0, joker, joker);
1432 } else {
1433 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1434 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1435 for (Index im = 0; im < nmodes; im++)
1436 reinterp(pha_mat_fou(find, joker, pdir, idir, ist1, ist2, im),
1437 Fou_int(joker, ist1, ist2),
1438 T_itw_lag,
1439 T_lag);
1440 }
1441 }
1442 }
1443 } else // derive Fourier component(s) on scattering elements own za_grid (using
1444 // given aa_grid), then 2D-interpolate inc and sca polar directions.
1445 // Do calcs only for required pha_mat components depending on stokes_dim.
1446 {
1447 Index nza = ssd.za_grid.nelem();
1448 Index naa = ssd.aa_grid.nelem();
1449 ConstVectorView za_datagrid = ssd.za_grid;
1450 ConstVectorView aa_datagrid = ssd.aa_grid;
1451 ARTS_ASSERT(aa_datagrid[0] == 0.);
1452 ARTS_ASSERT(aa_datagrid[naa - 1] == 180.);
1453 Vector daa(naa);
1454
1455 // Precalculate azimuth integration weights for this azimuthally randomly
1456 // oriented scat element (need to do this per scat element as ssd.aa_grid is
1457 // scat element specific (might change between elements) and need to do this
1458 // on actual grid instead of grid number since the grid can, at least
1459 // theoretically be non-equidistant).
1460 daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1461 for (Index iaa = 1; iaa < naa - 1; iaa++)
1462 daa[iaa] = (aa_datagrid[iaa + 1] - aa_datagrid[iaa - 1]) / 360.;
1463 daa[naa - 1] = (aa_datagrid[naa - 1] - aa_datagrid[naa - 2]) / 360.;
1464
1465 // Precalculate polar angle interpolation grid positions and weights
1466 ArrayOfGridPos pdir_za_gp, idir_za_gp;
1467 gridpos(pdir_za_gp, za_datagrid, pdir_array);
1468 gridpos(idir_za_gp, za_datagrid, idir_array);
1469 Tensor3 dir_itw(npDir, niDir, 4);
1470 interpweights(dir_itw, pdir_za_gp, idir_za_gp);
1471
1472 Tensor4 Fou_ssd(nza, nza, stokes_dim, stokes_dim);
1473 Tensor5 Fou_angint(nTin, npDir, niDir, stokes_dim, stokes_dim);
1474
1475 for (Index find = 0; find < nf; find++) {
1476 for (Index Tind = 0; Tind < nTin; Tind++) {
1477 // first, extract the phase matrix at the scatt elements own polar angle
1478 // grid and integrate over azimuth deriving their respective azimuthal
1479 // (Fourier series) 0-mode
1480 Fou_ssd = 0.;
1481 for (Index iza = 0; iza < nza; iza++)
1482 for (Index sza = 0; sza < nza; sza++) {
1483 for (Index iaa = 0; iaa < naa; iaa++) {
1484 // FIXME: are there any possible shortcuts for specific stokes
1485 // components (specifically, some where F0(ist1,ist2)==0?)
1486 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1487 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1488 Fou_ssd(sza, iza, ist1, ist2) +=
1489 daa[iaa] * ssd.pha_mat_data(find + f_start,
1490 Tind,
1491 sza,
1492 iaa,
1493 iza,
1494 0,
1495 ist1 * 4 + ist2);
1496 }
1497 }
1498
1499 // second, interpolate the extracted azimuthal mode to the stream directions
1500 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1501 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1502 // FIXME: do we need to apply any sign changes?
1503 interp(Fou_angint(Tind, joker, joker, ist1, ist2),
1504 dir_itw,
1505 Fou_ssd(joker, joker, ist1, ist2),
1506 pdir_za_gp,
1507 idir_za_gp);
1508 }
1509
1510 if (this_T_interp_order <
1511 0) // T only needs to be sorted into pha_mat_fou.
1512 {
1513 for (Index Tind = 0; Tind < nTout; Tind++)
1514 pha_mat_fou(find, Tind, joker, joker, joker, joker, 0) =
1515 Fou_angint(0, joker, joker, joker, joker);
1516 } else {
1517 for (Index pdir = 0; pdir < npDir; pdir++)
1518 for (Index idir = 0; idir < niDir; idir++)
1519 for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1520 for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1521 for (Index im = 0; im < nmodes; im++)
1522 reinterp(pha_mat_fou(find, joker, pdir, idir, ist1, ist2, im),
1523 Fou_angint(joker, pdir, idir, ist1, ist2),
1524 T_itw_lag,
1525 T_lag);
1526 }
1527 }
1528 }
1529}
1530
1532
1553void abs_vecTransform( //Output and Input
1554 StokesVector& abs_vec_lab,
1555 //Input
1556 ConstTensor3View abs_vec_data,
1557 ConstVectorView za_datagrid,
1558 ConstVectorView aa_datagrid _U_,
1559 const PType& ptype,
1560 const Numeric& za_sca _U_,
1561 const Numeric& aa_sca _U_,
1562 const Verbosity& verbosity) {
1563 const Index stokes_dim = abs_vec_lab.StokesDimensions();
1564 ARTS_ASSERT(abs_vec_lab.NumberOfFrequencies() == 1);
1565
1566 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1567 "The dimension of the stokes vector \n"
1568 "must be 1,2,3 or 4");
1569
1570 switch (ptype) {
1571 case PTYPE_GENERAL: {
1572 /*
1573 TO ANY DEVELOPER:
1574 current usage of coordinate systems in scattering solvers (RT and SSD
1575 extraction) and general radiative transfer is not consistent. Not an
1576 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1577 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1578 adding PTYPE_GENERAL support (see AUG appendix for more info).
1579 */
1580
1582 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1583 break;
1584 }
1585 case PTYPE_TOTAL_RND: {
1586 // The first element of the vector corresponds to the absorption
1587 // coefficient which is stored in the database, the others are 0.
1588
1589 abs_vec_lab.SetZero();
1590
1591 abs_vec_lab.Kjj()[0] = abs_vec_data(0, 0, 0);
1592 break;
1593 }
1594
1595 case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1596 {
1597 ARTS_ASSERT(abs_vec_data.ncols() == 2);
1598
1599 // In the case of azimuthally randomly oriented particles, only the first
1600 // two elements of the absorption coefficient vector are non-zero.
1601 // These values are dependent on the zenith angle of propagation.
1602
1603 // 1st interpolate data by za_sca
1604 GridPos gp;
1605 Vector itw(2);
1606
1607 gridpos(gp, za_datagrid, za_sca);
1608 interpweights(itw, gp);
1609
1610 abs_vec_lab.SetZero();
1611
1612 abs_vec_lab.Kjj()[0] = interp(itw, abs_vec_data(Range(joker), 0, 0), gp);
1613
1614 if (stokes_dim == 1) {
1615 break;
1616 }
1617 abs_vec_lab.K12()[0] = interp(itw, abs_vec_data(Range(joker), 0, 1), gp);
1618 break;
1619 }
1620 default: {
1622 out0 << "Not all ptype cases are implemented\n";
1623 }
1624 }
1625}
1626
1628
1649void ext_matTransform( //Output and Input
1650 PropagationMatrix& ext_mat_lab,
1651 //Input
1652 ConstTensor3View ext_mat_data,
1653 ConstVectorView za_datagrid,
1654 ConstVectorView aa_datagrid _U_,
1655 const PType& ptype,
1656 const Numeric& za_sca,
1657 const Numeric& aa_sca _U_,
1658 const Verbosity& verbosity) {
1659 const Index stokes_dim = ext_mat_lab.StokesDimensions();
1660 ARTS_ASSERT(ext_mat_lab.NumberOfFrequencies() == 1);
1661
1662 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1663 "The dimension of the stokes vector \n"
1664 "must be 1,2,3 or 4");
1665
1666 switch (ptype) {
1667 case PTYPE_GENERAL: {
1668 /*
1669 TO ANY DEVELOPER:
1670 current usage of coordinate systems in scattering solvers (RT and SSD
1671 extraction) and general radiative transfer is not consistent. Not an
1672 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1673 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1674 adding PTYPE_GENERAL support (see AUG appendix for more info).
1675 */
1676
1678 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1679 break;
1680 }
1681 case PTYPE_TOTAL_RND: {
1682 ARTS_ASSERT(ext_mat_data.ncols() == 1);
1683
1684 // In the case of randomly oriented particles the extinction matrix is
1685 // diagonal. The value of each element of the diagonal is the
1686 // extinction cross section, which is stored in the database.
1687
1688 ext_mat_lab.SetZero();
1689
1690 ext_mat_lab.Kjj() = ext_mat_data(0, 0, 0);
1691 break;
1692 }
1693
1694 case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1695 {
1696 ARTS_ASSERT(ext_mat_data.ncols() == 3);
1697
1698 // In the case of azimuthally randomly oriented particles, the extinction
1699 // matrix has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
1700 // These values are dependent on the zenith angle of propagation.
1701
1702 // 1st interpolate data by za_sca
1703 GridPos gp;
1704 Vector itw(2);
1705 Numeric Kjj;
1706 Numeric K12;
1707 Numeric K34;
1708
1709 gridpos(gp, za_datagrid, za_sca);
1710 interpweights(itw, gp);
1711
1712 ext_mat_lab.SetZero();
1713
1714 Kjj = interp(itw, ext_mat_data(Range(joker), 0, 0), gp);
1715 ext_mat_lab.Kjj() = Kjj;
1716
1717 if (stokes_dim < 2) {
1718 break;
1719 }
1720
1721 K12 = interp(itw, ext_mat_data(Range(joker), 0, 1), gp);
1722 ext_mat_lab.K12()[0] = K12;
1723
1724 if (stokes_dim < 4) {
1725 break;
1726 }
1727
1728 K34 = interp(itw, ext_mat_data(Range(joker), 0, 2), gp);
1729 ext_mat_lab.K34()[0] = K34;
1730 break;
1731 }
1732 default: {
1734 out0 << "Not all ptype cases are implemented\n";
1735 }
1736 }
1737}
1738
1740
1767void pha_matTransform( //Output
1768 MatrixView pha_mat_lab,
1769 //Input
1770 ConstTensor5View pha_mat_data,
1771 ConstVectorView za_datagrid,
1772 ConstVectorView aa_datagrid,
1773 const PType& ptype,
1774 const Index& za_sca_idx,
1775 const Index& aa_sca_idx,
1776 const Index& za_inc_idx,
1777 const Index& aa_inc_idx,
1778 ConstVectorView za_grid,
1779 ConstVectorView aa_grid,
1780 const Verbosity& verbosity) {
1781 const Index stokes_dim = pha_mat_lab.ncols();
1782
1783 Numeric za_sca = za_grid[za_sca_idx];
1784 Numeric aa_sca = aa_grid[aa_sca_idx];
1785 Numeric za_inc = za_grid[za_inc_idx];
1786 Numeric aa_inc = aa_grid[aa_inc_idx];
1787
1788 ARTS_USER_ERROR_IF (stokes_dim > 4 || stokes_dim < 1,
1789 "The dimension of the stokes vector \n"
1790 "must be 1,2,3 or 4");
1791
1792 switch (ptype) {
1793 case PTYPE_GENERAL: {
1794 /*
1795 TO ANY DEVELOPER:
1796 current usage of coordinate systems in scattering solvers (RT and SSD
1797 extraction) and general radiative transfer is not consistent. Not an
1798 issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1799 but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1800 adding PTYPE_GENERAL support (see AUG appendix for more info).
1801 */
1802
1804 out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1805 break;
1806 }
1807 case PTYPE_TOTAL_RND: {
1808 // Calculate the scattering angle and interpolate the data in it:
1809 Numeric theta_rad = scat_angle(za_sca, aa_sca, za_inc, aa_inc);
1810 const Numeric theta = RAD2DEG * theta_rad;
1811
1812 // Interpolation of the data on the scattering angle:
1813 Vector pha_mat_int(6);
1814 interpolate_scat_angle(pha_mat_int, pha_mat_data, za_datagrid, theta);
1815
1816 // Calculate the phase matrix in the laboratory frame:
1818 pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1819
1820 break;
1821 }
1822
1823 case PTYPE_AZIMUTH_RND: //Added by Cory Davis
1824 //Data is already stored in the laboratory frame,
1825 //but it is compressed a little. Details elsewhere.
1826 {
1827 ARTS_ASSERT(pha_mat_data.ncols() == 16);
1828 ARTS_ASSERT(pha_mat_data.npages() == za_datagrid.nelem());
1829 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
1830 (aa_sca - aa_inc > 180) *
1831 360; //delta_aa corresponds to the "books"
1832 //dimension of pha_mat_data
1833 GridPos za_sca_gp;
1834 GridPos delta_aa_gp;
1835 GridPos za_inc_gp;
1836 Vector itw(8);
1837
1838 gridpos(delta_aa_gp, aa_datagrid, abs(delta_aa));
1839 gridpos(za_inc_gp, za_datagrid, za_inc);
1840 gridpos(za_sca_gp, za_datagrid, za_sca);
1841
1842 interpweights(itw, za_sca_gp, delta_aa_gp, za_inc_gp);
1843
1844 pha_mat_lab(0, 0) =
1845 interp(itw,
1846 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 0),
1847 za_sca_gp,
1848 delta_aa_gp,
1849 za_inc_gp);
1850 if (stokes_dim == 1) {
1851 break;
1852 }
1853 pha_mat_lab(0, 1) =
1854 interp(itw,
1855 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 1),
1856 za_sca_gp,
1857 delta_aa_gp,
1858 za_inc_gp);
1859 pha_mat_lab(1, 0) =
1860 interp(itw,
1861 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 4),
1862 za_sca_gp,
1863 delta_aa_gp,
1864 za_inc_gp);
1865 pha_mat_lab(1, 1) =
1866 interp(itw,
1867 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 5),
1868 za_sca_gp,
1869 delta_aa_gp,
1870 za_inc_gp);
1871 if (stokes_dim == 2) {
1872 break;
1873 }
1874 if (delta_aa >= 0) {
1875 pha_mat_lab(0, 2) = interp(
1876 itw,
1877 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1878 za_sca_gp,
1879 delta_aa_gp,
1880 za_inc_gp);
1881 pha_mat_lab(1, 2) = interp(
1882 itw,
1883 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1884 za_sca_gp,
1885 delta_aa_gp,
1886 za_inc_gp);
1887 pha_mat_lab(2, 0) = interp(
1888 itw,
1889 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1890 za_sca_gp,
1891 delta_aa_gp,
1892 za_inc_gp);
1893 pha_mat_lab(2, 1) = interp(
1894 itw,
1895 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1896 za_sca_gp,
1897 delta_aa_gp,
1898 za_inc_gp);
1899 } else {
1900 pha_mat_lab(0, 2) = -interp(
1901 itw,
1902 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1903 za_sca_gp,
1904 delta_aa_gp,
1905 za_inc_gp);
1906 pha_mat_lab(1, 2) = -interp(
1907 itw,
1908 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1909 za_sca_gp,
1910 delta_aa_gp,
1911 za_inc_gp);
1912 pha_mat_lab(2, 0) = -interp(
1913 itw,
1914 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1915 za_sca_gp,
1916 delta_aa_gp,
1917 za_inc_gp);
1918 pha_mat_lab(2, 1) = -interp(
1919 itw,
1920 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1921 za_sca_gp,
1922 delta_aa_gp,
1923 za_inc_gp);
1924 }
1925 pha_mat_lab(2, 2) = interp(
1926 itw,
1927 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 10),
1928 za_sca_gp,
1929 delta_aa_gp,
1930 za_inc_gp);
1931 if (stokes_dim == 3) {
1932 break;
1933 }
1934 if (delta_aa >= 0) {
1935 pha_mat_lab(0, 3) = interp(
1936 itw,
1937 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1938 za_sca_gp,
1939 delta_aa_gp,
1940 za_inc_gp);
1941 pha_mat_lab(1, 3) = interp(
1942 itw,
1943 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1944 za_sca_gp,
1945 delta_aa_gp,
1946 za_inc_gp);
1947 pha_mat_lab(3, 0) = interp(
1948 itw,
1949 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1950 za_sca_gp,
1951 delta_aa_gp,
1952 za_inc_gp);
1953 pha_mat_lab(3, 1) = interp(
1954 itw,
1955 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1956 za_sca_gp,
1957 delta_aa_gp,
1958 za_inc_gp);
1959 } else {
1960 pha_mat_lab(0, 3) = -interp(
1961 itw,
1962 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1963 za_sca_gp,
1964 delta_aa_gp,
1965 za_inc_gp);
1966 pha_mat_lab(1, 3) = -interp(
1967 itw,
1968 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1969 za_sca_gp,
1970 delta_aa_gp,
1971 za_inc_gp);
1972 pha_mat_lab(3, 0) = -interp(
1973 itw,
1974 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1975 za_sca_gp,
1976 delta_aa_gp,
1977 za_inc_gp);
1978 pha_mat_lab(3, 1) = -interp(
1979 itw,
1980 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1981 za_sca_gp,
1982 delta_aa_gp,
1983 za_inc_gp);
1984 }
1985 pha_mat_lab(2, 3) = interp(
1986 itw,
1987 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 11),
1988 za_sca_gp,
1989 delta_aa_gp,
1990 za_inc_gp);
1991 pha_mat_lab(3, 2) = interp(
1992 itw,
1993 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 14),
1994 za_sca_gp,
1995 delta_aa_gp,
1996 za_inc_gp);
1997 pha_mat_lab(3, 3) = interp(
1998 itw,
1999 pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 15),
2000 za_sca_gp,
2001 delta_aa_gp,
2002 za_inc_gp);
2003 break;
2004 }
2005
2006 default: {
2008 out0 << "Not all ptype cases are implemented\n";
2009 }
2010 }
2011}
2012
2014
2043 MatrixView ext_mat,
2044 //Input
2045 ConstVectorView abs_vec,
2046 const Index& stokes_dim) {
2047 ARTS_ASSERT(stokes_dim >= 1 && stokes_dim <= 4);
2048 ARTS_ASSERT(ext_mat.nrows() == stokes_dim);
2049 ARTS_ASSERT(ext_mat.ncols() == stokes_dim);
2050 ARTS_ASSERT(abs_vec.nelem() == stokes_dim);
2051
2052 // first: diagonal elements
2053 for (Index is = 0; is < stokes_dim; is++) {
2054 ext_mat(is, is) += abs_vec[0];
2055 }
2056 // second: off-diagonal elements, namely first row and column
2057 for (Index is = 1; is < stokes_dim; is++) {
2058 ext_mat(0, is) += abs_vec[is];
2059 ext_mat(is, 0) += abs_vec[is];
2060 }
2061}
2062
2064
2078 const Numeric& aa_sca,
2079 const Numeric& za_inc,
2080 const Numeric& aa_inc) {
2081 Numeric theta_rad;
2082 Numeric ANG_TOL = 1e-7;
2083
2084 // CPD 5/10/03.
2085 // Two special cases are implemented here to avoid NaNs that can sometimes
2086 // occur in in the acos() formula in forward and backscattering cases.
2087 //
2088 // GH 2011-05-31
2089 // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
2090
2091 if ((abs(aa_sca - aa_inc) < ANG_TOL) ||
2092 (abs(abs(aa_sca - aa_inc) - 360) < ANG_TOL)) {
2093 theta_rad = Conversion::deg2rad(abs(za_sca - za_inc));
2094 } else if (abs(abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
2095 theta_rad = Conversion::deg2rad(za_sca + za_inc);
2096 if (theta_rad > PI) {
2097 theta_rad = 2 * PI - theta_rad;
2098 }
2099 } else {theta_rad =
2100 acos(Conversion::cosd(za_sca) * Conversion::cosd(za_inc) +
2101 Conversion::sind(za_sca) * Conversion::sind(za_inc) *
2102 Conversion::cosd(aa_sca - aa_inc));
2103 }
2104 return theta_rad;
2105}
2106
2108
2132 VectorView pha_mat_int,
2133 //Input:
2134 ConstTensor5View pha_mat_data,
2135 ConstVectorView za_datagrid,
2136 const Numeric theta) {
2137 GridPos thet_gp;
2138 gridpos(thet_gp, za_datagrid, theta);
2139
2140 Vector itw(2);
2141 interpweights(itw, thet_gp);
2142
2143 ARTS_ASSERT(pha_mat_data.ncols() == 6);
2144 for (Index i = 0; i < 6; i++) {
2145 pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), thet_gp);
2146 }
2147}
2148
2150
2175void pha_mat_labCalc( //Output:
2176 MatrixView pha_mat_lab,
2177 //Input:
2178 ConstVectorView pha_mat_int,
2179 const Numeric& za_sca,
2180 const Numeric& aa_sca,
2181 const Numeric& za_inc,
2182 const Numeric& aa_inc,
2183 const Numeric& theta_rad) {
2184 const Index stokes_dim = pha_mat_lab.ncols();
2185
2187 "NaN value(s) detected in *pha_mat_labCalc* (0,0). Could the "
2188 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2189 "input data are OK and you critically need the ongoing calculations, "
2190 "try to change the observation LOS slightly. If you can reproduce "
2191 "this error, please contact Patrick in order to help tracking down "
2192 "the reason to this problem. If you see this message occasionally "
2193 "when doing MC calculations, it should not be critical. This path "
2194 "sampling will be rejected and replaced with a new one.");
2195
2196 // For stokes_dim = 1, we only need Z11=F11:
2197 pha_mat_lab(0, 0) = F11;
2198
2199 if (stokes_dim > 1) {
2200 Numeric za_sca_rad = Conversion::deg2rad(za_sca);
2201 Numeric za_inc_rad = Conversion::deg2rad(za_inc);
2202 Numeric aa_sca_rad = Conversion::deg2rad(aa_sca);
2203 Numeric aa_inc_rad = Conversion::deg2rad(aa_inc);
2204
2205 const Numeric ANGTOL_RAD = 1e-6; //CPD: this constant is used to adjust
2206 //zenith angles close to 0 and PI. This is
2207 //also used to avoid float == float statements.
2208
2209 //
2210 // Several cases have to be considered:
2211 //
2212
2213 if ((abs(theta_rad) < ANGTOL_RAD) // forward scattering
2214 || (abs(theta_rad - Constant::pi) < ANGTOL_RAD) // backward scattering
2215 ||
2216 (abs(aa_inc_rad - aa_sca_rad) < ANGTOL_RAD) // inc and sca on meridian
2217 || (abs(abs(aa_inc_rad - aa_sca_rad) - Constant::two_pi) < ANGTOL_RAD) // "
2218 || (abs(abs(aa_inc_rad - aa_sca_rad) - Constant::pi) < ANGTOL_RAD) // "
2219 ) {
2220 pha_mat_lab(0, 1) = F12;
2221 pha_mat_lab(1, 0) = F12;
2222 pha_mat_lab(1, 1) = F22;
2223
2224 if (stokes_dim > 2) {
2225 pha_mat_lab(0, 2) = 0;
2226 pha_mat_lab(1, 2) = 0;
2227 pha_mat_lab(2, 0) = 0;
2228 pha_mat_lab(2, 1) = 0;
2229 pha_mat_lab(2, 2) = F33;
2230
2231 if (stokes_dim > 3) {
2232 pha_mat_lab(0, 3) = 0;
2233 pha_mat_lab(1, 3) = 0;
2234 pha_mat_lab(2, 3) = F34;
2235 pha_mat_lab(3, 0) = 0;
2236 pha_mat_lab(3, 1) = 0;
2237 pha_mat_lab(3, 2) = -F34;
2238 pha_mat_lab(3, 3) = F44;
2239 }
2240 }
2241 }
2242
2243 else {
2244 Numeric sigma1;
2245 Numeric sigma2;
2246
2247 Numeric s1, s2;
2248
2249 // In these cases we have to take limiting values.
2250
2251 if (za_inc_rad < ANGTOL_RAD) {
2252 sigma1 = pi + aa_sca_rad - aa_inc_rad;
2253 sigma2 = 0;
2254 } else if (za_inc_rad > pi - ANGTOL_RAD) {
2255 sigma1 = aa_sca_rad - aa_inc_rad;
2256 sigma2 = pi;
2257 } else if (za_sca_rad < ANGTOL_RAD) {
2258 sigma1 = 0;
2259 sigma2 = pi + aa_sca_rad - aa_inc_rad;
2260 } else if (za_sca_rad > pi - ANGTOL_RAD) {
2261 sigma1 = pi;
2262 sigma2 = aa_sca_rad - aa_inc_rad;
2263 } else {
2264 s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad)) /
2265 (sin(za_inc_rad) * sin(theta_rad));
2266 s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos(theta_rad)) /
2267 (sin(za_sca_rad) * sin(theta_rad));
2268
2269 sigma1 = acos(s1);
2270 sigma2 = acos(s2);
2271
2272 // Arccos is only defined in the range from -1 ... 1
2273 // Numerical problems can appear for values close to 1 or -1
2274 // this (also) catches the case when inc and sca are on one meridian
2275 if (std::isnan(sigma1) || std::isnan(sigma2)) {
2276 if (abs(s1 - 1) < ANGTOL_RAD) sigma1 = 0;
2277 if (abs(s1 + 1) < ANGTOL_RAD) sigma1 = pi;
2278 if (abs(s2 - 1) < ANGTOL_RAD) sigma2 = 0;
2279 if (abs(s2 + 1) < ANGTOL_RAD) sigma2 = pi;
2280 }
2281 }
2282
2283 const Numeric C1 = cos(2 * sigma1);
2284 const Numeric C2 = cos(2 * sigma2);
2285
2286 const Numeric S1 = sin(2 * sigma1);
2287 const Numeric S2 = sin(2 * sigma2);
2288
2289 pha_mat_lab(0, 1) = C1 * F12;
2290 pha_mat_lab(1, 0) = C2 * F12;
2291 pha_mat_lab(1, 1) = C1 * C2 * F22 - S1 * S2 * F33;
2292
2293 //ARTS_ASSERT(!std::isnan(pha_mat_lab(0,1)));
2294 //ARTS_ASSERT(!std::isnan(pha_mat_lab(1,0)));
2295 //ARTS_ASSERT(!std::isnan(pha_mat_lab(1,1)));
2296 ARTS_USER_ERROR_IF (std::isnan(pha_mat_lab(0, 1)) || std::isnan(pha_mat_lab(1, 0)) ||
2297 std::isnan(pha_mat_lab(1, 1)),
2298 "NaN value(s) detected in *pha_mat_labCalc* (0/1,1). Could the "
2299 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2300 "input data are OK and you critically need the ongoing calculations, "
2301 "try to change the observation LOS slightly. If you can reproduce "
2302 "this error, please contact Patrick in order to help tracking down "
2303 "the reason to this problem. If you see this message occasionally "
2304 "when doing MC calculations, it should not be critical. This path "
2305 "sampling will be rejected and replaced with a new one.");
2306
2307 if (stokes_dim > 2) {
2308 /*CPD: For skokes_dim > 2 some of the transformation formula
2309 for each element have a different sign depending on whether or
2310 not 0<aa_scat-aa_inc<180. For details see pages 94 and 95 of
2311 Mishchenkos chapter in :
2312 Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic
2313 scattering by nonspherical particles. In Exploring the Atmosphere
2314 by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag,
2315 Berlin, pp. 77-127.
2316 This is available at http://www.giss.nasa.gov/~crmim/publications/ */
2317 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
2318 (aa_sca - aa_inc > 180) * 360;
2319 if (delta_aa >= 0) {
2320 pha_mat_lab(0, 2) = S1 * F12;
2321 pha_mat_lab(1, 2) = S1 * C2 * F22 + C1 * S2 * F33;
2322 pha_mat_lab(2, 0) = -S2 * F12;
2323 pha_mat_lab(2, 1) = -C1 * S2 * F22 - S1 * C2 * F33;
2324 } else {
2325 pha_mat_lab(0, 2) = -S1 * F12;
2326 pha_mat_lab(1, 2) = -S1 * C2 * F22 - C1 * S2 * F33;
2327 pha_mat_lab(2, 0) = S2 * F12;
2328 pha_mat_lab(2, 1) = C1 * S2 * F22 + S1 * C2 * F33;
2329 }
2330 pha_mat_lab(2, 2) = -S1 * S2 * F22 + C1 * C2 * F33;
2331
2332 if (stokes_dim > 3) {
2333 if (delta_aa >= 0) {
2334 pha_mat_lab(1, 3) = S2 * F34;
2335 pha_mat_lab(3, 1) = S1 * F34;
2336 } else {
2337 pha_mat_lab(1, 3) = -S2 * F34;
2338 pha_mat_lab(3, 1) = -S1 * F34;
2339 }
2340 pha_mat_lab(0, 3) = 0;
2341 pha_mat_lab(2, 3) = C2 * F34;
2342 pha_mat_lab(3, 0) = 0;
2343 pha_mat_lab(3, 2) = -C1 * F34;
2344 pha_mat_lab(3, 3) = F44;
2345 }
2346 }
2347 }
2348 }
2349}
2350
2351ostream& operator<<(ostream& os, const SingleScatteringData& /*ssd*/) {
2352 os << "SingleScatteringData: Output operator not implemented";
2353 return os;
2354}
2355
2356ostream& operator<<(ostream& os, const ScatteringMetaData& /*ssd*/) {
2357 os << "ScatteringMetaData: Output operator not implemented";
2358 return os;
2359}
2360
2362
2378 PropagationMatrix& ext_mat,
2379 StokesVector& abs_vec,
2380 //Input:
2381 const PropagationMatrix& propmat_clearsky) {
2382 const Index stokes_dim = propmat_clearsky.StokesDimensions();
2383 const Index freq_dim = propmat_clearsky.NumberOfFrequencies();
2384
2385 // old abs_vecInit
2386 abs_vec = StokesVector(freq_dim, stokes_dim);
2387 abs_vec.SetZero();
2388
2389 ext_mat = PropagationMatrix(freq_dim, stokes_dim);
2390 ext_mat.SetZero();
2391
2392 // old ext_matAddGas and abs_vecAddGas for 0 vector and matrix
2393 abs_vec += propmat_clearsky;
2394 ext_mat += propmat_clearsky;
2395}
2396
2398
2408PType PTypeFromString(const String& ptype_string) {
2409 PType ptype;
2410 if (ptype_string == "general")
2411 ptype = PTYPE_GENERAL;
2412 else if (ptype_string == "totally_random")
2413 ptype = PTYPE_TOTAL_RND;
2414 else if (ptype_string == "azimuthally_random")
2415 ptype = PTYPE_AZIMUTH_RND;
2416 else {
2418 "Unknown ptype: ", ptype_string, "\n"
2419 "Valid types are: general, totally_random and "
2420 "azimuthally_random.")
2421 }
2422
2423 return ptype;
2424}
2425
2427
2437PType PType2FromString(const String& ptype_string) {
2438 PType ptype;
2439 if (ptype_string == "general")
2440 ptype = PTYPE_GENERAL;
2441 else if (ptype_string == "macroscopically_isotropic")
2442 ptype = PTYPE_TOTAL_RND;
2443 else if (ptype_string == "horizontally_aligned")
2444 ptype = PTYPE_AZIMUTH_RND;
2445 else {
2447 "Unknown ptype: ", ptype_string, "\n"
2448 "Valid types are: general, macroscopically_isotropic and "
2449 "horizontally_aligned.")
2450 }
2451
2452 return ptype;
2453}
2454
2456
2465 String ptype_string;
2466
2467 switch (ptype) {
2468 case PTYPE_GENERAL:
2469 ptype_string = "general";
2470 break;
2471 case PTYPE_TOTAL_RND:
2472 ptype_string = "totally_random";
2473 break;
2474 case PTYPE_AZIMUTH_RND:
2475 ptype_string = "azimuthally_random";
2476 break;
2477 default:
2479 "Internal error: Cannot map PType enum value ",
2480 ptype, " to String.")
2481 break;
2482 }
2483
2484 return ptype_string;
2485}
2486
2488
2496 // First check that input fulfills requirements on older data formats:
2497 // 1) Is za_grid symmetric and includes 90deg?
2498 Index nza = ssd.za_grid.nelem();
2499 for (Index i = 0; i < nza / 2; i++) {
2501 180. - ssd.za_grid[nza - 1 - i], ssd.za_grid[i], 2 * DBL_EPSILON),
2502 "Zenith grid of azimuthally_random single scattering data\n"
2503 "is not symmetric with respect to 90degree.")
2504 }
2505 ARTS_USER_ERROR_IF (!is_same_within_epsilon(ssd.za_grid[nza / 2], 90., 2 * DBL_EPSILON),
2506 "Zenith grid of azimuthally_random single scattering data\n"
2507 "does not contain 90 degree grid point.")
2508
2509 // 2) Are data sizes correct?
2510 ostringstream os_pha_mat;
2511 os_pha_mat << "pha_mat ";
2512 ostringstream os_ext_mat;
2513 os_ext_mat << "ext_mat ";
2514 ostringstream os_abs_vec;
2515 os_abs_vec << "abs_vec ";
2516 chk_size(os_pha_mat.str(),
2517 ssd.pha_mat_data,
2518 ssd.f_grid.nelem(),
2519 ssd.T_grid.nelem(),
2520 ssd.za_grid.nelem(),
2521 ssd.aa_grid.nelem(),
2522 ssd.za_grid.nelem() / 2 + 1,
2523 1,
2524 16);
2525
2526 chk_size(os_ext_mat.str(),
2527 ssd.ext_mat_data,
2528 ssd.f_grid.nelem(),
2529 ssd.T_grid.nelem(),
2530 ssd.za_grid.nelem() / 2 + 1,
2531 1,
2532 3);
2533
2534 chk_size(os_abs_vec.str(),
2535 ssd.abs_vec_data,
2536 ssd.f_grid.nelem(),
2537 ssd.T_grid.nelem(),
2538 ssd.za_grid.nelem() / 2 + 1,
2539 1,
2540 2);
2541
2542 // Now that we are sure that za_grid is properly symmetric, we just need to
2543 // copy over the data (ie no interpolation).
2544 Tensor5 tmpT5 = ssd.abs_vec_data;
2545 ssd.abs_vec_data.resize(tmpT5.nshelves(),
2546 tmpT5.nbooks(),
2547 ssd.za_grid.nelem(),
2548 tmpT5.nrows(),
2549 tmpT5.ncols());
2550 ssd.abs_vec_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2551 for (Index i = 0; i < nza / 2; i++) {
2552 ssd.abs_vec_data(joker, joker, nza - 1 - i, joker, joker) =
2553 tmpT5(joker, joker, i, joker, joker);
2554 }
2555
2556 tmpT5 = ssd.ext_mat_data;
2557 ssd.ext_mat_data.resize(tmpT5.nshelves(),
2558 tmpT5.nbooks(),
2559 ssd.za_grid.nelem(),
2560 tmpT5.nrows(),
2561 tmpT5.ncols());
2562 ssd.ext_mat_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2563 for (Index i = 0; i < nza / 2; i++) {
2564 ssd.ext_mat_data(joker, joker, nza - 1 - i, joker, joker) =
2565 tmpT5(joker, joker, i, joker, joker);
2566 }
2567
2568 Tensor7 tmpT7 = ssd.pha_mat_data;
2569 ssd.pha_mat_data.resize(tmpT7.nlibraries(),
2570 tmpT7.nvitrines(),
2571 tmpT7.nshelves(),
2572 tmpT7.nbooks(),
2573 ssd.za_grid.nelem(),
2574 tmpT7.nrows(),
2575 tmpT7.ncols());
2576 ssd.pha_mat_data(
2577 joker, joker, joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT7;
2578
2579 // scatt. matrix elements 13,23,31,32 and 14,24,41,42 (=elements 2,6,8,9 and
2580 // 3,7,12,13 in ARTS' flattened format, respectively) change sign.
2581 tmpT7(joker, joker, joker, joker, joker, joker, Range(2, 2)) *= -1.;
2582 tmpT7(joker, joker, joker, joker, joker, joker, Range(6, 4)) *= -1.;
2583 tmpT7(joker, joker, joker, joker, joker, joker, Range(12, 2)) *= -1.;
2584
2585 // For second half of incident polar angles (>90deg), we need to mirror the
2586 // original data in both incident and scattered polar angle around 90deg "planes".
2587 for (Index i = 0; i < nza / 2; i++)
2588 for (Index j = 0; j < nza; j++)
2589 ssd.pha_mat_data(
2590 joker, joker, nza - 1 - j, joker, nza - 1 - i, joker, joker) =
2591 tmpT7(joker, joker, j, joker, i, joker, joker);
2592}
2593
2595
2604 const String& particle_ssdmethod_string) {
2605 ParticleSSDMethod particle_ssdmethod;
2606 if (particle_ssdmethod_string == "tmatrix")
2607 particle_ssdmethod = PARTICLE_SSDMETHOD_TMATRIX;
2608 else {
2610 "Unknown particle SSD method: ",
2611 particle_ssdmethod_string, "\n"
2612 "Valid methods: tmatrix")
2613 }
2614
2615 return particle_ssdmethod;
2616}
2617
2619
2627String PTypeToString(const ParticleSSDMethod& particle_ssdmethod) {
2628 String particle_ssdmethod_string;
2629
2630 switch (particle_ssdmethod) {
2632 particle_ssdmethod_string = "tmatrix";
2633 break;
2634 default:
2636 "Internal error: Cannot map ParticleSSDMethod enum value ",
2637 particle_ssdmethod, " to String.")
2638 break;
2639 }
2640
2641 return particle_ssdmethod_string;
2642}
2643
2644
2646
2674 MatrixView abs_data,
2675 Tensor3View pfun_data,
2676 const ArrayOfSingleScatteringData& scat_data,
2677 const Index& iss,
2678 ConstMatrixView pnd_data,
2679 ArrayOfIndex& cloudbox_limits,
2680 ConstVectorView T_grid,
2681 ConstVectorView sa_grid)
2682{
2683 // Sizes
2684 const Index nse = scat_data.nelem();
2685 const Index nf = scat_data[0].f_grid.nelem();
2686 [[maybe_unused]] const Index nt = T_grid.nelem();
2687 const Index nsa = sa_grid.nelem();
2688 const Index ncl = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2689
2690 ARTS_ASSERT(ext_data.nrows() == nf);
2691 ARTS_ASSERT(ext_data.ncols() == nt);
2692 ARTS_ASSERT(abs_data.nrows() == nf);
2693 ARTS_ASSERT(abs_data.ncols() == nt);
2694 ARTS_ASSERT(pfun_data.npages() == nf);
2695 ARTS_ASSERT(pfun_data.nrows() == nt);
2696 ARTS_ASSERT(pfun_data.ncols() == nsa);
2697 ARTS_ASSERT(cloudbox_limits.nelem() == 2);
2698 ARTS_ASSERT(pnd_data.nrows() == nse);
2699 ARTS_ASSERT(pnd_data.ncols() == ncl);
2700
2701 // Check that all data is TRO
2702 {
2703 bool all_totrand = true;
2704 for (Index ie = 0; ie < nse; ie++) {
2705 if (scat_data[ie].ptype != PTYPE_TOTAL_RND)
2706 all_totrand = false;
2707 }
2708 ARTS_USER_ERROR_IF (!all_totrand,
2709 "This method demands that all scat_data are TRO");
2710 }
2711
2712 // Help variables to hold non-zero data inside the cloudbox
2713 const Index cl_start = cloudbox_limits[0];
2714 Vector T_values(ncl);
2715 ArrayOfIndex cboxlayer(ncl);
2716
2717 // Loop scattering elements
2718 for (Index ie = 0; ie < nse; ie++) {
2719 // Allowed temperature range
2720 const Index last = scat_data[ie].T_grid.nelem() - 1;
2721 const Numeric tmin = 1.5*scat_data[ie].T_grid[0] -
2722 0.5*scat_data[ie].T_grid[1];
2723 const Numeric tmax = 1.5*scat_data[ie].T_grid[last] -
2724 0.5*scat_data[ie].T_grid[last-1];
2725
2726 Index nvals = 0;
2727 for(Index icl=0; icl<ncl; ++icl){
2728 // Nothing to do if PND is zero
2729 if (abs(pnd_data(ie,icl)) > 1e-3) {
2730 const Numeric Tthis = T_grid[cl_start + icl];
2731 ARTS_USER_ERROR_IF(Tthis < tmin || Tthis > tmax,
2732 "Temperature interpolation error for scattering element ", ie,
2733 " of species ", iss, "\nYour temperature: ", Tthis, " K\n"
2734 "Allowed range of temperatures: ", tmin, " - ", tmax, " K");
2735 T_values[nvals] = Tthis;
2736 cboxlayer[nvals] = icl;
2737 nvals++;
2738 }
2739 }
2740
2741 if (nvals > 0) {
2742 // Temperature-only interpolation weights
2743 ArrayOfGridPos gp_t(nvals);
2744 gridpos(gp_t, scat_data[ie].T_grid, T_values[Range(0,nvals)]);
2745 Matrix itw1(nvals, 2);
2746 interpweights(itw1, gp_t);
2747
2748 // Temperature + scattering angle interpolation weights
2749 ArrayOfGridPos gp_sa(nsa);
2750 gridpos(gp_sa, scat_data[ie].za_grid, sa_grid);
2751 Tensor3 itw2(nvals, nsa, 4);
2752 interpweights(itw2, gp_t, gp_sa);
2753
2754 // Loop frequencies
2755 for (Index iv = 0; iv < nf; iv++) {
2756 // Interpolate
2757 Vector ext1(nvals), abs1(nvals);
2758 Matrix pfu1(nvals,nsa);
2759 interp(ext1, itw1, scat_data[ie].ext_mat_data(iv,joker,0,0,0), gp_t);
2760 interp(abs1, itw1, scat_data[ie].abs_vec_data(iv,joker,0,0,0), gp_t);
2761 interp(pfu1, itw2, scat_data[ie].pha_mat_data(iv,joker,joker,0,0,0,0),
2762 gp_t, gp_sa);
2763
2764 // Add to container variables
2765 for (Index i = 0; i < nvals; i++) {
2766 const Index ic = cboxlayer[i];
2767 const Index it = cl_start + ic;
2768 ext_data(iv,it) += pnd_data(ie,ic) * ext1[i];
2769 abs_data(iv,it) += pnd_data(ie,ic) * abs1[i];
2770 for (Index ia = 0; ia < nsa; ia++) {
2771 pfun_data(iv,it,ia) += pnd_data(ie,ic) * pfu1(i,ia);
2772 }
2773 }
2774 }
2775 }
2776 }
2777}
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:241
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:906
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:421
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Tensor3.
Definition: matpackIII.h:130
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index ncols() const noexcept
Definition: matpackIV.h:142
Index nrows() const noexcept
Definition: matpackIV.h:141
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
A constant view of a Tensor5.
Definition: matpackV.h:141
Index nrows() const noexcept
Definition: matpackV.h:152
Index ncols() const noexcept
Definition: matpackV.h:153
Index npages() const noexcept
Definition: matpackV.h:151
Index nbooks() const noexcept
Definition: matpackV.h:150
Index nshelves() const noexcept
Definition: matpackV.h:149
Index nbooks() const noexcept
Definition: matpackVI.h:157
Index nvitrines() const noexcept
Definition: matpackVI.h:155
Index ncols() const noexcept
Definition: matpackVI.h:160
Index npages() const noexcept
Definition: matpackVI.h:158
Index nshelves() const noexcept
Definition: matpackVI.h:156
Index nrows() const noexcept
Definition: matpackVI.h:159
Index ncols() const noexcept
Definition: matpackVII.h:159
Index npages() const noexcept
Definition: matpackVII.h:157
Index nrows() const noexcept
Definition: matpackVII.h:158
Index nlibraries() const noexcept
Definition: matpackVII.h:153
Index nvitrines() const noexcept
Definition: matpackVII.h:154
Index nshelves() const noexcept
Definition: matpackVII.h:155
Index nbooks() const noexcept
Definition: matpackVII.h:156
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
void SetZero()
Sets all data to zero.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
The range class.
Definition: matpackI.h:159
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3View class.
Definition: matpackIII.h:246
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4View class.
Definition: matpackIV.h:292
The Tensor4 class.
Definition: matpackIV.h:429
The Tensor5View class.
Definition: matpackV.h:343
The Tensor5 class.
Definition: matpackV.h:516
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1741
The Tensor6View class.
Definition: matpackVI.h:632
The Tensor6 class.
Definition: matpackVI.h:1099
The Tensor7View class.
Definition: matpackVII.h:1303
The Tensor7 class.
Definition: matpackVII.h:2399
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5488
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
#define _U_
Definition: config.h:180
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#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 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.
#define abs(x)
#define min(a, b)
#define max(a, b)
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:354
Header file for logic.cc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:161
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
Declarations having to do with the four output streams.
#define CREATE_OUT0
Definition: messages.h:203
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric two_pi
Two times pi.
constexpr Numeric e
Elementary charge convenience name [C].
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.
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:53
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.
void FouComp_1ScatElem(Tensor7View pha_mat_fou, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Vector &pdir_array, const Vector &idir_array, const Index &f_start, const Index &t_interp_order, const Index &naa_totran)
Preparing phase matrix fourier series components for one scattering element.
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:52
@ PTYPE_GENERAL
Definition: optproperties.h:53
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
ParticleSSDMethod
An attribute to classify the method to be used for SingleScatteringData.
Definition: optproperties.h:63
@ PARTICLE_SSDMETHOD_TMATRIX
Definition: optproperties.h:65
Structure to store a grid position.
Definition: interpolation.h:73
This file contains basic functions to handle XML data files.