ARTS  2.4.0(git:4fb77825)
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"
45 #include "logic.h"
46 #include "math_funcs.h"
47 #include "matpackVII.h"
48 #include "messages.h"
49 #include "xml_io.h"
50 
51 extern const Numeric DEG2RAD;
52 extern const Numeric RAD2DEG;
53 extern const Numeric PI;
54 
55 #define F11 pha_mat_int[0]
56 #define F12 pha_mat_int[1]
57 #define F22 pha_mat_int[2]
58 #define F33 pha_mat_int[3]
59 #define F34 pha_mat_int[4]
60 #define F44 pha_mat_int[5]
61 
63 
72 /*
73 void methodname(//Output
74  type& name,
75  //Input
76  const type& name)
77 {
78 }
79 */
80 
82 
102 void opt_prop_Bulk( //Output
103  Tensor5& ext_mat, // (nf,nT,ndir,nst,nst)
104  Tensor4& abs_vec, // (nf,nT,ndir,nst)
105  Index& ptype,
106  //Input
107  const ArrayOfTensor5& ext_mat_ss, // [nss](nf,nT,ndir,nst,nst)
108  const ArrayOfTensor4& abs_vec_ss, // [nss](nf,nT,ndir,nst)
109  const ArrayOfIndex& ptypes_ss) {
110  assert(ext_mat_ss.nelem() == abs_vec_ss.nelem());
111 
112  ext_mat = ext_mat_ss[0];
113  abs_vec = abs_vec_ss[0];
114 
115  for (Index i_ss = 1; i_ss < ext_mat_ss.nelem(); i_ss++) {
116  ext_mat += ext_mat_ss[i_ss];
117  abs_vec += abs_vec_ss[i_ss];
118  }
119  ptype = max(ptypes_ss);
120 }
121 
123 
148 void opt_prop_ScatSpecBulk( //Output
149  ArrayOfTensor5& ext_mat, // [nss](nf,nT,ndir,nst,nst)
150  ArrayOfTensor4& abs_vec, // [nss](nf,nT,ndir,nst)
151  ArrayOfIndex& ptype,
152  //Input
153  const ArrayOfArrayOfTensor5& ext_mat_se, // [nss][nse](nf,nT,ndir,nst,nst)
154  const ArrayOfArrayOfTensor4& abs_vec_se, // [nss][nse](nf,nT,ndir,nst)
155  const ArrayOfArrayOfIndex& ptypes_se,
156  ConstMatrixView pnds,
157  ConstMatrixView t_ok) {
158  assert(t_ok.nrows() == pnds.nrows());
159  assert(t_ok.ncols() == pnds.ncols());
160  assert(TotalNumberOfElements(ext_mat_se) == pnds.nrows());
161  assert(TotalNumberOfElements(abs_vec_se) == pnds.nrows());
162  assert(ext_mat_se.nelem() == abs_vec_se.nelem());
163 
164  const Index nT = pnds.ncols();
165  const Index nf = abs_vec_se[0][0].nbooks();
166  const Index nDir = abs_vec_se[0][0].nrows();
167  const Index stokes_dim = abs_vec_se[0][0].ncols();
168 
169  const Index nss = ext_mat_se.nelem();
170  ext_mat.resize(nss);
171  abs_vec.resize(nss);
172  ptype.resize(nss);
173  Tensor4 ext_tmp;
174  Tensor3 abs_tmp;
175 
176  Index i_se_flat = 0;
177 
178  for (Index i_ss = 0; i_ss < nss; i_ss++) {
179  assert(ext_mat_se[i_ss].nelem() == abs_vec_se[i_ss].nelem());
180  assert(nT == ext_mat_se[i_ss][0].nbooks());
181  assert(nT == abs_vec_se[i_ss][0].npages());
182 
183  ext_mat[i_ss].resize(nf, nT, nDir, stokes_dim, stokes_dim);
184  ext_mat[i_ss] = 0.;
185  abs_vec[i_ss].resize(nf, nT, nDir, stokes_dim);
186  abs_vec[i_ss] = 0.;
187 
188  for (Index i_se = 0; i_se < ext_mat_se[i_ss].nelem(); i_se++) {
189  assert(nT == ext_mat_se[i_ss][i_se].nbooks());
190  assert(nT == abs_vec_se[i_ss][i_se].npages());
191 
192  for (Index Tind = 0; Tind < nT; Tind++) {
193  if (pnds(i_se_flat, Tind) != 0.) {
194  if (t_ok(i_se_flat, Tind) > 0.) {
195  ext_tmp = ext_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker);
196  ext_tmp *= pnds(i_se_flat, Tind);
197  ext_mat[i_ss](joker, Tind, joker, joker, joker) += ext_tmp;
198 
199  abs_tmp = abs_vec_se[i_ss][i_se](joker, Tind, joker, joker);
200  abs_tmp *= pnds(i_se_flat, Tind);
201  abs_vec[i_ss](joker, Tind, joker, joker) += abs_tmp;
202  } else {
203  ostringstream os;
204  os << "Interpolation error for (flat-array) scattering element #"
205  << i_se_flat << "\n"
206  << "at location/temperature point #" << Tind << "\n";
207  throw runtime_error(os.str());
208  }
209  }
210  }
211  i_se_flat++;
212  }
213  ptype[i_ss] = max(ptypes_se[i_ss]);
214  }
215 }
216 
218 
248 void opt_prop_NScatElems( //Output
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
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  assert(i_se_flat == Nse_all);
309 }
310 
312 
335 void opt_prop_1ScatElem( //Output
336  Tensor5View ext_mat, // nf, nT, ndir, nst, nst
337  Tensor4View abs_vec, // nf, nT, ndir, nst
338  Index& ptype,
339  VectorView t_ok,
340  //Input
341  const SingleScatteringData& ssd,
342  const Vector& T_array,
343  const Matrix& dir_array,
344  const Index& f_start,
345  const Index& t_interp_order) {
346  // FIXME: this is prob best done in scat_data_checkedCalc (or
347  // cloudbox_checkedCalc) to have it done once and for all. Here at max assert.
348 
349  // At very first check validity of the scatt elements ptype (so far we only
350  // handle PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND).
351  /*
352  if( ssd.ptype != PTYPE_TOTAL_RND and ssd.ptype != PTYPE_AZIMUTH_RND )
353  {
354  ostringstream os;
355  os << "Only ptypes " << PTYPE_TOTAL_RND << " and " << PTYPE_AZIMUTH_RND
356  << " can be handled.\n"
357  << "Encountered scattering element with ptype " << ssd.ptype
358  << ", though.";
359  throw runtime_error( os.str() );
360  }
361  */
362 
363  assert(ssd.ptype == PTYPE_TOTAL_RND or ssd.ptype == PTYPE_AZIMUTH_RND);
364 
365  const Index nf = ext_mat.nshelves();
366  assert(abs_vec.nbooks() == nf);
367  if (nf > 1) {
368  assert(nf == ssd.f_grid.nelem());
369  }
370 
371  const Index nTout = T_array.nelem();
372  assert(ext_mat.nbooks() == nTout);
373  assert(abs_vec.npages() == nTout);
374  assert(t_ok.nelem() == nTout);
375 
376  const Index nDir = dir_array.nrows();
377  assert(ext_mat.npages() == nDir);
378  assert(abs_vec.nrows() == nDir);
379 
380  const Index stokes_dim = abs_vec.ncols();
381  assert(ext_mat.nrows() == stokes_dim);
382  assert(ext_mat.ncols() == stokes_dim);
383 
384  ptype = ssd.ptype;
385 
386  // Determine T-interpol order as well as interpol positions and weights (they
387  // are the same for all directions (and freqs), ie it is sufficient to
388  // calculate them once).
389  const Index nTin = ssd.T_grid.nelem();
390  Index this_T_interp_order;
391  Matrix T_itw;
392  ArrayOfGridPosPoly T_gp(nTout);
394  this_T_interp_order,
395  T_gp,
396  T_itw,
397  ssd.T_grid,
398  T_array,
399  t_interp_order);
400 
401  // Now loop over requested directions (and apply simultaneously for all freqs):
402  // 1) extract/interpolate direction (but not for tot.random)
403  // 2) apply T-interpol
404  if (ptype == PTYPE_TOTAL_RND)
405  if (this_T_interp_order < 0) // just extract (and unpack) ext and abs data
406  // for the fs, the Tin, and the dirin and sort
407  // (copy) into the output fs, Ts, and dirs.
408  {
409  Tensor3 ext_mat_tmp(nf, stokes_dim, stokes_dim);
410  Matrix abs_vec_tmp(nf, stokes_dim);
411  for (Index find = 0; find < nf; find++) {
412  ext_mat_SSD2Stokes(ext_mat_tmp(find, joker, joker),
413  ssd.ext_mat_data(find + f_start, 0, 0, 0, joker),
414  stokes_dim,
415  ptype);
416  abs_vec_SSD2Stokes(abs_vec_tmp(find, joker),
417  ssd.abs_vec_data(find + f_start, 0, 0, 0, joker),
418  stokes_dim,
419  ptype);
420  }
421  for (Index Tind = 0; Tind < nTout; Tind++)
422  for (Index dind = 0; dind < nDir; dind++) {
423  ext_mat(joker, Tind, dind, joker, joker) = ext_mat_tmp;
424  abs_vec(joker, Tind, dind, joker) = abs_vec_tmp;
425  }
426  } else // T-interpolation required (but not dir). To be done on the compact
427  // ssd format.
428  {
429  Tensor4 ext_mat_tmp(nf, nTout, stokes_dim, stokes_dim);
430  Tensor3 abs_vec_tmp(nf, nTout, stokes_dim);
431  Matrix ext_mat_tmp_ssd(nTout, ssd.ext_mat_data.ncols());
432  Matrix abs_vec_tmp_ssd(nTout, ssd.abs_vec_data.ncols());
433  for (Index find = 0; find < nf; find++) {
434  for (Index nst = 0; nst < ext_mat_tmp_ssd.ncols(); nst++)
435  interp(ext_mat_tmp_ssd(joker, nst),
436  T_itw,
437  ssd.ext_mat_data(find + f_start, joker, 0, 0, nst),
438  T_gp);
439  for (Index Tind = 0; Tind < nTout; Tind++)
440  if (t_ok[Tind] > 0.)
441  ext_mat_SSD2Stokes(ext_mat_tmp(find, Tind, joker, joker),
442  ext_mat_tmp_ssd(Tind, joker),
443  stokes_dim,
444  ptype);
445  else
446  ext_mat_tmp(find, Tind, joker, joker) = 0.;
447 
448  for (Index nst = 0; nst < abs_vec_tmp_ssd.ncols(); nst++)
449  interp(abs_vec_tmp_ssd(joker, nst),
450  T_itw,
451  ssd.abs_vec_data(find + f_start, joker, 0, 0, nst),
452  T_gp);
453  for (Index Tind = 0; Tind < nTout; Tind++)
454  if (t_ok[Tind] > 0.)
455  abs_vec_SSD2Stokes(abs_vec_tmp(find, Tind, joker),
456  abs_vec_tmp_ssd(Tind, joker),
457  stokes_dim,
458  ptype);
459  else
460  abs_vec_tmp(find, Tind, joker) = 0.;
461  }
462 
463  for (Index dind = 0; dind < nDir; dind++) {
464  ext_mat(joker, joker, dind, joker, joker) = ext_mat_tmp;
465  abs_vec(joker, joker, dind, joker) = abs_vec_tmp;
466  }
467  }
468 
469  else // dir-interpolation for non-tot-random particles
470  {
471  // as we don't allow other than az.random here, ext and abs will only vary
472  // with za, not with aa. Hence, we could be smart here and calc data only
473  // for unique za (and copy the rest). however, this smartness might cost as
474  // well. so for now, we leave that and blindly loop over the given direction
475  // array, and leave smartness in setting up directional array to the calling
476  // methods.
477 
478  // derive the direction interpolation weights.
479  ArrayOfGridPos dir_gp(nDir);
480  gridpos(dir_gp, ssd.za_grid, dir_array(joker, 0));
481  Matrix dir_itw(nDir, 2); // only interpolating in za, ie 1D linear interpol
482  interpweights(dir_itw, dir_gp);
483 
484  const Index next = ssd.ext_mat_data.ncols();
485  const Index nabs = ssd.abs_vec_data.ncols();
486 
487  if (this_T_interp_order < 0) // T only needs to be extracted.
488  {
489  Matrix ext_mat_tmp_ssd(nDir, next);
490  Matrix abs_vec_tmp_ssd(nDir, nabs);
491  Tensor4 ext_mat_tmp(nf, nDir, stokes_dim, stokes_dim);
492  Tensor3 abs_vec_tmp(nf, nDir, stokes_dim);
493  for (Index find = 0; find < nf; find++) {
494  for (Index nst = 0; nst < next; nst++)
495  interp(ext_mat_tmp_ssd(joker, nst),
496  dir_itw,
497  ssd.ext_mat_data(find + f_start, 0, joker, 0, nst),
498  dir_gp);
499  for (Index Dind = 0; Dind < nDir; Dind++)
500  ext_mat_SSD2Stokes(ext_mat_tmp(find, Dind, joker, joker),
501  ext_mat_tmp_ssd(Dind, joker),
502  stokes_dim,
503  ptype);
504 
505  for (Index nst = 0; nst < nabs; nst++)
506  interp(abs_vec_tmp_ssd(joker, nst),
507  dir_itw,
508  ssd.abs_vec_data(find + f_start, 0, joker, 0, nst),
509  dir_gp);
510  for (Index Dind = 0; Dind < nDir; Dind++)
511  abs_vec_SSD2Stokes(abs_vec_tmp(find, Dind, joker),
512  abs_vec_tmp_ssd(Dind, joker),
513  stokes_dim,
514  ptype);
515  }
516 
517  for (Index Tind = 0; Tind < nTout; Tind++) {
518  ext_mat(joker, Tind, joker, joker, joker) = ext_mat_tmp;
519  abs_vec(joker, Tind, joker, joker) = abs_vec_tmp;
520  }
521  } else // T- and dir-interpolation required. To be done on the compact ssd
522  // format.
523  {
524  Tensor3 ext_mat_tmp_ssd(nTin, nDir, next);
525  Tensor3 abs_vec_tmp_ssd(nTin, nDir, nabs);
526  Matrix ext_mat_tmp(nTout, next);
527  Matrix abs_vec_tmp(nTout, nabs);
528 
529  for (Index find = 0; find < nf; find++) {
530  for (Index Tind = 0; Tind < nTin; Tind++) {
531  for (Index nst = 0; nst < next; nst++)
532  interp(ext_mat_tmp_ssd(Tind, joker, nst),
533  dir_itw,
534  ssd.ext_mat_data(find + f_start, Tind, joker, 0, nst),
535  dir_gp);
536  for (Index nst = 0; nst < nabs; nst++)
537  interp(abs_vec_tmp_ssd(Tind, joker, nst),
538  dir_itw,
539  ssd.abs_vec_data(find + f_start, Tind, joker, 0, nst),
540  dir_gp);
541  }
542 
543  for (Index Dind = 0; Dind < nDir; Dind++) {
544  for (Index nst = 0; nst < next; nst++)
545  interp(ext_mat_tmp(joker, nst),
546  T_itw,
547  ext_mat_tmp_ssd(joker, Dind, nst),
548  T_gp);
549  for (Index Tind = 0; Tind < nTout; Tind++)
550  ext_mat_SSD2Stokes(ext_mat(find, Tind, Dind, joker, joker),
551  ext_mat_tmp(Tind, joker),
552  stokes_dim,
553  ptype);
554 
555  for (Index nst = 0; nst < nabs; nst++)
556  interp(abs_vec_tmp(joker, nst),
557  T_itw,
558  abs_vec_tmp_ssd(joker, Dind, nst),
559  T_gp);
560  for (Index Tind = 0; Tind < nTout; Tind++)
561  abs_vec_SSD2Stokes(abs_vec(find, Tind, Dind, joker),
562  abs_vec_tmp(Tind, joker),
563  stokes_dim,
564  ptype);
565  }
566  }
567  }
568  }
569 }
570 
572 
585 void ext_mat_SSD2Stokes( //Output
586  MatrixView ext_mat_stokes,
587  //Input
588  ConstVectorView ext_mat_ssd,
589  const Index& stokes_dim,
590  const Index& ptype) {
591  // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
592  // calling methods, though.
593  assert(ptype <= PTYPE_AZIMUTH_RND);
594 
595  ext_mat_stokes = 0.;
596 
597  for (Index ist = 0; ist < stokes_dim; ist++) {
598  ext_mat_stokes(ist, ist) = ext_mat_ssd[0];
599  }
600 
601  if (ptype > PTYPE_TOTAL_RND) {
602  switch (stokes_dim) {
603  case 4: {
604  ext_mat_stokes(2, 3) = ext_mat_ssd[2];
605  ext_mat_stokes(3, 2) = -ext_mat_ssd[2];
606  } /* FALLTHROUGH */
607  case 3: {
608  // nothing to be done here. but we need this for executing the below
609  // also in case of stokes_dim==3.
610  }
611  case 2: {
612  ext_mat_stokes(0, 1) = ext_mat_ssd[1];
613  ext_mat_stokes(1, 0) = ext_mat_ssd[1];
614  }
615  }
616  }
617 }
618 
620 
633 void abs_vec_SSD2Stokes( //Output
634  VectorView abs_vec_stokes,
635  //Input
636  ConstVectorView abs_vec_ssd,
637  const Index& stokes_dim,
638  const Index& ptype) {
639  // for now, no handling of PTYPE_GENERAL. should be ensured somewhere in the
640  // calling methods, though.
641  assert(ptype <= PTYPE_AZIMUTH_RND);
642 
643  abs_vec_stokes = 0.;
644 
645  abs_vec_stokes[0] = abs_vec_ssd[0];
646 
647  if (ptype > PTYPE_TOTAL_RND and stokes_dim > 1) {
648  abs_vec_stokes[1] = abs_vec_ssd[1];
649  }
650 }
651 
653 
670 void pha_mat_Bulk( //Output
671  Tensor6& pha_mat, // (nf,nT,npdir,nidir,nst,nst)
672  Index& ptype,
673  //Input
674  const ArrayOfTensor6& pha_mat_ss, // [nss](nf,nT,npdir,nidir,nst,nst)
675  const ArrayOfIndex& ptypes_ss) {
676  pha_mat = pha_mat_ss[0];
677 
678  for (Index i_ss = 1; i_ss < pha_mat_ss.nelem(); i_ss++)
679  pha_mat += pha_mat_ss[i_ss];
680 
681  ptype = max(ptypes_ss);
682 }
683 
685 
708 void pha_mat_ScatSpecBulk( //Output
709  ArrayOfTensor6& pha_mat, // [nss](nf,nT,npdir,nidir,nst,nst)
710  ArrayOfIndex& ptype,
711  //Input
712  const ArrayOfArrayOfTensor6&
713  pha_mat_se, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
714  const ArrayOfArrayOfIndex& ptypes_se,
715  ConstMatrixView pnds,
716  ConstMatrixView t_ok) {
717  assert(t_ok.nrows() == pnds.nrows());
718  assert(t_ok.ncols() == pnds.ncols());
719  assert(TotalNumberOfElements(pha_mat_se) == pnds.nrows());
720 
721  const Index nT = pnds.ncols();
722  const Index nf = pha_mat_se[0][0].nvitrines();
723  const Index npDir = pha_mat_se[0][0].nbooks();
724  const Index niDir = pha_mat_se[0][0].npages();
725  const Index stokes_dim = pha_mat_se[0][0].ncols();
726 
727  const Index nss = pha_mat_se.nelem();
728  pha_mat.resize(nss);
729  ptype.resize(nss);
730  Tensor5 pha_tmp;
731 
732  Index i_se_flat = 0;
733 
734  for (Index i_ss = 0; i_ss < nss; i_ss++) {
735  assert(nT == pha_mat_se[i_ss][0].nshelves());
736 
737  pha_mat[i_ss].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
738  pha_mat[i_ss] = 0.;
739 
740  for (Index i_se = 0; i_se < pha_mat_se[i_ss].nelem(); i_se++) {
741  assert(nT == pha_mat_se[i_ss][i_se].nshelves());
742 
743  for (Index Tind = 0; Tind < nT; Tind++) {
744  if (pnds(i_se_flat, Tind) != 0.) {
745  if (t_ok(i_se_flat, Tind) > 0.) {
746  pha_tmp =
747  pha_mat_se[i_ss][i_se](joker, Tind, joker, joker, joker, joker);
748  pha_tmp *= pnds(i_se_flat, Tind);
749  pha_mat[i_ss](joker, Tind, joker, joker, joker, joker) += pha_tmp;
750  } else {
751  ostringstream os;
752  os << "Interpolation error for (flat-array) scattering element #"
753  << i_se_flat << "\n"
754  << "at location/temperature point #" << Tind << "\n";
755  throw runtime_error(os.str());
756  }
757  }
758  }
759  i_se_flat++;
760  }
761  ptype[i_ss] = max(ptypes_se[i_ss]);
762  }
763 }
764 
766 
796 void pha_mat_NScatElems( //Output
797  ArrayOfArrayOfTensor6& pha_mat, // [nss][nse](nf,nT,npdir,nidir,nst,nst)
798  ArrayOfArrayOfIndex& ptypes,
799  Matrix& t_ok,
800  //Input
802  const Index& stokes_dim,
803  const Vector& T_array,
804  const Matrix& pdir_array,
805  const Matrix& idir_array,
806  const Index& f_index,
807  const Index& t_interp_order) {
808  Index f_start, nf;
809  if (f_index < 0) {
810  nf = scat_data[0][0].pha_mat_data.nlibraries();
811  f_start = 0;
812  //f_end = f_start+nf;
813  } else {
814  nf = 1;
815  if (scat_data[0][0].pha_mat_data.nlibraries() == 1)
816  f_start = 0;
817  else
818  f_start = f_index;
819  //f_end = f_start+nf;
820  }
821 
822  const Index nT = T_array.nelem();
823  const Index npDir = pdir_array.nrows();
824  const Index niDir = idir_array.nrows();
825 
826  const Index nss = scat_data.nelem();
827  pha_mat.resize(nss);
828  ptypes.resize(nss);
829 
830  const Index Nse_all = TotalNumberOfElements(scat_data);
831  t_ok.resize(Nse_all, nT);
832  Index i_se_flat = 0;
833 
834  for (Index i_ss = 0; i_ss < nss; i_ss++) {
835  Index nse = scat_data[i_ss].nelem();
836  pha_mat[i_ss].resize(nse);
837  ptypes[i_ss].resize(nse);
838 
839  for (Index i_se = 0; i_se < nse; i_se++) {
840  pha_mat[i_ss][i_se].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
841 
842  pha_mat_1ScatElem(pha_mat[i_ss][i_se],
843  ptypes[i_ss][i_se],
844  t_ok(i_se_flat, joker),
845  scat_data[i_ss][i_se],
846  T_array,
847  pdir_array,
848  idir_array,
849  f_start,
850  t_interp_order);
851  i_se_flat++;
852  }
853  }
854  assert(i_se_flat == Nse_all);
855 }
856 
858 
881 void pha_mat_1ScatElem( //Output
882  Tensor6View pha_mat, // nf, nT, npdir, nidir, nst, nst
883  Index& ptype,
884  VectorView t_ok,
885  //Input
886  const SingleScatteringData& ssd,
887  const Vector& T_array,
888  const Matrix& pdir_array,
889  const Matrix& idir_array,
890  const Index& f_start,
891  const Index& t_interp_order) {
892  assert(ssd.ptype == PTYPE_TOTAL_RND or ssd.ptype == PTYPE_AZIMUTH_RND);
893 
894  const Index nf = pha_mat.nvitrines();
895  if (nf > 1) {
896  assert(nf == ssd.f_grid.nelem());
897  }
898 
899  const Index nTout = T_array.nelem();
900  assert(pha_mat.nshelves() == nTout);
901  assert(t_ok.nelem() == nTout);
902 
903  const Index npDir = pdir_array.nrows();
904  assert(pha_mat.nbooks() == npDir);
905 
906  const Index niDir = idir_array.nrows();
907  assert(pha_mat.npages() == niDir);
908 
909  const Index stokes_dim = pha_mat.ncols();
910  assert(pha_mat.nrows() == stokes_dim);
911 
912  ptype = ssd.ptype;
913 
914  // Determine T-interpol order as well as interpol positions and weights (they
915  // are the same for all directions (and freqs), ie it is sufficient to
916  // calculate them once).
917  const Index nTin = ssd.T_grid.nelem();
918  Index this_T_interp_order;
919  Matrix T_itw;
920  ArrayOfGridPosPoly T_gp(nTout);
922  this_T_interp_order,
923  T_gp,
924  T_itw,
925  ssd.T_grid,
926  T_array,
927  t_interp_order);
928 
929  // Now loop over requested directions (and apply simultaneously for all freqs):
930  // 1) interpolate direction
931  // 2) apply T-interpol
932  if (ptype == PTYPE_TOTAL_RND) {
933  // determine how many of the compact stokes elements we will need.
934  // restrict interpolations to those.
935  Index npha;
936  if (stokes_dim == 1)
937  npha = 1;
938  else if (stokes_dim < 4) // stokes_dim==2 || stokes_dim==3
939  npha = 4;
940  else
941  npha = 6;
942  if (this_T_interp_order < 0) // just extract (and unpack) pha data for the
943  // fs and Tin, and interpolate in sca-angs, and
944  // sort (copy) into the output fs, Ts, and dirs.
945  {
946  for (Index pdir = 0; pdir < npDir; pdir++)
947  for (Index idir = 0; idir < niDir; idir++) {
948  // calc scat ang theta from incident and prop dirs
949  Numeric theta = scat_angle(pdir_array(pdir, 0),
950  pdir_array(pdir, 1),
951  idir_array(idir, 0),
952  idir_array(idir, 1));
953 
954  // get scat angle interpolation weights
955  GridPos dir_gp;
956  gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
957  Vector dir_itw(2);
958  interpweights(dir_itw, dir_gp);
959 
960  Vector pha_mat_int(npha, 0.);
961  Matrix pha_mat_tmp(stokes_dim, stokes_dim);
962  for (Index find = 0; find < nf; find++) {
963  // perform the scat angle interpolation
964  for (Index nst = 0; nst < npha; nst++)
965  pha_mat_int[nst] = interp(
966  dir_itw,
967  ssd.pha_mat_data(find + f_start, 0, joker, 0, 0, 0, nst),
968  dir_gp);
969 
970  // convert from scat to lab frame
971  pha_mat_labCalc(pha_mat_tmp(joker, joker),
972  pha_mat_int,
973  pdir_array(pdir, 0),
974  pdir_array(pdir, 1),
975  idir_array(idir, 0),
976  idir_array(idir, 1),
977  theta);
978 
979  for (Index Tind = 0; Tind < nTout; Tind++)
980  pha_mat(find, Tind, pdir, idir, joker, joker) = pha_mat_tmp;
981  }
982  }
983  } else // T-interpolation required. To be done on the compact ssd format.
984  {
985  for (Index pdir = 0; pdir < npDir; pdir++)
986  for (Index idir = 0; idir < niDir; idir++) {
987  // calc scat ang theta from incident and prop dirs
988  Numeric theta = scat_angle(pdir_array(pdir, 0),
989  pdir_array(pdir, 1),
990  idir_array(idir, 0),
991  idir_array(idir, 1));
992 
993  // get scat angle interpolation weights
994  GridPos dir_gp;
995  gridpos(dir_gp, ssd.za_grid, theta * RAD2DEG);
996  Vector dir_itw(2);
997  interpweights(dir_itw, dir_gp);
998 
999  Matrix pha_mat_int(nTin, npha, 0.);
1000  Matrix pha_mat_tmp(nTout, npha, 0.);
1001  for (Index find = 0; find < nf; find++) {
1002  for (Index Tind = 0; Tind < nTin; Tind++)
1003  // perform the scat angle interpolation
1004  for (Index nst = 0; nst < npha; nst++) {
1005  pha_mat_int(Tind, nst) = interp(
1006  dir_itw,
1007  ssd.pha_mat_data(find + f_start, Tind, joker, 0, 0, 0, nst),
1008  dir_gp);
1009  }
1010  // perform the T-interpolation
1011  for (Index nst = 0; nst < npha; nst++) {
1012  interp(pha_mat_tmp(joker, nst),
1013  T_itw,
1014  pha_mat_int(joker, nst),
1015  T_gp);
1016  }
1017  // FIXME: it's probably better to do the frame conversion first,
1018  // then the T-interpolation (which is better depends on how many T-
1019  // aka vertical grid points we have. for single point, T-interpol
1020  // first is better, for a full vertical grid, frame conversion first
1021  // should be faster.)
1022  for (Index Tind = 0; Tind < nTout; Tind++) {
1023  // convert from scat to lab frame
1024  pha_mat_labCalc(pha_mat(find, Tind, pdir, idir, joker, joker),
1025  pha_mat_tmp(Tind, joker),
1026  pdir_array(pdir, 0),
1027  pdir_array(pdir, 1),
1028  idir_array(idir, 0),
1029  idir_array(idir, 1),
1030  theta);
1031  }
1032  }
1033  }
1034  }
1035  } else // dir-interpolation for non-tot-random particles
1036  // Data is already stored in the laboratory frame,
1037  // but it is compressed a little. Details elsewhere.
1038  {
1039  Index nDir = npDir * niDir;
1040  Vector adelta_aa(nDir);
1041  Matrix delta_aa(npDir, niDir);
1042  ArrayOfGridPos daa_gp(nDir), pza_gp(nDir), iza_gp(nDir);
1043  ArrayOfGridPos pza_gp_tmp(npDir), iza_gp_tmp(niDir);
1044 
1045  gridpos(pza_gp_tmp, ssd.za_grid, pdir_array(joker, 0));
1046  gridpos(iza_gp_tmp, ssd.za_grid, idir_array(joker, 0));
1047 
1048  Index j = 0;
1049  for (Index pdir = 0; pdir < npDir; pdir++) {
1050  for (Index idir = 0; idir < niDir; idir++) {
1051  delta_aa(pdir, idir) =
1052  pdir_array(pdir, 1) - idir_array(idir, 1) +
1053  (pdir_array(pdir, 1) - idir_array(idir, 1) < -180) * 360 -
1054  (pdir_array(pdir, 1) - idir_array(idir, 1) > 180) * 360;
1055  adelta_aa[j] = abs(delta_aa(pdir, idir));
1056  pza_gp[j] = pza_gp_tmp[pdir];
1057  iza_gp[j] = iza_gp_tmp[idir];
1058  j++;
1059  }
1060  }
1061 
1062  gridpos(daa_gp, ssd.aa_grid, adelta_aa);
1063 
1064  Matrix dir_itw(nDir, 8);
1065  interpweights(dir_itw, pza_gp, daa_gp, iza_gp);
1066 
1067  if (this_T_interp_order < 0) // T only needs to be extracted.
1068  {
1069  Tensor3 pha_mat_int(nDir, stokes_dim, stokes_dim, 0.);
1070  Tensor4 pha_mat_tmp(npDir, niDir, stokes_dim, stokes_dim);
1071 
1072  for (Index find = 0; find < nf; find++) {
1073  // perform the (tri-linear) angle interpolation. but only for the
1074  // pha_mat elements that we actually need.
1075 
1076  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1077  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1078  interp(
1079  pha_mat_int(joker, ist1, ist2),
1080  dir_itw,
1081  ssd.pha_mat_data(
1082  find + f_start, 0, joker, joker, joker, 0, ist1 * 4 + ist2),
1083  pza_gp,
1084  daa_gp,
1085  iza_gp);
1086 
1087  // sort direction-combined 1D-array back into prop and incident
1088  // direction matrix
1089  Index i = 0;
1090  for (Index pdir = 0; pdir < npDir; pdir++)
1091  for (Index idir = 0; idir < niDir; idir++) {
1092  pha_mat_tmp(pdir, idir, joker, joker) =
1093  pha_mat_int(i, joker, joker);
1094  i++;
1095  }
1096 
1097  if (stokes_dim > 2) {
1098  for (Index pdir = 0; pdir < npDir; pdir++)
1099  for (Index idir = 0; idir < niDir; idir++)
1100  if (delta_aa(pdir, idir) < 0.) {
1101  pha_mat_tmp(pdir, idir, 0, 2) *= -1;
1102  pha_mat_tmp(pdir, idir, 1, 2) *= -1;
1103  pha_mat_tmp(pdir, idir, 2, 0) *= -1;
1104  pha_mat_tmp(pdir, idir, 2, 1) *= -1;
1105  }
1106  }
1107 
1108  if (stokes_dim > 3) {
1109  for (Index pdir = 0; pdir < npDir; pdir++)
1110  for (Index idir = 0; idir < niDir; idir++)
1111  if (delta_aa(pdir, idir) < 0.) {
1112  pha_mat_tmp(pdir, idir, 0, 3) *= -1;
1113  pha_mat_tmp(pdir, idir, 1, 3) *= -1;
1114  pha_mat_tmp(pdir, idir, 3, 0) *= -1;
1115  pha_mat_tmp(pdir, idir, 3, 1) *= -1;
1116  }
1117  }
1118 
1119  for (Index Tind = 0; Tind < nTout; Tind++)
1120  pha_mat(find, Tind, joker, joker, joker, joker) = pha_mat_tmp;
1121  }
1122  }
1123 
1124  else // T- and dir-interpolation required. To be done on the compact ssd
1125  // format.
1126  {
1127  Tensor4 pha_mat_int(nTin, nDir, stokes_dim, stokes_dim, 0.);
1128 
1129  for (Index find = 0; find < nf; find++) {
1130  // perform the (tri-linear) angle interpolation. but only for the
1131  // pha_mat elements that we actually need.
1132  for (Index Tind = 0; Tind < nTin; Tind++) {
1133  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1134  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1135  interp(pha_mat_int(Tind, joker, ist1, ist2),
1136  dir_itw,
1137  ssd.pha_mat_data(find + f_start,
1138  Tind,
1139  joker,
1140  joker,
1141  joker,
1142  0,
1143  ist1 * 4 + ist2),
1144  pza_gp,
1145  daa_gp,
1146  iza_gp);
1147  }
1148 
1149  // perform the T-interpolation and simultaneously sort
1150  // direction-combined 1D-array back into prop and incident direction
1151  // matrix.
1152  Index i = 0;
1153  for (Index pdir = 0; pdir < npDir; pdir++)
1154  for (Index idir = 0; idir < niDir; idir++) {
1155  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1156  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1157  interp(pha_mat(find, joker, pdir, idir, ist1, ist2),
1158  T_itw,
1159  pha_mat_int(joker, i, ist1, ist2),
1160  T_gp);
1161  i++;
1162  }
1163 
1164  if (stokes_dim > 2) {
1165  for (Index pdir = 0; pdir < npDir; pdir++)
1166  for (Index idir = 0; idir < niDir; idir++)
1167  if (delta_aa(pdir, idir) < 0.) {
1168  pha_mat(find, joker, pdir, idir, 0, 2) *= -1;
1169  pha_mat(find, joker, pdir, idir, 1, 2) *= -1;
1170  pha_mat(find, joker, pdir, idir, 2, 0) *= -1;
1171  pha_mat(find, joker, pdir, idir, 2, 1) *= -1;
1172  }
1173  }
1174 
1175  if (stokes_dim > 3) {
1176  for (Index pdir = 0; pdir < npDir; pdir++)
1177  for (Index idir = 0; idir < niDir; idir++)
1178  if (delta_aa(pdir, idir) < 0.) {
1179  pha_mat(find, joker, pdir, idir, 0, 3) *= -1;
1180  pha_mat(find, joker, pdir, idir, 1, 3) *= -1;
1181  pha_mat(find, joker, pdir, idir, 3, 0) *= -1;
1182  pha_mat(find, joker, pdir, idir, 3, 1) *= -1;
1183  }
1184  }
1185  }
1186  }
1187  }
1188 }
1189 
1191 
1218 void FouComp_1ScatElem( //Output
1219  Tensor7View pha_mat_fou, // nf, nT, npdir, nidir, nst, nst, m
1220  Index& ptype,
1221  VectorView t_ok,
1222  //Input
1223  const SingleScatteringData& ssd,
1224  const Vector& T_array,
1225  const Vector& pdir_array,
1226  const Vector& idir_array,
1227  const Index& f_start,
1228  const Index& t_interp_order,
1229  const Index& naa_totran) {
1230  assert(ssd.ptype == PTYPE_TOTAL_RND or ssd.ptype == PTYPE_AZIMUTH_RND);
1231 
1232  const Index nf = pha_mat_fou.nlibraries();
1233  if (nf > 1) {
1234  assert(nf == ssd.f_grid.nelem());
1235  }
1236 
1237  const Index nTout = T_array.nelem();
1238  assert(pha_mat_fou.nvitrines() == nTout);
1239  assert(t_ok.nelem() == nTout);
1240 
1241  const Index npDir = pdir_array.nelem();
1242  assert(pha_mat_fou.nshelves() == npDir);
1243  const Index niDir = idir_array.nelem();
1244  assert(pha_mat_fou.nbooks() == niDir);
1245 
1246  const Index stokes_dim = pha_mat_fou.nrows();
1247  assert(pha_mat_fou.npages() == stokes_dim);
1248  // currently code is only prepared for stokes_dim up to 2 (nothing else needed in
1249  // RT4 and generally in azimuth-symmetrical system)
1250  assert(stokes_dim < 3);
1251 
1252  const Index nmodes = pha_mat_fou.ncols();
1253  // currently code is only prepared for fourier mode 0 (nothing else needed in
1254  // RT4 and generally in azimuth-symmetrical system)
1255  assert(nmodes == 0);
1256 
1257  ptype = ssd.ptype;
1258 
1259  // Determine T-interpol order as well as interpol positions and weights (they
1260  // are the same for all directions (and freqs), ie it is sufficient to
1261  // calculate them once).
1262  const Index nTin = ssd.T_grid.nelem();
1263  Index this_T_interp_order;
1264  Matrix T_itw;
1265  ArrayOfGridPosPoly T_gp(nTout);
1267  this_T_interp_order,
1268  T_gp,
1269  T_itw,
1270  ssd.T_grid,
1271  T_array,
1272  t_interp_order);
1273 
1274  // 1) derive Fourier component(s)
1275  // 2) apply T-interpol
1276  if (ptype == PTYPE_TOTAL_RND) {
1277  // DCalculate azimuth angles and their integration weights for Fourier
1278  // component derivation (they are only determined by naa_totran).
1279  if (naa_totran < 3) {
1280  ostringstream os;
1281  os << "Azimuth grid size for scatt matrix extraction"
1282  << " (*naa_totran*) must be >3.\n"
1283  << "Yours is " << naa_totran << ".\n";
1284  throw runtime_error(os.str());
1285  }
1286  Vector aa_grid;
1287  nlinspace(aa_grid, 0, 180, naa_totran);
1288  Numeric daa_totran =
1289  1. / float(naa_totran - 1); // 2*180./360./(naa_totran-1)
1290  Vector theta(naa_totran);
1291  ArrayOfGridPos theta_gp;
1292  Matrix theta_itw(naa_totran, 2);
1293 
1294  Index npha;
1295  if (stokes_dim == 1)
1296  npha = 1;
1297  else if (stokes_dim < 4) // stokes_dim==2 || stokes_dim==3
1298  npha = 4;
1299  else
1300  npha = 6;
1301 
1303  Matrix pha_mat_angint(naa_totran, npha);
1304  Tensor3 Fou_int(nTin, stokes_dim, stokes_dim);
1305 
1306  for (Index idir = 0; idir < niDir; idir++)
1307  for (Index pdir = 0; pdir < npDir; pdir++) {
1308  for (Index iaa = 0; iaa < naa_totran; iaa++)
1309  // calc scat ang theta from incident and prop dirs
1310  theta[iaa] =
1311  scat_angle(pdir_array[pdir], aa_grid[iaa], idir_array[idir], 0.);
1312 
1313  // get scat angle interpolation weights
1314  gridpos(theta_gp, ssd.za_grid, theta * RAD2DEG);
1315  interpweights(theta_itw, theta_gp);
1316 
1317  for (Index find = 0; find < nf; find++) {
1318  Fou_int = 0.;
1319  for (Index Tind = 0; Tind < nTin; Tind++) {
1320  // perform the scat angle interpolation
1321  for (Index nst = 0; nst < npha; nst++)
1322  interp(
1323  pha_mat_angint(joker, nst),
1324  theta_itw,
1325  ssd.pha_mat_data(find + f_start, Tind, joker, 0, 0, 0, nst),
1326  theta_gp);
1327  for (Index iaa = 0; iaa < naa_totran; iaa++) {
1328  // convert from scat to lab frame
1330  pha_mat_angint(iaa, joker),
1331  pdir_array[pdir],
1332  aa_grid[iaa],
1333  idir_array[idir],
1334  0.,
1335  theta[iaa]);
1336  // and sum up/integrate
1337  // FIXME: can the integration probably be done in lab frame?
1338  // test!
1339  if (iaa == 0 || iaa == naa_totran - 1)
1340  pha_mat *= (daa_totran / 2.);
1341  else
1342  pha_mat *= daa_totran;
1343  Fou_int(Tind, joker, joker) += pha_mat;
1344  }
1345  }
1346 
1347  if (this_T_interp_order <
1348  0) // T only needs to be sorted into pha_mat_fou.
1349  {
1350  for (Index Tind = 0; Tind < nTout; Tind++)
1351  pha_mat_fou(find, Tind, pdir, idir, joker, joker, 0) =
1352  Fou_int(0, joker, joker);
1353  } else {
1354  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1355  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1356  for (Index im = 0; im < nmodes; im++)
1357  interp(pha_mat_fou(find, joker, pdir, idir, ist1, ist2, im),
1358  T_itw,
1359  Fou_int(joker, ist1, ist2),
1360  T_gp);
1361  }
1362  }
1363  }
1364  } else // derive Fourier component(s) on scattering elements own za_grid (using
1365  // given aa_grid), then 2D-interpolate inc and sca polar directions.
1366  // Do calcs only for required pha_mat components depending on stokes_dim.
1367  {
1368  Index nza = ssd.za_grid.nelem();
1369  Index naa = ssd.aa_grid.nelem();
1370  ConstVectorView za_datagrid = ssd.za_grid;
1371  ConstVectorView aa_datagrid = ssd.aa_grid;
1372  assert(aa_datagrid[0] == 0.);
1373  assert(aa_datagrid[naa - 1] == 180.);
1374  Vector daa(naa);
1375 
1376  // Precalculate azimuth integration weights for this azimuthally randomly
1377  // oriented scat element (need to do this per scat element as ssd.aa_grid is
1378  // scat element specific (might change between elements) and need to do this
1379  // on actual grid instead of grid number since the grid can, at least
1380  // theoretically be non-equidistant).
1381  daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1382  for (Index iaa = 1; iaa < naa - 1; iaa++)
1383  daa[iaa] = (aa_datagrid[iaa + 1] - aa_datagrid[iaa - 1]) / 360.;
1384  daa[naa - 1] = (aa_datagrid[naa - 1] - aa_datagrid[naa - 2]) / 360.;
1385 
1386  // Precalculate polar angle interpolation grid positions and weights
1387  ArrayOfGridPos pdir_za_gp, idir_za_gp;
1388  gridpos(pdir_za_gp, za_datagrid, pdir_array);
1389  gridpos(idir_za_gp, za_datagrid, idir_array);
1390  Tensor3 dir_itw(npDir, niDir, 4);
1391  interpweights(dir_itw, pdir_za_gp, idir_za_gp);
1392 
1393  Tensor4 Fou_ssd(nza, nza, stokes_dim, stokes_dim);
1394  Tensor5 Fou_angint(nTin, npDir, niDir, stokes_dim, stokes_dim);
1395 
1396  for (Index find = 0; find < nf; find++) {
1397  for (Index Tind = 0; Tind < nTin; Tind++) {
1398  // first, extract the phase matrix at the scatt elements own polar angle
1399  // grid and integrate over azimuth deriving their respective azimuthal
1400  // (Fourier series) 0-mode
1401  Fou_ssd = 0.;
1402  for (Index iza = 0; iza < nza; iza++)
1403  for (Index sza = 0; sza < nza; sza++) {
1404  for (Index iaa = 0; iaa < naa; iaa++) {
1405  // FIXME: are there any possible shortcuts for specific stokes
1406  // components (specifically, some where F0(ist1,ist2)==0?)
1407  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1408  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1409  Fou_ssd(sza, iza, ist1, ist2) +=
1410  daa[iaa] * ssd.pha_mat_data(find + f_start,
1411  Tind,
1412  sza,
1413  iaa,
1414  iza,
1415  0,
1416  ist1 * 4 + ist2);
1417  }
1418  }
1419 
1420  // second, interpolate the extracted azimuthal mode to the stream directions
1421  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1422  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1423  // FIXME: do we need to apply any sign changes?
1424  interp(Fou_angint(Tind, joker, joker, ist1, ist2),
1425  dir_itw,
1426  Fou_ssd(joker, joker, ist1, ist2),
1427  pdir_za_gp,
1428  idir_za_gp);
1429  }
1430 
1431  if (this_T_interp_order <
1432  0) // T only needs to be sorted into pha_mat_fou.
1433  {
1434  for (Index Tind = 0; Tind < nTout; Tind++)
1435  pha_mat_fou(find, Tind, joker, joker, joker, joker, 0) =
1436  Fou_angint(0, joker, joker, joker, joker);
1437  } else {
1438  for (Index pdir = 0; pdir < npDir; pdir++)
1439  for (Index idir = 0; idir < niDir; idir++)
1440  for (Index ist1 = 0; ist1 < stokes_dim; ist1++)
1441  for (Index ist2 = 0; ist2 < stokes_dim; ist2++)
1442  for (Index im = 0; im < nmodes; im++)
1443  interp(pha_mat_fou(find, joker, pdir, idir, ist1, ist2, im),
1444  T_itw,
1445  Fou_angint(joker, pdir, idir, ist1, ist2),
1446  T_gp);
1447  }
1448  }
1449  }
1450 }
1451 
1453 
1474 void abs_vecTransform( //Output and Input
1475  StokesVector& abs_vec_lab,
1476  //Input
1477  ConstTensor3View abs_vec_data,
1478  ConstVectorView za_datagrid,
1479  ConstVectorView aa_datagrid _U_,
1480  const PType& ptype,
1481  const Numeric& za_sca _U_,
1482  const Numeric& aa_sca _U_,
1483  const Verbosity& verbosity) {
1484  const Index stokes_dim = abs_vec_lab.StokesDimensions();
1485  assert(abs_vec_lab.NumberOfFrequencies() == 1);
1486 
1487  if (stokes_dim > 4 || stokes_dim < 1) {
1488  throw runtime_error(
1489  "The dimension of the stokes vector \n"
1490  "must be 1,2,3 or 4");
1491  }
1492 
1493  switch (ptype) {
1494  case PTYPE_GENERAL: {
1495  /*
1496  TO ANY DEVELOPER:
1497  current usage of coordinate systems in scattering solvers (RT and SSD
1498  extraction) and general radiative transfer is not consistent. Not an
1499  issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1500  but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1501  adding PTYPE_GENERAL support (see AUG appendix for more info).
1502  */
1503 
1504  CREATE_OUT0;
1505  out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1506  break;
1507  }
1508  case PTYPE_TOTAL_RND: {
1509  // The first element of the vector corresponds to the absorption
1510  // coefficient which is stored in the database, the others are 0.
1511 
1512  abs_vec_lab.SetZero();
1513 
1514  abs_vec_lab.Kjj()[0] = abs_vec_data(0, 0, 0);
1515  break;
1516  }
1517 
1518  case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1519  {
1520  assert(abs_vec_data.ncols() == 2);
1521 
1522  // In the case of azimuthally randomly oriented particles, only the first
1523  // two elements of the absorption coefficient vector are non-zero.
1524  // These values are dependent on the zenith angle of propagation.
1525 
1526  // 1st interpolate data by za_sca
1527  GridPos gp;
1528  Vector itw(2);
1529 
1530  gridpos(gp, za_datagrid, za_sca);
1531  interpweights(itw, gp);
1532 
1533  abs_vec_lab.SetZero();
1534 
1535  abs_vec_lab.Kjj()[0] = interp(itw, abs_vec_data(Range(joker), 0, 0), gp);
1536 
1537  if (stokes_dim == 1) {
1538  break;
1539  }
1540  abs_vec_lab.K12()[0] = interp(itw, abs_vec_data(Range(joker), 0, 1), gp);
1541  break;
1542  }
1543  default: {
1544  CREATE_OUT0;
1545  out0 << "Not all ptype cases are implemented\n";
1546  }
1547  }
1548 }
1549 
1551 
1572 void ext_matTransform( //Output and Input
1573  PropagationMatrix& ext_mat_lab,
1574  //Input
1575  ConstTensor3View ext_mat_data,
1576  ConstVectorView za_datagrid,
1577  ConstVectorView aa_datagrid _U_,
1578  const PType& ptype,
1579  const Numeric& za_sca,
1580  const Numeric& aa_sca _U_,
1581  const Verbosity& verbosity) {
1582  const Index stokes_dim = ext_mat_lab.StokesDimensions();
1583  assert(ext_mat_lab.NumberOfFrequencies() == 1);
1584 
1585  if (stokes_dim > 4 || stokes_dim < 1) {
1586  throw runtime_error(
1587  "The dimension of the stokes vector \n"
1588  "must be 1,2,3 or 4");
1589  }
1590 
1591  switch (ptype) {
1592  case PTYPE_GENERAL: {
1593  /*
1594  TO ANY DEVELOPER:
1595  current usage of coordinate systems in scattering solvers (RT and SSD
1596  extraction) and general radiative transfer is not consistent. Not an
1597  issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1598  but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1599  adding PTYPE_GENERAL support (see AUG appendix for more info).
1600  */
1601 
1602  CREATE_OUT0;
1603  out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1604  break;
1605  }
1606  case PTYPE_TOTAL_RND: {
1607  assert(ext_mat_data.ncols() == 1);
1608 
1609  // In the case of randomly oriented particles the extinction matrix is
1610  // diagonal. The value of each element of the diagonal is the
1611  // extinction cross section, which is stored in the database.
1612 
1613  ext_mat_lab.SetZero();
1614 
1615  ext_mat_lab.Kjj() = ext_mat_data(0, 0, 0);
1616  break;
1617  }
1618 
1619  case PTYPE_AZIMUTH_RND: //Added by Cory Davis 9/12/03
1620  {
1621  assert(ext_mat_data.ncols() == 3);
1622 
1623  // In the case of azimuthally randomly oriented particles, the extinction
1624  // matrix has only 3 independent non-zero elements Kjj, K12=K21, and K34=-K43.
1625  // These values are dependent on the zenith angle of propagation.
1626 
1627  // 1st interpolate data by za_sca
1628  GridPos gp;
1629  Vector itw(2);
1630  Numeric Kjj;
1631  Numeric K12;
1632  Numeric K34;
1633 
1634  gridpos(gp, za_datagrid, za_sca);
1635  interpweights(itw, gp);
1636 
1637  ext_mat_lab.SetZero();
1638 
1639  Kjj = interp(itw, ext_mat_data(Range(joker), 0, 0), gp);
1640  ext_mat_lab.Kjj() = Kjj;
1641 
1642  if (stokes_dim < 2) {
1643  break;
1644  }
1645 
1646  K12 = interp(itw, ext_mat_data(Range(joker), 0, 1), gp);
1647  ext_mat_lab.K12()[0] = K12;
1648 
1649  if (stokes_dim < 4) {
1650  break;
1651  }
1652 
1653  K34 = interp(itw, ext_mat_data(Range(joker), 0, 2), gp);
1654  ext_mat_lab.K34()[0] = K34;
1655  break;
1656  }
1657  default: {
1658  CREATE_OUT0;
1659  out0 << "Not all ptype cases are implemented\n";
1660  }
1661  }
1662 }
1663 
1665 
1692 void pha_matTransform( //Output
1693  MatrixView pha_mat_lab,
1694  //Input
1695  ConstTensor5View pha_mat_data,
1696  ConstVectorView za_datagrid,
1697  ConstVectorView aa_datagrid,
1698  const PType& ptype,
1699  const Index& za_sca_idx,
1700  const Index& aa_sca_idx,
1701  const Index& za_inc_idx,
1702  const Index& aa_inc_idx,
1705  const Verbosity& verbosity) {
1706  const Index stokes_dim = pha_mat_lab.ncols();
1707 
1708  Numeric za_sca = za_grid[za_sca_idx];
1709  Numeric aa_sca = aa_grid[aa_sca_idx];
1710  Numeric za_inc = za_grid[za_inc_idx];
1711  Numeric aa_inc = aa_grid[aa_inc_idx];
1712 
1713  if (stokes_dim > 4 || stokes_dim < 1) {
1714  throw runtime_error(
1715  "The dimension of the stokes vector \n"
1716  "must be 1,2,3 or 4");
1717  }
1718 
1719  switch (ptype) {
1720  case PTYPE_GENERAL: {
1721  /*
1722  TO ANY DEVELOPER:
1723  current usage of coordinate systems in scattering solvers (RT and SSD
1724  extraction) and general radiative transfer is not consistent. Not an
1725  issue as long as only PTYPE_TOTAL_RND and PTYPE_AZIMUTH_RND are used,
1726  but will be a problem for PTYPE_GENERAL, ie needs to be fixed BEFORE
1727  adding PTYPE_GENERAL support (see AUG appendix for more info).
1728  */
1729 
1730  CREATE_OUT0;
1731  out0 << "Case PTYPE_GENERAL not yet implemented. \n";
1732  break;
1733  }
1734  case PTYPE_TOTAL_RND: {
1735  // Calculate the scattering angle and interpolate the data in it:
1736  Numeric theta_rad = scat_angle(za_sca, aa_sca, za_inc, aa_inc);
1737  const Numeric theta = RAD2DEG * theta_rad;
1738 
1739  // Interpolation of the data on the scattering angle:
1740  Vector pha_mat_int(6);
1741  interpolate_scat_angle(pha_mat_int, pha_mat_data, za_datagrid, theta);
1742 
1743  // Calculate the phase matrix in the laboratory frame:
1745  pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1746 
1747  break;
1748  }
1749 
1750  case PTYPE_AZIMUTH_RND: //Added by Cory Davis
1751  //Data is already stored in the laboratory frame,
1752  //but it is compressed a little. Details elsewhere.
1753  {
1754  assert(pha_mat_data.ncols() == 16);
1755  assert(pha_mat_data.npages() == za_datagrid.nelem());
1756  Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
1757  (aa_sca - aa_inc > 180) *
1758  360; //delta_aa corresponds to the "books"
1759  //dimension of pha_mat_data
1760  GridPos za_sca_gp;
1761  GridPos delta_aa_gp;
1762  GridPos za_inc_gp;
1763  Vector itw(8);
1764 
1765  gridpos(delta_aa_gp, aa_datagrid, abs(delta_aa));
1766  gridpos(za_inc_gp, za_datagrid, za_inc);
1767  gridpos(za_sca_gp, za_datagrid, za_sca);
1768 
1769  interpweights(itw, za_sca_gp, delta_aa_gp, za_inc_gp);
1770 
1771  pha_mat_lab(0, 0) =
1772  interp(itw,
1773  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 0),
1774  za_sca_gp,
1775  delta_aa_gp,
1776  za_inc_gp);
1777  if (stokes_dim == 1) {
1778  break;
1779  }
1780  pha_mat_lab(0, 1) =
1781  interp(itw,
1782  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 1),
1783  za_sca_gp,
1784  delta_aa_gp,
1785  za_inc_gp);
1786  pha_mat_lab(1, 0) =
1787  interp(itw,
1788  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 4),
1789  za_sca_gp,
1790  delta_aa_gp,
1791  za_inc_gp);
1792  pha_mat_lab(1, 1) =
1793  interp(itw,
1794  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 5),
1795  za_sca_gp,
1796  delta_aa_gp,
1797  za_inc_gp);
1798  if (stokes_dim == 2) {
1799  break;
1800  }
1801  if (delta_aa >= 0) {
1802  pha_mat_lab(0, 2) = interp(
1803  itw,
1804  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1805  za_sca_gp,
1806  delta_aa_gp,
1807  za_inc_gp);
1808  pha_mat_lab(1, 2) = interp(
1809  itw,
1810  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1811  za_sca_gp,
1812  delta_aa_gp,
1813  za_inc_gp);
1814  pha_mat_lab(2, 0) = interp(
1815  itw,
1816  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1817  za_sca_gp,
1818  delta_aa_gp,
1819  za_inc_gp);
1820  pha_mat_lab(2, 1) = interp(
1821  itw,
1822  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1823  za_sca_gp,
1824  delta_aa_gp,
1825  za_inc_gp);
1826  } else {
1827  pha_mat_lab(0, 2) = -interp(
1828  itw,
1829  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 2),
1830  za_sca_gp,
1831  delta_aa_gp,
1832  za_inc_gp);
1833  pha_mat_lab(1, 2) = -interp(
1834  itw,
1835  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 6),
1836  za_sca_gp,
1837  delta_aa_gp,
1838  za_inc_gp);
1839  pha_mat_lab(2, 0) = -interp(
1840  itw,
1841  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 8),
1842  za_sca_gp,
1843  delta_aa_gp,
1844  za_inc_gp);
1845  pha_mat_lab(2, 1) = -interp(
1846  itw,
1847  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 9),
1848  za_sca_gp,
1849  delta_aa_gp,
1850  za_inc_gp);
1851  }
1852  pha_mat_lab(2, 2) = interp(
1853  itw,
1854  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 10),
1855  za_sca_gp,
1856  delta_aa_gp,
1857  za_inc_gp);
1858  if (stokes_dim == 3) {
1859  break;
1860  }
1861  if (delta_aa >= 0) {
1862  pha_mat_lab(0, 3) = interp(
1863  itw,
1864  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1865  za_sca_gp,
1866  delta_aa_gp,
1867  za_inc_gp);
1868  pha_mat_lab(1, 3) = interp(
1869  itw,
1870  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1871  za_sca_gp,
1872  delta_aa_gp,
1873  za_inc_gp);
1874  pha_mat_lab(3, 0) = interp(
1875  itw,
1876  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1877  za_sca_gp,
1878  delta_aa_gp,
1879  za_inc_gp);
1880  pha_mat_lab(3, 1) = interp(
1881  itw,
1882  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1883  za_sca_gp,
1884  delta_aa_gp,
1885  za_inc_gp);
1886  } else {
1887  pha_mat_lab(0, 3) = -interp(
1888  itw,
1889  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 3),
1890  za_sca_gp,
1891  delta_aa_gp,
1892  za_inc_gp);
1893  pha_mat_lab(1, 3) = -interp(
1894  itw,
1895  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 7),
1896  za_sca_gp,
1897  delta_aa_gp,
1898  za_inc_gp);
1899  pha_mat_lab(3, 0) = -interp(
1900  itw,
1901  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 12),
1902  za_sca_gp,
1903  delta_aa_gp,
1904  za_inc_gp);
1905  pha_mat_lab(3, 1) = -interp(
1906  itw,
1907  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 13),
1908  za_sca_gp,
1909  delta_aa_gp,
1910  za_inc_gp);
1911  }
1912  pha_mat_lab(2, 3) = interp(
1913  itw,
1914  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 11),
1915  za_sca_gp,
1916  delta_aa_gp,
1917  za_inc_gp);
1918  pha_mat_lab(3, 2) = interp(
1919  itw,
1920  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 14),
1921  za_sca_gp,
1922  delta_aa_gp,
1923  za_inc_gp);
1924  pha_mat_lab(3, 3) = interp(
1925  itw,
1926  pha_mat_data(Range(joker), Range(joker), Range(joker), 0, 15),
1927  za_sca_gp,
1928  delta_aa_gp,
1929  za_inc_gp);
1930  break;
1931  }
1932 
1933  default: {
1934  CREATE_OUT0;
1935  out0 << "Not all ptype cases are implemented\n";
1936  }
1937  }
1938 }
1939 
1941 
1969 void ext_matFromabs_vec( //Output
1971  //Input
1973  const Index& stokes_dim) {
1974  assert(stokes_dim >= 1 && stokes_dim <= 4);
1975  assert(ext_mat.nrows() == stokes_dim);
1976  assert(ext_mat.ncols() == stokes_dim);
1977  assert(abs_vec.nelem() == stokes_dim);
1978 
1979  // first: diagonal elements
1980  for (Index is = 0; is < stokes_dim; is++) {
1981  ext_mat(is, is) += abs_vec[0];
1982  }
1983  // second: off-diagonal elements, namely first row and column
1984  for (Index is = 1; is < stokes_dim; is++) {
1985  ext_mat(0, is) += abs_vec[is];
1986  ext_mat(is, 0) += abs_vec[is];
1987  }
1988 }
1989 
1991 
2009  VectorView t_ok,
2010  Index& this_T_interp_order,
2011  ArrayOfGridPosPoly& T_gp,
2012  Matrix& T_itw,
2013  //Input
2014  ConstVectorView T_grid,
2015  const Vector& T_array,
2016  const Index& t_interp_order) {
2017  const Index nTin = T_grid.nelem();
2018  const Index nTout = T_array.nelem();
2019 
2020  this_T_interp_order = -1;
2021 
2022  if (nTin > 1) {
2023  this_T_interp_order = min(t_interp_order, nTin - 1);
2024  T_itw.resize(nTout, this_T_interp_order + 1);
2025 
2026  // we need to check T-grid exceedance. and catch these cases (because T
2027  // is assumed to correspond to a location and T-exceedance is ok when pnd=0
2028  // there. however, here we do not know pnd.) and report them to have the
2029  // calling method dealing with it.
2030 
2031  // we set the extrapolfax explicitly and use it here as well as in
2032  // gridpos_poly below.
2033  const Numeric extrapolfac = 0.5;
2034  const Numeric lowlim = T_grid[0] - extrapolfac * (T_grid[1] - T_grid[0]);
2035  const Numeric uplim =
2036  T_grid[nTin - 1] + extrapolfac * (T_grid[nTin - 1] - T_grid[nTin - 2]);
2037 
2038  bool any_T_exceed = false;
2039  for (Index Tind = 0; Tind < nTout; Tind++) {
2040  if (T_array[Tind] < lowlim || T_array[Tind] > uplim) {
2041  t_ok[Tind] = -1.;
2042  any_T_exceed = true;
2043  } else
2044  t_ok[Tind] = 1.;
2045  }
2046 
2047  if (any_T_exceed) {
2048  GridPosPoly dummy_gp;
2049  dummy_gp.idx.resize(this_T_interp_order + 1);
2050  dummy_gp.w.resize(this_T_interp_order + 1);
2051  for (Index i = 0; i <= this_T_interp_order; ++i) dummy_gp.idx[i] = i;
2052  dummy_gp.w = 0.;
2053  dummy_gp.w[0] = 1.;
2054  bool grid_unchecked = true;
2055 
2056  for (Index iT = 0; iT < nTout; iT++) {
2057  if (t_ok[iT] < 0)
2058  // setup T-exceed points with dummy grid positions
2059  T_gp[iT] = dummy_gp;
2060  else {
2061  if (grid_unchecked) {
2063  "Temperature interpolation in pha_mat_1ScatElem",
2064  T_grid,
2065  T_array[Range(iT, 1)],
2066  this_T_interp_order);
2067  grid_unchecked = false;
2068  }
2069  gridpos_poly(
2070  T_gp[iT], T_grid, T_array[iT], this_T_interp_order, extrapolfac);
2071  }
2072  }
2073  } else
2074  gridpos_poly(T_gp, T_grid, T_array, this_T_interp_order, extrapolfac);
2075  interpweights(T_itw, T_gp);
2076  } else {
2077  t_ok = 1.;
2078  }
2079 }
2080 
2082 
2096  const Numeric& aa_sca,
2097  const Numeric& za_inc,
2098  const Numeric& aa_inc) {
2099  Numeric theta_rad;
2100  Numeric ANG_TOL = 1e-7;
2101 
2102  // CPD 5/10/03.
2103  // Two special cases are implemented here to avoid NaNs that can sometimes
2104  // occur in in the acos() formula in forward and backscattering cases.
2105  //
2106  // GH 2011-05-31
2107  // Consider not only aa_sca-aa_inc ~= 0, but also aa_sca-aa_inc ~= 360.
2108 
2109  if ((abs(aa_sca - aa_inc) < ANG_TOL) ||
2110  (abs(abs(aa_sca - aa_inc) - 360) < ANG_TOL)) {
2111  theta_rad = DEG2RAD * abs(za_sca - za_inc);
2112  } else if (abs(abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
2113  theta_rad = DEG2RAD * (za_sca + za_inc);
2114  if (theta_rad > PI) {
2115  theta_rad = 2 * PI - theta_rad;
2116  }
2117  } else {
2118  const Numeric za_sca_rad = za_sca * DEG2RAD;
2119  const Numeric za_inc_rad = za_inc * DEG2RAD;
2120  const Numeric aa_sca_rad = aa_sca * DEG2RAD;
2121  const Numeric aa_inc_rad = aa_inc * DEG2RAD;
2122 
2123  theta_rad =
2124  acos(cos(za_sca_rad) * cos(za_inc_rad) +
2125  sin(za_sca_rad) * sin(za_inc_rad) * cos(aa_sca_rad - aa_inc_rad));
2126  }
2127  return theta_rad;
2128 }
2129 
2131 
2154 void interpolate_scat_angle( //Output:
2155  VectorView pha_mat_int,
2156  //Input:
2157  ConstTensor5View pha_mat_data,
2158  ConstVectorView za_datagrid,
2159  const Numeric theta) {
2160  GridPos thet_gp;
2161  gridpos(thet_gp, za_datagrid, theta);
2162 
2163  Vector itw(2);
2164  interpweights(itw, thet_gp);
2165 
2166  assert(pha_mat_data.ncols() == 6);
2167  for (Index i = 0; i < 6; i++) {
2168  pha_mat_int[i] = interp(itw, pha_mat_data(joker, 0, 0, 0, i), thet_gp);
2169  }
2170 }
2171 
2173 
2198 void pha_mat_labCalc( //Output:
2199  MatrixView pha_mat_lab,
2200  //Input:
2201  ConstVectorView pha_mat_int,
2202  const Numeric& za_sca,
2203  const Numeric& aa_sca,
2204  const Numeric& za_inc,
2205  const Numeric& aa_inc,
2206  const Numeric& theta_rad) {
2207  const Index stokes_dim = pha_mat_lab.ncols();
2208 
2209  if (std::isnan(F11)) {
2210  throw runtime_error(
2211  "NaN value(s) detected in *pha_mat_labCalc* (0,0). Could the "
2212  "input data contain NaNs? Please check with *scat_dataCheck*. If "
2213  "input data are OK and you critically need the ongoing calculations, "
2214  "try to change the observation LOS slightly. If you can reproduce "
2215  "this error, please contact Patrick in order to help tracking down "
2216  "the reason to this problem. If you see this message occasionally "
2217  "when doing MC calculations, it should not be critical. This path "
2218  "sampling will be rejected and replaced with a new one.");
2219  }
2220 
2221  // For stokes_dim = 1, we only need Z11=F11:
2222  pha_mat_lab(0, 0) = F11;
2223 
2224  if (stokes_dim > 1) {
2225  Numeric za_sca_rad = za_sca * DEG2RAD;
2226  Numeric za_inc_rad = za_inc * DEG2RAD;
2227  Numeric aa_sca_rad = aa_sca * DEG2RAD;
2228  Numeric aa_inc_rad = aa_inc * DEG2RAD;
2229 
2230  const Numeric ANGTOL_RAD = 1e-6; //CPD: this constant is used to adjust
2231  //zenith angles close to 0 and PI. This is
2232  //also used to avoid float == float statements.
2233 
2234  //
2235  // Several cases have to be considered:
2236  //
2237 
2238  if ((abs(theta_rad) < ANGTOL_RAD) // forward scattering
2239  || (abs(theta_rad - PI) < ANGTOL_RAD) // backward scattering
2240  ||
2241  (abs(aa_inc_rad - aa_sca_rad) < ANGTOL_RAD) // inc and sca on meridian
2242  || (abs(abs(aa_inc_rad - aa_sca_rad) - 360.) < ANGTOL_RAD) // "
2243  || (abs(abs(aa_inc_rad - aa_sca_rad) - 180.) < ANGTOL_RAD) // "
2244  ) {
2245  pha_mat_lab(0, 1) = F12;
2246  pha_mat_lab(1, 0) = F12;
2247  pha_mat_lab(1, 1) = F22;
2248 
2249  if (stokes_dim > 2) {
2250  pha_mat_lab(0, 2) = 0;
2251  pha_mat_lab(1, 2) = 0;
2252  pha_mat_lab(2, 0) = 0;
2253  pha_mat_lab(2, 1) = 0;
2254  pha_mat_lab(2, 2) = F33;
2255 
2256  if (stokes_dim > 3) {
2257  pha_mat_lab(0, 3) = 0;
2258  pha_mat_lab(1, 3) = 0;
2259  pha_mat_lab(2, 3) = F34;
2260  pha_mat_lab(3, 0) = 0;
2261  pha_mat_lab(3, 1) = 0;
2262  pha_mat_lab(3, 2) = -F34;
2263  pha_mat_lab(3, 3) = F44;
2264  }
2265  }
2266  }
2267 
2268  else {
2269  Numeric sigma1;
2270  Numeric sigma2;
2271 
2272  Numeric s1, s2;
2273 
2274  // In these cases we have to take limiting values.
2275 
2276  if (za_inc_rad < ANGTOL_RAD) {
2277  sigma1 = PI + aa_sca_rad - aa_inc_rad;
2278  sigma2 = 0;
2279  } else if (za_inc_rad > PI - ANGTOL_RAD) {
2280  sigma1 = aa_sca_rad - aa_inc_rad;
2281  sigma2 = PI;
2282  } else if (za_sca_rad < ANGTOL_RAD) {
2283  sigma1 = 0;
2284  sigma2 = PI + aa_sca_rad - aa_inc_rad;
2285  } else if (za_sca_rad > PI - ANGTOL_RAD) {
2286  sigma1 = PI;
2287  sigma2 = aa_sca_rad - aa_inc_rad;
2288  } else {
2289  s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad)) /
2290  (sin(za_inc_rad) * sin(theta_rad));
2291  s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos(theta_rad)) /
2292  (sin(za_sca_rad) * sin(theta_rad));
2293 
2294  sigma1 = acos(s1);
2295  sigma2 = acos(s2);
2296 
2297  // Arccos is only defined in the range from -1 ... 1
2298  // Numerical problems can appear for values close to 1 or -1
2299  // this (also) catches the case when inc and sca are on one meridian
2300  if (std::isnan(sigma1) || std::isnan(sigma2)) {
2301  if (abs(s1 - 1) < ANGTOL_RAD) sigma1 = 0;
2302  if (abs(s1 + 1) < ANGTOL_RAD) sigma1 = PI;
2303  if (abs(s2 - 1) < ANGTOL_RAD) sigma2 = 0;
2304  if (abs(s2 + 1) < ANGTOL_RAD) sigma2 = PI;
2305  }
2306  }
2307 
2308  const Numeric C1 = cos(2 * sigma1);
2309  const Numeric C2 = cos(2 * sigma2);
2310 
2311  const Numeric S1 = sin(2 * sigma1);
2312  const Numeric S2 = sin(2 * sigma2);
2313 
2314  pha_mat_lab(0, 1) = C1 * F12;
2315  pha_mat_lab(1, 0) = C2 * F12;
2316  pha_mat_lab(1, 1) = C1 * C2 * F22 - S1 * S2 * F33;
2317 
2318  //assert(!std::isnan(pha_mat_lab(0,1)));
2319  //assert(!std::isnan(pha_mat_lab(1,0)));
2320  //assert(!std::isnan(pha_mat_lab(1,1)));
2321  if (std::isnan(pha_mat_lab(0, 1)) || std::isnan(pha_mat_lab(1, 0)) ||
2322  std::isnan(pha_mat_lab(1, 1))) {
2323  throw runtime_error(
2324  "NaN value(s) detected in *pha_mat_labCalc* (0/1,1). Could the "
2325  "input data contain NaNs? Please check with *scat_dataCheck*. If "
2326  "input data are OK and you critically need the ongoing calculations, "
2327  "try to change the observation LOS slightly. If you can reproduce "
2328  "this error, please contact Patrick in order to help tracking down "
2329  "the reason to this problem. If you see this message occasionally "
2330  "when doing MC calculations, it should not be critical. This path "
2331  "sampling will be rejected and replaced with a new one.");
2332  }
2333 
2334  if (stokes_dim > 2) {
2335  /*CPD: For skokes_dim > 2 some of the transformation formula
2336  for each element have a different sign depending on whether or
2337  not 0<aa_scat-aa_inc<180. For details see pages 94 and 95 of
2338  Mishchenkos chapter in :
2339  Mishchenko, M. I., and L. D. Travis, 2003: Electromagnetic
2340  scattering by nonspherical particles. In Exploring the Atmosphere
2341  by Remote Sensing Techniques (R. Guzzi, Ed.), Springer-Verlag,
2342  Berlin, pp. 77-127.
2343  This is available at http://www.giss.nasa.gov/~crmim/publications/ */
2344  Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
2345  (aa_sca - aa_inc > 180) * 360;
2346  if (delta_aa >= 0) {
2347  pha_mat_lab(0, 2) = S1 * F12;
2348  pha_mat_lab(1, 2) = S1 * C2 * F22 + C1 * S2 * F33;
2349  pha_mat_lab(2, 0) = -S2 * F12;
2350  pha_mat_lab(2, 1) = -C1 * S2 * F22 - S1 * C2 * F33;
2351  } else {
2352  pha_mat_lab(0, 2) = -S1 * F12;
2353  pha_mat_lab(1, 2) = -S1 * C2 * F22 - C1 * S2 * F33;
2354  pha_mat_lab(2, 0) = S2 * F12;
2355  pha_mat_lab(2, 1) = C1 * S2 * F22 + S1 * C2 * F33;
2356  }
2357  pha_mat_lab(2, 2) = -S1 * S2 * F22 + C1 * C2 * F33;
2358 
2359  if (stokes_dim > 3) {
2360  if (delta_aa >= 0) {
2361  pha_mat_lab(1, 3) = S2 * F34;
2362  pha_mat_lab(3, 1) = S1 * F34;
2363  } else {
2364  pha_mat_lab(1, 3) = -S2 * F34;
2365  pha_mat_lab(3, 1) = -S1 * F34;
2366  }
2367  pha_mat_lab(0, 3) = 0;
2368  pha_mat_lab(2, 3) = C2 * F34;
2369  pha_mat_lab(3, 0) = 0;
2370  pha_mat_lab(3, 2) = -C1 * F34;
2371  pha_mat_lab(3, 3) = F44;
2372  }
2373  }
2374  }
2375  }
2376 }
2377 
2378 ostream& operator<<(ostream& os, const SingleScatteringData& /*ssd*/) {
2379  os << "SingleScatteringData: Output operator not implemented";
2380  return os;
2381 }
2382 
2383 ostream& operator<<(ostream& os, const ArrayOfSingleScatteringData& /*assd*/) {
2384  os << "ArrayOfSingleScatteringData: Output operator not implemented";
2385  return os;
2386 }
2387 
2388 ostream& operator<<(ostream& os, const ScatteringMetaData& /*ssd*/) {
2389  os << "ScatteringMetaData: Output operator not implemented";
2390  return os;
2391 }
2392 
2393 ostream& operator<<(ostream& os, const ArrayOfScatteringMetaData& /*assd*/) {
2394  os << "ArrayOfScatteringMetaData: Output operator not implemented";
2395  return os;
2396 }
2397 
2399 
2417  //Input:
2419  Index stokes_dim, freq_dim;
2420 
2421  if (propmat_clearsky.nelem() > 1) {
2422  stokes_dim = propmat_clearsky[0].StokesDimensions();
2423  freq_dim = propmat_clearsky[0].NumberOfFrequencies();
2424  } else {
2425  stokes_dim = 0;
2426  freq_dim = 0;
2427  }
2428 
2429  // old abs_vecInit
2430  abs_vec = StokesVector(freq_dim, stokes_dim);
2431  abs_vec.SetZero();
2432 
2433  ext_mat = PropagationMatrix(freq_dim, stokes_dim);
2434  ext_mat.SetZero();
2435 
2436  // old ext_matAddGas and abs_vecAddGas for 0 vector and matrix
2437  for (auto& pm : propmat_clearsky) {
2438  abs_vec += pm;
2439  ext_mat += pm;
2440  }
2441 }
2442 
2444 
2454 PType PTypeFromString(const String& ptype_string) {
2455  PType ptype;
2456  if (ptype_string == "general")
2457  ptype = PTYPE_GENERAL;
2458  else if (ptype_string == "totally_random")
2459  ptype = PTYPE_TOTAL_RND;
2460  else if (ptype_string == "azimuthally_random")
2461  ptype = PTYPE_AZIMUTH_RND;
2462  else {
2463  ostringstream os;
2464  os << "Unknown ptype: " << ptype_string << endl
2465  << "Valid types are: general, totally_random and "
2466  << "azimuthally_random.";
2467  throw std::runtime_error(os.str());
2468  }
2469 
2470  return ptype;
2471 }
2472 
2474 
2484 PType PType2FromString(const String& ptype_string) {
2485  PType ptype;
2486  if (ptype_string == "general")
2487  ptype = PTYPE_GENERAL;
2488  else if (ptype_string == "macroscopically_isotropic")
2489  ptype = PTYPE_TOTAL_RND;
2490  else if (ptype_string == "horizontally_aligned")
2491  ptype = PTYPE_AZIMUTH_RND;
2492  else {
2493  ostringstream os;
2494  os << "Unknown ptype: " << ptype_string << endl
2495  << "Valid types are: general, macroscopically_isotropic and "
2496  << "horizontally_aligned.";
2497  throw std::runtime_error(os.str());
2498  }
2499 
2500  return ptype;
2501 }
2502 
2504 
2512 String PTypeToString(const PType& ptype) {
2513  String ptype_string;
2514 
2515  switch (ptype) {
2516  case PTYPE_GENERAL:
2517  ptype_string = "general";
2518  break;
2519  case PTYPE_TOTAL_RND:
2520  ptype_string = "totally_random";
2521  break;
2522  case PTYPE_AZIMUTH_RND:
2523  ptype_string = "azimuthally_random";
2524  break;
2525  default:
2526  ostringstream os;
2527  os << "Internal error: Cannot map PType enum value " << ptype
2528  << " to String.";
2529  throw std::runtime_error(os.str());
2530  break;
2531  }
2532 
2533  return ptype_string;
2534 }
2535 
2537 
2545  // First check that input fulfills requirements on older data formats:
2546  // 1) Is za_grid symmetric and includes 90deg?
2547  Index nza = ssd.za_grid.nelem();
2548  for (Index i = 0; i < nza / 2; i++) {
2550  180. - ssd.za_grid[nza - 1 - i], ssd.za_grid[i], 2 * DBL_EPSILON)) {
2551  ostringstream os;
2552  os << "Zenith grid of azimuthally_random single scattering data\n"
2553  << "is not symmetric with respect to 90degree.";
2554  throw std::runtime_error(os.str());
2555  }
2556  }
2557  if (!is_same_within_epsilon(ssd.za_grid[nza / 2], 90., 2 * DBL_EPSILON)) {
2558  ostringstream os;
2559  os << "Zenith grid of azimuthally_random single scattering data\n"
2560  << "does not contain 90 degree grid point.";
2561  throw std::runtime_error(os.str());
2562  }
2563 
2564  // 2) Are data sizes correct?
2565  ostringstream os_pha_mat;
2566  os_pha_mat << "pha_mat ";
2567  ostringstream os_ext_mat;
2568  os_ext_mat << "ext_mat ";
2569  ostringstream os_abs_vec;
2570  os_abs_vec << "abs_vec ";
2571  chk_size(os_pha_mat.str(),
2572  ssd.pha_mat_data,
2573  ssd.f_grid.nelem(),
2574  ssd.T_grid.nelem(),
2575  ssd.za_grid.nelem(),
2576  ssd.aa_grid.nelem(),
2577  ssd.za_grid.nelem() / 2 + 1,
2578  1,
2579  16);
2580 
2581  chk_size(os_ext_mat.str(),
2582  ssd.ext_mat_data,
2583  ssd.f_grid.nelem(),
2584  ssd.T_grid.nelem(),
2585  ssd.za_grid.nelem() / 2 + 1,
2586  1,
2587  3);
2588 
2589  chk_size(os_abs_vec.str(),
2590  ssd.abs_vec_data,
2591  ssd.f_grid.nelem(),
2592  ssd.T_grid.nelem(),
2593  ssd.za_grid.nelem() / 2 + 1,
2594  1,
2595  2);
2596 
2597  // Now that we are sure that za_grid is properly symmetric, we just need to
2598  // copy over the data (ie no interpolation).
2599  Tensor5 tmpT5 = ssd.abs_vec_data;
2600  ssd.abs_vec_data.resize(tmpT5.nshelves(),
2601  tmpT5.nbooks(),
2602  ssd.za_grid.nelem(),
2603  tmpT5.nrows(),
2604  tmpT5.ncols());
2605  ssd.abs_vec_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2606  for (Index i = 0; i < nza / 2; i++) {
2607  ssd.abs_vec_data(joker, joker, nza - 1 - i, joker, joker) =
2608  tmpT5(joker, joker, i, joker, joker);
2609  }
2610 
2611  tmpT5 = ssd.ext_mat_data;
2612  ssd.ext_mat_data.resize(tmpT5.nshelves(),
2613  tmpT5.nbooks(),
2614  ssd.za_grid.nelem(),
2615  tmpT5.nrows(),
2616  tmpT5.ncols());
2617  ssd.ext_mat_data(joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT5;
2618  for (Index i = 0; i < nza / 2; i++) {
2619  ssd.ext_mat_data(joker, joker, nza - 1 - i, joker, joker) =
2620  tmpT5(joker, joker, i, joker, joker);
2621  }
2622 
2623  Tensor7 tmpT7 = ssd.pha_mat_data;
2624  ssd.pha_mat_data.resize(tmpT7.nlibraries(),
2625  tmpT7.nvitrines(),
2626  tmpT7.nshelves(),
2627  tmpT7.nbooks(),
2628  ssd.za_grid.nelem(),
2629  tmpT7.nrows(),
2630  tmpT7.ncols());
2631  ssd.pha_mat_data(
2632  joker, joker, joker, joker, Range(0, nza / 2 + 1), joker, joker) = tmpT7;
2633 
2634  // scatt. matrix elements 13,23,31,32 and 14,24,41,42 (=elements 2,6,8,9 and
2635  // 3,7,12,13 in ARTS' flattened format, respectively) change sign.
2636  tmpT7(joker, joker, joker, joker, joker, joker, Range(2, 2)) *= -1.;
2637  tmpT7(joker, joker, joker, joker, joker, joker, Range(6, 4)) *= -1.;
2638  tmpT7(joker, joker, joker, joker, joker, joker, Range(12, 2)) *= -1.;
2639 
2640  // For second half of incident polar angles (>90deg), we need to mirror the
2641  // original data in both incident and scattered polar angle around 90deg "planes".
2642  for (Index i = 0; i < nza / 2; i++)
2643  for (Index j = 0; j < nza; j++)
2644  ssd.pha_mat_data(
2645  joker, joker, nza - 1 - j, joker, nza - 1 - i, joker, joker) =
2646  tmpT7(joker, joker, j, joker, i, joker, joker);
2647 }
2648 
2650 
2659  const String& particle_ssdmethod_string) {
2660  ParticleSSDMethod particle_ssdmethod;
2661  if (particle_ssdmethod_string == "tmatrix")
2662  particle_ssdmethod = PARTICLE_SSDMETHOD_TMATRIX;
2663  else {
2664  ostringstream os;
2665  os << "Unknown particle SSD method: " << particle_ssdmethod_string << endl
2666  << "Valid methods: tmatrix";
2667  throw std::runtime_error(os.str());
2668  }
2669 
2670  return particle_ssdmethod;
2671 }
2672 
2674 
2682 String PTypeToString(const ParticleSSDMethod& particle_ssdmethod) {
2683  String particle_ssdmethod_string;
2684 
2685  switch (particle_ssdmethod) {
2687  particle_ssdmethod_string = "tmatrix";
2688  break;
2689  default:
2690  ostringstream os;
2691  os << "Internal error: Cannot map ParticleSSDMethod enum value "
2692  << particle_ssdmethod << " to String.";
2693  throw std::runtime_error(os.str());
2694  break;
2695  }
2696 
2697  return particle_ssdmethod_string;
2698 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::abs_vec
StokesVector abs_vec(Workspace &ws) noexcept
Definition: autoarts.h:2239
Tensor7View
The Tensor7View class.
Definition: matpackVII.h:1286
is_same_within_epsilon
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
pha_mat_1ScatElem
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.
Definition: optproperties.cc:881
RAD2DEG
const Numeric RAD2DEG
F11
#define F11
Definition: optproperties.cc:55
ConstTensor7View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:48
ARTS::Group::PropagationMatrix
PropagationMatrix PropagationMatrix
Definition: autoarts.h:93
ConstTensor5View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackV.cc:47
SingleScatteringData::za_grid
Vector za_grid
Definition: optproperties.h:85
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
SingleScatteringData::f_grid
Vector f_grid
Definition: optproperties.h:83
PTYPE_GENERAL
@ PTYPE_GENERAL
Definition: optproperties.h:54
PTypeFromString
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
Definition: optproperties.cc:2454
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
ConvertAzimuthallyRandomSingleScatteringData
void ConvertAzimuthallyRandomSingleScatteringData(SingleScatteringData &ssd)
Convert azimuthally-random oriented SingleScatteringData to latest version.
Definition: optproperties.cc:2544
F22
#define F22
Definition: optproperties.cc:57
ConstTensor5View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:56
ConstTensor7View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:60
PARTICLE_SSDMETHOD_TMATRIX
@ PARTICLE_SSDMETHOD_TMATRIX
Definition: optproperties.h:66
joker
const Joker joker
abs_vecTransform
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.
Definition: optproperties.cc:1474
interpolation.h
Header file for interpolation.cc.
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
ext_matTransform
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.
Definition: optproperties.cc:1572
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
F34
#define F34
Definition: optproperties.cc:59
chk_interpolation_grids
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:989
ARTS::Var::f_index
Index f_index(Workspace &ws) noexcept
Definition: autoarts.h:3464
SingleScatteringData::aa_grid
Vector aa_grid
Definition: optproperties.h:86
PropagationMatrix
Definition: propagationmatrix.h:87
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
opt_prop_ScatSpecBulk
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.
Definition: optproperties.cc:148
FouComp_1ScatElem
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.
Definition: optproperties.cc:1218
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
ARTS::Var::nshelves
Index nshelves(Workspace &ws) noexcept
Definition: autoarts.h:4689
PType2FromString
PType PType2FromString(const String &ptype_string)
Convert ptype name to enum value.
Definition: optproperties.cc:2484
ext_mat_SSD2Stokes
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.
Definition: optproperties.cc:585
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
array.h
This file contains the definition of Array.
ConstTensor5View::npages
Index npages() const
Returns the number of pages.
Definition: matpackV.cc:50
scat_angle
Numeric scat_angle(const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc)
Calculates the scattering angle.
Definition: optproperties.cc:2095
SingleScatteringData
Definition: optproperties.h:80
ARTS::Var::pha_mat
Tensor4 pha_mat(Workspace &ws) noexcept
Definition: autoarts.h:4881
DEG2RAD
const Numeric DEG2RAD
ConstTensor7View::nlibraries
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:42
opt_prop_NScatElems
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.
Definition: optproperties.cc:248
Array
This can be used to make arrays out of anything.
Definition: array.h:108
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
PTypeToString
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
Definition: optproperties.cc:2512
abs_vec_SSD2Stokes
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.
Definition: optproperties.cc:633
min
#define min(a, b)
Definition: legacy_continua.cc:20628
Tensor7::resize
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5484
SingleScatteringData::ptype
PType ptype
Definition: optproperties.h:81
PType
PType
An attribute to classify the particle type (ptype) of a SingleScatteringData.
Definition: optproperties.h:53
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:204
Tensor5::resize
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1743
pha_mat_Bulk
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
Definition: optproperties.cc:670
messages.h
Declarations having to do with the four output streams.
PropagationMatrix::K12
VectorView K12(const Index iz=0, const Index ia=0)
Vector view to K(0, 1) elements.
Definition: propagationmatrix.h:809
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
GridPosPoly
Structure to store a grid position for higher order interpolation.
Definition: interpolation_poly.h:64
_U_
#define _U_
Definition: config.h:183
my_basic_string< char >
ScatteringMetaData
Definition: optproperties.h:109
SingleScatteringData::ext_mat_data
Tensor5 ext_mat_data
Definition: optproperties.h:88
PropagationMatrix::K34
VectorView K34(const Index iz=0, const Index ia=0)
Vector view to K(2, 3) elements.
Definition: propagationmatrix.h:859
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
ConstTensor7View::nvitrines
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:45
PropagationMatrix::NumberOfFrequencies
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Definition: propagationmatrix.h:204
optproperties.h
Scattering database structure and functions.
GridPosPoly::idx
ArrayOfIndex idx
Definition: interpolation_poly.h:67
ConstTensor7View::npages
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:54
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
interpolate_scat_angle
void interpolate_scat_angle(VectorView pha_mat_int, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, const Numeric theta)
Interpolate data on the scattering angle.
Definition: optproperties.cc:2154
Tensor5View
The Tensor5View class.
Definition: matpackV.h:333
pha_mat_NScatElems
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.
Definition: optproperties.cc:796
math_funcs.h
ARTS::Var::npages
Index npages(Workspace &ws) noexcept
Definition: autoarts.h:4675
ARTS::Var::ext_mat
PropagationMatrix ext_mat(Workspace &ws) noexcept
Definition: autoarts.h:3390
max
#define max(a, b)
Definition: legacy_continua.cc:20629
opt_prop_1ScatElem
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.
Definition: optproperties.cc:335
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
ssd_tinterp_parameters
void ssd_tinterp_parameters(VectorView t_ok, Index &this_T_interp_order, ArrayOfGridPosPoly &T_gp, Matrix &T_itw, ConstVectorView T_grid, const Vector &T_array, const Index &t_interp_order)
Determine T-interpol parameters for a specific scattering element.
Definition: optproperties.cc:2008
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
PTYPE_AZIMUTH_RND
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:55
pha_matTransform
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.
Definition: optproperties.cc:1692
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
ParticleSSDMethod
ParticleSSDMethod
An attribute to classify the method to be used for SingleScatteringData.
Definition: optproperties.h:64
PropagationMatrix::SetZero
void SetZero()
Sets all data to zero.
Definition: propagationmatrix.h:934
ext_matFromabs_vec
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
Definition: optproperties.cc:1969
Range
The range class.
Definition: matpackI.h:160
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ParticleSSDMethodFromString
ParticleSSDMethod ParticleSSDMethodFromString(const String &particle_ssdmethod_string)
Convert particle ssd method name to enum value.
Definition: optproperties.cc:2658
ConstTensor7View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:51
ConstTensor5View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackV.cc:53
opt_prop_sum_propmat_clearsky
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const ArrayOfPropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
Definition: optproperties.cc:2414
logic.h
Header file for logic.cc.
opt_prop_Bulk
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
Definition: optproperties.cc:102
pha_mat_labCalc
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.
Definition: optproperties.cc:2198
F33
#define F33
Definition: optproperties.cc:58
ConstTensor5View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackV.cc:44
PropagationMatrix::StokesDimensions
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Definition: propagationmatrix.h:201
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:799
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:284
SingleScatteringData::abs_vec_data
Tensor5 abs_vec_data
Definition: optproperties.h:89
operator<<
ostream & operator<<(ostream &os, const SingleScatteringData &)
Definition: optproperties.cc:2378
F12
#define F12
Definition: optproperties.cc:56
GridPosPoly::w
Vector w
Definition: interpolation_poly.h:70
F44
#define F44
Definition: optproperties.cc:60
chk_size
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:466
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstTensor7View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:57
PI
const Numeric PI
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
check_input.h
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
Vector
The Vector class.
Definition: matpackI.h:860
SingleScatteringData::pha_mat_data
Tensor7 pha_mat_data
Definition: optproperties.h:87
ARTS::Var::propmat_clearsky
ArrayOfPropagationMatrix propmat_clearsky(Workspace &ws) noexcept
Definition: autoarts.h:5398
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
PTYPE_TOTAL_RND
@ PTYPE_TOTAL_RND
Definition: optproperties.h:56
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
SingleScatteringData::T_grid
Vector T_grid
Definition: optproperties.h:84
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Group::StokesVector
StokesVector StokesVector
Definition: autoarts.h:101
ARTS::Var::nbooks
Index nbooks(Workspace &ws) noexcept
Definition: autoarts.h:4539
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:143
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2382
matpackVII.h
Tensor6View
The Tensor6View class.
Definition: matpackVI.h:621
pha_mat_ScatSpecBulk
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.
Definition: optproperties.cc:708
arts.h
The global header file for ARTS.
xml_io.h
This file contains basic functions to handle XML data files.