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