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