72 const Numeric& rtp_planck_value,
73 const bool& trans_is_precalc) {
100 if (!trans_is_precalc) {
102 trans_mat, lstep, ext_mat_av, 0);
106 if (stokes_dim == 1) {
107 stokes_vec[0] = stokes_vec[0] * trans_mat(0, 0) +
108 (abs_vec_av.
Kjj()[0] * rtp_planck_value + sca_vec_av[0]) /
109 ext_mat_av.
Kjj()[0] * (1 - trans_mat(0, 0));
135 Matrix invK(stokes_dim, stokes_dim);
139 source *= rtp_planck_value;
141 for (
Index i = 0; i < stokes_dim; i++)
142 source[i] += sca_vec_av[i];
146 mult(x, invK, source);
151 Matrix ImT(stokes_dim, stokes_dim);
158 mult(term1, trans_mat, stokes_vec);
160 for (
Index i = 0; i < stokes_dim; i++)
161 stokes_vec[i] = term1[i] + term2[i];
170 const Agenda& spt_calc_agenda,
171 const Index& za_index,
172 const Index& aa_index,
182 out3 <<
"Calculate scattering properties in cloudbox \n";
184 const Index atmosphere_dim = cloudbox_limits.
nelem() / 2;
186 const Index stokes_dim = ext_mat_field.
ncols();
188 ARTS_ASSERT(atmosphere_dim == 1 || atmosphere_dim == 3);
190 ext_mat_field.
ncols() == abs_vec_field.
ncols());
192 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
195 Index Nlat_cloud = 1;
196 Index Nlon_cloud = 1;
198 if (atmosphere_dim == 3) {
199 Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
200 Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
207 for (
auto& av : abs_vec_spt_local) {
212 for (
auto& pm : ext_mat_spt_local) {
232 for (
Index scat_p_index_local = 0; scat_p_index_local < Np_cloud;
233 scat_p_index_local++) {
234 for (
Index scat_lat_index_local = 0; scat_lat_index_local < Nlat_cloud;
235 scat_lat_index_local++) {
236 for (
Index scat_lon_index_local = 0; scat_lon_index_local < Nlon_cloud;
237 scat_lon_index_local++) {
238 if (atmosphere_dim == 1)
239 rtp_temperature_local =
240 t_field(scat_p_index_local + cloudbox_limits[0], 0, 0);
242 rtp_temperature_local =
243 t_field(scat_p_index_local + cloudbox_limits[0],
244 scat_lat_index_local + cloudbox_limits[2],
245 scat_lon_index_local + cloudbox_limits[4]);
253 scat_lat_index_local,
254 scat_lon_index_local,
255 rtp_temperature_local,
283 scat_lat_index_local,
284 scat_lon_index_local,
288 abs_vec_field(scat_p_index_local,
289 scat_lat_index_local,
290 scat_lon_index_local,
294 scat_lat_index_local,
295 scat_lon_index_local,
307 const Index& p_index,
308 const Index& za_index,
313 const Agenda& propmat_clearsky_agenda,
316 const Agenda& ppath_step_agenda,
318 const Numeric& ppath_lraytrace,
325 const Index& f_index,
329 const Agenda& surface_rtprop_agenda,
331 const Index& scat_za_interp,
343 ppath_step.
pos(0, 0) = z_field(p_index, 0, 0);
344 ppath_step.
r[0] = refellipsoid[0] + z_field(p_index, 0, 0);
347 ppath_step.
los(0, 0) = za_grid[za_index];
350 ppath_step.
gp_p[0].idx = p_index;
351 ppath_step.
gp_p[0].fd[0] = 0;
352 ppath_step.
gp_p[0].fd[1] = 1;
359 Vector(1, f_grid[f_index]),
367 if ((cloudbox_limits[0] <= ppath_step.
gp_p[1].idx &&
368 cloudbox_limits[1] > ppath_step.
gp_p[1].idx) ||
369 (cloudbox_limits[1] == ppath_step.
gp_p[1].idx &&
370 abs(ppath_step.
gp_p[1].fd[0]) < 1e-6)) {
372 const Index stokes_dim = cloudbox_field_mono.
ncols();
383 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.
np, 0.);
384 Matrix abs_vec_int(stokes_dim, ppath_step.
np, 0.);
385 Matrix sca_vec_int(stokes_dim, ppath_step.
np, 0.);
386 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.
np, 0.);
388 Matrix vmr_list_int(N_species, ppath_step.
np, 0.);
394 cloudbox_field_mono_int,
422 propmat_clearsky_agenda,
429 cloudbox_field_mono_int,
446 surface_rtprop_agenda,
463 const Index& p_index,
464 const Index& za_index,
470 const Agenda& propmat_clearsky_agenda,
473 const Agenda& ppath_step_agenda,
475 const Numeric& ppath_lraytrace,
483 const Index& f_index,
487 const Agenda& surface_rtprop_agenda,
488 const Index& scat_za_interp,
500 ppath_step.
pos(0, 0) = z_field(p_index, 0, 0);
501 ppath_step.
r[0] = refellipsoid[0] + z_field(p_index, 0, 0);
504 ppath_step.
los(0, 0) = za_grid[za_index];
507 ppath_step.
gp_p[0].idx = p_index;
508 ppath_step.
gp_p[0].fd[0] = 0;
509 ppath_step.
gp_p[0].fd[1] = 1;
516 Vector(1, f_grid[f_index]),
524 if ((cloudbox_limits[0] <= ppath_step.
gp_p[1].idx &&
525 cloudbox_limits[1] > ppath_step.
gp_p[1].idx) ||
526 (cloudbox_limits[1] == ppath_step.
gp_p[1].idx &&
527 abs(ppath_step.
gp_p[1].fd[0]) < 1e-6)) {
529 const Index stokes_dim = cloudbox_field_mono.
ncols();
540 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.
np, 0.);
541 Matrix abs_vec_int(stokes_dim, ppath_step.
np, 0.);
542 Matrix sca_vec_int(stokes_dim, ppath_step.
np, 0.);
543 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.
np, 0.);
545 Matrix vmr_list_int(N_species, ppath_step.
np, 0.);
551 cloudbox_field_mono_int,
558 cloudbox_field_mono_old,
583 propmat_clearsky_agenda,
590 cloudbox_field_mono_int,
605 surface_rtprop_agenda,
619 const Index& p_index,
620 const Index& za_index,
625 const Agenda& propmat_clearsky_agenda,
633 const Index& f_index,
641 const Index stokes_dim = cloudbox_field_mono.
ncols();
642 const Index atmosphere_dim = 1;
646 Matrix matrix_tmp(stokes_dim, stokes_dim);
647 Vector vector_tmp(stokes_dim);
648 Vector rtp_vmr(N_species, 0.);
650 Vector sca_vec_av(stokes_dim, 0);
655 Vector stokes_vec(stokes_dim);
657 if ((p_index == 0) && (za_grid[za_index] > 90)) {
665 if (za_grid[za_index] <= 90.0) {
666 stokes_vec = cloudbox_field_mono(
667 p_index - cloudbox_limits[0] + 1, 0, 0, za_index, 0,
joker);
668 Numeric z_field_above = z_field(p_index + 1, 0, 0);
669 Numeric z_field_0 = z_field(p_index, 0, 0);
672 if (za_grid[za_index] == 90.0) {
678 cos_theta = cos(za_grid[za_index] *
PI / 180.0);
680 Numeric dz = (z_field_above - z_field_0);
682 lstep =
abs(dz / cos_theta);
686 0.5 * (t_field(p_index, 0, 0) + t_field(p_index + 1, 0, 0));
689 Numeric rtp_pressure = 0.5 * (p_grid[p_index] + p_grid[p_index + 1]);
692 for (
Index j = 0; j < N_species; j++)
693 rtp_vmr[j] = 0.5 * (vmr_field(j, p_index, 0, 0) +
694 vmr_field(j, p_index + 1, 0, 0));
700 const Vector rtp_mag_dummy(3, 0);
701 const Vector ppath_los_dummy;
713 f_grid[
Range(f_index, 1)],
720 propmat_clearsky_agenda);
724 for (
Index k = 0; k < stokes_dim; k++) {
731 p_index - cloudbox_limits[0], 0, 0, za_index, 0, k) +
733 p_index - cloudbox_limits[0] + 1, 0, 0, za_index, 0, k));
740 abs_vec_field(p_index - cloudbox_limits[0], 0, 0,
joker),
741 abs_vec_field(p_index - cloudbox_limits[0] + 1, 0, 0,
joker));
746 ext_mat_field(p_index - cloudbox_limits[0], 0, 0,
joker,
joker),
747 ext_mat_field(p_index - cloudbox_limits[0] + 1, 0, 0,
joker,
joker));
757 if (out3.sufficient_priority()) {
760 out3 <<
"-----------------------------------------\n";
761 out3 <<
"Input for radiative transfer step \n"
762 <<
"calculation inside"
765 out3 <<
"Stokes vector at intersection point: \n" << stokes_vec <<
"\n";
766 out3 <<
"lstep: ..." << lstep <<
"\n";
767 out3 <<
"------------------------------------------\n";
768 out3 <<
"Averaged coefficients: \n";
769 out3 <<
"Planck function: " << rte_planck_value <<
"\n";
770 out3 <<
"Scattering vector: " << sca_vec_av <<
"\n";
771 out3 <<
"Absorption vector: " << vector_tmp <<
"\n";
772 out3 <<
"Extinction matrix: " << matrix_tmp <<
"\n";
780 Matrix(stokes_dim, stokes_dim),
789 p_index - cloudbox_limits[0], 0, 0, za_index, 0,
joker) =
792 if (za_grid[za_index] > 90) {
793 stokes_vec = cloudbox_field_mono(
794 p_index - cloudbox_limits[0] - 1, 0, 0, za_index, 0,
joker);
795 Numeric z_field_below = z_field(p_index - 1, 0, 0);
796 Numeric z_field_0 = z_field(p_index, 0, 0);
799 if (za_grid[za_index] == 90.0) {
800 cos_theta = cos((za_grid[za_index] - 1) *
PI / 180.);
804 cos_theta = cos(za_grid[za_index] *
PI / 180.0);
806 Numeric dz = (z_field_0 - z_field_below);
807 lstep =
abs(dz / cos_theta);
811 0.5 * (t_field(p_index, 0, 0) + t_field(p_index - 1, 0, 0));
814 Numeric rtp_pressure = 0.5 * (p_grid[p_index] + p_grid[p_index - 1]);
818 for (
Index k = 0; k < N_species; k++)
819 rtp_vmr[k] = 0.5 * (vmr_field(k, p_index, 0, 0) +
820 vmr_field(k, p_index - 1, 0, 0));
826 const Vector rtp_mag_dummy(3, 0);
827 const Vector ppath_los_dummy;
838 f_grid[
Range(f_index, 1)],
845 propmat_clearsky_agenda);
855 for (
Index k = 0; k < stokes_dim; k++) {
862 p_index - cloudbox_limits[0], 0, 0, za_index, 0, k) +
864 p_index - cloudbox_limits[0] - 1, 0, 0, za_index, 0, k));
870 abs_vec_field(p_index - cloudbox_limits[0], 0, 0,
joker),
871 abs_vec_field(p_index - cloudbox_limits[0] - 1, 0, 0,
joker));
876 ext_mat_field(p_index - cloudbox_limits[0], 0, 0,
joker,
joker),
877 ext_mat_field(p_index - cloudbox_limits[0] + 1, 0, 0,
joker,
joker));
887 if (out3.sufficient_priority()) {
890 out3 <<
"-----------------------------------------\n";
891 out3 <<
"Input for radiative transfer step \n"
892 <<
"calculation inside"
895 out3 <<
"Stokes vector at intersection point: \n" << stokes_vec <<
"\n";
896 out3 <<
"lstep: ..." << lstep <<
"\n";
897 out3 <<
"------------------------------------------\n";
898 out3 <<
"Averaged coefficients: \n";
899 out3 <<
"Planck function: " << rte_planck_value <<
"\n";
900 out3 <<
"Scattering vector: " << sca_vec_av <<
"\n";
901 out3 <<
"Absorption vector: " << vector_tmp <<
"\n";
902 out3 <<
"Extinction matrix: " << matrix_tmp <<
"\n";
910 Matrix(stokes_dim, stokes_dim),
919 p_index - cloudbox_limits[0], 0, 0, za_index, 0,
joker) =
926 else if (bkgr == 2) {
930 Vector rte_pos(atmosphere_dim);
931 Numeric z_field_0 = z_field(0, 0, 0);
935 rte_los = za_grid[za_index];
947 "Surface reflections inside cloud box not yet handled.");
1109 const Index& p_index,
1110 const Index& lat_index,
1111 const Index& lon_index,
1112 const Index& za_index,
1113 const Index& aa_index,
1119 const Agenda& propmat_clearsky_agenda,
1122 const Agenda& ppath_step_agenda,
1124 const Numeric& ppath_lraytrace,
1133 const Index& f_index,
1142 const Index stokes_dim = cloudbox_field_mono.
ncols();
1144 Vector sca_vec_av(stokes_dim, 0);
1148 aa_g[i] = aa_grid[i] - 180.;
1160 ppath_step.
pos(0, 2) = lon_grid[lon_index];
1161 ppath_step.
pos(0, 1) = lat_grid[lat_index];
1162 ppath_step.
pos(0, 0) = z_field(p_index, lat_index, lon_index);
1165 refell2r(refellipsoid, ppath_step.
pos(0, 1)) + ppath_step.
pos(0, 0);
1168 ppath_step.
los(0, 0) = za_grid[za_index];
1169 ppath_step.
los(0, 1) = aa_g[aa_index];
1172 ppath_step.
gp_p[0].idx = p_index;
1173 ppath_step.
gp_p[0].fd[0] = 0.;
1174 ppath_step.
gp_p[0].fd[1] = 1.;
1176 ppath_step.
gp_lat[0].idx = lat_index;
1177 ppath_step.
gp_lat[0].fd[0] = 0.;
1178 ppath_step.
gp_lat[0].fd[1] = 1.;
1180 ppath_step.
gp_lon[0].idx = lon_index;
1181 ppath_step.
gp_lon[0].fd[0] = 0.;
1182 ppath_step.
gp_lon[0].fd[1] = 1.;
1189 Vector(1, f_grid[f_index]),
1206 for (
Index i = 0; i < ppath_step.
np; i++) {
1207 cloud_gp_p[i].idx -= cloudbox_limits[0];
1208 cloud_gp_lat[i].idx -= cloudbox_limits[2];
1209 cloud_gp_lon[i].idx -= cloudbox_limits[4];
1211 const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
1212 const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
1213 const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
1234 for (
Index i = 0; i < los_grid_aa.
nelem(); i++)
1235 los_grid_aa[i] = los_grid_aa[i] + 180.;
1238 gridpos(gp_za, za_grid, los_grid_za);
1241 gridpos(gp_aa, aa_grid, los_grid_aa);
1243 Matrix itw_p_za(ppath_step.
np, 32);
1245 itw_p_za, cloud_gp_p, cloud_gp_lat, cloud_gp_lon, gp_za, gp_aa);
1253 Tensor3 ext_mat_int(stokes_dim, stokes_dim, ppath_step.
np);
1254 Matrix abs_vec_int(stokes_dim, ppath_step.
np);
1255 Matrix sca_vec_int(stokes_dim, ppath_step.
np, 0.);
1256 Matrix cloudbox_field_mono_int(stokes_dim, ppath_step.
np, 0.);
1260 Vector stokes_vec(stokes_dim);
1268 for (
Index i = 0; i < stokes_dim; i++) {
1271 out3 <<
"Interpolate ext_mat:\n";
1272 for (
Index j = 0; j < stokes_dim; j++) {
1298 out3 <<
"Interpolate doit_scat_field:\n";
1307 out3 <<
"Interpolate cloudbox_field_mono:\n";
1322 out3 <<
"Interpolate temperature field\n";
1339 Matrix vmr_list_int(N_species, ppath_step.
np);
1341 for (
Index i = 0; i < N_species; i++) {
1342 out3 <<
"Interpolate vmr field\n";
1350 vmr_list_int(i,
joker) = vmr_int;
1354 itw2p(p_int, p_grid, ppath_step.
gp_p, itw_p);
1356 out3 <<
"Calculate radiative transfer inside cloudbox.\n";
1358 cloudbox_field_mono,
1359 propmat_clearsky_agenda,
1366 cloudbox_field_mono_int,
1384 const Agenda& propmat_clearsky_agenda,
1385 const Ppath& ppath_step,
1395 const Index& f_index,
1396 const Index& p_index,
1397 const Index& lat_index,
1398 const Index& lon_index,
1399 const Index& za_index,
1400 const Index& aa_index,
1404 const Index N_species = vmr_list_int.
nrows();
1405 const Index stokes_dim = cloudbox_field_mono.
ncols();
1406 const Index atmosphere_dim = cloudbox_limits.
nelem() / 2;
1408 Vector sca_vec_av(stokes_dim, 0);
1409 Vector stokes_vec(stokes_dim, 0.);
1411 Vector rtp_vmr_local(N_species, 0.);
1419 Matrix matrix_tmp(stokes_dim, stokes_dim);
1420 Vector vector_tmp(stokes_dim);
1423 stokes_vec = cloudbox_field_mono_int(
joker, ppath_step.
np - 1);
1425 for (
Index k = ppath_step.
np - 1; k >= 0; k--) {
1428 swap(cur_propmat_clearsky, prev_propmat_clearsky);
1433 const Vector rtp_mag_dummy(3, 0);
1434 const Vector ppath_los_dummy;
1440 cur_propmat_clearsky,
1446 f_grid[
Range(f_index, 1)],
1452 vmr_list_int(
joker, k),
1453 propmat_clearsky_agenda);
1457 if (k == ppath_step.
np - 1)
continue;
1460 prev_propmat_clearsky += cur_propmat_clearsky;
1461 prev_propmat_clearsky *= 0.5;
1464 ext_mat_local, abs_vec_local, prev_propmat_clearsky);
1466 for (
Index i = 0; i < stokes_dim; i++) {
1470 sca_vec_av[i] = 0.5 * (sca_vec_int(i, k) + sca_vec_int(i, k + 1));
1477 abs_vec_int(
joker, k + 1));
1490 Numeric rte_planck_value =
planck(f, 0.5 * (t_int[k] + t_int[k + 1]));
1496 if (out3.sufficient_priority()) {
1499 out3 <<
"-----------------------------------------\n";
1500 out3 <<
"Input for radiative transfer step \n"
1501 <<
"calculation inside"
1504 out3 <<
"Stokes vector at intersection point: \n" << stokes_vec <<
"\n";
1505 out3 <<
"lstep: ..." << lstep <<
"\n";
1506 out3 <<
"------------------------------------------\n";
1507 out3 <<
"Averaged coefficients: \n";
1508 out3 <<
"Planck function: " << rte_planck_value <<
"\n";
1509 out3 <<
"Scattering vector: " << sca_vec_av <<
"\n";
1510 out3 <<
"Absorption vector: " << vector_tmp <<
"\n";
1511 out3 <<
"Extinction matrix: " << matrix_tmp <<
"\n";
1519 Matrix(stokes_dim, stokes_dim),
1528 if (atmosphere_dim == 1)
1529 cloudbox_field_mono(
1530 p_index - cloudbox_limits[0], 0, 0, za_index, 0,
joker) =
1532 else if (atmosphere_dim == 3)
1533 cloudbox_field_mono(p_index - cloudbox_limits[0],
1534 lat_index - cloudbox_limits[2],
1535 lon_index - cloudbox_limits[4],
1538 joker) = stokes_vec;
1545 const Agenda& surface_rtprop_agenda,
1547 const Index& f_index,
1548 const Index& stokes_dim,
1549 const Ppath& ppath_step,
1552 const Index& za_index) {
1553 chk_not_empty(
"surface_rtprop_agenda", surface_rtprop_agenda);
1569 rte_pos = ppath_step.
pos(np - 1,
Range(0, ppath_step.
dim));
1573 rte_los = ppath_step.
los(np - 1,
joker);
1582 Vector(1, f_grid[f_index]),
1585 surface_rtprop_agenda);
1587 iy = surface_emission;
1594 for (
Index ilos = 0; ilos < nlos; ilos++) {
1601 cloudbox_field_mono(cloudbox_limits[0],
1604 (za_grid.
nelem() - 1 - za_index),
1607 iy(0,
joker) += rtmp;
1610 cloudbox_field_mono(cloudbox_limits[0], 0, 0, za_index, 0,
joker) =
1616 const Index& accelerated,
1622 for (
Index i = 0; i < accelerated; ++i) {
1657 for (
Index p_index = 0; p_index < N_p; ++p_index) {
1658 for (
Index za_index = 0; za_index < N_za; ++za_index) {
1659 A1 += Q1(p_index, za_index) * Q1(p_index, za_index) *
1660 J(p_index, za_index);
1661 A2B1 += Q2(p_index, za_index) * Q1(p_index, za_index) *
1662 J(p_index, za_index);
1663 B2 += Q2(p_index, za_index) * Q2(p_index, za_index) *
1664 J(p_index, za_index);
1665 C1 += Q1(p_index, za_index) * Q3(p_index, za_index) *
1666 J(p_index, za_index);
1667 C2 += Q2(p_index, za_index) * Q3(p_index, za_index) *
1668 J(p_index, za_index);
1672 NGA = (C1 * B2 - C2 * A2B1) / (A1 * B2 - A2B1 * A2B1);
1673 NGB = (C2 * A1 - C1 * A2B1) / (A1 * B2 - A2B1 * A2B1);
1675 if (!std::isnan(NGB) && !std::isnan(NGA)) {
1677 for (
Index p_index = 0; p_index < N_p; ++p_index) {
1678 for (
Index za_index = 0; za_index < N_za; ++za_index) {
1679 Q1(p_index, za_index) = (1 - NGA - NGB) * S4(p_index, za_index) +
1680 NGA * S3(p_index, za_index) +
1681 NGB * S2(p_index, za_index);
1684 cloudbox_field_mono(
joker, 0, 0,
joker, 0, i) = Q1;
1705 const Ppath& ppath_step,
1708 const Index& scat_za_interp,
1713 const Index stokes_dim = cloudbox_field_mono.
ncols();
1721 for (
Index i = 0; i < ppath_step.
np; i++)
1722 cloud_gp_p[i].idx -= cloudbox_limits[0];
1725 const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
1738 gridpos(gp_za, za_grid, los_grid);
1747 for (
Index i = 0; i < stokes_dim; i++) {
1750 out3 <<
"Interpolate ext_mat:\n";
1751 for (
Index j = 0; j < stokes_dim; j++) {
1757 ext_mat_field(
joker, 0, 0, i, j),
1764 out3 <<
"Interpolate abs_vec:\n";
1766 abs_vec_int(i,
joker), itw, abs_vec_field(
joker, 0, 0, i), cloud_gp_p);
1772 out3 <<
"Interpolate doit_scat_field and cloudbox_field_mono:\n";
1773 if (scat_za_interp == 0)
1782 cloudbox_field_mono(
joker, 0, 0,
joker, 0, i),
1785 }
else if (scat_za_interp == 1)
1791 stokes_dim, ppath_step.
np, za_grid.
nelem(), 0.);
1792 Tensor3 cloudbox_field_mono_int_za(
1793 stokes_dim, ppath_step.
np, za_grid.
nelem(), 0.);
1794 for (
Index za = 0; za < za_grid.
nelem(); za++) {
1797 doit_scat_field(
joker, 0, 0, za, 0, i),
1799 out3 <<
"Interpolate cloudbox_field_mono:\n";
1802 cloudbox_field_mono(
joker, 0, 0, za, 0, i),
1805 for (
Index ip = 0; ip < ppath_step.
np; ip++) {
1807 sca_vec_int_za(i, ip,
joker),
1810 cloudbox_field_mono_int(i, ip) =
1812 cloudbox_field_mono_int_za(i, ip,
joker),
1823 out3 <<
"Interpolate temperature field\n";
1836 for (
Index i_sp = 0; i_sp < N_species; i_sp++) {
1837 out3 <<
"Interpolate vmr field\n";
1839 vmr_list_int(i_sp,
joker) = vmr_int;
1844 itw2p(p_int, p_grid, ppath_step.
gp_p, itw);
1849 Matrix& cloudbox_field_opt,
1854 const Index& scat_za_interp) {
1861 Vector i_approx_interp(N_za);
1866 idx.push_back(N_za - 1);
1881 while (max_diff > acc) {
1889 for (
Index i_p = 0; i_p < N_p; i_p++) {
1890 for (
Index i_za_red = 0; i_za_red < idx.
nelem(); i_za_red++) {
1891 za_reduced[i_za_red] = za_grid_fine[idx[i_za_red]];
1892 cloudbox_field_opt(i_p, i_za_red) =
1893 cloudbox_field_mono(i_p, 0, 0, idx[i_za_red], 0, 0);
1896 gridpos(gp_za, za_reduced, za_grid_fine);
1898 if (scat_za_interp == 0 || idx.
nelem() < 3) {
1900 interp(i_approx_interp, itw, cloudbox_field_opt(i_p,
joker), gp_za);
1901 }
else if (scat_za_interp == 1) {
1902 for (
Index i_za = 0; i_za < N_za; i_za++) {
1904 cloudbox_field_opt(i_p,
joker),
1914 for (
Index i_za = 0; i_za < N_za; i_za++) {
1915 diff_vec[i_za] =
abs(cloudbox_field_mono(i_p, 0, 0, i_za, 0, 0) -
1916 i_approx_interp[i_za]);
1917 if (diff_vec[i_za] > max_diff_za[i_p]) {
1918 max_diff_za[i_p] = diff_vec[i_za];
1923 if (max_diff_za[i_p] > max_diff_p) {
1924 max_diff_p = max_diff_za[i_p];
1931 max_diff_p / cloudbox_field_mono(ind_p, 0, 0, ind_za[ind_p], 0, 0) * 100.;
1933 idx.push_back(ind_za[ind_p]);
1936 i_sort.resize(idx_unsorted.
nelem());
1939 for (
Index i = 0; i < idx_unsorted.
nelem(); i++)
1940 idx[i] = idx_unsorted[i_sort[i]];
1948 za_grid_opt[i] = za_grid_fine[idx[i]];
1949 cloudbox_field_opt(
joker, i) = cloudbox_field_mono(
joker, 0, 0, idx[i], 0, 0);
1955 const Tensor6& cloudbox_field_mono,
1957 const Agenda& spt_calc_agenda,
1958 const Index& atmosphere_dim,
1963 const Numeric& norm_error_threshold,
1964 const Index& norm_debug,
1967 "Only 1D is supported here for now");
1976 "The range of *za_grid* must [0 180].");
1982 "The range of *aa_grid* must [0 360].");
1985 const Index stokes_dim = doit_scat_field.
ncols();
1990 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1997 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
2003 doit_scat_field.
nbooks(),
2004 doit_scat_field.
npages(),
2005 doit_scat_field.
nrows(),
2008 Index aa_index_local = 0;
2011 for (
Index za_index_local = 0; za_index_local < Nza;
2027 for (
Index p_index = 0;
2028 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
2034 for (
Index i = 0; i < stokes_dim; i++) {
2035 doit_scat_ext_field(p_index, 0, 0, za_index_local, 0) +=
2036 cloudbox_field_mono(p_index, 0, 0, za_index_local, 0, i) *
2037 (ext_mat_field(p_index, 0, 0, 0, i) -
2038 abs_vec_field(p_index, 0, 0, i));
2044 Index corr_max_p_index = -1;
2046 for (
Index p_index = 0; p_index < Np; p_index++) {
2049 doit_scat_field(p_index, 0, 0,
joker, 0, 0), za_grid);
2052 doit_scat_ext_field(p_index, 0, 0,
joker, 0), za_grid);
2056 const Numeric corr_factor = scat_ext_int / scat_int;
2060 if (!std::isnan(corr_factor) && !std::isinf(corr_factor)) {
2061 if (
abs(corr_factor) >
abs(corr_max)) {
2062 corr_max = corr_factor;
2063 corr_max_p_index = p_index;
2066 out0 <<
" DOIT corr_factor: " << 1. - corr_factor
2067 <<
" ( scat_ext_int: " << scat_ext_int
2068 <<
", scat_int: " << scat_int <<
")"
2069 <<
" at p_index " << p_index <<
"\n";
2072 "ERROR: DOIT correction factor exceeds threshold (=",
2073 norm_error_threshold,
"): ", setprecision(4),
2074 1. - corr_factor,
" at p_index ", p_index,
"\n")
2075 if (
abs(1. - corr_factor) > norm_error_threshold / 2.) {
2076 out0 <<
" WARNING: DOIT correction factor above threshold/2: "
2077 << 1. - corr_factor <<
" at p_index " << p_index <<
"\n";
2081 doit_scat_field(p_index, 0, 0,
joker, 0,
joker) *= corr_factor;
2082 }
else if (norm_debug) {
2083 out0 <<
" DOIT corr_factor ignored: " << 1. - corr_factor
2084 <<
" ( scat_ext_int: " << scat_ext_int <<
", scat_int: " << scat_int
2086 <<
" at p_index " << p_index <<
"\n";
2091 if (corr_max_p_index != -1) {
2092 os <<
" Max. DOIT correction factor in this iteration: " << 1. - corr_max
2093 <<
" at p_index " << corr_max_p_index <<
"\n";
2095 os <<
" No DOIT correction performed in this iteration.\n";
Declarations for agendas.
This file contains the definition of Array.
Constants of physical expressions as constexpr.
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
void spt_calc_agendaExecute(Workspace &ws, ArrayOfPropagationMatrix &ext_mat_spt, ArrayOfStokesVector &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Numeric rtp_temperature, const Index za_index, const Index aa_index, const Agenda &input_agenda)
void propmat_clearsky_agendaExecute(Workspace &ws, PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfSpeciesTag &select_abs_species, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
void surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
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.
A constant view of a Tensor4.
Index ncols() const noexcept
Index nbooks() const noexcept
A constant view of a Tensor5.
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Tensor6.
Index nbooks() const noexcept
Index nvitrines() const noexcept
Index ncols() const noexcept
Index npages() const noexcept
Index nshelves() const noexcept
Index nrows() 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.
void MatrixAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Sets the dense matrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
void AddAverageAtPosition(const ConstMatrixView &mat1, const ConstMatrixView &mat2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of the two input at position.
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
Stokes vector is as Propagation matrix but only has 4 possible values.
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
void AddAverageAtPosition(const ConstVectorView &vec1, const ConstVectorView &vec2, const Index iv=0, const Index iz=0, const Index ia=0)
Add the average of both inputs at position.
void resize(Index n)
Resize function.
bool is_inside_cloudbox(const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, const bool include_boundaries)
Internal cloudbox functions.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
void za_gridOpt(Vector &za_grid_opt, Matrix &cloudbox_field_opt, ConstVectorView za_grid_fine, ConstTensor6View cloudbox_field_mono, const Numeric &acc, const Index &scat_za_interp)
Optimize the zenith angle grid.
void doit_scat_fieldNormalize(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &cloudbox_field_mono, const ArrayOfIndex &cloudbox_limits, const Agenda &spt_calc_agenda, const Index &atmosphere_dim, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Tensor3 &t_field, const Numeric &norm_error_threshold, const Index &norm_debug, const Verbosity &verbosity)
Normalization of scattered field.
constexpr Numeric RAD2DEG
void cloud_fieldsCalc(Workspace &ws, Tensor5View ext_mat_field, Tensor4View abs_vec_field, const Agenda &spt_calc_agenda, const Index &za_index, const Index &aa_index, const ArrayOfIndex &cloudbox_limits, ConstTensor3View t_field, ConstTensor4View pnd_field, const Verbosity &verbosity)
Calculate ext_mat, abs_vec for all points inside the cloudbox for one.
void cloud_ppath_update1D_noseq(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View cloudbox_field_mono_old, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Agenda &surface_rtprop_agenda, const Index &scat_za_interp, const Verbosity &verbosity)
Calculation of radiation field along a propagation path step for specified zenith direction and press...
void cloud_RT_no_background(Workspace &ws, Tensor6View cloudbox_field_mono, const Agenda &propmat_clearsky_agenda, const Ppath &ppath_step, ConstVectorView t_int, ConstMatrixView vmr_list_int, ConstTensor3View ext_mat_int, ConstMatrixView abs_vec_int, ConstMatrixView sca_vec_int, ConstMatrixView cloudbox_field_mono_int, ConstVectorView p_int, const ArrayOfIndex &cloudbox_limits, ConstVectorView f_grid, const Index &f_index, const Index &p_index, const Index &lat_index, const Index &lon_index, const Index &za_index, const Index &aa_index, const Verbosity &verbosity)
Calculates RT in the cloudbox (1D)
void cloudbox_field_ngAcceleration(Tensor6 &cloudbox_field_mono, const ArrayOfTensor6 &acceleration_input, const Index &accelerated, const Verbosity &)
Convergence acceleration.
void interp_cloud_coeff1D(Tensor3View ext_mat_int, MatrixView abs_vec_int, MatrixView sca_vec_int, MatrixView cloudbox_field_mono_int, VectorView t_int, MatrixView vmr_list_int, VectorView p_int, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, ConstTensor6View doit_scat_field, ConstTensor6View cloudbox_field_mono, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView p_grid, const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, ConstVectorView za_grid, const Index &scat_za_interp, const Verbosity &verbosity)
Interpolate all inputs of the VRTE on a propagation path step.
void rte_step_doit_replacement(VectorView stokes_vec, MatrixView trans_mat, const PropagationMatrix &ext_mat_av, const StokesVector &abs_vec_av, ConstVectorView sca_vec_av, const Numeric &lstep, const Numeric &rtp_planck_value, const bool &trans_is_precalc)
Solves monochromatic VRTE for an atmospheric slab with constant conditions.
void cloud_RT_surface(Workspace &ws, Tensor6View cloudbox_field_mono, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, const Index &f_index, const Index &stokes_dim, const Ppath &ppath_step, const ArrayOfIndex &cloudbox_limits, ConstVectorView za_grid, const Index &za_index)
Calculates RT in the cloudbox.
void cloud_ppath_update3D(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &lat_index, const Index &lon_index, const Index &za_index, const Index &aa_index, ConstVectorView za_grid, ConstVectorView aa_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Index &, const Verbosity &verbosity)
Radiative transfer calculation along a path inside the cloudbox (3D).
void cloud_ppath_update1D_planeparallel(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, ConstVectorView p_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Verbosity &verbosity)
Radiative transfer calculation inside cloudbox for planeparallel case.
void cloud_ppath_update1D(Workspace &ws, Tensor6View cloudbox_field_mono, const Index &p_index, const Index &za_index, ConstVectorView za_grid, const ArrayOfIndex &cloudbox_limits, ConstTensor6View doit_scat_field, const Agenda &propmat_clearsky_agenda, ConstTensor4View vmr_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, ConstVectorView p_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstTensor3View t_field, ConstVectorView f_grid, const Index &f_index, ConstTensor5View ext_mat_field, ConstTensor4View abs_vec_field, const Agenda &surface_rtprop_agenda, const Index &scat_za_interp, const Verbosity &verbosity)
Calculates radiation field along a propagation path step for specified zenith direction and pressure ...
Radiative transfer in cloudbox.
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
void id_mat(MatrixView I)
Identity Matrix.
Linear algebra functions.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
bool is_singular(ConstMatrixView A)
Checks if a square matrix is singular.
Header file for logic.cc.
void opt_prop_bulkCalc(PropagationMatrix &ext_mat, StokesVector &abs_vec, const ArrayOfPropagationMatrix &ext_mat_spt, const ArrayOfStokesVector &abs_vec_spt, const Tensor4 &pnd_field, const Index &scat_p_index, const Index &scat_lat_index, const Index &scat_lon_index, const Verbosity &)
WORKSPACE METHOD: opt_prop_bulkCalc.
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
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.
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
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 auto rad2deg(auto x) noexcept
Converts radians to degrees.
void opt_prop_sum_propmat_clearsky(PropagationMatrix &ext_mat, StokesVector &abs_vec, const PropagationMatrix &propmat_clearsky)
Get optical properties from propmat_clearsky.
Numeric planck(const Numeric &f, const Numeric &t)
planck
This file contains declerations of functions of physical character.
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
Initiates a Ppath structure to hold the given number of points.
Propagation path structure and functions.
void compute_transmission_matrix_from_averaged_matrix_at_frequency(MatrixView T, const Numeric &r, const PropagationMatrix &averaged_propagation_matrix, const Index iv, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Stuff related to the propagation matrix.
Declaration of functions in rte.cc.
Contains sorting routines.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Header file for special_interp.cc.
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Vector lstep
The length between ppath points.
Vector r
Radius of each ppath point.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Index dim
Atmospheric dimensionality.
This file contains basic functions to handle XML data files.