57#define F11 pha_mat_int[0]
58#define F12 pha_mat_int[1]
59#define F22 pha_mat_int[2]
60#define F33 pha_mat_int[3]
61#define F34 pha_mat_int[4]
62#define F44 pha_mat_int[5]
114 ext_mat = ext_mat_ss[0];
115 abs_vec = abs_vec_ss[0];
117 for (
Index i_ss = 1; i_ss < ext_mat_ss.
nelem(); i_ss++) {
118 ext_mat += ext_mat_ss[i_ss];
119 abs_vec += abs_vec_ss[i_ss];
121 ptype =
max(ptypes_ss);
167 const Index nf = abs_vec_se[0][0].nbooks();
168 const Index nDir = abs_vec_se[0][0].nrows();
169 const Index stokes_dim = abs_vec_se[0][0].ncols();
180 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
185 ext_mat[i_ss].resize(nf, nT, nDir, stokes_dim, stokes_dim);
187 abs_vec[i_ss].resize(nf, nT, nDir, stokes_dim);
190 for (
Index i_se = 0; i_se < ext_mat_se[i_ss].
nelem(); i_se++) {
191 ARTS_ASSERT(nT == ext_mat_se[i_ss][i_se].nbooks());
192 ARTS_ASSERT(nT == abs_vec_se[i_ss][i_se].npages());
194 for (
Index Tind = 0; Tind < nT; Tind++) {
195 if (pnds(i_se_flat, Tind) != 0.) {
196 if (t_ok(i_se_flat, Tind) > 0.) {
198 ext_tmp *= pnds(i_se_flat, Tind);
202 abs_tmp *= pnds(i_se_flat, Tind);
206 "Interpolation error for (flat-array) scattering element #",
208 "at location/temperature point #", Tind,
"\n")
214 ptype[i_ss] =
max(ptypes_se[i_ss]);
256 const Index& stokes_dim,
259 const Index& f_index,
260 const Index& t_interp_order) {
263 nf = scat_data[0][0].ext_mat_data.nshelves();
268 if (scat_data[0][0].ext_mat_data.nshelves() == 1)
287 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
289 ext_mat[i_ss].resize(nse);
290 abs_vec[i_ss].resize(nse);
291 ptypes[i_ss].resize(nse);
293 for (
Index i_se = 0; i_se < nse; i_se++) {
294 ext_mat[i_ss][i_se].resize(nf, nT, nDir, stokes_dim, stokes_dim);
295 abs_vec[i_ss][i_se].resize(nf, nT, nDir, stokes_dim);
300 t_ok(i_se_flat,
joker),
301 scat_data[i_ss][i_se],
332 Index& this_T_interp_order,
336 const Index& t_interp_order) {
340 this_T_interp_order = -1;
344 this_T_interp_order =
min(t_interp_order, nTin - 1);
353 const Numeric extrapolfac = 0.5;
354 const Numeric lowlim = T_grid[0] - extrapolfac * (T_grid[1] - T_grid[0]);
356 T_grid[nTin - 1] + extrapolfac * (T_grid[nTin - 1] - T_grid[nTin - 2]);
358 bool any_T_exceed =
false;
359 for (
Index Tind = 0; Tind < nTout; Tind++) {
360 if (T_array[Tind] < lowlim || T_array[Tind] > uplim) {
370 T_lag.reserve(nTout);
372 bool grid_unchecked =
true;
374 for (
Index iT = 0; iT < nTout; iT++) {
376 T_lag.emplace_back();
378 if (grid_unchecked) {
380 "Temperature interpolation in pha_mat_1ScatElem",
382 T_array[
Range(iT, 1)],
383 this_T_interp_order);
384 grid_unchecked =
false;
386 T_lag.emplace_back(0, T_array[iT], T_grid, this_T_interp_order,
false);
432 const Index& f_start,
433 const Index& t_interp_order) {
478 Index this_T_interp_order;
490 if (this_T_interp_order < 0)
494 Tensor3 ext_mat_tmp(nf, stokes_dim, stokes_dim);
495 Matrix abs_vec_tmp(nf, stokes_dim);
496 for (
Index find = 0; find < nf; find++) {
506 for (
Index Tind = 0; Tind < nTout; Tind++)
507 for (
Index dind = 0; dind < nDir; dind++) {
509 abs_vec(
joker, Tind, dind,
joker) = abs_vec_tmp;
514 Tensor4 ext_mat_tmp(nf, nTout, stokes_dim, stokes_dim);
515 Tensor3 abs_vec_tmp(nf, nTout, stokes_dim);
518 for (
Index find = 0; find < nf; find++) {
519 for (
Index nst = 0; nst < ext_mat_tmp_ssd.
ncols(); nst++) {
524 for (
Index Tind = 0; Tind < nTout; Tind++)
527 ext_mat_tmp_ssd(Tind,
joker),
533 for (
Index nst = 0; nst < abs_vec_tmp_ssd.
ncols(); nst++) {
539 for (
Index Tind = 0; Tind < nTout; Tind++)
542 abs_vec_tmp_ssd(Tind,
joker),
546 abs_vec_tmp(find, Tind,
joker) = 0.;
549 for (
Index dind = 0; dind < nDir; dind++) {
573 if (this_T_interp_order < 0)
575 Matrix ext_mat_tmp_ssd(nDir, next);
576 Matrix abs_vec_tmp_ssd(nDir, nabs);
577 Tensor4 ext_mat_tmp(nf, nDir, stokes_dim, stokes_dim);
578 Tensor3 abs_vec_tmp(nf, nDir, stokes_dim);
579 for (
Index find = 0; find < nf; find++) {
580 for (
Index nst = 0; nst < next; nst++)
585 for (
Index Dind = 0; Dind < nDir; Dind++)
587 ext_mat_tmp_ssd(Dind,
joker),
591 for (
Index nst = 0; nst < nabs; nst++)
596 for (
Index Dind = 0; Dind < nDir; Dind++)
598 abs_vec_tmp_ssd(Dind,
joker),
603 for (
Index Tind = 0; Tind < nTout; Tind++) {
610 Tensor3 ext_mat_tmp_ssd(nTin, nDir, next);
611 Tensor3 abs_vec_tmp_ssd(nTin, nDir, nabs);
612 Matrix ext_mat_tmp(nTout, next);
613 Matrix abs_vec_tmp(nTout, nabs);
615 for (
Index find = 0; find < nf; find++) {
616 for (
Index Tind = 0; Tind < nTin; Tind++) {
617 for (
Index nst = 0; nst < next; nst++)
622 for (
Index nst = 0; nst < nabs; nst++)
629 for (
Index Dind = 0; Dind < nDir; Dind++) {
630 for (
Index nst = 0; nst < next; nst++) {
632 ext_mat_tmp_ssd(
joker, Dind, nst),
636 for (
Index Tind = 0; Tind < nTout; Tind++)
638 ext_mat_tmp(Tind,
joker),
642 for (
Index nst = 0; nst < nabs; nst++) {
644 abs_vec_tmp_ssd(
joker, Dind, nst),
648 for (
Index Tind = 0; Tind < nTout; Tind++)
650 abs_vec_tmp(Tind,
joker),
677 const Index& stokes_dim,
678 const Index& ptype) {
685 for (
Index ist = 0; ist < stokes_dim; ist++) {
686 ext_mat_stokes(ist, ist) = ext_mat_ssd[0];
690 switch (stokes_dim) {
692 ext_mat_stokes(2, 3) = ext_mat_ssd[2];
693 ext_mat_stokes(3, 2) = -ext_mat_ssd[2];
700 ext_mat_stokes(0, 1) = ext_mat_ssd[1];
701 ext_mat_stokes(1, 0) = ext_mat_ssd[1];
725 const Index& stokes_dim,
726 const Index& ptype) {
733 abs_vec_stokes[0] = abs_vec_ssd[0];
736 abs_vec_stokes[1] = abs_vec_ssd[1];
764 pha_mat = pha_mat_ss[0];
766 for (
Index i_ss = 1; i_ss < pha_mat_ss.
nelem(); i_ss++)
767 pha_mat += pha_mat_ss[i_ss];
769 ptype =
max(ptypes_ss);
810 const Index nf = pha_mat_se[0][0].nvitrines();
811 const Index npDir = pha_mat_se[0][0].nbooks();
812 const Index niDir = pha_mat_se[0][0].npages();
813 const Index stokes_dim = pha_mat_se[0][0].ncols();
822 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
825 pha_mat[i_ss].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
828 for (
Index i_se = 0; i_se < pha_mat_se[i_ss].
nelem(); i_se++) {
829 ARTS_ASSERT(nT == pha_mat_se[i_ss][i_se].nshelves());
831 for (
Index Tind = 0; Tind < nT; Tind++) {
832 if (pnds(i_se_flat, Tind) != 0.) {
833 if (t_ok(i_se_flat, Tind) > 0.) {
836 pha_tmp *= pnds(i_se_flat, Tind);
840 "Interpolation error for (flat-array) scattering element #",
842 "at location/temperature point #", Tind,
"\n")
848 ptype[i_ss] =
max(ptypes_se[i_ss]);
889 const Index& stokes_dim,
893 const Index& f_index,
894 const Index& t_interp_order) {
897 nf = scat_data[0][0].pha_mat_data.nlibraries();
902 if (scat_data[0][0].pha_mat_data.nlibraries() == 1)
921 for (
Index i_ss = 0; i_ss < nss; i_ss++) {
923 pha_mat[i_ss].resize(nse);
924 ptypes[i_ss].resize(nse);
926 for (
Index i_se = 0; i_se < nse; i_se++) {
927 pha_mat[i_ss][i_se].resize(nf, nT, npDir, niDir, stokes_dim, stokes_dim);
931 t_ok(i_se_flat,
joker),
932 scat_data[i_ss][i_se],
977 const Index& f_start,
978 const Index& t_interp_order) {
1005 Index this_T_interp_order;
1007 this_T_interp_order,
1020 if (stokes_dim == 1)
1022 else if (stokes_dim < 4)
1026 if (this_T_interp_order < 0)
1030 for (
Index pdir = 0; pdir < npDir; pdir++)
1031 for (
Index idir = 0; idir < niDir; idir++) {
1034 pdir_array(pdir, 1),
1035 idir_array(idir, 0),
1036 idir_array(idir, 1));
1044 Vector pha_mat_int(npha, 0.);
1045 Matrix pha_mat_tmp(stokes_dim, stokes_dim);
1046 for (
Index find = 0; find < nf; find++) {
1048 for (
Index nst = 0; nst < npha; nst++)
1049 pha_mat_int[nst] =
interp(
1057 pdir_array(pdir, 0),
1058 pdir_array(pdir, 1),
1059 idir_array(idir, 0),
1060 idir_array(idir, 1),
1063 for (
Index Tind = 0; Tind < nTout; Tind++)
1064 pha_mat(find, Tind, pdir, idir,
joker,
joker) = pha_mat_tmp;
1069 for (
Index pdir = 0; pdir < npDir; pdir++)
1070 for (
Index idir = 0; idir < niDir; idir++) {
1073 pdir_array(pdir, 1),
1074 idir_array(idir, 0),
1075 idir_array(idir, 1));
1083 Matrix pha_mat_int(nTin, npha, 0.);
1084 Matrix pha_mat_tmp(nTout, npha, 0.);
1085 for (
Index find = 0; find < nf; find++) {
1086 for (
Index Tind = 0; Tind < nTin; Tind++)
1088 for (
Index nst = 0; nst < npha; nst++) {
1089 pha_mat_int(Tind, nst) =
interp(
1095 for (
Index nst = 0; nst < npha; nst++) {
1097 pha_mat_int(
joker, nst),
1106 for (
Index Tind = 0; Tind < nTout; Tind++) {
1109 pha_mat_tmp(Tind,
joker),
1110 pdir_array(pdir, 0),
1111 pdir_array(pdir, 1),
1112 idir_array(idir, 0),
1113 idir_array(idir, 1),
1123 Index nDir = npDir * niDir;
1125 Matrix delta_aa(npDir, niDir);
1133 for (
Index pdir = 0; pdir < npDir; pdir++) {
1134 for (
Index idir = 0; idir < niDir; idir++) {
1135 delta_aa(pdir, idir) =
1136 pdir_array(pdir, 1) - idir_array(idir, 1) +
1137 (pdir_array(pdir, 1) - idir_array(idir, 1) < -180) * 360 -
1138 (pdir_array(pdir, 1) - idir_array(idir, 1) > 180) * 360;
1139 adelta_aa[j] =
abs(delta_aa(pdir, idir));
1140 pza_gp[j] = pza_gp_tmp[pdir];
1141 iza_gp[j] = iza_gp_tmp[idir];
1151 if (this_T_interp_order < 0)
1153 Tensor3 pha_mat_int(nDir, stokes_dim, stokes_dim, 0.);
1154 Tensor4 pha_mat_tmp(npDir, niDir, stokes_dim, stokes_dim);
1156 for (
Index find = 0; find < nf; find++) {
1160 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1161 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1163 pha_mat_int(
joker, ist1, ist2),
1174 for (
Index pdir = 0; pdir < npDir; pdir++)
1175 for (
Index idir = 0; idir < niDir; idir++) {
1181 if (stokes_dim > 2) {
1182 for (
Index pdir = 0; pdir < npDir; pdir++)
1183 for (
Index idir = 0; idir < niDir; idir++)
1184 if (delta_aa(pdir, idir) < 0.) {
1185 pha_mat_tmp(pdir, idir, 0, 2) *= -1;
1186 pha_mat_tmp(pdir, idir, 1, 2) *= -1;
1187 pha_mat_tmp(pdir, idir, 2, 0) *= -1;
1188 pha_mat_tmp(pdir, idir, 2, 1) *= -1;
1192 if (stokes_dim > 3) {
1193 for (
Index pdir = 0; pdir < npDir; pdir++)
1194 for (
Index idir = 0; idir < niDir; idir++)
1195 if (delta_aa(pdir, idir) < 0.) {
1196 pha_mat_tmp(pdir, idir, 0, 3) *= -1;
1197 pha_mat_tmp(pdir, idir, 1, 3) *= -1;
1198 pha_mat_tmp(pdir, idir, 3, 0) *= -1;
1199 pha_mat_tmp(pdir, idir, 3, 1) *= -1;
1203 for (
Index Tind = 0; Tind < nTout; Tind++)
1211 Tensor4 pha_mat_int(nTin, nDir, stokes_dim, stokes_dim, 0.);
1213 for (
Index find = 0; find < nf; find++) {
1216 for (
Index Tind = 0; Tind < nTin; Tind++) {
1217 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1218 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1237 for (
Index pdir = 0; pdir < npDir; pdir++)
1238 for (
Index idir = 0; idir < niDir; idir++) {
1239 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1240 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1242 pha_mat_int(
joker, i, ist1, ist2),
1248 if (stokes_dim > 2) {
1249 for (
Index pdir = 0; pdir < npDir; pdir++)
1250 for (
Index idir = 0; idir < niDir; idir++)
1251 if (delta_aa(pdir, idir) < 0.) {
1252 pha_mat(find,
joker, pdir, idir, 0, 2) *= -1;
1253 pha_mat(find,
joker, pdir, idir, 1, 2) *= -1;
1254 pha_mat(find,
joker, pdir, idir, 2, 0) *= -1;
1255 pha_mat(find,
joker, pdir, idir, 2, 1) *= -1;
1259 if (stokes_dim > 3) {
1260 for (
Index pdir = 0; pdir < npDir; pdir++)
1261 for (
Index idir = 0; idir < niDir; idir++)
1262 if (delta_aa(pdir, idir) < 0.) {
1263 pha_mat(find,
joker, pdir, idir, 0, 3) *= -1;
1264 pha_mat(find,
joker, pdir, idir, 1, 3) *= -1;
1265 pha_mat(find,
joker, pdir, idir, 3, 0) *= -1;
1266 pha_mat(find,
joker, pdir, idir, 3, 1) *= -1;
1309 const Vector& pdir_array,
1310 const Vector& idir_array,
1311 const Index& f_start,
1312 const Index& t_interp_order,
1313 const Index& naa_totran) {
1330 const Index stokes_dim = pha_mat_fou.
nrows();
1347 Index this_T_interp_order;
1349 this_T_interp_order,
1361 "Azimuth grid size for scatt matrix extraction"
1362 " (*naa_totran*) must be >3.\n"
1363 "Yours is ", naa_totran,
".\n")
1367 1. / float(naa_totran - 1);
1368 Vector theta(naa_totran);
1370 Matrix theta_itw(naa_totran, 2);
1373 if (stokes_dim == 1)
1375 else if (stokes_dim < 4)
1380 Matrix pha_mat(stokes_dim, stokes_dim);
1381 Matrix pha_mat_angint(naa_totran, npha);
1382 Tensor3 Fou_int(nTin, stokes_dim, stokes_dim);
1384 for (
Index idir = 0; idir < niDir; idir++)
1385 for (
Index pdir = 0; pdir < npDir; pdir++) {
1386 for (
Index iaa = 0; iaa < naa_totran; iaa++)
1389 scat_angle(pdir_array[pdir], aa_grid[iaa], idir_array[idir], 0.);
1395 for (
Index find = 0; find < nf; find++) {
1397 for (
Index Tind = 0; Tind < nTin; Tind++) {
1399 for (
Index nst = 0; nst < npha; nst++)
1401 pha_mat_angint(
joker, nst),
1405 for (
Index iaa = 0; iaa < naa_totran; iaa++) {
1408 pha_mat_angint(iaa,
joker),
1417 if (iaa == 0 || iaa == naa_totran - 1)
1418 pha_mat *= (daa_totran / 2.);
1420 pha_mat *= daa_totran;
1425 if (this_T_interp_order <
1428 for (
Index Tind = 0; Tind < nTout; Tind++)
1429 pha_mat_fou(find, Tind, pdir, idir,
joker,
joker, 0) =
1432 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1433 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1434 for (
Index im = 0; im < nmodes; im++)
1435 reinterp(pha_mat_fou(find,
joker, pdir, idir, ist1, ist2, im),
1436 Fou_int(
joker, ist1, ist2),
1459 daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1460 for (
Index iaa = 1; iaa < naa - 1; iaa++)
1461 daa[iaa] = (aa_datagrid[iaa + 1] - aa_datagrid[iaa - 1]) / 360.;
1462 daa[naa - 1] = (aa_datagrid[naa - 1] - aa_datagrid[naa - 2]) / 360.;
1466 gridpos(pdir_za_gp, za_datagrid, pdir_array);
1467 gridpos(idir_za_gp, za_datagrid, idir_array);
1468 Tensor3 dir_itw(npDir, niDir, 4);
1471 Tensor4 Fou_ssd(nza, nza, stokes_dim, stokes_dim);
1472 Tensor5 Fou_angint(nTin, npDir, niDir, stokes_dim, stokes_dim);
1474 for (
Index find = 0; find < nf; find++) {
1475 for (
Index Tind = 0; Tind < nTin; Tind++) {
1480 for (
Index iza = 0; iza < nza; iza++)
1481 for (
Index sza = 0; sza < nza; sza++) {
1482 for (
Index iaa = 0; iaa < naa; iaa++) {
1485 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1486 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1487 Fou_ssd(sza, iza, ist1, ist2) +=
1499 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1500 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1509 if (this_T_interp_order <
1512 for (
Index Tind = 0; Tind < nTout; Tind++)
1516 for (
Index pdir = 0; pdir < npDir; pdir++)
1517 for (
Index idir = 0; idir < niDir; idir++)
1518 for (
Index ist1 = 0; ist1 < stokes_dim; ist1++)
1519 for (
Index ist2 = 0; ist2 < stokes_dim; ist2++)
1520 for (
Index im = 0; im < nmodes; im++)
1521 reinterp(pha_mat_fou(find,
joker, pdir, idir, ist1, ist2, im),
1522 Fou_angint(
joker, pdir, idir, ist1, ist2),
1566 "The dimension of the stokes vector \n"
1567 "must be 1,2,3 or 4");
1581 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1590 abs_vec_lab.
Kjj()[0] = abs_vec_data(0, 0, 0);
1606 gridpos(gp, za_datagrid, za_sca);
1613 if (stokes_dim == 1) {
1621 out0 <<
"Not all ptype cases are implemented\n";
1662 "The dimension of the stokes vector \n"
1663 "must be 1,2,3 or 4");
1677 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1689 ext_mat_lab.
Kjj() = ext_mat_data(0, 0, 0);
1708 gridpos(gp, za_datagrid, za_sca);
1714 ext_mat_lab.
Kjj() = Kjj;
1716 if (stokes_dim < 2) {
1721 ext_mat_lab.
K12()[0] = K12;
1723 if (stokes_dim < 4) {
1728 ext_mat_lab.
K34()[0] = K34;
1733 out0 <<
"Not all ptype cases are implemented\n";
1773 const Index& za_sca_idx,
1774 const Index& aa_sca_idx,
1775 const Index& za_inc_idx,
1776 const Index& aa_inc_idx,
1780 const Index stokes_dim = pha_mat_lab.
ncols();
1782 Numeric za_sca = za_grid[za_sca_idx];
1783 Numeric aa_sca = aa_grid[aa_sca_idx];
1784 Numeric za_inc = za_grid[za_inc_idx];
1785 Numeric aa_inc = aa_grid[aa_inc_idx];
1788 "The dimension of the stokes vector \n"
1789 "must be 1,2,3 or 4");
1803 out0 <<
"Case PTYPE_GENERAL not yet implemented. \n";
1817 pha_mat_lab, pha_mat_int, za_sca, aa_sca, za_inc, aa_inc, theta_rad);
1828 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
1829 (aa_sca - aa_inc > 180) *
1837 gridpos(delta_aa_gp, aa_datagrid,
abs(delta_aa));
1838 gridpos(za_inc_gp, za_datagrid, za_inc);
1839 gridpos(za_sca_gp, za_datagrid, za_sca);
1849 if (stokes_dim == 1) {
1870 if (stokes_dim == 2) {
1873 if (delta_aa >= 0) {
1874 pha_mat_lab(0, 2) =
interp(
1880 pha_mat_lab(1, 2) =
interp(
1886 pha_mat_lab(2, 0) =
interp(
1892 pha_mat_lab(2, 1) =
interp(
1899 pha_mat_lab(0, 2) = -
interp(
1905 pha_mat_lab(1, 2) = -
interp(
1911 pha_mat_lab(2, 0) = -
interp(
1917 pha_mat_lab(2, 1) = -
interp(
1924 pha_mat_lab(2, 2) =
interp(
1930 if (stokes_dim == 3) {
1933 if (delta_aa >= 0) {
1934 pha_mat_lab(0, 3) =
interp(
1940 pha_mat_lab(1, 3) =
interp(
1946 pha_mat_lab(3, 0) =
interp(
1952 pha_mat_lab(3, 1) =
interp(
1959 pha_mat_lab(0, 3) = -
interp(
1965 pha_mat_lab(1, 3) = -
interp(
1971 pha_mat_lab(3, 0) = -
interp(
1977 pha_mat_lab(3, 1) = -
interp(
1984 pha_mat_lab(2, 3) =
interp(
1990 pha_mat_lab(3, 2) =
interp(
1996 pha_mat_lab(3, 3) =
interp(
2007 out0 <<
"Not all ptype cases are implemented\n";
2045 const Index& stokes_dim) {
2052 for (
Index is = 0; is < stokes_dim; is++) {
2053 ext_mat(is, is) += abs_vec[0];
2056 for (
Index is = 1; is < stokes_dim; is++) {
2057 ext_mat(0, is) += abs_vec[is];
2058 ext_mat(is, 0) += abs_vec[is];
2090 if ((
abs(aa_sca - aa_inc) < ANG_TOL) ||
2091 (
abs(
abs(aa_sca - aa_inc) - 360) < ANG_TOL)) {
2093 }
else if (
abs(
abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
2095 if (theta_rad >
PI) {
2096 theta_rad = 2 *
PI - theta_rad;
2137 gridpos(thet_gp, za_datagrid, theta);
2143 for (
Index i = 0; i < 6; i++) {
2144 pha_mat_int[i] =
interp(itw, pha_mat_data(
joker, 0, 0, 0, i), thet_gp);
2183 const Index stokes_dim = pha_mat_lab.
ncols();
2186 "NaN value(s) detected in *pha_mat_labCalc* (0,0). Could the "
2187 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2188 "input data are OK and you critically need the ongoing calculations, "
2189 "try to change the observation LOS slightly. If you can reproduce "
2190 "this error, please contact Patrick in order to help tracking down "
2191 "the reason to this problem. If you see this message occasionally "
2192 "when doing MC calculations, it should not be critical. This path "
2193 "sampling will be rejected and replaced with a new one.");
2196 pha_mat_lab(0, 0) =
F11;
2198 if (stokes_dim > 1) {
2204 const Numeric ANGTOL_RAD = 1e-6;
2212 if ((
abs(theta_rad) < ANGTOL_RAD)
2213 || (
abs(theta_rad -
PI) < ANGTOL_RAD)
2215 (
abs(aa_inc_rad - aa_sca_rad) < ANGTOL_RAD)
2216 || (
abs(
abs(aa_inc_rad - aa_sca_rad) - 360.) < ANGTOL_RAD)
2217 || (
abs(
abs(aa_inc_rad - aa_sca_rad) - 180.) < ANGTOL_RAD)
2219 pha_mat_lab(0, 1) =
F12;
2220 pha_mat_lab(1, 0) =
F12;
2221 pha_mat_lab(1, 1) =
F22;
2223 if (stokes_dim > 2) {
2224 pha_mat_lab(0, 2) = 0;
2225 pha_mat_lab(1, 2) = 0;
2226 pha_mat_lab(2, 0) = 0;
2227 pha_mat_lab(2, 1) = 0;
2228 pha_mat_lab(2, 2) =
F33;
2230 if (stokes_dim > 3) {
2231 pha_mat_lab(0, 3) = 0;
2232 pha_mat_lab(1, 3) = 0;
2233 pha_mat_lab(2, 3) =
F34;
2234 pha_mat_lab(3, 0) = 0;
2235 pha_mat_lab(3, 1) = 0;
2236 pha_mat_lab(3, 2) = -
F34;
2237 pha_mat_lab(3, 3) =
F44;
2250 if (za_inc_rad < ANGTOL_RAD) {
2251 sigma1 =
PI + aa_sca_rad - aa_inc_rad;
2253 }
else if (za_inc_rad >
PI - ANGTOL_RAD) {
2254 sigma1 = aa_sca_rad - aa_inc_rad;
2256 }
else if (za_sca_rad < ANGTOL_RAD) {
2258 sigma2 =
PI + aa_sca_rad - aa_inc_rad;
2259 }
else if (za_sca_rad >
PI - ANGTOL_RAD) {
2261 sigma2 = aa_sca_rad - aa_inc_rad;
2263 s1 = (cos(za_sca_rad) - cos(za_inc_rad) * cos(theta_rad)) /
2264 (sin(za_inc_rad) * sin(theta_rad));
2265 s2 = (cos(za_inc_rad) - cos(za_sca_rad) * cos(theta_rad)) /
2266 (sin(za_sca_rad) * sin(theta_rad));
2275 if (
abs(s1 - 1) < ANGTOL_RAD) sigma1 = 0;
2276 if (
abs(s1 + 1) < ANGTOL_RAD) sigma1 =
PI;
2277 if (
abs(s2 - 1) < ANGTOL_RAD) sigma2 = 0;
2278 if (
abs(s2 + 1) < ANGTOL_RAD) sigma2 =
PI;
2282 const Numeric C1 = cos(2 * sigma1);
2283 const Numeric C2 = cos(2 * sigma2);
2285 const Numeric S1 = sin(2 * sigma1);
2286 const Numeric S2 = sin(2 * sigma2);
2288 pha_mat_lab(0, 1) = C1 *
F12;
2289 pha_mat_lab(1, 0) = C2 *
F12;
2290 pha_mat_lab(1, 1) = C1 * C2 *
F22 - S1 * S2 *
F33;
2297 "NaN value(s) detected in *pha_mat_labCalc* (0/1,1). Could the "
2298 "input data contain NaNs? Please check with *scat_dataCheck*. If "
2299 "input data are OK and you critically need the ongoing calculations, "
2300 "try to change the observation LOS slightly. If you can reproduce "
2301 "this error, please contact Patrick in order to help tracking down "
2302 "the reason to this problem. If you see this message occasionally "
2303 "when doing MC calculations, it should not be critical. This path "
2304 "sampling will be rejected and replaced with a new one.");
2306 if (stokes_dim > 2) {
2316 Numeric delta_aa = aa_sca - aa_inc + (aa_sca - aa_inc < -180) * 360 -
2317 (aa_sca - aa_inc > 180) * 360;
2318 if (delta_aa >= 0) {
2319 pha_mat_lab(0, 2) = S1 *
F12;
2320 pha_mat_lab(1, 2) = S1 * C2 *
F22 + C1 * S2 *
F33;
2321 pha_mat_lab(2, 0) = -S2 *
F12;
2322 pha_mat_lab(2, 1) = -C1 * S2 *
F22 - S1 * C2 *
F33;
2324 pha_mat_lab(0, 2) = -S1 *
F12;
2325 pha_mat_lab(1, 2) = -S1 * C2 *
F22 - C1 * S2 *
F33;
2326 pha_mat_lab(2, 0) = S2 *
F12;
2327 pha_mat_lab(2, 1) = C1 * S2 *
F22 + S1 * C2 *
F33;
2329 pha_mat_lab(2, 2) = -S1 * S2 *
F22 + C1 * C2 *
F33;
2331 if (stokes_dim > 3) {
2332 if (delta_aa >= 0) {
2333 pha_mat_lab(1, 3) = S2 *
F34;
2334 pha_mat_lab(3, 1) = S1 *
F34;
2336 pha_mat_lab(1, 3) = -S2 *
F34;
2337 pha_mat_lab(3, 1) = -S1 *
F34;
2339 pha_mat_lab(0, 3) = 0;
2340 pha_mat_lab(2, 3) = C2 *
F34;
2341 pha_mat_lab(3, 0) = 0;
2342 pha_mat_lab(3, 2) = -C1 *
F34;
2343 pha_mat_lab(3, 3) =
F44;
2351 os <<
"SingleScatteringData: Output operator not implemented";
2356 os <<
"ArrayOfSingleScatteringData: Output operator not implemented";
2361 os <<
"ScatteringMetaData: Output operator not implemented";
2366 os <<
"ArrayOfScatteringMetaData: Output operator not implemented";
2402 abs_vec += propmat_clearsky;
2403 ext_mat += propmat_clearsky;
2419 if (ptype_string ==
"general")
2421 else if (ptype_string ==
"totally_random")
2423 else if (ptype_string ==
"azimuthally_random")
2427 "Unknown ptype: ", ptype_string,
"\n"
2428 "Valid types are: general, totally_random and "
2429 "azimuthally_random.")
2448 if (ptype_string ==
"general")
2450 else if (ptype_string ==
"macroscopically_isotropic")
2452 else if (ptype_string ==
"horizontally_aligned")
2456 "Unknown ptype: ", ptype_string,
"\n"
2457 "Valid types are: general, macroscopically_isotropic and "
2458 "horizontally_aligned.")
2478 ptype_string =
"general";
2481 ptype_string =
"totally_random";
2484 ptype_string =
"azimuthally_random";
2488 "Internal error: Cannot map PType enum value ",
2489 ptype,
" to String.")
2493 return ptype_string;
2508 for (
Index i = 0; i < nza / 2; i++) {
2510 180. - ssd.
za_grid[nza - 1 - i], ssd.
za_grid[i], 2 * DBL_EPSILON),
2511 "Zenith grid of azimuthally_random single scattering data\n"
2512 "is not symmetric with respect to 90degree.")
2515 "Zenith grid of azimuthally_random single scattering data\n"
2516 "does not contain 90 degree grid point.")
2519 ostringstream os_pha_mat;
2520 os_pha_mat <<
"pha_mat ";
2521 ostringstream os_ext_mat;
2522 os_ext_mat <<
"ext_mat ";
2523 ostringstream os_abs_vec;
2524 os_abs_vec <<
"abs_vec ";
2560 for (
Index i = 0; i < nza / 2; i++) {
2572 for (
Index i = 0; i < nza / 2; i++) {
2596 for (
Index i = 0; i < nza / 2; i++)
2597 for (
Index j = 0; j < nza; j++)
2613 const String& particle_ssdmethod_string) {
2615 if (particle_ssdmethod_string ==
"tmatrix")
2619 "Unknown particle SSD method: ",
2620 particle_ssdmethod_string,
"\n"
2621 "Valid methods: tmatrix")
2624 return particle_ssdmethod;
2637 String particle_ssdmethod_string;
2639 switch (particle_ssdmethod) {
2641 particle_ssdmethod_string =
"tmatrix";
2645 "Internal error: Cannot map ParticleSSDMethod enum value ",
2646 particle_ssdmethod,
" to String.")
2650 return particle_ssdmethod_string;
2694 const Index nf = scat_data[0].f_grid.
nelem();
2695 [[maybe_unused]]
const Index nt = T_grid.
nelem();
2697 const Index ncl = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2712 bool all_totrand =
true;
2713 for (
Index ie = 0; ie < nse; ie++) {
2715 all_totrand =
false;
2718 "This method demands that all scat_data are TRO");
2722 const Index cl_start = cloudbox_limits[0];
2727 for (
Index ie = 0; ie < nse; ie++) {
2730 const Numeric tmin = 1.5*scat_data[ie].T_grid[0] -
2731 0.5*scat_data[ie].T_grid[1];
2732 const Numeric tmax = 1.5*scat_data[ie].T_grid[
last] -
2733 0.5*scat_data[ie].T_grid[
last-1];
2736 for(
Index icl=0; icl<ncl; ++icl){
2738 if (
abs(pnd_data(ie,icl)) > 1e-3) {
2739 const Numeric Tthis = T_grid[cl_start + icl];
2741 "Temperature interpolation error for scattering element ", ie,
2742 " of species ", iss,
"\nYour temperature: ", Tthis,
" K\n"
2743 "Allowed range of temperatures: ", tmin,
" - ", tmax,
" K");
2744 T_values[nvals] = Tthis;
2745 cboxlayer[nvals] = icl;
2753 gridpos(gp_t, scat_data[ie].T_grid, T_values[
Range(0,nvals)]);
2759 gridpos(gp_sa, scat_data[ie].za_grid, sa_grid);
2764 for (
Index iv = 0; iv < nf; iv++) {
2766 Vector ext1(nvals), abs1(nvals);
2768 interp(ext1, itw1, scat_data[ie].ext_mat_data(iv,
joker,0,0,0), gp_t);
2769 interp(abs1, itw1, scat_data[ie].abs_vec_data(iv,
joker,0,0,0), gp_t);
2774 for (
Index i = 0; i < nvals; i++) {
2775 const Index ic = cboxlayer[i];
2776 const Index it = cl_start + ic;
2777 ext_data(iv,it) += pnd_data(ie,ic) * ext1[i];
2778 abs_data(iv,it) += pnd_data(ie,ic) * abs1[i];
2779 for (
Index ia = 0; ia < nsa; ia++) {
2780 pfun_data(iv,it,ia) += pnd_data(ie,ic) * pfu1(i,ia);
This file contains the definition of Array.
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
The global header file for ARTS.
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
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.
Constants of physical expressions as constexpr.
#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
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.
Index nelem(const Lines &l)
Number of lines.
auto sind(T x) noexcept -> decltype(std::sin(deg2rad(x)))
Returns the sine of the deg2rad of the input.
auto cosd(T x) noexcept -> decltype(std::cos(deg2rad(x)))
Returns the cosine of the deg2rad of the input.
constexpr auto deg2rad(T x) noexcept -> decltype(x *one_degree_in_radians)
Converts degrees to radians.
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
constexpr bool isnan(double d) noexcept
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.
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.
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.