58#define F11 pha_mat_int[0]
59#define F12 pha_mat_int[1]
60#define F22 pha_mat_int[2]
61#define F33 pha_mat_int[3]
62#define F34 pha_mat_int[4]
63#define F44 pha_mat_int[5]
115 ext_mat = ext_mat_ss[0];
116 abs_vec = abs_vec_ss[0];
118 for (
Index i_ss = 1; i_ss < ext_mat_ss.
nelem(); i_ss++) {
119 ext_mat += ext_mat_ss[i_ss];
120 abs_vec += abs_vec_ss[i_ss];
122 ptype =
max(ptypes_ss);
168 const Index nf = abs_vec_se[0][0].nbooks();
169 const Index nDir = abs_vec_se[0][0].nrows();
170 const Index stokes_dim = abs_vec_se[0][0].ncols();
181 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
182 ARTS_ASSERT(ext_mat_se[i_ss].nelem() == abs_vec_se[i_ss].nelem());
186 ext_mat[i_ss].resize(nf, nT, nDir, stokes_dim, stokes_dim);
188 abs_vec[i_ss].resize(nf, nT, nDir, stokes_dim);
191 for (
Index i_se = 0; i_se < ext_mat_se[i_ss].
nelem(); i_se++) {
192 ARTS_ASSERT(nT == ext_mat_se[i_ss][i_se].nbooks());
193 ARTS_ASSERT(nT == abs_vec_se[i_ss][i_se].npages());
195 for (
Index Tind = 0; Tind < nT; Tind++) {
196 if (pnds(i_se_flat, Tind) != 0.) {
197 if (t_ok(i_se_flat, Tind) > 0.) {
199 ext_tmp *= pnds(i_se_flat, Tind);
203 abs_tmp *= pnds(i_se_flat, Tind);
207 "Interpolation error for (flat-array) scattering element #",
209 "at location/temperature point #", Tind,
"\n")
215 ptype[i_ss] =
max(ptypes_se[i_ss]);
257 const Index& stokes_dim,
260 const Index& f_index,
261 const Index& t_interp_order) {
264 nf = scat_data[0][0].ext_mat_data.nshelves();
269 if (scat_data[0][0].ext_mat_data.nshelves() == 1)
288 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
290 ext_mat[i_ss].resize(nse);
291 abs_vec[i_ss].resize(nse);
292 ptypes[i_ss].resize(nse);
294 for (
Index i_se = 0; i_se < nse; i_se++) {
295 ext_mat[i_ss][i_se].resize(nf, nT, nDir, stokes_dim, stokes_dim);
296 abs_vec[i_ss][i_se].resize(nf, nT, nDir, stokes_dim);
301 t_ok(i_se_flat,
joker),
302 scat_data[i_ss][i_se],
333 Index& this_T_interp_order,
337 const Index& t_interp_order) {
341 this_T_interp_order = -1;
345 this_T_interp_order =
min(t_interp_order, nTin - 1);
354 const Numeric extrapolfac = 0.5;
355 const Numeric lowlim = T_grid[0] - extrapolfac * (T_grid[1] - T_grid[0]);
357 T_grid[nTin - 1] + extrapolfac * (T_grid[nTin - 1] - T_grid[nTin - 2]);
359 bool any_T_exceed =
false;
360 for (
Index Tind = 0; Tind < nTout; Tind++) {
361 if (T_array[Tind] < lowlim || T_array[Tind] > uplim) {
371 T_lag.reserve(nTout);
373 bool grid_unchecked =
true;
375 for (
Index iT = 0; iT < nTout; iT++) {
377 T_lag.emplace_back();
379 if (grid_unchecked) {
381 "Temperature interpolation in pha_mat_1ScatElem",
383 T_array[
Range(iT, 1)],
384 this_T_interp_order);
385 grid_unchecked =
false;
387 T_lag.emplace_back(0, T_array[iT], T_grid, this_T_interp_order,
false);
433 const Index& f_start,
434 const Index& t_interp_order) {
479 Index this_T_interp_order;
491 if (this_T_interp_order < 0)
495 Tensor3 ext_mat_tmp(nf, stokes_dim, stokes_dim);
496 Matrix abs_vec_tmp(nf, stokes_dim);
497 for (
Index find = 0; find < nf; find++) {
507 for (
Index Tind = 0; Tind < nTout; Tind++)
508 for (
Index dind = 0; dind < nDir; dind++) {
510 abs_vec(
joker, Tind, dind,
joker) = abs_vec_tmp;
515 Tensor4 ext_mat_tmp(nf, nTout, stokes_dim, stokes_dim);
516 Tensor3 abs_vec_tmp(nf, nTout, stokes_dim);
519 for (
Index find = 0; find < nf; find++) {
520 for (
Index nst = 0; nst < ext_mat_tmp_ssd.
ncols(); nst++) {
521 reinterp(ext_mat_tmp_ssd(
joker, nst),
525 for (
Index Tind = 0; Tind < nTout; Tind++)
528 ext_mat_tmp_ssd(Tind,
joker),
534 for (
Index nst = 0; nst < abs_vec_tmp_ssd.
ncols(); nst++) {
535 reinterp(abs_vec_tmp_ssd(
joker, nst),
540 for (
Index Tind = 0; Tind < nTout; Tind++)
543 abs_vec_tmp_ssd(Tind,
joker),
547 abs_vec_tmp(find, Tind,
joker) = 0.;
550 for (
Index dind = 0; dind < nDir; dind++) {
574 if (this_T_interp_order < 0)
576 Matrix ext_mat_tmp_ssd(nDir, next);
577 Matrix abs_vec_tmp_ssd(nDir, nabs);
578 Tensor4 ext_mat_tmp(nf, nDir, stokes_dim, stokes_dim);
579 Tensor3 abs_vec_tmp(nf, nDir, stokes_dim);
580 for (
Index find = 0; find < nf; find++) {
581 for (
Index nst = 0; nst < next; nst++)
586 for (
Index Dind = 0; Dind < nDir; Dind++)
588 ext_mat_tmp_ssd(Dind,
joker),
592 for (
Index nst = 0; nst < nabs; nst++)
597 for (
Index Dind = 0; Dind < nDir; Dind++)
599 abs_vec_tmp_ssd(Dind,
joker),
604 for (
Index Tind = 0; Tind < nTout; Tind++) {
611 Tensor3 ext_mat_tmp_ssd(nTin, nDir, next);
612 Tensor3 abs_vec_tmp_ssd(nTin, nDir, nabs);
613 Matrix ext_mat_tmp(nTout, next);
614 Matrix abs_vec_tmp(nTout, nabs);
616 for (
Index find = 0; find < nf; find++) {
617 for (
Index Tind = 0; Tind < nTin; Tind++) {
618 for (
Index nst = 0; nst < next; nst++)
623 for (
Index nst = 0; nst < nabs; nst++)
630 for (
Index Dind = 0; Dind < nDir; Dind++) {
631 for (
Index nst = 0; nst < next; nst++) {
632 reinterp(ext_mat_tmp(
joker, nst),
633 ext_mat_tmp_ssd(
joker, Dind, nst),
637 for (
Index Tind = 0; Tind < nTout; Tind++)
639 ext_mat_tmp(Tind,
joker),
643 for (
Index nst = 0; nst < nabs; nst++) {
644 reinterp(abs_vec_tmp(
joker, nst),
645 abs_vec_tmp_ssd(
joker, Dind, nst),
649 for (
Index Tind = 0; Tind < nTout; Tind++)
651 abs_vec_tmp(Tind,
joker),
678 const Index& stokes_dim,
679 const Index& ptype) {
686 for (
Index ist = 0; ist < stokes_dim; ist++) {
687 ext_mat_stokes(ist, ist) = ext_mat_ssd[0];
691 switch (stokes_dim) {
693 ext_mat_stokes(2, 3) = ext_mat_ssd[2];
694 ext_mat_stokes(3, 2) = -ext_mat_ssd[2];
701 ext_mat_stokes(0, 1) = ext_mat_ssd[1];
702 ext_mat_stokes(1, 0) = ext_mat_ssd[1];
726 const Index& stokes_dim,
727 const Index& ptype) {
734 abs_vec_stokes[0] = abs_vec_ssd[0];
737 abs_vec_stokes[1] = abs_vec_ssd[1];
765 pha_mat = pha_mat_ss[0];
767 for (
Index i_ss = 1; i_ss < pha_mat_ss.
nelem(); i_ss++)
768 pha_mat += pha_mat_ss[i_ss];
770 ptype =
max(ptypes_ss);
811 const Index nf = pha_mat_se[0][0].nvitrines();
812 const Index npDir = pha_mat_se[0][0].nbooks();
813 const Index niDir = pha_mat_se[0][0].npages();
814 const Index stokes_dim = pha_mat_se[0][0].ncols();
823 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
826 pha_mat[i_ss].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
829 for (
Index i_se = 0; i_se < pha_mat_se[i_ss].
nelem(); i_se++) {
830 ARTS_ASSERT(nT == pha_mat_se[i_ss][i_se].nshelves());
832 for (
Index Tind = 0; Tind < nT; Tind++) {
833 if (pnds(i_se_flat, Tind) != 0.) {
834 if (t_ok(i_se_flat, Tind) > 0.) {
837 pha_tmp *= pnds(i_se_flat, Tind);
841 "Interpolation error for (flat-array) scattering element #",
843 "at location/temperature point #", Tind,
"\n")
849 ptype[i_ss] =
max(ptypes_se[i_ss]);
890 const Index& stokes_dim,
894 const Index& f_index,
895 const Index& t_interp_order) {
898 nf = scat_data[0][0].pha_mat_data.nlibraries();
903 if (scat_data[0][0].pha_mat_data.nlibraries() == 1)
922 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
924 pha_mat[i_ss].resize(nse);
925 ptypes[i_ss].resize(nse);
927 for (
Index i_se = 0; i_se < nse; i_se++) {
928 pha_mat[i_ss][i_se].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
932 t_ok(i_se_flat,
joker),
933 scat_data[i_ss][i_se],
978 const Index& f_start,
979 const Index& t_interp_order) {
1006 Index this_T_interp_order;
1008 this_T_interp_order,
1021 if (stokes_dim == 1)
1023 else if (stokes_dim < 4)
1027 if (this_T_interp_order < 0)
1031 for (
Index pdir = 0; pdir < npDir; pdir++)
1032 for (
Index idir = 0; idir < niDir; idir++) {
1035 pdir_array(pdir, 1),
1036 idir_array(idir, 0),
1037 idir_array(idir, 1));
1045 Vector pha_mat_int(npha, 0.);
1046 Matrix pha_mat_tmp(stokes_dim, stokes_dim);
1047 for (
Index find = 0; find < nf; find++) {
1049 for (
Index nst = 0; nst < npha; nst++)
1050 pha_mat_int[nst] =
interp(
1058 pdir_array(pdir, 0),
1059 pdir_array(pdir, 1),
1060 idir_array(idir, 0),
1061 idir_array(idir, 1),
1064 for (
Index Tind = 0; Tind < nTout; Tind++)
1065 pha_mat(find, Tind, pdir, idir,
joker,
joker) = pha_mat_tmp;
1070 for (
Index pdir = 0; pdir < npDir; pdir++)
1071 for (
Index idir = 0; idir < niDir; idir++) {
1074 pdir_array(pdir, 1),
1075 idir_array(idir, 0),
1076 idir_array(idir, 1));
1084 Matrix pha_mat_int(nTin, npha, 0.);
1085 Matrix pha_mat_tmp(nTout, npha, 0.);
1086 for (
Index find = 0; find < nf; find++) {
1087 for (
Index Tind = 0; Tind < nTin; Tind++)
1089 for (
Index nst = 0; nst < npha; nst++) {
1090 pha_mat_int(Tind, nst) =
interp(
1096 for (
Index nst = 0; nst < npha; nst++) {
1097 reinterp(pha_mat_tmp(
joker, nst),
1098 pha_mat_int(
joker, nst),
1107 for (
Index Tind = 0; Tind < nTout; Tind++) {
1110 pha_mat_tmp(Tind,
joker),
1111 pdir_array(pdir, 0),
1112 pdir_array(pdir, 1),
1113 idir_array(idir, 0),
1114 idir_array(idir, 1),
1124 Index nDir = npDir * niDir;
1126 Matrix delta_aa(npDir, niDir);
1134 for (
Index pdir = 0; pdir < npDir; pdir++) {
1135 for (
Index idir = 0; idir < niDir; idir++) {
1136 delta_aa(pdir, idir) =
1137 pdir_array(pdir, 1) - idir_array(idir, 1) +
1138 (pdir_array(pdir, 1) - idir_array(idir, 1) < -180) * 360 -
1139 (pdir_array(pdir, 1) - idir_array(idir, 1) > 180) * 360;
1140 adelta_aa[j] =
abs(delta_aa(pdir, idir));
1141 pza_gp[j] = pza_gp_tmp[pdir];
1142 iza_gp[j] = iza_gp_tmp[idir];
1152 if (this_T_interp_order < 0)
1154 Tensor3 pha_mat_int(nDir, stokes_dim, stokes_dim, 0.);
1155 Tensor4 pha_mat_tmp(npDir, niDir, stokes_dim, stokes_dim);
1157 for (
Index find = 0; find < nf; find++) {
1161 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1162 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1164 pha_mat_int(
joker, ist1, ist2),
1175 for (
Index pdir = 0; pdir < npDir; pdir++)
1176 for (
Index idir = 0; idir < niDir; idir++) {
1182 if (stokes_dim > 2) {
1183 for (
Index pdir = 0; pdir < npDir; pdir++)
1184 for (
Index idir = 0; idir < niDir; idir++)
1185 if (delta_aa(pdir, idir) < 0.) {
1186 pha_mat_tmp(pdir, idir, 0, 2) *= -1;
1187 pha_mat_tmp(pdir, idir, 1, 2) *= -1;
1188 pha_mat_tmp(pdir, idir, 2, 0) *= -1;
1189 pha_mat_tmp(pdir, idir, 2, 1) *= -1;
1193 if (stokes_dim > 3) {
1194 for (
Index pdir = 0; pdir < npDir; pdir++)
1195 for (
Index idir = 0; idir < niDir; idir++)
1196 if (delta_aa(pdir, idir) < 0.) {
1197 pha_mat_tmp(pdir, idir, 0, 3) *= -1;
1198 pha_mat_tmp(pdir, idir, 1, 3) *= -1;
1199 pha_mat_tmp(pdir, idir, 3, 0) *= -1;
1200 pha_mat_tmp(pdir, idir, 3, 1) *= -1;
1204 for (
Index Tind = 0; Tind < nTout; Tind++)
1212 Tensor4 pha_mat_int(nTin, nDir, stokes_dim, stokes_dim, 0.);
1214 for (
Index find = 0; find < nf; find++) {
1217 for (
Index Tind = 0; Tind < nTin; Tind++) {
1218 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1219 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1238 for (
Index pdir = 0; pdir < npDir; pdir++)
1239 for (
Index idir = 0; idir < niDir; idir++) {
1240 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1241 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1242 reinterp(pha_mat(find,
joker, pdir, idir, ist1, ist2),
1243 pha_mat_int(
joker, i, ist1, ist2),
1249 if (stokes_dim > 2) {
1250 for (
Index pdir = 0; pdir < npDir; pdir++)
1251 for (
Index idir = 0; idir < niDir; idir++)
1252 if (delta_aa(pdir, idir) < 0.) {
1253 pha_mat(find,
joker, pdir, idir, 0, 2) *= -1;
1254 pha_mat(find,
joker, pdir, idir, 1, 2) *= -1;
1255 pha_mat(find,
joker, pdir, idir, 2, 0) *= -1;
1256 pha_mat(find,
joker, pdir, idir, 2, 1) *= -1;
1260 if (stokes_dim > 3) {
1261 for (
Index pdir = 0; pdir < npDir; pdir++)
1262 for (
Index idir = 0; idir < niDir; idir++)
1263 if (delta_aa(pdir, idir) < 0.) {
1264 pha_mat(find,
joker, pdir, idir, 0, 3) *= -1;
1265 pha_mat(find,
joker, pdir, idir, 1, 3) *= -1;
1266 pha_mat(find,
joker, pdir, idir, 3, 0) *= -1;
1267 pha_mat(find,
joker, pdir, idir, 3, 1) *= -1;
1310 const Vector& pdir_array,
1311 const Vector& idir_array,
1312 const Index& f_start,
1313 const Index& t_interp_order,
1314 const Index& naa_totran) {
1331 const Index stokes_dim = pha_mat_fou.
nrows();
1348 Index this_T_interp_order;
1350 this_T_interp_order,
1362 "Azimuth grid size for scatt matrix extraction"
1363 " (*naa_totran*) must be >3.\n"
1364 "Yours is ", naa_totran,
".\n")
1368 1. / float(naa_totran - 1);
1369 Vector theta(naa_totran);
1371 Matrix theta_itw(naa_totran, 2);
1374 if (stokes_dim == 1)
1376 else if (stokes_dim < 4)
1381 Matrix pha_mat(stokes_dim, stokes_dim);
1382 Matrix pha_mat_angint(naa_totran, npha);
1383 Tensor3 Fou_int(nTin, stokes_dim, stokes_dim);
1385 for (
Index idir = 0; idir < niDir; idir++)
1386 for (
Index pdir = 0; pdir < npDir; pdir++) {
1387 for (
Index iaa = 0; iaa < naa_totran; iaa++)
1390 scat_angle(pdir_array[pdir], aa_grid[iaa], idir_array[idir], 0.);
1396 for (
Index find = 0; find < nf; find++) {
1398 for (
Index Tind = 0; Tind < nTin; Tind++) {
1400 for (
Index nst = 0; nst < npha; nst++)
1402 pha_mat_angint(
joker, nst),
1406 for (
Index iaa = 0; iaa < naa_totran; iaa++) {
1409 pha_mat_angint(iaa,
joker),
1418 if (iaa == 0 || iaa == naa_totran - 1)
1419 pha_mat *= (daa_totran / 2.);
1421 pha_mat *= daa_totran;
1426 if (this_T_interp_order <
1429 for (
Index Tind = 0; Tind < nTout; Tind++)
1430 pha_mat_fou(find, Tind, pdir, idir,
joker,
joker, 0) =
1433 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1434 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1435 for (
Index im = 0; im < nmodes; im++)
1436 reinterp(pha_mat_fou(find,
joker, pdir, idir, ist1, ist2, im),
1437 Fou_int(
joker, ist1, ist2),
1460 daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1461 for (
Index iaa = 1; iaa < naa - 1; iaa++)
1462 daa[iaa] = (aa_datagrid[iaa + 1] - aa_datagrid[iaa - 1]) / 360.;
1463 daa[naa - 1] = (aa_datagrid[naa - 1] - aa_datagrid[naa - 2]) / 360.;
1467 gridpos(pdir_za_gp, za_datagrid, pdir_array);
1468 gridpos(idir_za_gp, za_datagrid, idir_array);
1469 Tensor3 dir_itw(npDir, niDir, 4);
1472 Tensor4 Fou_ssd(nza, nza, stokes_dim, stokes_dim);
1473 Tensor5 Fou_angint(nTin, npDir, niDir, stokes_dim, stokes_dim);
1475 for (
Index find = 0; find < nf; find++) {
1476 for (
Index Tind = 0; Tind < nTin; Tind++) {
1481 for (
Index iza = 0; iza < nza; iza++)
1482 for (
Index sza = 0; sza < nza; sza++) {
1483 for (
Index iaa = 0; iaa < naa; iaa++) {
1486 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1487 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1488 Fou_ssd(sza, iza, ist1, ist2) +=
1500 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1501 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1510 if (this_T_interp_order <
1513 for (
Index Tind = 0; Tind < nTout; Tind++)
1517 for (
Index pdir = 0; pdir < npDir; pdir++)
1518 for (
Index idir = 0; idir < niDir; idir++)
1519 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1520 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1521 for (
Index im = 0; im < nmodes; im++)
1522 reinterp(pha_mat_fou(find,
joker, pdir, idir, ist1, ist2, im),
1523 Fou_angint(
joker, pdir, idir, ist1, ist2),
1567 "The dimension of the stokes vector \n"
1568 "must be 1,2,3 or 4");
1582 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1591 abs_vec_lab.
Kjj()[0] = abs_vec_data(0, 0, 0);
1607 gridpos(gp, za_datagrid, za_sca);
1614 if (stokes_dim == 1) {
1622 out0 <<
"Not all ptype cases are implemented\n";
1663 "The dimension of the stokes vector \n"
1664 "must be 1,2,3 or 4");
1678 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1690 ext_mat_lab.
Kjj() = ext_mat_data(0, 0, 0);
1709 gridpos(gp, za_datagrid, za_sca);
1715 ext_mat_lab.
Kjj() = Kjj;
1717 if (stokes_dim < 2) {
1722 ext_mat_lab.
K12()[0] = K12;
1724 if (stokes_dim < 4) {
1729 ext_mat_lab.
K34()[0] = K34;
1734 out0 <<
"Not all ptype cases are implemented\n";
1774 const Index& za_sca_idx,
1775 const Index& aa_sca_idx,
1776 const Index& za_inc_idx,
1777 const Index& aa_inc_idx,
1781 const Index stokes_dim = pha_mat_lab.
ncols();
1783 Numeric za_sca = za_grid[za_sca_idx];
1784 Numeric aa_sca = aa_grid[aa_sca_idx];
1785 Numeric za_inc = za_grid[za_inc_idx];
1786 Numeric aa_inc = aa_grid[aa_inc_idx];
1789 "The dimension of the stokes vector \n"
1790 "must be 1,2,3 or 4");
1804 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1818 pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1829 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
1830 (aa_sca - aa_inc > 180) *
1838 gridpos(delta_aa_gp, aa_datagrid,
abs(delta_aa));
1839 gridpos(za_inc_gp, za_datagrid, za_inc);
1840 gridpos(za_sca_gp, za_datagrid, za_sca);
1850 if (stokes_dim == 1) {
1871 if (stokes_dim == 2) {
1874 if (delta_aa >= 0) {
1875 pha_mat_lab(0, 2) =
interp(
1881 pha_mat_lab(1, 2) =
interp(
1887 pha_mat_lab(2, 0) =
interp(
1893 pha_mat_lab(2, 1) =
interp(
1900 pha_mat_lab(0, 2) = -
interp(
1906 pha_mat_lab(1, 2) = -
interp(
1912 pha_mat_lab(2, 0) = -
interp(
1918 pha_mat_lab(2, 1) = -
interp(
1925 pha_mat_lab(2, 2) =
interp(
1931 if (stokes_dim == 3) {
1934 if (delta_aa >= 0) {
1935 pha_mat_lab(0, 3) =
interp(
1941 pha_mat_lab(1, 3) =
interp(
1947 pha_mat_lab(3, 0) =
interp(
1953 pha_mat_lab(3, 1) =
interp(
1960 pha_mat_lab(0, 3) = -
interp(
1966 pha_mat_lab(1, 3) = -
interp(
1972 pha_mat_lab(3, 0) = -
interp(
1978 pha_mat_lab(3, 1) = -
interp(
1985 pha_mat_lab(2, 3) =
interp(
1991 pha_mat_lab(3, 2) =
interp(
1997 pha_mat_lab(3, 3) =
interp(
2008 out0 <<
"Not all ptype cases are implemented\n";
2046 const Index& stokes_dim) {
2053 for (
Index is = 0; is < stokes_dim; is++) {
2054 ext_mat(is, is) += abs_vec[0];
2057 for (
Index is = 1; is < stokes_dim; is++) {
2058 ext_mat(0, is) += abs_vec[is];
2059 ext_mat(is, 0) += abs_vec[is];
2091 if ((
abs(aa_sca - aa_inc) < ANG_TOL) ||
2092 (
abs(
abs(aa_sca - aa_inc) - 360) < ANG_TOL)) {
2094 }
else if (
abs(
abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
2096 if (theta_rad >
PI) {
2097 theta_rad = 2 *
PI - theta_rad;
2138 gridpos(thet_gp, za_datagrid, theta);
2144 for (
Index i = 0; i < 6; i++) {
2145 pha_mat_int[i] =
interp(itw, pha_mat_data(
joker, 0, 0, 0, i), thet_gp);
2184 const Index stokes_dim = pha_mat_lab.
ncols();
2187 "NaN value(s) detected in *pha_mat_labCalc* (0,0). Could the "
2188 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2189 "input data are OK and you critically need the ongoing calculations, "
2190 "try to change the observation LOS slightly. If you can reproduce "
2191 "this error, please contact Patrick in order to help tracking down "
2192 "the reason to this problem. If you see this message occasionally "
2193 "when doing MC calculations, it should not be critical. This path "
2194 "sampling will be rejected and replaced with a new one.");
2197 pha_mat_lab(0, 0) =
F11;
2199 if (stokes_dim > 1) {
2205 const Numeric ANGTOL_RAD = 1e-6;
2213 if ((
abs(theta_rad) < ANGTOL_RAD)
2216 (
abs(aa_inc_rad - aa_sca_rad) < ANGTOL_RAD)
2220 pha_mat_lab(0, 1) =
F12;
2221 pha_mat_lab(1, 0) =
F12;
2222 pha_mat_lab(1, 1) =
F22;
2224 if (stokes_dim > 2) {
2225 pha_mat_lab(0, 2) = 0;
2226 pha_mat_lab(1, 2) = 0;
2227 pha_mat_lab(2, 0) = 0;
2228 pha_mat_lab(2, 1) = 0;
2229 pha_mat_lab(2, 2) =
F33;
2231 if (stokes_dim > 3) {
2232 pha_mat_lab(0, 3) = 0;
2233 pha_mat_lab(1, 3) = 0;
2234 pha_mat_lab(2, 3) =
F34;
2235 pha_mat_lab(3, 0) = 0;
2236 pha_mat_lab(3, 1) = 0;
2237 pha_mat_lab(3, 2) = -
F34;
2238 pha_mat_lab(3, 3) =
F44;
2251 if (za_inc_rad < ANGTOL_RAD) {
2252 sigma1 = pi + aa_sca_rad - aa_inc_rad;
2254 }
else if (za_inc_rad > pi - ANGTOL_RAD) {
2255 sigma1 = aa_sca_rad - aa_inc_rad;
2257 }
else if (za_sca_rad < ANGTOL_RAD) {
2259 sigma2 = pi + aa_sca_rad - aa_inc_rad;
2260 }
else if (za_sca_rad > pi - ANGTOL_RAD) {
2262 sigma2 = aa_sca_rad - aa_inc_rad;
2264 s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad)) /
2265 (sin(za_inc_rad) * sin(theta_rad));
2266 s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos(theta_rad)) /
2267 (sin(za_sca_rad) * sin(theta_rad));
2275 if (std::isnan(sigma1) || std::isnan(sigma2)) {
2276 if (
abs(s1 - 1) < ANGTOL_RAD) sigma1 = 0;
2277 if (
abs(s1 + 1) < ANGTOL_RAD) sigma1 = pi;
2278 if (
abs(s2 - 1) < ANGTOL_RAD) sigma2 = 0;
2279 if (
abs(s2 + 1) < ANGTOL_RAD) sigma2 = pi;
2283 const Numeric C1 = cos(2 * sigma1);
2284 const Numeric C2 = cos(2 * sigma2);
2286 const Numeric S1 = sin(2 * sigma1);
2287 const Numeric S2 = sin(2 * sigma2);
2289 pha_mat_lab(0, 1) = C1 *
F12;
2290 pha_mat_lab(1, 0) = C2 *
F12;
2291 pha_mat_lab(1, 1) = C1 * C2 *
F22 - S1 * S2 *
F33;
2296 ARTS_USER_ERROR_IF (std::isnan(pha_mat_lab(0, 1)) || std::isnan(pha_mat_lab(1, 0)) ||
2297 std::isnan(pha_mat_lab(1, 1)),
2298 "NaN value(s) detected in *pha_mat_labCalc* (0/1,1). Could the "
2299 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2300 "input data are OK and you critically need the ongoing calculations, "
2301 "try to change the observation LOS slightly. If you can reproduce "
2302 "this error, please contact Patrick in order to help tracking down "
2303 "the reason to this problem. If you see this message occasionally "
2304 "when doing MC calculations, it should not be critical. This path "
2305 "sampling will be rejected and replaced with a new one.");
2307 if (stokes_dim > 2) {
2317 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
2318 (aa_sca - aa_inc > 180) * 360;
2319 if (delta_aa >= 0) {
2320 pha_mat_lab(0, 2) = S1 *
F12;
2321 pha_mat_lab(1, 2) = S1 * C2 *
F22 + C1 * S2 *
F33;
2322 pha_mat_lab(2, 0) = -S2 *
F12;
2323 pha_mat_lab(2, 1) = -C1 * S2 *
F22 - S1 * C2 *
F33;
2325 pha_mat_lab(0, 2) = -S1 *
F12;
2326 pha_mat_lab(1, 2) = -S1 * C2 *
F22 - C1 * S2 *
F33;
2327 pha_mat_lab(2, 0) = S2 *
F12;
2328 pha_mat_lab(2, 1) = C1 * S2 *
F22 + S1 * C2 *
F33;
2330 pha_mat_lab(2, 2) = -S1 * S2 *
F22 + C1 * C2 *
F33;
2332 if (stokes_dim > 3) {
2333 if (delta_aa >= 0) {
2334 pha_mat_lab(1, 3) = S2 *
F34;
2335 pha_mat_lab(3, 1) = S1 *
F34;
2337 pha_mat_lab(1, 3) = -S2 *
F34;
2338 pha_mat_lab(3, 1) = -S1 *
F34;
2340 pha_mat_lab(0, 3) = 0;
2341 pha_mat_lab(2, 3) = C2 *
F34;
2342 pha_mat_lab(3, 0) = 0;
2343 pha_mat_lab(3, 2) = -C1 *
F34;
2344 pha_mat_lab(3, 3) =
F44;
2352 os <<
"SingleScatteringData: Output operator not implemented";
2357 os <<
"ScatteringMetaData: Output operator not implemented";
2393 abs_vec += propmat_clearsky;
2394 ext_mat += propmat_clearsky;
2410 if (ptype_string ==
"general")
2412 else if (ptype_string ==
"totally_random")
2414 else if (ptype_string ==
"azimuthally_random")
2418 "Unknown ptype: ", ptype_string,
"\n"
2419 "Valid types are: general, totally_random and "
2420 "azimuthally_random.")
2439 if (ptype_string ==
"general")
2441 else if (ptype_string ==
"macroscopically_isotropic")
2443 else if (ptype_string ==
"horizontally_aligned")
2447 "Unknown ptype: ", ptype_string,
"\n"
2448 "Valid types are: general, macroscopically_isotropic and "
2449 "horizontally_aligned.")
2469 ptype_string =
"general";
2472 ptype_string =
"totally_random";
2475 ptype_string =
"azimuthally_random";
2479 "Internal error: Cannot map PType enum value ",
2480 ptype,
" to String.")
2484 return ptype_string;
2499 for (
Index i = 0; i < nza / 2; i++) {
2501 180. - ssd.
za_grid[nza - 1 - i], ssd.
za_grid[i], 2 * DBL_EPSILON),
2502 "Zenith grid of azimuthally_random single scattering data\n"
2503 "is not symmetric with respect to 90degree.")
2506 "Zenith grid of azimuthally_random single scattering data\n"
2507 "does not contain 90 degree grid point.")
2510 ostringstream os_pha_mat;
2511 os_pha_mat <<
"pha_mat ";
2512 ostringstream os_ext_mat;
2513 os_ext_mat <<
"ext_mat ";
2514 ostringstream os_abs_vec;
2515 os_abs_vec <<
"abs_vec ";
2551 for (
Index i = 0; i < nza / 2; i++) {
2563 for (
Index i = 0; i < nza / 2; i++) {
2587 for (
Index i = 0; i < nza / 2; i++)
2588 for (
Index j = 0; j < nza; j++)
2604 const String& particle_ssdmethod_string) {
2606 if (particle_ssdmethod_string ==
"tmatrix")
2610 "Unknown particle SSD method: ",
2611 particle_ssdmethod_string,
"\n"
2612 "Valid methods: tmatrix")
2615 return particle_ssdmethod;
2628 String particle_ssdmethod_string;
2630 switch (particle_ssdmethod) {
2632 particle_ssdmethod_string =
"tmatrix";
2636 "Internal error: Cannot map ParticleSSDMethod enum value ",
2637 particle_ssdmethod,
" to String.")
2641 return particle_ssdmethod_string;
2685 const Index nf = scat_data[0].f_grid.
nelem();
2686 [[maybe_unused]]
const Index nt = T_grid.
nelem();
2688 const Index ncl = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2703 bool all_totrand =
true;
2704 for (
Index ie = 0; ie < nse; ie++) {
2706 all_totrand =
false;
2709 "This method demands that all scat_data are TRO");
2713 const Index cl_start = cloudbox_limits[0];
2718 for (
Index ie = 0; ie < nse; ie++) {
2721 const Numeric tmin = 1.5*scat_data[ie].T_grid[0] -
2722 0.5*scat_data[ie].T_grid[1];
2723 const Numeric tmax = 1.5*scat_data[ie].T_grid[
last] -
2724 0.5*scat_data[ie].T_grid[
last-1];
2727 for(
Index icl=0; icl<ncl; ++icl){
2729 if (
abs(pnd_data(ie,icl)) > 1e-3) {
2730 const Numeric Tthis = T_grid[cl_start + icl];
2732 "Temperature interpolation error for scattering element ", ie,
2733 " of species ", iss,
"\nYour temperature: ", Tthis,
" K\n"
2734 "Allowed range of temperatures: ", tmin,
" - ", tmax,
" K");
2735 T_values[nvals] = Tthis;
2736 cboxlayer[nvals] = icl;
2744 gridpos(gp_t, scat_data[ie].T_grid, T_values[
Range(0,nvals)]);
2750 gridpos(gp_sa, scat_data[ie].za_grid, sa_grid);
2755 for (
Index iv = 0; iv < nf; iv++) {
2757 Vector ext1(nvals), abs1(nvals);
2759 interp(ext1, itw1, scat_data[ie].ext_mat_data(iv,
joker,0,0,0), gp_t);
2760 interp(abs1, itw1, scat_data[ie].abs_vec_data(iv,
joker,0,0,0), gp_t);
2765 for (
Index i = 0; i < nvals; i++) {
2766 const Index ic = cboxlayer[i];
2767 const Index it = cl_start + ic;
2768 ext_data(iv,it) += pnd_data(ie,ic) * ext1[i];
2769 abs_data(iv,it) += pnd_data(ie,ic) * abs1[i];
2770 for (
Index ia = 0; ia < nsa; ia++) {
2771 pfun_data(iv,it,ia) += pnd_data(ie,ic) * pfu1(i,ia);
This file contains the definition of Array.
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
The global header file for ARTS.
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Tensor3.
Index npages() const
Returns the number of pages.
Index nrows() const
Returns the number of rows.
Index ncols() const
Returns the number of columns.
Index ncols() const noexcept
Index nrows() const noexcept
Index nbooks() const noexcept
Index npages() const noexcept
A constant view of a Tensor5.
Index nrows() const noexcept
Index ncols() const noexcept
Index npages() const noexcept
Index nbooks() const noexcept
Index nshelves() const noexcept
Index nbooks() const noexcept
Index nvitrines() const noexcept
Index ncols() const noexcept
Index npages() const noexcept
Index nshelves() const noexcept
Index nrows() const noexcept
Index ncols() const noexcept
Index npages() const noexcept
Index nrows() const noexcept
Index nlibraries() const noexcept
Index nvitrines() const noexcept
Index nshelves() const noexcept
Index nbooks() const noexcept
A constant view of a Vector.
Index nelem() const noexcept
Returns the number of elements.
void resize(Index r, Index c)
Resize function.
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.
Stokes vector is as Propagation matrix but only has 4 possible values.
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
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.
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Header file for logic.cc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Numeric last(ConstVectorView x)
last
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
Declarations having to do with the four output streams.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric two_pi
Two times pi.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto sind(auto x) noexcept
Returns the sine of the deg2rad of the input.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
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)
PType PType2FromString(const String &ptype_string)
Convert ptype name to enum value.
void ext_mat_SSD2Stokes(MatrixView ext_mat_stokes, ConstVectorView ext_mat_ssd, const Index &stokes_dim, const Index &ptype)
Extinction matrix scat_data to stokes format conversion.
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void FouComp_1ScatElem(Tensor7View pha_mat_fou, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Vector &pdir_array, const Vector &idir_array, const Index &f_start, const Index &t_interp_order, const Index &naa_totran)
Preparing phase matrix fourier series components for one scattering element.
constexpr Numeric DEG2RAD
void ext_matFromabs_vec(MatrixView ext_mat, ConstVectorView abs_vec, const Index &stokes_dim)
Derive extinction matrix from absorption vector.
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const PropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
void pha_mat_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
Numeric scat_angle(const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc)
Calculates the scattering angle.
ostream & operator<<(ostream &os, const SingleScatteringData &)
void pha_mat_labCalc(MatrixView pha_mat_lab, ConstVectorView pha_mat_int, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &theta_rad)
Calculate phase matrix in laboratory coordinate system.
void abs_vecTransform(StokesVector &abs_vec_lab, ConstTensor3View abs_vec_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of absorption vector.
ArrayOfLagrangeInterpolation ssd_tinterp_parameters(VectorView t_ok, Index &this_T_interp_order, ConstVectorView T_grid, const Vector &T_array, const Index &t_interp_order)
Determine T-interpol parameters for a specific scattering element.
constexpr Numeric RAD2DEG
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
void pha_mat_ScatSpecBulk(ArrayOfTensor6 &pha_mat, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor6 &pha_mat_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk phase matrices.
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
void ext_matTransform(PropagationMatrix &ext_mat_lab, ConstTensor3View ext_mat_data, ConstVectorView za_datagrid, ConstVectorView aa_datagrid, const PType &ptype, const Numeric &za_sca, const Numeric &aa_sca, const Verbosity &verbosity)
Transformation of extinction matrix.
void ConvertAzimuthallyRandomSingleScatteringData(SingleScatteringData &ssd)
Convert azimuthally-random oriented SingleScatteringData to latest version.
void abs_vec_SSD2Stokes(VectorView abs_vec_stokes, ConstVectorView abs_vec_ssd, const Index &stokes_dim, const Index &ptype)
Absorption vector scat_data to stokes format conversion.
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
ParticleSSDMethod ParticleSSDMethodFromString(const String &particle_ssdmethod_string)
Convert particle ssd method name to enum value.
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
void pha_mat_NScatElems(ArrayOfArrayOfTensor6 &pha_mat, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_index, const Index &t_interp_order)
Phase matrices from all scattering elements.
void interpolate_scat_angle(VectorView pha_mat_int, ConstTensor5View pha_mat_data, ConstVectorView za_datagrid, const Numeric theta)
Interpolate data on the scattering angle.
void opt_prop_1ScatElem(Tensor5View ext_mat, Tensor4View abs_vec, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &dir_array, const Index &f_start, const Index &t_interp_order)
Preparing extinction and absorption from one scattering element.
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.
ParticleSSDMethod
An attribute to classify the method to be used for SingleScatteringData.
@ PARTICLE_SSDMETHOD_TMATRIX
Structure to store a grid position.
This file contains basic functions to handle XML data files.