73 Index& doit_za_grid_size,
77 const Index& N_za_grid,
78 const Index& N_aa_grid,
79 const String& za_grid_opt_file,
87 "N_aa_grid must be > 0 (even for 1D / DISORT cases).")
95 "N_za_grid must be >= 0.")
96 doit_za_grid_size = N_za_grid;
98 if (za_grid_opt_file ==
"")
103 "N_za_grid must be >1 or =0 (the latter only allowed for RT4).")
112 Index& doit_conv_flag,
113 Index& doit_iteration_counter,
116 const Tensor6& cloudbox_field_mono_old,
119 const Index& max_iterations,
120 const Index& throw_nonconv_error,
127 "Convergence flag is non-zero, which means that this\n"
128 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
129 "be used only in *doit_conv_test_agenda*\n");
135 const Index N_aa = cloudbox_field_mono.
nrows();
136 const Index stokes_dim = cloudbox_field_mono.
ncols();
140 "You have to specify limiting values for the "
141 "convergence test for each Stokes component "
142 "separately. That means that *epsilon* must "
143 "have *stokes_dim* elements!");
147 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
148 "The fields (Tensor6) *cloudbox_field* and \n"
149 "*cloudbox_field_old* which are compared in the \n"
150 "convergence test do not have the same size.\n");
154 doit_iteration_counter += 1;
155 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
157 if (doit_iteration_counter > max_iterations) {
159 out <<
"Method does not converge (number of iterations \n"
160 <<
"is > " << max_iterations <<
"). Either the cloud "
161 <<
"particle number density \n"
162 <<
"is too large or the numerical setup for the DOIT \n"
163 <<
"calculation is not correct. In case of limb \n"
164 <<
"simulations please make sure that you use an \n"
165 <<
"optimized zenith angle grid. \n"
166 <<
"*cloudbox_field* might be wrong.\n";
167 if (throw_nonconv_error != 0) {
173 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
174 cloudbox_field_mono = NAN;
177 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
182 for (
Index p_index = 0; p_index < N_p; p_index++) {
183 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
184 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
185 for (
Index za_index = 0; za_index < N_za; za_index++) {
186 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
187 for (
Index stokes_index = 0; stokes_index < stokes_dim;
189 Numeric diff = (cloudbox_field_mono(p_index,
195 cloudbox_field_mono_old(p_index,
206 if (
abs(diff) > epsilon[stokes_index]) {
207 out1 <<
"difference: " << diff <<
"\n";
225 Index& doit_conv_flag,
226 Index& doit_iteration_counter,
229 const Tensor6& cloudbox_field_mono_old,
231 const Index& f_index,
234 const Index& max_iterations,
235 const Index& throw_nonconv_error,
243 "Convergence flag is non-zero, which means that this \n"
244 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
245 "be used only in *doit_conv_test_agenda*\n");
251 const Index N_aa = cloudbox_field_mono.
nrows();
252 const Index stokes_dim = cloudbox_field_mono.
ncols();
256 "You have to specify limiting values for the "
257 "convergence test for each Stokes component "
258 "separately. That means that *epsilon* must "
259 "have *stokes_dim* elements!");
263 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
264 "The fields (Tensor6) *cloudbox_field* and \n"
265 "*cloudbox_field_old* which are compared in the \n"
266 "convergence test do not have the same size.\n");
275 "*f_index* is greater than number of elements in the\n"
276 "frequency grid.\n");
280 doit_iteration_counter += 1;
281 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
284 if (doit_iteration_counter > max_iterations) {
286 out <<
"At frequency " << f_grid[f_index] <<
" GHz \n"
287 <<
"method does not converge (number of iterations \n"
288 <<
"is > " << max_iterations <<
"). Either the particle"
289 <<
" number density \n"
290 <<
"is too large or the numerical setup for the DOIT \n"
291 <<
"calculation is not correct. In case of limb \n"
292 <<
"simulations please make sure that you use an \n"
293 <<
"optimized zenith angle grid. \n"
294 <<
"*cloudbox_field* might be wrong.\n";
295 if (throw_nonconv_error != 0) {
301 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
302 cloudbox_field_mono = NAN;
305 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
310 for (
Index p_index = 0; p_index < N_p; p_index++) {
311 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
312 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
313 for (
Index za_index = 0; za_index < N_za; za_index++) {
314 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
315 for (
Index stokes_index = 0; stokes_index < stokes_dim;
317 Numeric diff = cloudbox_field_mono(p_index,
323 cloudbox_field_mono_old(p_index,
336 if (
abs(diff_bt) > epsilon[stokes_index]) {
337 out1 <<
"BT difference: " << diff_bt <<
" in stokes dim "
338 << stokes_index <<
"\n";
360 Index& doit_conv_flag,
361 Index& doit_iteration_counter,
364 const Tensor6& cloudbox_field_mono_old,
366 const Index& f_index,
369 const Index& max_iterations,
370 const Index& throw_nonconv_error,
378 "Convergence flag is non-zero, which means that this \n"
379 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
380 "be used only in *doit_conv_test_agenda*\n");
386 const Index N_aa = cloudbox_field_mono.
nrows();
387 const Index stokes_dim = cloudbox_field_mono.
ncols();
391 "You have to specify limiting values for the "
392 "convergence test for each Stokes component "
393 "separately. That means that *epsilon* must "
394 "have *stokes_dim* elements!");
398 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
399 "The fields (Tensor6) *cloudbox_field* and \n"
400 "*cloudbox_field_old* which are compared in the \n"
401 "convergence test do not have the same size.\n");
410 "*f_index* is greater than number of elements in the\n"
411 "frequency grid.\n");
415 doit_iteration_counter += 1;
416 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
418 if (doit_iteration_counter > max_iterations) {
420 out <<
"Method does not converge (number of iterations \n"
421 <<
"is > " << max_iterations <<
"). Either the"
422 <<
" particle number density \n"
423 <<
"is too large or the numerical setup for the DOIT \n"
424 <<
"calculation is not correct. In case of limb \n"
425 <<
"simulations please make sure that you use an \n"
426 <<
"optimized zenith angle grid. \n";
427 if (throw_nonconv_error != 0) {
433 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
434 cloudbox_field_mono = NAN;
437 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
447 for (
Index p_index = 0; p_index < N_p; p_index++) {
448 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
449 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
450 for (
Index za_index = 0; za_index < N_za; za_index++) {
451 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
454 p_index, lat_index, lon_index, za_index, aa_index, i) -
455 cloudbox_field_mono_old(p_index,
468 lqs[i] =
sqrt(lqs[i]);
469 lqs[i] /= (
Numeric)(N_p * N_lat * N_lon * N_za * N_aa);
474 if (lqs[i] >= epsilon[i]) doit_conv_flag = 0;
477 out1 <<
"lqs [I]: " << lqs[0] <<
"\n";
487 const Agenda& doit_scat_field_agenda,
488 const Agenda& doit_rte_agenda,
489 const Agenda& doit_conv_test_agenda,
490 const Index& accelerated,
497 chk_not_empty(
"doit_scat_field_agenda", doit_scat_field_agenda);
499 chk_not_empty(
"doit_conv_test_agenda", doit_conv_test_agenda);
504 for (
Index p = 0; p < cloudbox_field_mono.
npages(); p++)
505 for (
Index r = 0; r < cloudbox_field_mono.
nrows(); r++)
508 "*cloudbox_field_mono* contains at least one NaN value.\n"
509 "This indicates an improper initialization of *cloudbox_field*.");
516 Tensor6 cloudbox_field_mono_old_local;
517 Index doit_conv_flag_local;
518 Index doit_iteration_counter_local;
524 cloudbox_field_mono.
nbooks(),
525 cloudbox_field_mono.
npages(),
526 cloudbox_field_mono.
nrows(),
527 cloudbox_field_mono.
ncols(),
530 doit_conv_flag_local = 0;
531 doit_iteration_counter_local = 0;
535 acceleration_input.resize(4);
537 while (doit_conv_flag_local == 0) {
539 cloudbox_field_mono_old_local = cloudbox_field_mono;
544 out2 <<
" Execute doit_scat_field_agenda. \n";
546 ws, doit_scat_field_local, cloudbox_field_mono, doit_scat_field_agenda);
549 out2 <<
" Execute doit_rte_agenda. \n";
551 ws, cloudbox_field_mono, doit_scat_field_local, doit_rte_agenda);
555 doit_conv_flag_local,
556 doit_iteration_counter_local,
558 cloudbox_field_mono_old_local,
559 doit_conv_test_agenda);
562 if (accelerated > 0 && doit_conv_flag_local == 0) {
563 acceleration_input[(doit_iteration_counter_local - 1) % 4] =
566 if (doit_iteration_counter_local % 4 == 0) {
568 cloudbox_field_mono, acceleration_input, accelerated, verbosity);
580 const Tensor6& doit_scat_field,
583 const Agenda& propmat_clearsky_agenda,
586 const Agenda& spt_calc_agenda,
590 const Agenda& ppath_step_agenda,
592 const Numeric& ppath_lraytrace,
595 const Vector& refellipsoid,
599 const Index& f_index,
600 const Agenda& surface_rtprop_agenda,
601 const Index& doit_za_interp,
607 <<
" cloudbox_fieldUpdate1D: Radiative transfer calculation in cloudbox\n";
608 out2 <<
" ------------------------------------------------------------- \n";
617 "The cloudbox dimension is not 1D! \n"
618 "Do you really want to do a 1D calculation? \n"
619 "If not, use *cloudbox_fieldUpdateSeq3D*.\n");
625 "The range of *za_grid* must [0 180].");
628 "The length of *p_grid* must be >= 2.");
641 "*f_index* is greater than number of elements in the\n"
642 "frequency grid.\n");
645 "Interpolation method is not defined. Use \n"
646 "*doit_za_interpSet*.\n");
648 const Index stokes_dim = doit_scat_field.
ncols();
653 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
661 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
675 out3 <<
"Calculate optical properties of individual scattering elements\n";
685 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
692 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
694 Tensor6 cloudbox_field_old(cloudbox_field_mono);
697 Index aa_index_local = 0;
700 for (
Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
719 for (
Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
721 if ((p_index != 0) || (za_grid[za_index_local] <= 90.)) {
730 propmat_clearsky_agenda,
743 surface_rtprop_agenda,
760 const Agenda& propmat_clearsky_agenda,
763 const Agenda& spt_calc_agenda,
768 const Agenda& ppath_step_agenda,
770 const Numeric& ppath_lraytrace,
773 const Vector& refellipsoid,
777 const Index& f_index,
778 const Agenda& surface_rtprop_agenda,
779 const Index& doit_za_interp,
780 const Index& normalize,
781 const Numeric& norm_error_threshold,
782 const Index& norm_debug,
788 <<
" cloudbox_fieldUpdateSeq1D: Radiative transfer calculation in cloudbox\n";
789 out2 <<
" ------------------------------------------------------------- \n";
794 chk_not_empty(
"propmat_clearsky_agenda", propmat_clearsky_agenda);
799 "The cloudbox dimension is not 1D! \n"
800 "Do you really want to do a 1D calculation? \n"
801 "For 3D use *cloudbox_fieldUpdateSeq3D*.\n");
807 "The range of *za_grid* must [0 180].");
810 "The length of *p_grid* must be >= 2.");
823 "*f_index* is greater than number of elements in the\n"
824 "frequency grid.\n");
827 "Interpolation method is not defined. Use \n"
828 "*doit_za_interpSet*.\n");
830 const Index stokes_dim = doit_scat_field.
ncols();
835 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
843 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
857 out3 <<
"Calculate optical properties of individual scattering elements\n";
867 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
874 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
879 180. - asin((refellipsoid[0] + z_field(cloudbox_limits[0], 0, 0)) /
880 (refellipsoid[0] + z_field(cloudbox_limits[1], 0, 0))) *
890 Matrix cloudbox_field_limb;
893 Index aa_index_local = 0;
907 norm_error_threshold,
913 for (
Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
933 if (za_grid[za_index_local] <= 90.) {
939 for (
Index p_index = cloudbox_limits[1] - 1;
940 p_index >= cloudbox_limits[0];
949 propmat_clearsky_agenda,
962 surface_rtprop_agenda,
966 }
else if (za_grid[za_index_local] >= theta_lim) {
970 for (
Index p_index = cloudbox_limits[0] + 1;
971 p_index <= cloudbox_limits[1];
980 propmat_clearsky_agenda,
993 surface_rtprop_agenda,
1007 bool conv_flag =
false;
1009 while (!conv_flag && limb_it < 10) {
1011 cloudbox_field_limb =
1012 cloudbox_field_mono(
joker, 0, 0, za_index_local, 0,
joker);
1013 for (
Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
1021 cloudbox_field_mono,
1027 propmat_clearsky_agenda,
1040 surface_rtprop_agenda,
1047 for (
Index p_index = 0;
1048 conv_flag && p_index < cloudbox_field_mono.
nvitrines();
1050 for (
Index stokes_index = 0; conv_flag && stokes_index < stokes_dim;
1052 Numeric diff = cloudbox_field_mono(
1053 p_index, 0, 0, za_index_local, 0, stokes_index) -
1054 cloudbox_field_limb(p_index, stokes_index);
1060 if (
abs(diff_bt) > epsilon[stokes_index]) {
1061 out2 <<
"Limb BT difference: " << diff_bt <<
" in stokes dim "
1062 << stokes_index <<
"\n";
1068 out2 <<
"Limb iterations: " << limb_it <<
"\n";
1079 const Tensor6& doit_scat_field,
1082 const Agenda& propmat_clearsky_agenda,
1085 const Agenda& spt_calc_agenda,
1090 const Agenda& ppath_step_agenda,
1092 const Numeric& ppath_lraytrace,
1097 const Vector& refellipsoid,
1101 const Index& f_index,
1102 const Index& doit_za_interp,
1108 <<
" cloudbox_fieldUpdateSeq3D: Radiative transfer calculatiuon in cloudbox.\n";
1109 out2 <<
" ------------------------------------------------------------- \n";
1114 chk_not_empty(
"propmat_clearsky_agenda", propmat_clearsky_agenda);
1119 "The cloudbox dimension is not 3D! \n"
1120 "Do you really want to do a 3D calculation? \n"
1121 "For 1D use *cloudbox_fieldUpdateSeq1D*.\n");
1127 "The range of *za_grid* must [0 180].");
1133 "The range of *za_grid* must [0 360].");
1140 "z_field", z_field, p_grid.
nelem(), lat_grid.
nelem(), lon_grid.
nelem());
1142 "t_field", t_field, p_grid.
nelem(), lat_grid.
nelem(), lon_grid.
nelem());
1151 "*f_index* is greater than number of elements in the\n"
1152 "frequency grid.\n");
1155 "Interpolation method is not defined. Use \n"
1156 "*doit_za_interpSet*.\n");
1158 const Index stokes_dim = doit_scat_field.
ncols();
1163 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1164 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1165 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1171 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1172 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1173 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1185 out3 <<
"Calculate optical properties of individual scattering elements\n";
1195 const Index p_low = cloudbox_limits[0];
1196 const Index p_up = cloudbox_limits[1];
1197 const Index lat_low = cloudbox_limits[2];
1198 const Index lat_up = cloudbox_limits[3];
1199 const Index lon_low = cloudbox_limits[4];
1200 const Index lon_up = cloudbox_limits[5];
1204 Tensor5 ext_mat_field(p_up - p_low + 1,
1205 lat_up - lat_low + 1,
1206 lon_up - lon_low + 1,
1210 Tensor4 abs_vec_field(p_up - p_low + 1,
1211 lat_up - lat_low + 1,
1212 lon_up - lon_low + 1,
1217 for (
Index za_index = 0; za_index < N_scat_za; za_index++) {
1220 for (
Index aa_index = 1; aa_index < N_scat_aa; aa_index++) {
1239 Vector stokes_vec(stokes_dim, 0.);
1241 Numeric theta_lim = 180. - asin((refellipsoid[0] + z_field(p_low, 0, 0)) /
1242 (refellipsoid[0] + z_field(p_up, 0, 0))) *
1246 if (za_grid[za_index] <= 90.) {
1252 for (
Index p_index = p_up - 1; p_index >= p_low; p_index--) {
1253 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1254 for (
Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1256 cloudbox_field_mono,
1266 propmat_clearsky_agenda,
1287 else if (za_grid[za_index] > theta_lim) {
1291 for (
Index p_index = p_low + 1; p_index <= p_up; p_index++) {
1292 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1293 for (
Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1295 cloudbox_field_mono,
1305 propmat_clearsky_agenda,
1334 else if (za_grid[za_index] > 90. && za_grid[za_index] < theta_lim) {
1335 for (
Index p_index = p_low; p_index <= p_up; p_index++) {
1340 if (!(p_index == 0 && za_grid[za_index] > 90.)) {
1341 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1342 for (
Index lon_index = lon_low; lon_index <= lon_up;
1345 cloudbox_field_mono,
1355 propmat_clearsky_agenda,
1392 const Tensor6& doit_scat_field,
1395 const Agenda& propmat_clearsky_agenda,
1398 const Agenda& spt_calc_agenda,
1407 const Index& f_index,
1413 <<
" cloudbox_fieldUpdateSeq1DPP: Radiative transfer calculation in cloudbox.\n";
1415 <<
" --------------------------------------------------------------------- \n";
1417 const Index stokes_dim = doit_scat_field.
ncols();
1423 "The dimension of stokes vector must be"
1427 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1435 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1453 out3 <<
"Calculate optical properties of individual scattering elements\n";
1464 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1471 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
1474 for (za_index = 0; za_index < N_scat_za; za_index++) {
1493 Vector stokes_vec(stokes_dim, 0.);
1495 if (za_grid[za_index] <= 90) {
1504 for (
Index p_index = cloudbox_limits[1] - 1;
1505 p_index >= cloudbox_limits[0];
1508 cloudbox_field_mono,
1514 propmat_clearsky_agenda,
1525 }
else if (za_grid[za_index] > 90) {
1529 for (
Index p_index = cloudbox_limits[0] + 1;
1530 p_index <= cloudbox_limits[1];
1533 cloudbox_field_mono,
1539 propmat_clearsky_agenda,
1559 Index& doit_is_initialized,
1561 const Index& stokes_dim,
1562 const Index& atmosphere_dim,
1566 const Index& doit_za_grid_size,
1567 const Index& cloudbox_on,
1572 doit_is_initialized = 0;
1573 out0 <<
" Cloudbox is off, DOIT calculation will be skipped.\n";
1582 "The dimension of stokes vector must be"
1594 "For accurate results, za_grid must have "
1595 "more than 15 elements.");
1596 if (N_scat_za > 100) {
1598 out1 <<
"Warning: za_grid is very large, which means that the \n"
1599 <<
"calculation will be very slow.\n";
1603 "The range of *za_grid* must [0 180].");
1606 "*za_grid* must be increasing.");
1612 "For accurate results, aa_grid must have "
1613 "more than 5 elements.");
1614 if (N_scat_aa > 100) {
1616 out1 <<
"Warning: aa_grid is very large which means that the \n"
1617 <<
"calculation will be very slow.\n";
1621 "The range of *aa_grid* must [0 360].");
1624 "*doit_za_grid_size* must be greater than 15 for accurate results");
1625 if (doit_za_grid_size > 100) {
1627 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
1628 <<
"calculation will be very slow.\n";
1632 "*cloudbox_limits* is a vector which contains the"
1633 "upper and lower limit of the cloud for all "
1634 "atmospheric dimensions. So its dimension must"
1635 "be 2 x *atmosphere_dim*");
1640 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1642 const Index Ns = stokes_dim;
1645 if (atmosphere_dim == 1) {
1646 cloudbox_field.
resize(Nf, Np_cloud, 1, 1, Nza, 1, Ns);
1647 doit_scat_field.
resize(Np_cloud, 1, 1, Nza, 1, Ns);
1648 }
else if (atmosphere_dim == 3) {
1649 const Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1650 const Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
1653 cloudbox_field.
resize(Nf, Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1654 doit_scat_field.
resize(Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1657 "Scattering calculations are not possible for a 2D"
1658 "atmosphere. If you want to do scattering calculations"
1659 "*atmosphere_dim* has to be either 1 or 3");
1662 cloudbox_field = NAN;
1663 doit_scat_field = NAN;
1664 doit_is_initialized = 1;
1670 const Vector& p_grid_orig,
1674 Tensor6 cloudbox_field_mono_opt(
1675 p_grid_orig.
nelem(), 1, 1, cloudbox_field_mono.
npages(), 1, 1);
1681 cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1683 os <<
"There is a problem with the pressure grid interpolation";
1687 gridpos(Z_gp, p_grid_cloudbox, p_grid_orig);
1691 for (
Index i = 0; i < cloudbox_field_mono.
npages(); i++) {
1692 interp(cloudbox_field_mono_opt(
joker, 0, 0, i, 0, 0),
1694 cloudbox_field_mono(
joker, 0, 0, i, 0, 0),
1697 cloudbox_field_mono = cloudbox_field_mono_opt;
1715 const Index& f_index,
1716 const Agenda& propmat_clearsky_agenda,
1719 const Index& cloudbox_size_max,
1725 " This method currently only works for unpolarized radiation "
1726 "( stokes_dim = 1)");
1730 pnd_field.
nbooks() != cloudbox_limits[1] - cloudbox_limits[0] + 1,
1731 " ScatSpeciesMerge has to be applied in order to use this method");
1733 bool was_too_much =
false;
1734 Numeric tau_scat_max_internal = tau_scat_max;
1736 p_grid_orig = p_grid;
1737 vector<Numeric> z_grid_new;
1738 Vector ext_mat(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1739 Vector abs_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1740 Vector scat_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1742 z_grid_new.reserve(1000);
1743 for (
Index i = 0; i < cloudbox_limits[0]; i++)
1744 z_grid_new.push_back(z_field(i, 0, 0));
1750 const Vector rtp_mag_dummy(3, 0);
1751 const Vector ppath_los_dummy;
1758 Index scat_data_insert_offset = 0;
1759 Vector single_scattering_albedo(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1761 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1] + 1; ++k) {
1762 Index cloudbox_index = k - cloudbox_limits[0];
1765 ext_mat[cloudbox_index] =
1766 scat_data_local[cloudbox_index].ext_mat_data(0, 0, 0, 0, 0);
1767 abs_vec[cloudbox_index] =
1768 scat_data_local[cloudbox_index].abs_vec_data(0, 0, 0, 0, 0);
1769 scat_vec[cloudbox_index] =
1770 ext_mat[cloudbox_index] - abs_vec[cloudbox_index];
1773 cur_propmat_clearsky,
1779 f_grid[
Range(f_index, 1)],
1785 vmr_field(
joker, k, 0, 0),
1786 propmat_clearsky_agenda);
1787 abs_coeff += cur_propmat_clearsky.
Kjj()[0];
1789 single_scattering_albedo[cloudbox_index] =
1790 scat_vec[cloudbox_index] / (ext_mat[cloudbox_index] + abs_coeff);
1796 scat_data_insert_offset = 0;
1797 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1798 Index cloudbox_index = k - cloudbox_limits[0];
1799 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1800 single_scattering_albedo[cloudbox_index + 1]) /
1803 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1804 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1805 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1806 Index factor = (
Index)ceil(scat_opt_thk / tau_scat_max_internal);
1807 for (
Index j = 1; j < factor; j++) {
1808 scat_data_insert_offset++;
1813 if (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] + 1 >
1814 cloudbox_size_max) {
1815 tau_scat_max_internal += 0.01;
1816 was_too_much =
true;
1818 }
while (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] +
1821 scat_data_insert_offset = 0;
1825 <<
"Warning: Pressure grid optimization with the thresholds tau_scat_max = "
1826 << tau_scat_max <<
" and sgl_slb_max = " << sgl_alb_max
1827 <<
" has lead to an enhancement larger than the value of"
1828 <<
" cloudbox_size_max = " << cloudbox_size_max
1829 <<
". This is why I changed tau_scat_max to " << tau_scat_max_internal
1830 <<
". This may lead to errors of a too coarse grid! \n";
1836 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1837 Index cloudbox_index = k - cloudbox_limits[0];
1838 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1839 single_scattering_albedo[cloudbox_index + 1]) /
1842 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1843 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1844 z_grid_new.push_back(z_field(k, 0, 0));
1846 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1847 Index factor = (
Index)ceil(scat_opt_thk / tau_scat_max_internal);
1849 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) / (
Numeric)factor;
1851 scat_data_local[cloudbox_index + scat_data_insert_offset + 1];
1853 scat_data_local[cloudbox_index + scat_data_insert_offset];
1855 for (
Index j = 1; j < factor; j++) {
1856 z_grid_new.push_back(z_field(k, 0, 0) + (
Numeric)j * step);
1864 weightednextLayerPhamat *= weight;
1865 weightednextLayerExtmat *= weight;
1866 weightednextLayerAbsvec *= weight;
1877 scat_data_local.insert(scat_data_local.begin() + cloudbox_index +
1878 scat_data_insert_offset + 1,
1879 std::move(newLayer));
1881 scat_data_insert_offset++;
1886 cloudbox_limits_opt[0] = cloudbox_limits[0];
1887 cloudbox_limits_opt[1] = scat_data_insert_offset + cloudbox_limits[1];
1888 const Index cloudbox_opt_size =
1889 cloudbox_limits_opt[1] - cloudbox_limits_opt[0] + 1;
1890 out3 <<
"Frequency: " << f_grid[f_index]
1891 <<
" old: " << cloudbox_limits[1] - cloudbox_limits[0] + 1
1892 <<
" new: " << cloudbox_opt_size <<
" Factor: "
1893 << (
Numeric)cloudbox_opt_size /
1894 (
Numeric)(cloudbox_limits[1] - cloudbox_limits[0] + 1)
1897 for (
Index i = cloudbox_limits[1]; i < z_field.
npages(); i++)
1898 z_grid_new.push_back(z_field(i, 0, 0));
1900 Vector z_grid(z_grid_new.size());
1901 for (
Index i = 0; i < z_grid.
nelem(); i++) z_grid[i] = z_grid_new[i];
1902 p_grid_orig = p_grid[
Range(cloudbox_limits[0],
1903 cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1910 Tensor6 cloudbox_field_mono_opt(
1911 cloudbox_opt_size, 1, 1, cloudbox_field_mono.
npages(), 1, 1);
1912 Tensor7 pha_mat_doit_opt(cloudbox_opt_size,
1922 os <<
"At the current frequency " << f_grid[f_index]
1923 <<
" there was an error while interpolating the fields to the new z_field";
1934 for (
Index k = cloudbox_limits_opt[0]; k < cloudbox_limits_opt[1]; k++) {
1935 Index i = k - cloudbox_limits[0];
1936 scat_data_local[i].T_grid = t_field_new(i, 0, 0);
1940 interp(p_grid_opt, itw_z, p_grid, Z_gp);
1946 vmr_field_opt(i,
joker, 0, 0), itw_z, vmr_field(i,
joker, 0, 0), Z_gp);
1952 Range(cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1);
1953 Range r2 =
Range(cloudbox_limits_opt[0], cloudbox_opt_size);
1956 z_field(
Range(cloudbox_limits[0],
1957 cloudbox_limits[1] - cloudbox_limits[0] + 1),
1960 z_grid[
Range(cloudbox_limits_opt[0], cloudbox_opt_size)]);
1963 for (
Index i = 0; i < cloudbox_field_mono.
npages(); i++) {
1964 interp(cloudbox_field_mono_opt(
joker, 0, 0, i, 0, 0),
1966 cloudbox_field_mono(
joker, 0, 0, i, 0, 0),
1970 for (
Index j = 0; j < pha_mat_doit.
nbooks(); j++) {
1971 for (
Index k = 0; k < pha_mat_doit.
npages(); k++) {
1974 pha_mat_doit(
joker, i, 0, j, k, 0, 0),
1981 pnd_field.
resize(cloudbox_opt_size, cloudbox_opt_size, 1, 1);
1983 for (
Index i = 0; i < cloudbox_opt_size; i++) pnd_field(i, i, 0, 0) = 1.;
1986 p_grid = p_grid_opt;
1987 t_field = t_field_new;
1988 cloudbox_limits = cloudbox_limits_opt;
1989 cloudbox_field_mono = cloudbox_field_mono_opt;
1990 pha_mat_doit = pha_mat_doit_opt;
1992 z_field(
joker, 0, 0) = z_grid;
1993 vmr_field = vmr_field_opt;
1998 const Index& doit_iteration_counter,
1999 const Tensor6& cloudbox_field_mono,
2000 const Index& f_index,
2005 if (!frequencies.
nelem() || !iterations.
nelem())
return;
2011 os <<
"doit_iteration_f" << f_index <<
"_i" << doit_iteration_counter;
2014 if (frequencies[0] == -1 && iterations[0] == -1) {
2016 os.str() +
".xml", cloudbox_field_mono,
FILE_TYPE_ASCII, 0, verbosity);
2019 for (
Index j = 0; j < frequencies.
nelem(); j++) {
2020 if (f_index == frequencies[j] || (!j && frequencies[j] == -1)) {
2022 if (iterations[0] == -1) {
2024 cloudbox_field_mono,
2032 for (
Index i = 0; i < iterations.
nelem(); i++) {
2033 if (doit_iteration_counter == iterations[i])
2035 cloudbox_field_mono,
2050 const Agenda& pha_mat_spt_agenda,
2051 const Tensor6& cloudbox_field_mono,
2054 const Index& atmosphere_dim,
2058 const Index& doit_za_grid_size,
2075 "The range of *za_grid* must [0 180].");
2081 "The range of *aa_grid* must [0 360].");
2084 const Index stokes_dim = doit_scat_field.
ncols();
2095 if (atmosphere_dim == 1) {
2097 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2104 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2110 }
else if (atmosphere_dim == 3) {
2112 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2113 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2114 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2119 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2120 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2121 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2127 "The atmospheric dimension must be 1D or 3D \n"
2128 "for scattering calculations using the DOIT\n"
2129 "module, but it is not. The value of *atmosphere_dim*\n"
2130 "is ", atmosphere_dim,
".")
2134 "*cloudbox_limits* is a vector which contains the"
2135 "upper and lower limit of the cloud for all "
2136 "atmospheric dimensions. So its dimension must"
2137 "be 2 x *atmosphere_dim*");
2142 "The zenith angle grids for the computation of\n"
2143 "the scattering integral and the RT part must \n"
2144 "be equal. Check definitions in \n"
2145 "*DOAngularGridsSet*. The keyword \n"
2146 "'za_grid_opt_file' should be empty. \n");
2152 doit_za_grid_size, aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
2163 grid_stepsize[0] = 180. / (
Numeric)(doit_za_grid_size - 1);
2165 grid_stepsize[1] = 360. / (
Numeric)(Naa - 1);
2168 Tensor3 product_field(Nza, Naa, stokes_dim, 0);
2170 out2 <<
" Calculate the scattered field\n";
2172 if (atmosphere_dim == 1) {
2175 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2178 for (
Index za_index_local = 0; za_index_local < Nza; za_index_local++) {
2180 out3 <<
"Multiplication of phase matrix with incoming"
2181 <<
" intensities \n";
2187 for (
Index za_in = 0; za_in < Nza; ++za_in) {
2188 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2192 for (
Index i = 0; i < stokes_dim; i++) {
2193 for (
Index j = 0; j < stokes_dim; j++) {
2194 product_field(za_in, aa_in, i) +=
2196 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2197 cloudbox_field_mono(p_index, 0, 0, za_in, 0, j);
2206 for (
Index i = 0; i < stokes_dim; i++) {
2207 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2212 for (
Index i = 0; i < stokes_dim; i++) {
2213 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2225 else if (atmosphere_dim == 3) {
2230 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2232 for (
Index lat_index = 0;
2233 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2235 for (
Index lon_index = 0;
2236 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2238 Numeric rtp_temperature_local =
2239 t_field(p_index + cloudbox_limits[0],
2240 lat_index + cloudbox_limits[2],
2241 lon_index + cloudbox_limits[4]);
2243 for (
Index aa_index_local = 1; aa_index_local < Naa;
2245 for (
Index za_index_local = 0; za_index_local < Nza;
2247 out3 <<
"Calculate phase matrix \n";
2255 rtp_temperature_local,
2256 pha_mat_spt_agenda);
2271 for (
Index za_in = 0; za_in < Nza; ++za_in) {
2272 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2275 for (
Index i = 0; i < stokes_dim; i++) {
2276 for (
Index j = 0; j < stokes_dim; j++) {
2277 product_field(za_in, aa_in, i) +=
2278 pha_mat_local(za_in, aa_in, i, j) *
2279 cloudbox_field_mono(p_index,
2293 for (
Index i = 0; i < stokes_dim; i++) {
2294 doit_scat_field(p_index,
2321 const Agenda& pha_mat_spt_agenda,
2322 const Tensor6& cloudbox_field_mono,
2325 const Index& atmosphere_dim,
2329 const Index& doit_za_grid_size,
2330 const Index& doit_za_interp,
2345 "The range of *za_grid* must [0 180].");
2351 "The range of *aa_grid* must [0 360].");
2354 const Index stokes_dim = doit_scat_field.
ncols();
2365 if (atmosphere_dim == 1) {
2367 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2374 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2380 }
else if (atmosphere_dim == 3) {
2382 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2383 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2384 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2389 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2390 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2391 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2397 "The atmospheric dimension must be 1D or 3D \n"
2398 "for scattering calculations using the DOIT\n"
2399 "module, but it is not. The value of *atmosphere_dim*\n"
2400 "is ", atmosphere_dim,
".")
2404 "Interpolation method is not defined. Use \n"
2405 "*doit_za_interpSet*.\n");
2408 "*cloudbox_limits* is a vector which contains the"
2409 "upper and lower limit of the cloud for all "
2410 "atmospheric dimensions. So its dimension must"
2411 "be 2 x *atmosphere_dim*");
2414 "*doit_za_grid_size* must be greater than 15 for"
2415 "accurate results");
2416 if (doit_za_grid_size > 100) {
2418 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
2419 <<
"calculation will be very slow. The recommended value is 19.\n";
2426 doit_za_grid_size, aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
2437 nlinspace(za_g, 0, 180, doit_za_grid_size);
2442 gridpos(gp_za_i, za_grid, za_g);
2444 Matrix itw_za_i(doit_za_grid_size, 2);
2448 Matrix cloudbox_field_int(doit_za_grid_size, stokes_dim, 0);
2453 gridpos(gp_za, za_g, za_grid);
2459 Matrix doit_scat_field_org(doit_za_grid_size, stokes_dim, 0);
2464 grid_stepsize[0] = 180. / (
Numeric)(doit_za_grid_size - 1);
2466 grid_stepsize[1] = 360. / (
Numeric)(Naa - 1);
2469 Tensor3 product_field(doit_za_grid_size, Naa, stokes_dim, 0);
2471 if (atmosphere_dim == 1) {
2474 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2477 for (
Index i = 0; i < stokes_dim; i++) {
2478 if (doit_za_interp == 0) {
2481 cloudbox_field_mono(p_index, 0, 0,
joker, 0, i),
2483 }
else if (doit_za_interp == 1) {
2485 for (
Index za = 0; za < za_g.
nelem(); za++) {
2486 cloudbox_field_int(za, i) =
2488 cloudbox_field_mono(p_index, 0, 0,
joker, 0, i),
2499 for (
Index za_index_local = 0; za_index_local < doit_za_grid_size;
2501 out3 <<
"Multiplication of phase matrix with incoming"
2502 <<
" intensities \n";
2508 for (
Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2509 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2513 for (
Index i = 0; i < stokes_dim; i++) {
2514 for (
Index j = 0; j < stokes_dim; j++) {
2515 product_field(za_in, aa_in, i) +=
2517 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2518 cloudbox_field_int(za_in, j);
2525 out3 <<
"Compute integral. \n";
2527 for (
Index i = 0; i < stokes_dim; i++) {
2528 doit_scat_field_org(za_index_local, i) =
2533 for (
Index i = 0; i < stokes_dim; i++) {
2534 doit_scat_field_org(za_index_local, i) =
2546 for (
Index i = 0; i < stokes_dim; i++) {
2547 if (doit_za_interp == 0)
2549 interp(doit_scat_field(p_index, 0, 0,
joker, 0, i),
2551 doit_scat_field_org(
joker, i),
2555 for (
Index za = 0; za < za_grid.
nelem(); za++) {
2556 doit_scat_field(p_index, 0, 0, za, 0, i) =
interp_poly(
2557 za_g, doit_scat_field_org(
joker, i), za_grid[za], gp_za[za]);
2565 else if (atmosphere_dim == 3) {
2567 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2569 for (
Index lat_index = 0;
2570 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2572 for (
Index lon_index = 0;
2573 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2575 Numeric rtp_temperature_local =
2576 t_field(p_index + cloudbox_limits[0],
2577 lat_index + cloudbox_limits[2],
2578 lon_index + cloudbox_limits[4]);
2581 for (
Index aa_index_local = 1; aa_index_local < Naa;
2584 for (
Index i = 0; i < stokes_dim; i++) {
2586 cloudbox_field_int(
joker, i),
2588 cloudbox_field_mono(
2589 p_index, lat_index, lon_index,
joker, aa_index_local, i),
2593 for (
Index za_index_local = 0; za_index_local < doit_za_grid_size;
2595 out3 <<
"Calculate phase matrix \n";
2603 rtp_temperature_local,
2604 pha_mat_spt_agenda);
2619 out3 <<
"Multiplication of phase matrix with"
2620 <<
"incoming intensity \n";
2622 for (
Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2623 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2626 for (
Index i = 0; i < stokes_dim; i++) {
2627 for (
Index j = 0; j < stokes_dim; j++) {
2628 product_field(za_in, aa_in, i) +=
2629 pha_mat_local(za_in, aa_in, i, j) *
2630 cloudbox_field_int(za_in, j);
2636 out3 <<
"Compute the integral \n";
2638 for (
Index i = 0; i < stokes_dim; i++) {
2639 doit_scat_field_org(za_index_local, i) =
2648 for (
Index i = 0; i < stokes_dim; i++) {
2651 p_index, lat_index, lon_index,
joker, aa_index_local, i),
2653 doit_scat_field_org(
joker, i),
2663 out2 <<
" Finished scattered field.\n";
2668 Vector& doit_za_grid_opt,
2670 const Tensor6& cloudbox_field_mono,
2672 const Index& doit_za_interp,
2683 cloudbox_field_mono,
2689 cloudbox_field_mono.
ncols());
2692 "The last dimension of *cloudbox_field* corresponds\n"
2693 "to the Stokes dimension, therefore the number of\n"
2694 "columns in *cloudbox_field* must be a number between\n"
2695 "1 and 4, but it is not!");
2698 "Interpolation method is not defined. Use \n"
2699 "*doit_za_interpSet*.\n");
2701 if (za_grid.
nelem() < 500) {
2702 out1 <<
"Warning: The fine grid (*za_grid*) has less than\n"
2703 <<
"500 grid points which is likely not sufficient for\n"
2704 <<
"grid_optimization\n";
2713 Matrix cloudbox_field_opt_mat;
2714 cloudbox_field_opt_mat = 0.;
2718 cloudbox_field_opt_mat,
2720 cloudbox_field_mono,
2727 const Index& atmosphere_dim,
2734 "Polynomial interpolation is only implemented for\n"
2735 "1D DOIT calculations as \n"
2736 "in 3D there can be numerical problems.\n"
2737 "Please use 'linear' interpolation method.");
2739 if (method ==
"linear")
2741 else if (method ==
"polynomial")
2745 "Possible interpolation methods are 'linear' "
2746 "and 'polynomial'.\n");
2753 const Index& atmfields_checked,
2754 const Index& atmgeom_checked,
2755 const Index& cloudbox_checked,
2756 const Index& scat_data_checked,
2757 const Index& cloudbox_on,
2759 const Agenda& doit_mono_agenda,
2760 const Index& doit_is_initialized,
2768 out0 <<
" Cloudbox is off, DOIT calculation will be skipped.\n";
2777 "The atmospheric fields must be flagged to have "
2778 "passed a consistency check (atmfields_checked=1).");
2780 "The atmospheric geometry must be flagged to have "
2781 "passed a consistency check (atmgeom_checked=1).");
2783 "The cloudbox must be flagged to have "
2784 "passed a consistency check (cloudbox_checked=1).");
2787 if (!cloudbox_on)
return;
2790 "The scattering data must be flagged to have "
2791 "passed a consistency check (scat_data_checked=1).");
2802 "Initialization method *DoitInit* has to be called before "
2812 bool failed =
false;
2816#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) firstprivate(wss)
2817 for (
Index f_index = 0; f_index < nf; f_index++) {
2825 os <<
"Frequency: " << f_grid[f_index] / 1e9 <<
" GHz \n";
2828 Tensor6 cloudbox_field_mono_local =
2831 cloudbox_field_mono_local,
2836 cloudbox_field_mono_local;
2837 }
catch (
const std::exception& e) {
2840 os <<
"Error for f_index = " << f_index <<
" (" << f_grid[f_index]
2843#pragma omp critical(DoitCalc_fail)
2846 fail_msg = os.str();
2859 const Index& atmfields_checked,
2860 const Index& atmgeom_checked,
2861 const Index& cloudbox_checked,
2862 const Index& doit_is_initialized,
2863 const Agenda& iy_main_agenda,
2864 const Index& atmosphere_dim,
2869 const Index& cloudbox_on,
2872 const Index& stokes_dim,
2875 const Index& rigorous,
2880 "The atmospheric fields must be flagged to have "
2881 "passed a consistency check (atmfields_checked=1).");
2883 "The atmospheric geometry must be flagged to have "
2884 "passed a consistency check (atmgeom_checked=1).");
2886 "The cloudbox must be flagged to have "
2887 "passed a consistency check (cloudbox_checked=1).");
2890 if (!cloudbox_on)
return;
2894 "Initialization method *DoitInit* has to be called before "
2895 "*DoitGetIncoming*.")
2898 const String iy_unit =
"1";
2901 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2910 Tensor3 iy_transmittance(0, 0, 0);
2917 "The atmospheric dimensionality must be 1 or 3.");
2919 "*za_grid* must include 0 and 180 degrees as endpoints.");
2922 if (atmosphere_dim == 1) {
2928 for (
Index boundary = 0; boundary <= 1; boundary++) {
2929 const Index boundary_index =
2930 boundary ? cloudbox_field.
nvitrines() - 1 : 0;
2931 pos[0] = z_field(cloudbox_limits[boundary], 0, 0);
2935 los[0] = za_grid[0];
2957 cloudbox_field(
joker, boundary_index, 0, 0, 0, 0,
joker) = iy;
2959 for (
Index za_index = 1; za_index < Nza; za_index++) {
2960 los[0] = za_grid[za_index];
2982 cloudbox_field(
joker, boundary_index, 0, 0, za_index, 0,
joker) = iy;
2985 for (
Index fi = 0; fi < Nf; fi++) {
2988 fi, boundary_index, 0, 0, za_index, 0, 0) >
2990 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
2992 fi, boundary_index, 0, 0, za_index, 0, 0) <
2994 "ERROR: Radiance difference between "
2995 "interpolation points is too large (factor ", maxratio,
2997 "to safely interpolate. This might be due to "
2998 "za_grid being too coarse or the radiance field "
2999 "being a step-like function.\n"
3000 "Happens at boundary ", boundary_index,
3001 " between zenith angles ", za_grid[za_index - 1],
3002 " and ", za_grid[za_index],
"deg\n"
3003 "for frequency #", fi,
", where radiances are ",
3004 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0),
3006 cloudbox_field(fi, boundary_index, 0, 0, za_index, 0, 0),
3019 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3021 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3022 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3028 for (
Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3035 for (
Index boundary = 0; boundary <= 1; boundary++) {
3036 const Index boundary_index =
3037 boundary ? cloudbox_field.
nvitrines() - 1 : 0;
3038 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3039 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3040 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3041 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3042 pos[0] = z_field(cloudbox_limits[boundary],
3043 lat_index + cloudbox_limits[2],
3044 lon_index + cloudbox_limits[4]);
3046 for (
Index za_index = 0; za_index < Nza; za_index++) {
3047 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3048 los[0] = za_grid[za_index];
3049 los[1] = aa_g[aa_index];
3054 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3076 cloudbox_field(
joker,
3090 for (
Index boundary = 0; boundary <= 1; boundary++) {
3091 const Index boundary_index = boundary ? cloudbox_field.
nshelves() - 1 : 0;
3092 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3093 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3094 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3095 pos[1] = lat_grid[cloudbox_limits[boundary + 2]];
3096 pos[0] = z_field(p_index + cloudbox_limits[0],
3097 cloudbox_limits[boundary + 2],
3098 lon_index + cloudbox_limits[4]);
3100 for (
Index za_index = 0; za_index < Nza; za_index++) {
3101 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3102 los[0] = za_grid[za_index];
3103 los[1] = aa_g[aa_index];
3107 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3129 cloudbox_field(
joker,
3143 for (
Index boundary = 0; boundary <= 1; boundary++) {
3144 const Index boundary_index = boundary ? cloudbox_field.
nbooks() - 1 : 0;
3145 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3146 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3147 pos[2] = lon_grid[cloudbox_limits[boundary + 4]];
3148 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3149 pos[0] = z_field(p_index + cloudbox_limits[0],
3150 lat_index + cloudbox_limits[2],
3151 cloudbox_limits[boundary + 4]);
3153 for (
Index za_index = 0; za_index < Nza; za_index++) {
3154 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3155 los[0] = za_grid[za_index];
3156 los[1] = aa_g[aa_index];
3160 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3182 cloudbox_field(
joker,
3201 const Index& atmfields_checked,
3202 const Index& atmgeom_checked,
3203 const Index& cloudbox_checked,
3204 const Index& doit_is_initialized,
3205 const Agenda& iy_main_agenda,
3206 const Index& atmosphere_dim,
3213 const Index& stokes_dim,
3219 "The atmospheric fields must be flagged to have "
3220 "passed a consistency check (atmfields_checked=1).");
3222 "The atmospheric geometry must be flagged to have "
3223 "passed a consistency check (atmgeom_checked=1).");
3225 "The cloudbox must be flagged to have "
3226 "passed a consistency check (cloudbox_checked=1).");
3229 if (!cloudbox_on)
return;
3233 "Initialization method *DoitInit* has to be called before "
3234 "*DoitGetIncoming1DAtm*.")
3237 const String iy_unit =
"1";
3239 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3240 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3241 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3251 Tensor3 iy_transmittance(0, 0, 0);
3257 "The atmospheric dimensionality must be 3.");
3259 "*za_grid* must include 0 and 180 degrees as endpoints.");
3261 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3273 for (
Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3283 pos[1] = lat_grid[cloudbox_limits[2]];
3284 pos[2] = lon_grid[cloudbox_limits[4]];
3285 los[1] = aa_g[aa_index];
3288 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3290 cloudbox_limits[0] + p_index, cloudbox_limits[2], cloudbox_limits[4]);
3292 for (
Index za_index = 0; za_index < Nza; za_index++) {
3293 los[0] = za_grid[za_index];
3315 for (
Index aa = 0; aa < Naa; aa++) {
3318 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3319 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3320 cloudbox_field(
joker, 0, lat, lon, za_index, aa,
joker) = iy;
3325 else if (p_index == Np_cloud - 1)
3326 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3327 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3328 cloudbox_field(
joker,
3339 for (
Index lat = 0; lat < 2; lat++) {
3340 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3341 const Index boundary_index =
3342 lat ? cloudbox_field.
nshelves() - 1 : 0;
3344 joker, p_index, boundary_index, lon, za_index, aa,
joker) = iy;
3349 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3350 for (
Index lon = 0; lon < 2; lon++) {
3351 const Index boundary_index = lon ? cloudbox_field.
nbooks() - 1 : 0;
3353 joker, p_index, lat, boundary_index, za_index, aa,
joker) = iy;
3366 const Index& atmosphere_dim,
3367 const Index& stokes_dim,
3369 const Index& doit_is_initialized,
3370 const Tensor7& cloudbox_field_precalc,
3375 "This method is currently only implemented for 1D atmospheres!\n")
3379 "Initialization method *DoitInit* has to be called before "
3380 "*cloudbox_fieldSetFromPrecalc*")
3385 Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3388 "cloudbox_field_precalc has wrong size in frequency dimension.\n",
3389 nf,
" frequency points are expected, but cloudbox_field_precalc "
3390 "contains ", cloudbox_field_precalc.
nlibraries(),
3391 "frequency points.\n")
3393 "cloudbox_field_precalc has wrong size in pressure level dimension.\n",
3394 np,
" pressure levels expected, but cloudbox_field_precalc "
3395 "contains ", cloudbox_field_precalc.
nvitrines(),
3396 "pressure levels.\n")
3398 "cloudbox_field_precalc has wrong size in polar angle dimension.\n",
3399 nza,
" angles expected, but cloudbox_field_precalc "
3400 "contains ", cloudbox_field_precalc.
npages(),
"angles.\n")
3402 "cloudbox_field_precalc has wrong stokes dimension.\n"
3403 "Dimension ", stokes_dim,
3404 " expected, but cloudbox_field_precalc is dimesnion ",
3405 cloudbox_field_precalc.
ncols(),
".\n")
3413 Index first_upwell = 0;
3416 while (za_grid[first_upwell] < 90.) first_upwell++;
3418 Range downwell(0, first_upwell);
3419 Range upwell(first_upwell, za_grid.
nelem() - first_upwell);
3430 cloudbox_field(
joker, np - 1, 0, 0, upwell, 0,
joker) =
3431 cloudbox_field_precalc(
joker, np - 1, 0, 0, upwell, 0,
joker);
3433 if (cloudbox_limits[0] != 0)
3436 cloudbox_field(
joker, 0, 0, 0, downwell, 0,
joker) =
3437 cloudbox_field_precalc(
joker, 0, 0, 0, downwell, 0,
joker);
3453 const Index& atmosphere_dim,
3454 const Index& cloudbox_on,
3455 const Index& doit_is_initialized,
3456 const Index& all_frequencies,
3461 if (!cloudbox_on)
return;
3465 "Initialization method *DoitInit* has to be called before "
3466 "*cloudbox_fieldSetClearsky*.")
3469 <<
" Interpolate boundary clearsky field to obtain the initial field.\n";
3474 if (atmosphere_dim == 1) {
3477 for (
Index f_index = 0; f_index < nf; f_index++) {
3488 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3491 p_grid[
Range(cloudbox_limits[0],
3493 (cloudbox_limits[1] - cloudbox_limits[0]))],
3494 p_grid[
Range(cloudbox_limits[0],
3495 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3497 Matrix itw((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3500 Tensor6 scat_i_p(2, 1, 1, N_za, 1, N_i);
3504 cloudbox_field(f_index,
3512 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3513 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3514 for (
Index i = 0; i < N_i; ++i) {
3516 f_index,
Range(
joker), 0, 0, za_index, aa_index, i);
3519 scat_i_p(
Range(
joker), 0, 0, za_index, aa_index, i);
3521 interp(target_field, itw, source_field, p_gp);
3526 }
else if (atmosphere_dim == 3) {
3528 "Error in cloudbox_fieldSetClearsky: For 3D "
3529 "all_frequencies option is not implemented \n");
3531 for (
Index f_index = 0; f_index < cloudbox_field.
nvitrines(); f_index++) {
3539 Tensor6 scat_i_p(2, N_lat, N_lon, N_za, N_aa, N_i);
3545 Tensor6 scat_i_lat(N_p, 2, N_lon, N_za, N_aa, N_i);
3551 Tensor6 scat_i_lon(N_p, N_lat, 2, N_za, N_aa, N_i);
3559 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3560 ArrayOfGridPos lat_gp((cloudbox_limits[3] - cloudbox_limits[2]) + 1);
3561 ArrayOfGridPos lon_gp((cloudbox_limits[5] - cloudbox_limits[4]) + 1);
3568 p_grid[
Range(cloudbox_limits[0],
3570 (cloudbox_limits[1] - cloudbox_limits[0]))],
3571 p_grid[
Range(cloudbox_limits[0],
3572 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3574 lat_grid[
Range(cloudbox_limits[2],
3576 (cloudbox_limits[3] - cloudbox_limits[2]))],
3577 lat_grid[
Range(cloudbox_limits[2],
3578 (cloudbox_limits[3] - cloudbox_limits[2]) + 1)]);
3580 lon_grid[
Range(cloudbox_limits[4],
3582 (cloudbox_limits[5] - cloudbox_limits[4]))],
3583 lon_grid[
Range(cloudbox_limits[4],
3584 (cloudbox_limits[5] - cloudbox_limits[4]) + 1)]);
3589 Matrix itw_p((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3590 Matrix itw_lat((cloudbox_limits[3] - cloudbox_limits[2]) + 1, 2);
3591 Matrix itw_lon((cloudbox_limits[5] - cloudbox_limits[4]) + 1, 2);
3598 for (
Index lat_index = 0;
3599 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3601 for (
Index lon_index = 0;
3602 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3604 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3605 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3606 for (
Index i = 0; i < N_i; ++i) {
3607 VectorView target_field = cloudbox_field(f_index,
3616 Range(
joker), lat_index, lon_index, za_index, aa_index, i);
3618 interp(target_field, itw_p, source_field, p_gp);
3625 for (
Index p_index = 0;
3626 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3628 for (
Index lon_index = 0;
3629 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3631 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3632 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3633 for (
Index i = 0; i < N_i; ++i) {
3634 VectorView target_field = cloudbox_field(f_index,
3643 p_index,
Range(
joker), lon_index, za_index, aa_index, i);
3645 interp(target_field, itw_lat, source_field, lat_gp);
3652 for (
Index p_index = 0;
3653 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3655 for (
Index lat_index = 0;
3656 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3658 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3659 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3660 for (
Index i = 0; i < N_i; ++i) {
3661 VectorView target_field = cloudbox_field(f_index,
3670 p_index, lat_index,
Range(
joker), za_index, aa_index, i);
3672 interp(target_field, itw_lon, source_field, lon_gp);
3690 const Index& atmosphere_dim,
3691 const Index& stokes_dim,
3693 const Vector& cloudbox_field_values,
3701 cloudbox_field.
nrows(),
3702 cloudbox_field.
ncols());
3704 for (
Index f_index = 0; f_index < cloudbox_field.
nlibraries(); f_index++) {
3705 cloudbox_field_mono =
3715 cloudbox_field_values,
3719 cloudbox_field_mono;
3731 const Index& atmosphere_dim,
3732 const Index& stokes_dim,
3734 const Matrix& cloudbox_field_values,
3739 "Number of rows in *cloudbox_field_values* has to be equal"
3740 " the frequency dimension of *cloudbox_field*.");
3746 cloudbox_field.
nrows(),
3747 cloudbox_field.
ncols());
3749 for (
Index f_index = 0; f_index < cloudbox_field.
nlibraries(); f_index++) {
3750 cloudbox_field_mono =
3760 cloudbox_field_values(f_index,
joker),
3764 cloudbox_field_mono;
3776 const Index& atmosphere_dim,
3777 const Index& stokes_dim,
3779 const Vector& cloudbox_field_values,
3784 out2 <<
" Set initial field to constant values: " << cloudbox_field_values
3794 "The dimension of stokes vector must be"
3798 "Length of *cloudbox_field_values* has to be equal"
3801 "*cloudbox_limits* is a vector which contains the"
3802 "upper and lower limit of the cloud for all "
3803 "atmospheric dimensions. So its dimension must"
3804 "be 2 x *atmosphere_dim*.");
3806 for (
Index i = 0; i < stokes_dim; i++) {
3808 cloudbox_field_values[i];
Declarations for agendas.
This file contains the definition of Array.
The global header file for ARTS.
Constants of physical expressions as constexpr.
void doit_rte_agendaExecute(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const Agenda &input_agenda)
void doit_scat_field_agendaExecute(Workspace &ws, Tensor6 &doit_scat_field, const Tensor6 &cloudbox_field_mono, const Agenda &input_agenda)
void doit_mono_agendaExecute(Workspace &ws, Tensor6 &cloudbox_field_mono, const Vector &f_grid, const Index f_index, const Agenda &input_agenda)
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, Vector &geo_pos, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const String &iy_unit, const Index cloudbox_on, const Index jacobian_do, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, 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 doit_conv_test_agendaExecute(Workspace &ws, Index &doit_conv_flag, Index &doit_iteration_counter, const Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Agenda &input_agenda)
void pha_mat_spt_agendaExecute(Workspace &ws, Tensor5 &pha_mat_spt, const Index za_index, const Index scat_lat_index, const Index scat_lon_index, const Index scat_p_index, const Index aa_index, const Numeric rtp_temperature, const Agenda &input_agenda)
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
Index nrows() const noexcept
Index npages() const
Returns the number of pages.
Index nbooks() 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.
bool empty() const noexcept
Returns true if variable size is zero.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Stokes vector is as Propagation matrix but only has 4 possible values.
void resize(Index p, Index r, Index c)
Resize function.
void resize(Index b, Index p, Index r, Index c)
Resize function.
void resize(Index v, 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.
void resize(Index n)
Resize function.
#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.
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 cloudbox_field_ngAcceleration(Tensor6 &cloudbox_field_mono, const ArrayOfTensor6 &acceleration_input, const Index &accelerated, const Verbosity &)
Convergence acceleration.
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.
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_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
Linear algebra functions.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Header file for logic.cc.
void DoitCalc(Workspace &ws, Tensor7 &cloudbox_field, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &cloudbox_on, const Vector &f_grid, const Agenda &doit_mono_agenda, const Index &doit_is_initialized, const Verbosity &verbosity)
WORKSPACE METHOD: DoitCalc.
void cloudbox_fieldUpdateSeq3D(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Index &doit_za_interp, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq3D.
void doit_za_grid_optCalc(Vector &doit_za_grid_opt, const Tensor6 &cloudbox_field_mono, const Vector &za_grid, const Index &doit_za_interp, const Numeric &acc, const Verbosity &verbosity)
WORKSPACE METHOD: doit_za_grid_optCalc.
void DoitGetIncoming1DAtm(Workspace &ws, Tensor7 &cloudbox_field, Index &cloudbox_on, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &doit_is_initialized, const Agenda &iy_main_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const EnergyLevelMap &nlte_field, const ArrayOfIndex &cloudbox_limits, const Vector &f_grid, const Index &stokes_dim, const Vector &za_grid, const Vector &aa_grid, const Verbosity &)
WORKSPACE METHOD: DoitGetIncoming1DAtm.
void DoitWriteIterationFields(const Index &doit_iteration_counter, const Tensor6 &cloudbox_field_mono, const Index &f_index, const ArrayOfIndex &iterations, const ArrayOfIndex &frequencies, const Verbosity &verbosity)
WORKSPACE METHOD: DoitWriteIterationFields.
void doit_conv_flagLsq(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &f_grid, const Index &f_index, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagLsq.
void cloudbox_fieldUpdate1D(Workspace &ws, Tensor6 &cloudbox_field_mono, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Agenda &surface_rtprop_agenda, const Index &doit_za_interp, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdate1D.
void doit_conv_flagAbs(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagAbs.
void cloudbox_field_monoOptimizeReverse(Tensor6 &cloudbox_field_mono, const Vector &p_grid_orig, const Vector &p_grid, const ArrayOfIndex &cloudbox_limits, const Verbosity &)
WORKSPACE METHOD: cloudbox_field_monoOptimizeReverse.
void doit_za_interpSet(Index &doit_za_interp, const Index &atmosphere_dim, const String &method, const Verbosity &)
WORKSPACE METHOD: doit_za_interpSet.
constexpr Numeric RAD2DEG
void doit_scat_fieldCalcLimb(Workspace &ws, Tensor6 &doit_scat_field, const Agenda &pha_mat_spt_agenda, const Tensor6 &cloudbox_field_mono, const Tensor4 &pnd_field, const Tensor3 &t_field, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Index &doit_za_interp, const Tensor7 &pha_mat_doit, const Verbosity &verbosity)
WORKSPACE METHOD: doit_scat_fieldCalcLimb.
void cloudbox_fieldSetConstPerFreq(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Matrix &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetConstPerFreq.
void cloudbox_field_monoIterate(Workspace &ws, Tensor6 &cloudbox_field_mono, const Agenda &doit_scat_field_agenda, const Agenda &doit_rte_agenda, const Agenda &doit_conv_test_agenda, const Index &accelerated, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_field_monoIterate.
void cloudbox_fieldSetFromPrecalc(Tensor7 &cloudbox_field, const Vector &za_grid, const Vector &f_grid, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfIndex &cloudbox_limits, const Index &doit_is_initialized, const Tensor7 &cloudbox_field_precalc, const Verbosity &)
WORKSPACE METHOD: cloudbox_fieldSetFromPrecalc.
void DoitGetIncoming(Workspace &ws, Tensor7 &cloudbox_field, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &doit_is_initialized, const Agenda &iy_main_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &f_grid, const Index &stokes_dim, const Vector &za_grid, const Vector &aa_grid, const Index &rigorous, const Numeric &maxratio, const Verbosity &)
WORKSPACE METHOD: DoitGetIncoming.
void DOAngularGridsSet(Index &doit_za_grid_size, Vector &aa_grid, Vector &za_grid, const Index &N_za_grid, const Index &N_aa_grid, const String &za_grid_opt_file, const Verbosity &verbosity)
WORKSPACE METHOD: DOAngularGridsSet.
void cloudbox_field_monoSetConst(Tensor6 &cloudbox_field_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_field_monoSetConst.
void cloudbox_fieldSetConst(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &stokes_dim, const Vector &cloudbox_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetConst.
void cloudbox_fieldSetClearsky(Tensor7 &cloudbox_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &cloudbox_on, const Index &doit_is_initialized, const Index &all_frequencies, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldSetClearsky.
void OptimizeDoitPressureGrid(Workspace &ws, Vector &p_grid, Tensor4 &pnd_field, Tensor3 &t_field, ArrayOfArrayOfSingleScatteringData &scat_data_mono, Tensor3 &z_field, ArrayOfIndex &cloudbox_limits, Tensor6 &cloudbox_field_mono, Tensor7 &pha_mat_doit, Tensor4 &vmr_field, Vector &p_grid_orig, const Vector &f_grid, const Index &f_index, const Agenda &propmat_clearsky_agenda, const Numeric &tau_scat_max, const Numeric &sgl_alb_max, const Index &cloudbox_size_max, const Verbosity &verbosity)
WORKSPACE METHOD: OptimizeDoitPressureGrid.
void doit_conv_flagAbsBT(Index &doit_conv_flag, Index &doit_iteration_counter, Tensor6 &cloudbox_field_mono, const Tensor6 &cloudbox_field_mono_old, const Vector &f_grid, const Index &f_index, const Vector &epsilon, const Index &max_iterations, const Index &throw_nonconv_error, const Verbosity &verbosity)
WORKSPACE METHOD: doit_conv_flagAbsBT.
void cloudbox_fieldUpdateSeq1D(Workspace &ws, Tensor6 &cloudbox_field_mono, Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Vector &aa_grid, const Tensor4 &pnd_field, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Vector &p_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Agenda &surface_rtprop_agenda, const Index &doit_za_interp, const Index &normalize, const Numeric &norm_error_threshold, const Index &norm_debug, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq1D.
void cloudbox_fieldUpdateSeq1DPP(Workspace &ws, Tensor6 &cloudbox_field_mono, Index &za_index, const Tensor6 &doit_scat_field, const ArrayOfIndex &cloudbox_limits, const Agenda &propmat_clearsky_agenda, const Tensor4 &vmr_field, const Agenda &spt_calc_agenda, const Vector &za_grid, const Tensor4 &pnd_field, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Vector &f_grid, const Index &f_index, const Verbosity &verbosity)
WORKSPACE METHOD: cloudbox_fieldUpdateSeq1DPP.
void DoitInit(Tensor6 &doit_scat_field, Tensor7 &cloudbox_field, Index &doit_is_initialized, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &f_grid, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Verbosity &verbosity)
WORKSPACE METHOD: DoitInit.
void doit_scat_fieldCalc(Workspace &ws, Tensor6 &doit_scat_field, const Agenda &pha_mat_spt_agenda, const Tensor6 &cloudbox_field_mono, const Tensor4 &pnd_field, const Tensor3 &t_field, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Vector &za_grid, const Vector &aa_grid, const Index &doit_za_grid_size, const Tensor7 &pha_mat_doit, const Verbosity &verbosity)
WORKSPACE METHOD: doit_scat_fieldCalc.
Template functions for general supergeneric ws methods.
void pha_matCalc(Tensor4 &pha_mat, const Tensor5 &pha_mat_spt, const Tensor4 &pnd_field, const Index &atmosphere_dim, const Index &scat_p_index, const Index &scat_lat_index, const Index &scat_lon_index, const Verbosity &)
WORKSPACE METHOD: pha_matCalc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
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 auto rad2deg(auto x) noexcept
Converts radians to degrees.
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
This file contains declerations of functions of physical character.
Propagation path structure and functions.
Numeric pow(const Rational base, Numeric exp)
Power of.
Numeric sqrt(const Rational r)
Square root.
Declaration of functions in rte.cc.
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Header file for special_interp.cc.
The structure to describe a propagation path and releated quantities.
Auxiliary header stuff related to workspace variable groups.
This file contains basic functions to handle XML data files.
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.