71 Index& doit_za_grid_size,
75 const Index& N_za_grid,
76 const Index& N_aa_grid,
77 const String& za_grid_opt_file,
85 "N_aa_grid must be > 0 (even for 1D / DISORT cases).")
93 "N_za_grid must be >= 0.")
94 doit_za_grid_size = N_za_grid;
96 if (za_grid_opt_file ==
"")
101 "N_za_grid must be >1 or =0 (the latter only allowed for RT4).")
110 Index& doit_conv_flag,
111 Index& doit_iteration_counter,
114 const Tensor6& cloudbox_field_mono_old,
117 const Index& max_iterations,
118 const Index& throw_nonconv_error,
125 "Convergence flag is non-zero, which means that this\n"
126 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
127 "be used only in *doit_conv_test_agenda*\n");
133 const Index N_aa = cloudbox_field_mono.
nrows();
134 const Index stokes_dim = cloudbox_field_mono.
ncols();
138 "You have to specify limiting values for the "
139 "convergence test for each Stokes component "
140 "separately. That means that *epsilon* must "
141 "have *stokes_dim* elements!");
145 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
146 "The fields (Tensor6) *cloudbox_field* and \n"
147 "*cloudbox_field_old* which are compared in the \n"
148 "convergence test do not have the same size.\n");
152 doit_iteration_counter += 1;
153 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
155 if (doit_iteration_counter > max_iterations) {
157 out <<
"Method does not converge (number of iterations \n"
158 <<
"is > " << max_iterations <<
"). Either the cloud "
159 <<
"particle number density \n"
160 <<
"is too large or the numerical setup for the DOIT \n"
161 <<
"calculation is not correct. In case of limb \n"
162 <<
"simulations please make sure that you use an \n"
163 <<
"optimized zenith angle grid. \n"
164 <<
"*cloudbox_field* might be wrong.\n";
165 if (throw_nonconv_error != 0) {
171 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
172 cloudbox_field_mono = NAN;
175 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
180 for (
Index p_index = 0; p_index < N_p; p_index++) {
181 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
182 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
183 for (
Index za_index = 0; za_index < N_za; za_index++) {
184 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
185 for (
Index stokes_index = 0; stokes_index < stokes_dim;
187 Numeric diff = (cloudbox_field_mono(p_index,
193 cloudbox_field_mono_old(p_index,
204 if (
abs(diff) > epsilon[stokes_index]) {
205 out1 <<
"difference: " << diff <<
"\n";
223 Index& doit_conv_flag,
224 Index& doit_iteration_counter,
227 const Tensor6& cloudbox_field_mono_old,
229 const Index& f_index,
232 const Index& max_iterations,
233 const Index& throw_nonconv_error,
241 "Convergence flag is non-zero, which means that this \n"
242 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
243 "be used only in *doit_conv_test_agenda*\n");
249 const Index N_aa = cloudbox_field_mono.
nrows();
250 const Index stokes_dim = cloudbox_field_mono.
ncols();
254 "You have to specify limiting values for the "
255 "convergence test for each Stokes component "
256 "separately. That means that *epsilon* must "
257 "have *stokes_dim* elements!");
261 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
262 "The fields (Tensor6) *cloudbox_field* and \n"
263 "*cloudbox_field_old* which are compared in the \n"
264 "convergence test do not have the same size.\n");
273 "*f_index* is greater than number of elements in the\n"
274 "frequency grid.\n");
278 doit_iteration_counter += 1;
279 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
282 if (doit_iteration_counter > max_iterations) {
284 out <<
"At frequency " << f_grid[f_index] <<
" GHz \n"
285 <<
"method does not converge (number of iterations \n"
286 <<
"is > " << max_iterations <<
"). Either the particle"
287 <<
" number density \n"
288 <<
"is too large or the numerical setup for the DOIT \n"
289 <<
"calculation is not correct. In case of limb \n"
290 <<
"simulations please make sure that you use an \n"
291 <<
"optimized zenith angle grid. \n"
292 <<
"*cloudbox_field* might be wrong.\n";
293 if (throw_nonconv_error != 0) {
299 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
300 cloudbox_field_mono = NAN;
303 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
308 for (
Index p_index = 0; p_index < N_p; p_index++) {
309 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
310 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
311 for (
Index za_index = 0; za_index < N_za; za_index++) {
312 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
313 for (
Index stokes_index = 0; stokes_index < stokes_dim;
315 Numeric diff = cloudbox_field_mono(p_index,
321 cloudbox_field_mono_old(p_index,
334 if (
abs(diff_bt) > epsilon[stokes_index]) {
335 out1 <<
"BT difference: " << diff_bt <<
" in stokes dim "
336 << stokes_index <<
"\n";
358 Index& doit_conv_flag,
359 Index& doit_iteration_counter,
362 const Tensor6& cloudbox_field_mono_old,
364 const Index& f_index,
367 const Index& max_iterations,
368 const Index& throw_nonconv_error,
376 "Convergence flag is non-zero, which means that this \n"
377 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
378 "be used only in *doit_conv_test_agenda*\n");
384 const Index N_aa = cloudbox_field_mono.
nrows();
385 const Index stokes_dim = cloudbox_field_mono.
ncols();
389 "You have to specify limiting values for the "
390 "convergence test for each Stokes component "
391 "separately. That means that *epsilon* must "
392 "have *stokes_dim* elements!");
396 cloudbox_field_mono_old, N_p, N_lat, N_lon, N_za, N_aa, stokes_dim),
397 "The fields (Tensor6) *cloudbox_field* and \n"
398 "*cloudbox_field_old* which are compared in the \n"
399 "convergence test do not have the same size.\n");
408 "*f_index* is greater than number of elements in the\n"
409 "frequency grid.\n");
413 doit_iteration_counter += 1;
414 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
416 if (doit_iteration_counter > max_iterations) {
418 out <<
"Method does not converge (number of iterations \n"
419 <<
"is > " << max_iterations <<
"). Either the"
420 <<
" particle number density \n"
421 <<
"is too large or the numerical setup for the DOIT \n"
422 <<
"calculation is not correct. In case of limb \n"
423 <<
"simulations please make sure that you use an \n"
424 <<
"optimized zenith angle grid. \n";
425 if (throw_nonconv_error != 0) {
431 out1 <<
"Warning in DOIT calculation (output set to NaN):\n" << out.str();
432 cloudbox_field_mono = NAN;
435 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
445 for (
Index p_index = 0; p_index < N_p; p_index++) {
446 for (
Index lat_index = 0; lat_index < N_lat; lat_index++) {
447 for (
Index lon_index = 0; lon_index < N_lon; lon_index++) {
448 for (
Index za_index = 0; za_index < N_za; za_index++) {
449 for (
Index aa_index = 0; aa_index < N_aa; aa_index++) {
452 p_index, lat_index, lon_index, za_index, aa_index, i) -
453 cloudbox_field_mono_old(p_index,
466 lqs[i] =
sqrt(lqs[i]);
467 lqs[i] /= (
Numeric)(N_p * N_lat * N_lon * N_za * N_aa);
472 if (lqs[i] >= epsilon[i]) doit_conv_flag = 0;
475 out1 <<
"lqs [I]: " << lqs[0] <<
"\n";
485 const Agenda& doit_scat_field_agenda,
486 const Agenda& doit_rte_agenda,
487 const Agenda& doit_conv_test_agenda,
488 const Index& accelerated,
495 chk_not_empty(
"doit_scat_field_agenda", doit_scat_field_agenda);
497 chk_not_empty(
"doit_conv_test_agenda", doit_conv_test_agenda);
502 for (
Index p = 0; p < cloudbox_field_mono.
npages(); p++)
503 for (
Index r = 0; r < cloudbox_field_mono.
nrows(); r++)
506 "*cloudbox_field_mono* contains at least one NaN value.\n"
507 "This indicates an improper initialization of *cloudbox_field*.");
514 Tensor6 cloudbox_field_mono_old_local;
515 Index doit_conv_flag_local;
516 Index doit_iteration_counter_local;
522 cloudbox_field_mono.
nbooks(),
523 cloudbox_field_mono.
npages(),
524 cloudbox_field_mono.
nrows(),
525 cloudbox_field_mono.
ncols(),
528 doit_conv_flag_local = 0;
529 doit_iteration_counter_local = 0;
533 acceleration_input.resize(4);
535 while (doit_conv_flag_local == 0) {
537 cloudbox_field_mono_old_local = cloudbox_field_mono;
542 out2 <<
" Execute doit_scat_field_agenda. \n";
544 ws, doit_scat_field_local, cloudbox_field_mono, doit_scat_field_agenda);
547 out2 <<
" Execute doit_rte_agenda. \n";
549 ws, cloudbox_field_mono, doit_scat_field_local, doit_rte_agenda);
553 doit_conv_flag_local,
554 doit_iteration_counter_local,
556 cloudbox_field_mono_old_local,
557 doit_conv_test_agenda);
560 if (accelerated > 0 && doit_conv_flag_local == 0) {
561 acceleration_input[(doit_iteration_counter_local - 1) % 4] =
564 if (doit_iteration_counter_local % 4 == 0) {
566 cloudbox_field_mono, acceleration_input, accelerated, verbosity);
578 const Tensor6& doit_scat_field,
581 const Agenda& propmat_clearsky_agenda,
584 const Agenda& spt_calc_agenda,
588 const Agenda& ppath_step_agenda,
590 const Numeric& ppath_lraytrace,
593 const Vector& refellipsoid,
597 const Index& f_index,
598 const Agenda& surface_rtprop_agenda,
599 const Index& doit_za_interp,
605 <<
" cloudbox_fieldUpdate1D: Radiative transfer calculation in cloudbox\n";
606 out2 <<
" ------------------------------------------------------------- \n";
615 "The cloudbox dimension is not 1D! \n"
616 "Do you really want to do a 1D calculation? \n"
617 "If not, use *cloudbox_fieldUpdateSeq3D*.\n");
623 "The range of *za_grid* must [0 180].");
626 "The length of *p_grid* must be >= 2.");
639 "*f_index* is greater than number of elements in the\n"
640 "frequency grid.\n");
643 "Interpolation method is not defined. Use \n"
644 "*doit_za_interpSet*.\n");
646 const Index stokes_dim = doit_scat_field.
ncols();
651 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
659 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
673 out3 <<
"Calculate optical properties of individual scattering elements\n";
683 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
690 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
692 Tensor6 cloudbox_field_old(cloudbox_field_mono);
695 Index aa_index_local = 0;
698 for (
Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
717 for (
Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
719 if ((p_index != 0) || (za_grid[za_index_local] <= 90.)) {
728 propmat_clearsky_agenda,
741 surface_rtprop_agenda,
758 const Agenda& propmat_clearsky_agenda,
761 const Agenda& spt_calc_agenda,
766 const Agenda& ppath_step_agenda,
768 const Numeric& ppath_lraytrace,
771 const Vector& refellipsoid,
775 const Index& f_index,
776 const Agenda& surface_rtprop_agenda,
777 const Index& doit_za_interp,
778 const Index& normalize,
779 const Numeric& norm_error_threshold,
780 const Index& norm_debug,
786 <<
" cloudbox_fieldUpdateSeq1D: Radiative transfer calculation in cloudbox\n";
787 out2 <<
" ------------------------------------------------------------- \n";
792 chk_not_empty(
"propmat_clearsky_agenda", propmat_clearsky_agenda);
797 "The cloudbox dimension is not 1D! \n"
798 "Do you really want to do a 1D calculation? \n"
799 "For 3D use *cloudbox_fieldUpdateSeq3D*.\n");
805 "The range of *za_grid* must [0 180].");
808 "The length of *p_grid* must be >= 2.");
821 "*f_index* is greater than number of elements in the\n"
822 "frequency grid.\n");
825 "Interpolation method is not defined. Use \n"
826 "*doit_za_interpSet*.\n");
828 const Index stokes_dim = doit_scat_field.
ncols();
833 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
841 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
855 out3 <<
"Calculate optical properties of individual scattering elements\n";
865 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
872 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
877 180. - asin((refellipsoid[0] + z_field(cloudbox_limits[0], 0, 0)) /
878 (refellipsoid[0] + z_field(cloudbox_limits[1], 0, 0))) *
888 Matrix cloudbox_field_limb;
891 Index aa_index_local = 0;
905 norm_error_threshold,
911 for (
Index za_index_local = 0; za_index_local < N_scat_za; za_index_local++) {
931 if (za_grid[za_index_local] <= 90.) {
937 for (
Index p_index = cloudbox_limits[1] - 1;
938 p_index >= cloudbox_limits[0];
947 propmat_clearsky_agenda,
960 surface_rtprop_agenda,
964 }
else if (za_grid[za_index_local] >= theta_lim) {
968 for (
Index p_index = cloudbox_limits[0] + 1;
969 p_index <= cloudbox_limits[1];
978 propmat_clearsky_agenda,
991 surface_rtprop_agenda,
1005 bool conv_flag =
false;
1007 while (!conv_flag && limb_it < 10) {
1009 cloudbox_field_limb =
1010 cloudbox_field_mono(
joker, 0, 0, za_index_local, 0,
joker);
1011 for (
Index p_index = cloudbox_limits[0]; p_index <= cloudbox_limits[1];
1019 cloudbox_field_mono,
1025 propmat_clearsky_agenda,
1038 surface_rtprop_agenda,
1045 for (
Index p_index = 0;
1046 conv_flag && p_index < cloudbox_field_mono.
nvitrines();
1048 for (
Index stokes_index = 0; conv_flag && stokes_index < stokes_dim;
1050 Numeric diff = cloudbox_field_mono(
1051 p_index, 0, 0, za_index_local, 0, stokes_index) -
1052 cloudbox_field_limb(p_index, stokes_index);
1058 if (
abs(diff_bt) > epsilon[stokes_index]) {
1059 out2 <<
"Limb BT difference: " << diff_bt <<
" in stokes dim "
1060 << stokes_index <<
"\n";
1066 out2 <<
"Limb iterations: " << limb_it <<
"\n";
1077 const Tensor6& doit_scat_field,
1080 const Agenda& propmat_clearsky_agenda,
1083 const Agenda& spt_calc_agenda,
1088 const Agenda& ppath_step_agenda,
1090 const Numeric& ppath_lraytrace,
1095 const Vector& refellipsoid,
1099 const Index& f_index,
1100 const Index& doit_za_interp,
1106 <<
" cloudbox_fieldUpdateSeq3D: Radiative transfer calculatiuon in cloudbox.\n";
1107 out2 <<
" ------------------------------------------------------------- \n";
1112 chk_not_empty(
"propmat_clearsky_agenda", propmat_clearsky_agenda);
1117 "The cloudbox dimension is not 3D! \n"
1118 "Do you really want to do a 3D calculation? \n"
1119 "For 1D use *cloudbox_fieldUpdateSeq1D*.\n");
1125 "The range of *za_grid* must [0 180].");
1131 "The range of *za_grid* must [0 360].");
1138 "z_field", z_field, p_grid.
nelem(), lat_grid.
nelem(), lon_grid.
nelem());
1140 "t_field", t_field, p_grid.
nelem(), lat_grid.
nelem(), lon_grid.
nelem());
1149 "*f_index* is greater than number of elements in the\n"
1150 "frequency grid.\n");
1153 "Interpolation method is not defined. Use \n"
1154 "*doit_za_interpSet*.\n");
1156 const Index stokes_dim = doit_scat_field.
ncols();
1161 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1162 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1163 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1169 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1170 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1171 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1183 out3 <<
"Calculate optical properties of individual scattering elements\n";
1193 const Index p_low = cloudbox_limits[0];
1194 const Index p_up = cloudbox_limits[1];
1195 const Index lat_low = cloudbox_limits[2];
1196 const Index lat_up = cloudbox_limits[3];
1197 const Index lon_low = cloudbox_limits[4];
1198 const Index lon_up = cloudbox_limits[5];
1202 Tensor5 ext_mat_field(p_up - p_low + 1,
1203 lat_up - lat_low + 1,
1204 lon_up - lon_low + 1,
1208 Tensor4 abs_vec_field(p_up - p_low + 1,
1209 lat_up - lat_low + 1,
1210 lon_up - lon_low + 1,
1215 for (
Index za_index = 0; za_index < N_scat_za; za_index++) {
1218 for (
Index aa_index = 1; aa_index < N_scat_aa; aa_index++) {
1237 Vector stokes_vec(stokes_dim, 0.);
1239 Numeric theta_lim = 180. - asin((refellipsoid[0] + z_field(p_low, 0, 0)) /
1240 (refellipsoid[0] + z_field(p_up, 0, 0))) *
1244 if (za_grid[za_index] <= 90.) {
1250 for (
Index p_index = p_up - 1; p_index >= p_low; p_index--) {
1251 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1252 for (
Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1254 cloudbox_field_mono,
1264 propmat_clearsky_agenda,
1285 else if (za_grid[za_index] > theta_lim) {
1289 for (
Index p_index = p_low + 1; p_index <= p_up; p_index++) {
1290 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1291 for (
Index lon_index = lon_low; lon_index <= lon_up; lon_index++) {
1293 cloudbox_field_mono,
1303 propmat_clearsky_agenda,
1332 else if (za_grid[za_index] > 90. && za_grid[za_index] < theta_lim) {
1333 for (
Index p_index = p_low; p_index <= p_up; p_index++) {
1338 if (!(p_index == 0 && za_grid[za_index] > 90.)) {
1339 for (
Index lat_index = lat_low; lat_index <= lat_up; lat_index++) {
1340 for (
Index lon_index = lon_low; lon_index <= lon_up;
1343 cloudbox_field_mono,
1353 propmat_clearsky_agenda,
1390 const Tensor6& doit_scat_field,
1393 const Agenda& propmat_clearsky_agenda,
1396 const Agenda& spt_calc_agenda,
1405 const Index& f_index,
1411 <<
" cloudbox_fieldUpdateSeq1DPP: Radiative transfer calculation in cloudbox.\n";
1413 <<
" --------------------------------------------------------------------- \n";
1415 const Index stokes_dim = doit_scat_field.
ncols();
1421 "The dimension of stokes vector must be"
1425 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1433 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1451 out3 <<
"Calculate optical properties of individual scattering elements\n";
1462 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1469 cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1, stokes_dim, 0.);
1472 for (za_index = 0; za_index < N_scat_za; za_index++) {
1491 Vector stokes_vec(stokes_dim, 0.);
1493 if (za_grid[za_index] <= 90) {
1502 for (
Index p_index = cloudbox_limits[1] - 1;
1503 p_index >= cloudbox_limits[0];
1506 cloudbox_field_mono,
1512 propmat_clearsky_agenda,
1523 }
else if (za_grid[za_index] > 90) {
1527 for (
Index p_index = cloudbox_limits[0] + 1;
1528 p_index <= cloudbox_limits[1];
1531 cloudbox_field_mono,
1537 propmat_clearsky_agenda,
1557 Index& doit_is_initialized,
1559 const Index& stokes_dim,
1560 const Index& atmosphere_dim,
1564 const Index& doit_za_grid_size,
1565 const Index& cloudbox_on,
1570 doit_is_initialized = 0;
1571 out0 <<
" Cloudbox is off, DOIT calculation will be skipped.\n";
1580 "The dimension of stokes vector must be"
1592 "For accurate results, za_grid must have "
1593 "more than 15 elements.");
1594 if (N_scat_za > 100) {
1596 out1 <<
"Warning: za_grid is very large, which means that the \n"
1597 <<
"calculation will be very slow.\n";
1601 "The range of *za_grid* must [0 180].");
1604 "*za_grid* must be increasing.");
1610 "For accurate results, aa_grid must have "
1611 "more than 5 elements.");
1612 if (N_scat_aa > 100) {
1614 out1 <<
"Warning: aa_grid is very large which means that the \n"
1615 <<
"calculation will be very slow.\n";
1619 "The range of *aa_grid* must [0 360].");
1622 "*doit_za_grid_size* must be greater than 15 for accurate results");
1623 if (doit_za_grid_size > 100) {
1625 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
1626 <<
"calculation will be very slow.\n";
1630 "*cloudbox_limits* is a vector which contains the"
1631 "upper and lower limit of the cloud for all "
1632 "atmospheric dimensions. So its dimension must"
1633 "be 2 x *atmosphere_dim*");
1638 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1640 const Index Ns = stokes_dim;
1643 if (atmosphere_dim == 1) {
1644 cloudbox_field.
resize(Nf, Np_cloud, 1, 1, Nza, 1, Ns);
1645 doit_scat_field.
resize(Np_cloud, 1, 1, Nza, 1, Ns);
1646 }
else if (atmosphere_dim == 3) {
1647 const Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1648 const Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
1651 cloudbox_field.
resize(Nf, Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1652 doit_scat_field.
resize(Np_cloud, Nlat_cloud, Nlon_cloud, Nza, Naa, Ns);
1655 "Scattering calculations are not possible for a 2D"
1656 "atmosphere. If you want to do scattering calculations"
1657 "*atmosphere_dim* has to be either 1 or 3");
1660 cloudbox_field = NAN;
1661 doit_scat_field = NAN;
1662 doit_is_initialized = 1;
1668 const Vector& p_grid_orig,
1672 Tensor6 cloudbox_field_mono_opt(
1673 p_grid_orig.
nelem(), 1, 1, cloudbox_field_mono.
npages(), 1, 1);
1679 cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1681 os <<
"There is a problem with the pressure grid interpolation";
1685 gridpos(Z_gp, p_grid_cloudbox, p_grid_orig);
1689 for (
Index i = 0; i < cloudbox_field_mono.
npages(); i++) {
1690 interp(cloudbox_field_mono_opt(
joker, 0, 0, i, 0, 0),
1692 cloudbox_field_mono(
joker, 0, 0, i, 0, 0),
1695 cloudbox_field_mono = cloudbox_field_mono_opt;
1713 const Index& f_index,
1714 const Agenda& propmat_clearsky_agenda,
1717 const Index& cloudbox_size_max,
1723 " This method currently only works for unpolarized radiation "
1724 "( stokes_dim = 1)");
1728 pnd_field.
nbooks() != cloudbox_limits[1] - cloudbox_limits[0] + 1,
1729 " ScatSpeciesMerge has to be applied in order to use this method");
1731 bool was_too_much =
false;
1732 Numeric tau_scat_max_internal = tau_scat_max;
1734 p_grid_orig = p_grid;
1735 vector<Numeric> z_grid_new;
1736 Vector ext_mat(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1737 Vector abs_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1738 Vector scat_vec(cloudbox_limits[1] - cloudbox_limits[0] + 1);
1740 z_grid_new.reserve(1000);
1741 for (
Index i = 0; i < cloudbox_limits[0]; i++)
1742 z_grid_new.push_back(z_field(i, 0, 0));
1748 const Vector rtp_mag_dummy(3, 0);
1749 const Vector ppath_los_dummy;
1756 Index scat_data_insert_offset = 0;
1757 Vector single_scattering_albedo(cloudbox_limits[1] - cloudbox_limits[0] + 1,
1759 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1] + 1; ++k) {
1760 Index cloudbox_index = k - cloudbox_limits[0];
1763 ext_mat[cloudbox_index] =
1764 scat_data_local[cloudbox_index].ext_mat_data(0, 0, 0, 0, 0);
1765 abs_vec[cloudbox_index] =
1766 scat_data_local[cloudbox_index].abs_vec_data(0, 0, 0, 0, 0);
1767 scat_vec[cloudbox_index] =
1768 ext_mat[cloudbox_index] - abs_vec[cloudbox_index];
1771 cur_propmat_clearsky,
1777 f_grid[
Range(f_index, 1)],
1783 vmr_field(
joker, k, 0, 0),
1784 propmat_clearsky_agenda);
1785 abs_coeff += cur_propmat_clearsky.
Kjj()[0];
1787 single_scattering_albedo[cloudbox_index] =
1788 scat_vec[cloudbox_index] / (ext_mat[cloudbox_index] + abs_coeff);
1794 scat_data_insert_offset = 0;
1795 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1796 Index cloudbox_index = k - cloudbox_limits[0];
1797 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1798 single_scattering_albedo[cloudbox_index + 1]) /
1801 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1802 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1803 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1804 Index factor = (
Index)ceil(scat_opt_thk / tau_scat_max_internal);
1805 for (
Index j = 1; j < factor; j++) {
1806 scat_data_insert_offset++;
1811 if (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] + 1 >
1812 cloudbox_size_max) {
1813 tau_scat_max_internal += 0.01;
1814 was_too_much =
true;
1816 }
while (scat_data_insert_offset + cloudbox_limits[1] - cloudbox_limits[0] +
1819 scat_data_insert_offset = 0;
1823 <<
"Warning: Pressure grid optimization with the thresholds tau_scat_max = "
1824 << tau_scat_max <<
" and sgl_slb_max = " << sgl_alb_max
1825 <<
" has lead to an enhancement larger than the value of"
1826 <<
" cloudbox_size_max = " << cloudbox_size_max
1827 <<
". This is why I changed tau_scat_max to " << tau_scat_max_internal
1828 <<
". This may lead to errors of a too coarse grid! \n";
1834 for (
Index k = cloudbox_limits[0]; k < cloudbox_limits[1]; ++k) {
1835 Index cloudbox_index = k - cloudbox_limits[0];
1836 const Numeric sgl_alb = (single_scattering_albedo[cloudbox_index] +
1837 single_scattering_albedo[cloudbox_index + 1]) /
1840 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) *
1841 ((scat_vec[cloudbox_index] + scat_vec[cloudbox_index + 1]) / 2);
1842 z_grid_new.push_back(z_field(k, 0, 0));
1844 if (scat_opt_thk > tau_scat_max_internal && sgl_alb > sgl_alb_max) {
1845 Index factor = (
Index)ceil(scat_opt_thk / tau_scat_max_internal);
1847 (z_field(k + 1, 0, 0) - z_field(k, 0, 0)) / (
Numeric)factor;
1849 scat_data_local[cloudbox_index + scat_data_insert_offset + 1];
1851 scat_data_local[cloudbox_index + scat_data_insert_offset];
1853 for (
Index j = 1; j < factor; j++) {
1854 z_grid_new.push_back(z_field(k, 0, 0) + (
Numeric)j * step);
1862 weightednextLayerPhamat *= weight;
1863 weightednextLayerExtmat *= weight;
1864 weightednextLayerAbsvec *= weight;
1875 scat_data_local.insert(scat_data_local.begin() + cloudbox_index +
1876 scat_data_insert_offset + 1,
1877 std::move(newLayer));
1879 scat_data_insert_offset++;
1884 cloudbox_limits_opt[0] = cloudbox_limits[0];
1885 cloudbox_limits_opt[1] = scat_data_insert_offset + cloudbox_limits[1];
1886 const Index cloudbox_opt_size =
1887 cloudbox_limits_opt[1] - cloudbox_limits_opt[0] + 1;
1888 out3 <<
"Frequency: " << f_grid[f_index]
1889 <<
" old: " << cloudbox_limits[1] - cloudbox_limits[0] + 1
1890 <<
" new: " << cloudbox_opt_size <<
" Factor: "
1891 << (
Numeric)cloudbox_opt_size /
1892 (
Numeric)(cloudbox_limits[1] - cloudbox_limits[0] + 1)
1895 for (
Index i = cloudbox_limits[1]; i < z_field.
npages(); i++)
1896 z_grid_new.push_back(z_field(i, 0, 0));
1898 Vector z_grid(z_grid_new.size());
1899 for (
Index i = 0; i < z_grid.
nelem(); i++) z_grid[i] = z_grid_new[i];
1900 p_grid_orig = p_grid[
Range(cloudbox_limits[0],
1901 cloudbox_limits[1] - cloudbox_limits[0] + 1)];
1908 Tensor6 cloudbox_field_mono_opt(
1909 cloudbox_opt_size, 1, 1, cloudbox_field_mono.
npages(), 1, 1);
1910 Tensor7 pha_mat_doit_opt(cloudbox_opt_size,
1920 os <<
"At the current frequency " << f_grid[f_index]
1921 <<
" there was an error while interpolating the fields to the new z_field";
1932 for (
Index k = cloudbox_limits_opt[0]; k < cloudbox_limits_opt[1]; k++) {
1933 Index i = k - cloudbox_limits[0];
1934 scat_data_local[i].T_grid = t_field_new(i, 0, 0);
1938 interp(p_grid_opt, itw_z, p_grid, Z_gp);
1944 vmr_field_opt(i,
joker, 0, 0), itw_z, vmr_field(i,
joker, 0, 0), Z_gp);
1950 Range(cloudbox_limits[0], cloudbox_limits[1] - cloudbox_limits[0] + 1);
1951 Range r2 =
Range(cloudbox_limits_opt[0], cloudbox_opt_size);
1954 z_field(
Range(cloudbox_limits[0],
1955 cloudbox_limits[1] - cloudbox_limits[0] + 1),
1958 z_grid[
Range(cloudbox_limits_opt[0], cloudbox_opt_size)]);
1961 for (
Index i = 0; i < cloudbox_field_mono.
npages(); i++) {
1962 interp(cloudbox_field_mono_opt(
joker, 0, 0, i, 0, 0),
1964 cloudbox_field_mono(
joker, 0, 0, i, 0, 0),
1968 for (
Index j = 0; j < pha_mat_doit.
nbooks(); j++) {
1969 for (
Index k = 0; k < pha_mat_doit.
npages(); k++) {
1972 pha_mat_doit(
joker, i, 0, j, k, 0, 0),
1979 pnd_field.
resize(cloudbox_opt_size, cloudbox_opt_size, 1, 1);
1981 for (
Index i = 0; i < cloudbox_opt_size; i++) pnd_field(i, i, 0, 0) = 1.;
1984 p_grid = p_grid_opt;
1985 t_field = t_field_new;
1986 cloudbox_limits = cloudbox_limits_opt;
1987 cloudbox_field_mono = cloudbox_field_mono_opt;
1988 pha_mat_doit = pha_mat_doit_opt;
1990 z_field(
joker, 0, 0) = z_grid;
1991 vmr_field = vmr_field_opt;
1996 const Index& doit_iteration_counter,
1997 const Tensor6& cloudbox_field_mono,
1998 const Index& f_index,
2003 if (!frequencies.
nelem() || !iterations.
nelem())
return;
2009 os <<
"doit_iteration_f" << f_index <<
"_i" << doit_iteration_counter;
2012 if (frequencies[0] == -1 && iterations[0] == -1) {
2014 os.str() +
".xml", cloudbox_field_mono,
FILE_TYPE_ASCII, 0, verbosity);
2017 for (
Index j = 0; j < frequencies.
nelem(); j++) {
2018 if (f_index == frequencies[j] || (!j && frequencies[j] == -1)) {
2020 if (iterations[0] == -1) {
2022 cloudbox_field_mono,
2030 for (
Index i = 0; i < iterations.
nelem(); i++) {
2031 if (doit_iteration_counter == iterations[i])
2033 cloudbox_field_mono,
2048 const Agenda& pha_mat_spt_agenda,
2049 const Tensor6& cloudbox_field_mono,
2052 const Index& atmosphere_dim,
2056 const Index& doit_za_grid_size,
2073 "The range of *za_grid* must [0 180].");
2079 "The range of *aa_grid* must [0 360].");
2082 const Index stokes_dim = doit_scat_field.
ncols();
2093 if (atmosphere_dim == 1) {
2095 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2102 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2108 }
else if (atmosphere_dim == 3) {
2110 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2111 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2112 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2117 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2118 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2119 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2125 "The atmospheric dimension must be 1D or 3D \n"
2126 "for scattering calculations using the DOIT\n"
2127 "module, but it is not. The value of *atmosphere_dim*\n"
2128 "is ", atmosphere_dim,
".")
2132 "*cloudbox_limits* is a vector which contains the"
2133 "upper and lower limit of the cloud for all "
2134 "atmospheric dimensions. So its dimension must"
2135 "be 2 x *atmosphere_dim*");
2140 "The zenith angle grids for the computation of\n"
2141 "the scattering integral and the RT part must \n"
2142 "be equal. Check definitions in \n"
2143 "*DOAngularGridsSet*. The keyword \n"
2144 "'za_grid_opt_file' should be empty. \n");
2150 doit_za_grid_size, aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
2161 grid_stepsize[0] = 180. / (
Numeric)(doit_za_grid_size - 1);
2163 grid_stepsize[1] = 360. / (
Numeric)(Naa - 1);
2166 Tensor3 product_field(Nza, Naa, stokes_dim, 0);
2168 out2 <<
" Calculate the scattered field\n";
2170 if (atmosphere_dim == 1) {
2173 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2176 for (
Index za_index_local = 0; za_index_local < Nza; za_index_local++) {
2178 out3 <<
"Multiplication of phase matrix with incoming"
2179 <<
" intensities \n";
2185 for (
Index za_in = 0; za_in < Nza; ++za_in) {
2186 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2190 for (
Index i = 0; i < stokes_dim; i++) {
2191 for (
Index j = 0; j < stokes_dim; j++) {
2192 product_field(za_in, aa_in, i) +=
2194 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2195 cloudbox_field_mono(p_index, 0, 0, za_in, 0, j);
2204 for (
Index i = 0; i < stokes_dim; i++) {
2205 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2210 for (
Index i = 0; i < stokes_dim; i++) {
2211 doit_scat_field(p_index, 0, 0, za_index_local, 0, i) =
2223 else if (atmosphere_dim == 3) {
2228 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2230 for (
Index lat_index = 0;
2231 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2233 for (
Index lon_index = 0;
2234 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2236 Numeric rtp_temperature_local =
2237 t_field(p_index + cloudbox_limits[0],
2238 lat_index + cloudbox_limits[2],
2239 lon_index + cloudbox_limits[4]);
2241 for (
Index aa_index_local = 1; aa_index_local < Naa;
2243 for (
Index za_index_local = 0; za_index_local < Nza;
2245 out3 <<
"Calculate phase matrix \n";
2253 rtp_temperature_local,
2254 pha_mat_spt_agenda);
2269 for (
Index za_in = 0; za_in < Nza; ++za_in) {
2270 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2273 for (
Index i = 0; i < stokes_dim; i++) {
2274 for (
Index j = 0; j < stokes_dim; j++) {
2275 product_field(za_in, aa_in, i) +=
2276 pha_mat_local(za_in, aa_in, i, j) *
2277 cloudbox_field_mono(p_index,
2291 for (
Index i = 0; i < stokes_dim; i++) {
2292 doit_scat_field(p_index,
2319 const Agenda& pha_mat_spt_agenda,
2320 const Tensor6& cloudbox_field_mono,
2323 const Index& atmosphere_dim,
2327 const Index& doit_za_grid_size,
2328 const Index& doit_za_interp,
2343 "The range of *za_grid* must [0 180].");
2349 "The range of *aa_grid* must [0 360].");
2352 const Index stokes_dim = doit_scat_field.
ncols();
2363 if (atmosphere_dim == 1) {
2365 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2372 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2378 }
else if (atmosphere_dim == 3) {
2380 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2381 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2382 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2387 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2388 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
2389 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
2395 "The atmospheric dimension must be 1D or 3D \n"
2396 "for scattering calculations using the DOIT\n"
2397 "module, but it is not. The value of *atmosphere_dim*\n"
2398 "is ", atmosphere_dim,
".")
2402 "Interpolation method is not defined. Use \n"
2403 "*doit_za_interpSet*.\n");
2406 "*cloudbox_limits* is a vector which contains the"
2407 "upper and lower limit of the cloud for all "
2408 "atmospheric dimensions. So its dimension must"
2409 "be 2 x *atmosphere_dim*");
2412 "*doit_za_grid_size* must be greater than 15 for"
2413 "accurate results");
2414 if (doit_za_grid_size > 100) {
2416 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
2417 <<
"calculation will be very slow. The recommended value is 19.\n";
2424 doit_za_grid_size, aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
2435 nlinspace(za_g, 0, 180, doit_za_grid_size);
2440 gridpos(gp_za_i, za_grid, za_g);
2442 Matrix itw_za_i(doit_za_grid_size, 2);
2446 Matrix cloudbox_field_int(doit_za_grid_size, stokes_dim, 0);
2451 gridpos(gp_za, za_g, za_grid);
2457 Matrix doit_scat_field_org(doit_za_grid_size, stokes_dim, 0);
2462 grid_stepsize[0] = 180. / (
Numeric)(doit_za_grid_size - 1);
2464 grid_stepsize[1] = 360. / (
Numeric)(Naa - 1);
2467 Tensor3 product_field(doit_za_grid_size, Naa, stokes_dim, 0);
2469 if (atmosphere_dim == 1) {
2472 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2475 for (
Index i = 0; i < stokes_dim; i++) {
2476 if (doit_za_interp == 0) {
2479 cloudbox_field_mono(p_index, 0, 0,
joker, 0, i),
2481 }
else if (doit_za_interp == 1) {
2483 for (
Index za = 0; za < za_g.
nelem(); za++) {
2484 cloudbox_field_int(za, i) =
2486 cloudbox_field_mono(p_index, 0, 0,
joker, 0, i),
2497 for (
Index za_index_local = 0; za_index_local < doit_za_grid_size;
2499 out3 <<
"Multiplication of phase matrix with incoming"
2500 <<
" intensities \n";
2506 for (
Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2507 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2511 for (
Index i = 0; i < stokes_dim; i++) {
2512 for (
Index j = 0; j < stokes_dim; j++) {
2513 product_field(za_in, aa_in, i) +=
2515 p_index, za_index_local, 0, za_in, aa_in, i, j) *
2516 cloudbox_field_int(za_in, j);
2523 out3 <<
"Compute integral. \n";
2525 for (
Index i = 0; i < stokes_dim; i++) {
2526 doit_scat_field_org(za_index_local, i) =
2531 for (
Index i = 0; i < stokes_dim; i++) {
2532 doit_scat_field_org(za_index_local, i) =
2544 for (
Index i = 0; i < stokes_dim; i++) {
2545 if (doit_za_interp == 0)
2547 interp(doit_scat_field(p_index, 0, 0,
joker, 0, i),
2549 doit_scat_field_org(
joker, i),
2553 for (
Index za = 0; za < za_grid.
nelem(); za++) {
2554 doit_scat_field(p_index, 0, 0, za, 0, i) =
interp_poly(
2555 za_g, doit_scat_field_org(
joker, i), za_grid[za], gp_za[za]);
2563 else if (atmosphere_dim == 3) {
2565 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2567 for (
Index lat_index = 0;
2568 lat_index <= cloudbox_limits[3] - cloudbox_limits[2];
2570 for (
Index lon_index = 0;
2571 lon_index <= cloudbox_limits[5] - cloudbox_limits[4];
2573 Numeric rtp_temperature_local =
2574 t_field(p_index + cloudbox_limits[0],
2575 lat_index + cloudbox_limits[2],
2576 lon_index + cloudbox_limits[4]);
2579 for (
Index aa_index_local = 1; aa_index_local < Naa;
2582 for (
Index i = 0; i < stokes_dim; i++) {
2584 cloudbox_field_int(
joker, i),
2586 cloudbox_field_mono(
2587 p_index, lat_index, lon_index,
joker, aa_index_local, i),
2591 for (
Index za_index_local = 0; za_index_local < doit_za_grid_size;
2593 out3 <<
"Calculate phase matrix \n";
2601 rtp_temperature_local,
2602 pha_mat_spt_agenda);
2617 out3 <<
"Multiplication of phase matrix with"
2618 <<
"incoming intensity \n";
2620 for (
Index za_in = 0; za_in < doit_za_grid_size; za_in++) {
2621 for (
Index aa_in = 0; aa_in < Naa; ++aa_in) {
2624 for (
Index i = 0; i < stokes_dim; i++) {
2625 for (
Index j = 0; j < stokes_dim; j++) {
2626 product_field(za_in, aa_in, i) +=
2627 pha_mat_local(za_in, aa_in, i, j) *
2628 cloudbox_field_int(za_in, j);
2634 out3 <<
"Compute the integral \n";
2636 for (
Index i = 0; i < stokes_dim; i++) {
2637 doit_scat_field_org(za_index_local, i) =
2646 for (
Index i = 0; i < stokes_dim; i++) {
2649 p_index, lat_index, lon_index,
joker, aa_index_local, i),
2651 doit_scat_field_org(
joker, i),
2661 out2 <<
" Finished scattered field.\n";
2666 Vector& doit_za_grid_opt,
2668 const Tensor6& cloudbox_field_mono,
2670 const Index& doit_za_interp,
2681 cloudbox_field_mono,
2687 cloudbox_field_mono.
ncols());
2690 "The last dimension of *cloudbox_field* corresponds\n"
2691 "to the Stokes dimension, therefore the number of\n"
2692 "columns in *cloudbox_field* must be a number between\n"
2693 "1 and 4, but it is not!");
2696 "Interpolation method is not defined. Use \n"
2697 "*doit_za_interpSet*.\n");
2699 if (za_grid.
nelem() < 500) {
2700 out1 <<
"Warning: The fine grid (*za_grid*) has less than\n"
2701 <<
"500 grid points which is likely not sufficient for\n"
2702 <<
"grid_optimization\n";
2711 Matrix cloudbox_field_opt_mat;
2712 cloudbox_field_opt_mat = 0.;
2716 cloudbox_field_opt_mat,
2718 cloudbox_field_mono,
2725 const Index& atmosphere_dim,
2732 "Polynomial interpolation is only implemented for\n"
2733 "1D DOIT calculations as \n"
2734 "in 3D there can be numerical problems.\n"
2735 "Please use 'linear' interpolation method.");
2737 if (method ==
"linear")
2739 else if (method ==
"polynomial")
2743 "Possible interpolation methods are 'linear' "
2744 "and 'polynomial'.\n");
2751 const Index& atmfields_checked,
2752 const Index& atmgeom_checked,
2753 const Index& cloudbox_checked,
2754 const Index& scat_data_checked,
2755 const Index& cloudbox_on,
2757 const Agenda& doit_mono_agenda,
2758 const Index& doit_is_initialized,
2766 out0 <<
" Cloudbox is off, DOIT calculation will be skipped.\n";
2775 "The atmospheric fields must be flagged to have "
2776 "passed a consistency check (atmfields_checked=1).");
2778 "The atmospheric geometry must be flagged to have "
2779 "passed a consistency check (atmgeom_checked=1).");
2781 "The cloudbox must be flagged to have "
2782 "passed a consistency check (cloudbox_checked=1).");
2785 if (!cloudbox_on)
return;
2788 "The scattering data must be flagged to have "
2789 "passed a consistency check (scat_data_checked=1).");
2800 "Initialization method *DoitInit* has to be called before "
2810 bool failed =
false;
2814#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) firstprivate(wss)
2815 for (
Index f_index = 0; f_index < nf; f_index++) {
2823 os <<
"Frequency: " << f_grid[f_index] / 1e9 <<
" GHz \n";
2826 Tensor6 cloudbox_field_mono_local =
2829 cloudbox_field_mono_local,
2834 cloudbox_field_mono_local;
2835 }
catch (
const std::exception& e) {
2838 os <<
"Error for f_index = " << f_index <<
" (" << f_grid[f_index]
2841#pragma omp critical(DoitCalc_fail)
2844 fail_msg = os.str();
2857 const Index& atmfields_checked,
2858 const Index& atmgeom_checked,
2859 const Index& cloudbox_checked,
2860 const Index& doit_is_initialized,
2861 const Agenda& iy_main_agenda,
2862 const Index& atmosphere_dim,
2867 const Index& cloudbox_on,
2870 const Index& stokes_dim,
2873 const Index& rigorous,
2878 "The atmospheric fields must be flagged to have "
2879 "passed a consistency check (atmfields_checked=1).");
2881 "The atmospheric geometry must be flagged to have "
2882 "passed a consistency check (atmgeom_checked=1).");
2884 "The cloudbox must be flagged to have "
2885 "passed a consistency check (cloudbox_checked=1).");
2888 if (!cloudbox_on)
return;
2892 "Initialization method *DoitInit* has to be called before "
2893 "*DoitGetIncoming*.")
2896 const String iy_unit =
"1";
2899 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2908 Tensor3 iy_transmittance(0, 0, 0);
2915 "The atmospheric dimensionality must be 1 or 3.");
2917 "*za_grid* must include 0 and 180 degrees as endpoints.");
2920 if (atmosphere_dim == 1) {
2926 for (
Index boundary = 0; boundary <= 1; boundary++) {
2927 const Index boundary_index =
2928 boundary ? cloudbox_field.
nvitrines() - 1 : 0;
2929 pos[0] = z_field(cloudbox_limits[boundary], 0, 0);
2933 los[0] = za_grid[0];
2955 cloudbox_field(
joker, boundary_index, 0, 0, 0, 0,
joker) = iy;
2957 for (
Index za_index = 1; za_index < Nza; za_index++) {
2958 los[0] = za_grid[za_index];
2980 cloudbox_field(
joker, boundary_index, 0, 0, za_index, 0,
joker) = iy;
2983 for (
Index fi = 0; fi < Nf; fi++) {
2986 fi, boundary_index, 0, 0, za_index, 0, 0) >
2988 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0) /
2990 fi, boundary_index, 0, 0, za_index, 0, 0) <
2992 "ERROR: Radiance difference between "
2993 "interpolation points is too large (factor ", maxratio,
2995 "to safely interpolate. This might be due to "
2996 "za_grid being too coarse or the radiance field "
2997 "being a step-like function.\n"
2998 "Happens at boundary ", boundary_index,
2999 " between zenith angles ", za_grid[za_index - 1],
3000 " and ", za_grid[za_index],
"deg\n"
3001 "for frequency #", fi,
", where radiances are ",
3002 cloudbox_field(fi, boundary_index, 0, 0, za_index - 1, 0, 0),
3004 cloudbox_field(fi, boundary_index, 0, 0, za_index, 0, 0),
3017 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3019 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3020 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3026 for (
Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3033 for (
Index boundary = 0; boundary <= 1; boundary++) {
3034 const Index boundary_index =
3035 boundary ? cloudbox_field.
nvitrines() - 1 : 0;
3036 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3037 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3038 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3039 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3040 pos[0] = z_field(cloudbox_limits[boundary],
3041 lat_index + cloudbox_limits[2],
3042 lon_index + cloudbox_limits[4]);
3044 for (
Index za_index = 0; za_index < Nza; za_index++) {
3045 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3046 los[0] = za_grid[za_index];
3047 los[1] = aa_g[aa_index];
3052 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3074 cloudbox_field(
joker,
3088 for (
Index boundary = 0; boundary <= 1; boundary++) {
3089 const Index boundary_index = boundary ? cloudbox_field.
nshelves() - 1 : 0;
3090 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3091 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++) {
3092 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
3093 pos[1] = lat_grid[cloudbox_limits[boundary + 2]];
3094 pos[0] = z_field(p_index + cloudbox_limits[0],
3095 cloudbox_limits[boundary + 2],
3096 lon_index + cloudbox_limits[4]);
3098 for (
Index za_index = 0; za_index < Nza; za_index++) {
3099 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3100 los[0] = za_grid[za_index];
3101 los[1] = aa_g[aa_index];
3105 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3127 cloudbox_field(
joker,
3141 for (
Index boundary = 0; boundary <= 1; boundary++) {
3142 const Index boundary_index = boundary ? cloudbox_field.
nbooks() - 1 : 0;
3143 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3144 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++) {
3145 pos[2] = lon_grid[cloudbox_limits[boundary + 4]];
3146 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3147 pos[0] = z_field(p_index + cloudbox_limits[0],
3148 lat_index + cloudbox_limits[2],
3149 cloudbox_limits[boundary + 4]);
3151 for (
Index za_index = 0; za_index < Nza; za_index++) {
3152 for (
Index aa_index = 0; aa_index < Naa; aa_index++) {
3153 los[0] = za_grid[za_index];
3154 los[1] = aa_g[aa_index];
3158 if ((za_index != 0 && za_index != (Nza - 1)) || aa_index == 0) {
3180 cloudbox_field(
joker,
3199 const Index& atmfields_checked,
3200 const Index& atmgeom_checked,
3201 const Index& cloudbox_checked,
3202 const Index& doit_is_initialized,
3203 const Agenda& iy_main_agenda,
3204 const Index& atmosphere_dim,
3211 const Index& stokes_dim,
3217 "The atmospheric fields must be flagged to have "
3218 "passed a consistency check (atmfields_checked=1).");
3220 "The atmospheric geometry must be flagged to have "
3221 "passed a consistency check (atmgeom_checked=1).");
3223 "The cloudbox must be flagged to have "
3224 "passed a consistency check (cloudbox_checked=1).");
3227 if (!cloudbox_on)
return;
3231 "Initialization method *DoitInit* has to be called before "
3232 "*DoitGetIncoming1DAtm*.")
3235 const String iy_unit =
"1";
3237 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3238 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3239 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3249 Tensor3 iy_transmittance(0, 0, 0);
3255 "The atmospheric dimensionality must be 3.");
3257 "*za_grid* must include 0 and 180 degrees as endpoints.");
3259 "*aa_grid* must include 0 and 360 degrees as endpoints.");
3271 for (
Index i = 0; i < Naa; i++) aa_g[i] = aa_grid[i] - 180;
3281 pos[1] = lat_grid[cloudbox_limits[2]];
3282 pos[2] = lon_grid[cloudbox_limits[4]];
3283 los[1] = aa_g[aa_index];
3286 for (
Index p_index = 0; p_index < Np_cloud; p_index++) {
3288 cloudbox_limits[0] + p_index, cloudbox_limits[2], cloudbox_limits[4]);
3290 for (
Index za_index = 0; za_index < Nza; za_index++) {
3291 los[0] = za_grid[za_index];
3313 for (
Index aa = 0; aa < Naa; aa++) {
3316 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3317 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3318 cloudbox_field(
joker, 0, lat, lon, za_index, aa,
joker) = iy;
3323 else if (p_index == Np_cloud - 1)
3324 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3325 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3326 cloudbox_field(
joker,
3337 for (
Index lat = 0; lat < 2; lat++) {
3338 for (
Index lon = 0; lon < Nlon_cloud; lon++) {
3339 const Index boundary_index =
3340 lat ? cloudbox_field.
nshelves() - 1 : 0;
3342 joker, p_index, boundary_index, lon, za_index, aa,
joker) = iy;
3347 for (
Index lat = 0; lat < Nlat_cloud; lat++) {
3348 for (
Index lon = 0; lon < 2; lon++) {
3349 const Index boundary_index = lon ? cloudbox_field.
nbooks() - 1 : 0;
3351 joker, p_index, lat, boundary_index, za_index, aa,
joker) = iy;
3364 const Index& atmosphere_dim,
3365 const Index& stokes_dim,
3367 const Index& doit_is_initialized,
3368 const Tensor7& cloudbox_field_precalc,
3373 "This method is currently only implemented for 1D atmospheres!\n")
3377 "Initialization method *DoitInit* has to be called before "
3378 "*cloudbox_fieldSetFromPrecalc*")
3383 Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3386 "cloudbox_field_precalc has wrong size in frequency dimension.\n",
3387 nf,
" frequency points are expected, but cloudbox_field_precalc "
3388 "contains ", cloudbox_field_precalc.
nlibraries(),
3389 "frequency points.\n")
3391 "cloudbox_field_precalc has wrong size in pressure level dimension.\n",
3392 np,
" pressure levels expected, but cloudbox_field_precalc "
3393 "contains ", cloudbox_field_precalc.
nvitrines(),
3394 "pressure levels.\n")
3396 "cloudbox_field_precalc has wrong size in polar angle dimension.\n",
3397 nza,
" angles expected, but cloudbox_field_precalc "
3398 "contains ", cloudbox_field_precalc.
npages(),
"angles.\n")
3400 "cloudbox_field_precalc has wrong stokes dimension.\n"
3401 "Dimension ", stokes_dim,
3402 " expected, but cloudbox_field_precalc is dimesnion ",
3403 cloudbox_field_precalc.
ncols(),
".\n")
3411 Index first_upwell = 0;
3414 while (za_grid[first_upwell] < 90.) first_upwell++;
3416 Range downwell(0, first_upwell);
3417 Range upwell(first_upwell, za_grid.
nelem() - first_upwell);
3428 cloudbox_field(
joker, np - 1, 0, 0, upwell, 0,
joker) =
3429 cloudbox_field_precalc(
joker, np - 1, 0, 0, upwell, 0,
joker);
3431 if (cloudbox_limits[0] != 0)
3434 cloudbox_field(
joker, 0, 0, 0, downwell, 0,
joker) =
3435 cloudbox_field_precalc(
joker, 0, 0, 0, downwell, 0,
joker);
3451 const Index& atmosphere_dim,
3452 const Index& cloudbox_on,
3453 const Index& doit_is_initialized,
3454 const Index& all_frequencies,
3459 if (!cloudbox_on)
return;
3463 "Initialization method *DoitInit* has to be called before "
3464 "*cloudbox_fieldSetClearsky*.")
3467 <<
" Interpolate boundary clearsky field to obtain the initial field.\n";
3472 if (atmosphere_dim == 1) {
3475 for (
Index f_index = 0; f_index < nf; f_index++) {
3486 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3489 p_grid[
Range(cloudbox_limits[0],
3491 (cloudbox_limits[1] - cloudbox_limits[0]))],
3492 p_grid[
Range(cloudbox_limits[0],
3493 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3495 Matrix itw((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3498 Tensor6 scat_i_p(2, 1, 1, N_za, 1, N_i);
3502 cloudbox_field(f_index,
3510 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3511 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3512 for (
Index i = 0; i < N_i; ++i) {
3514 f_index,
Range(
joker), 0, 0, za_index, aa_index, i);
3517 scat_i_p(
Range(
joker), 0, 0, za_index, aa_index, i);
3519 interp(target_field, itw, source_field, p_gp);
3524 }
else if (atmosphere_dim == 3) {
3526 "Error in cloudbox_fieldSetClearsky: For 3D "
3527 "all_frequencies option is not implemented \n");
3529 for (
Index f_index = 0; f_index < cloudbox_field.
nvitrines(); f_index++) {
3537 Tensor6 scat_i_p(2, N_lat, N_lon, N_za, N_aa, N_i);
3543 Tensor6 scat_i_lat(N_p, 2, N_lon, N_za, N_aa, N_i);
3549 Tensor6 scat_i_lon(N_p, N_lat, 2, N_za, N_aa, N_i);
3557 ArrayOfGridPos p_gp((cloudbox_limits[1] - cloudbox_limits[0]) + 1);
3558 ArrayOfGridPos lat_gp((cloudbox_limits[3] - cloudbox_limits[2]) + 1);
3559 ArrayOfGridPos lon_gp((cloudbox_limits[5] - cloudbox_limits[4]) + 1);
3566 p_grid[
Range(cloudbox_limits[0],
3568 (cloudbox_limits[1] - cloudbox_limits[0]))],
3569 p_grid[
Range(cloudbox_limits[0],
3570 (cloudbox_limits[1] - cloudbox_limits[0]) + 1)]);
3572 lat_grid[
Range(cloudbox_limits[2],
3574 (cloudbox_limits[3] - cloudbox_limits[2]))],
3575 lat_grid[
Range(cloudbox_limits[2],
3576 (cloudbox_limits[3] - cloudbox_limits[2]) + 1)]);
3578 lon_grid[
Range(cloudbox_limits[4],
3580 (cloudbox_limits[5] - cloudbox_limits[4]))],
3581 lon_grid[
Range(cloudbox_limits[4],
3582 (cloudbox_limits[5] - cloudbox_limits[4]) + 1)]);
3587 Matrix itw_p((cloudbox_limits[1] - cloudbox_limits[0]) + 1, 2);
3588 Matrix itw_lat((cloudbox_limits[3] - cloudbox_limits[2]) + 1, 2);
3589 Matrix itw_lon((cloudbox_limits[5] - cloudbox_limits[4]) + 1, 2);
3596 for (
Index lat_index = 0;
3597 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3599 for (
Index lon_index = 0;
3600 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3602 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3603 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3604 for (
Index i = 0; i < N_i; ++i) {
3605 VectorView target_field = cloudbox_field(f_index,
3614 Range(
joker), lat_index, lon_index, za_index, aa_index, i);
3616 interp(target_field, itw_p, source_field, p_gp);
3623 for (
Index p_index = 0;
3624 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3626 for (
Index lon_index = 0;
3627 lon_index <= (cloudbox_limits[5] - cloudbox_limits[4]);
3629 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3630 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3631 for (
Index i = 0; i < N_i; ++i) {
3632 VectorView target_field = cloudbox_field(f_index,
3641 p_index,
Range(
joker), lon_index, za_index, aa_index, i);
3643 interp(target_field, itw_lat, source_field, lat_gp);
3650 for (
Index p_index = 0;
3651 p_index <= (cloudbox_limits[1] - cloudbox_limits[0]);
3653 for (
Index lat_index = 0;
3654 lat_index <= (cloudbox_limits[3] - cloudbox_limits[2]);
3656 for (
Index za_index = 0; za_index < N_za; ++za_index) {
3657 for (
Index aa_index = 0; aa_index < N_aa; ++aa_index) {
3658 for (
Index i = 0; i < N_i; ++i) {
3659 VectorView target_field = cloudbox_field(f_index,
3668 p_index, lat_index,
Range(
joker), za_index, aa_index, i);
3670 interp(target_field, itw_lon, source_field, lon_gp);
3688 const Index& atmosphere_dim,
3689 const Index& stokes_dim,
3691 const Vector& cloudbox_field_values,
3699 cloudbox_field.
nrows(),
3700 cloudbox_field.
ncols());
3702 for (
Index f_index = 0; f_index < cloudbox_field.
nlibraries(); f_index++) {
3703 cloudbox_field_mono =
3713 cloudbox_field_values,
3717 cloudbox_field_mono;
3729 const Index& atmosphere_dim,
3730 const Index& stokes_dim,
3732 const Matrix& cloudbox_field_values,
3737 "Number of rows in *cloudbox_field_values* has to be equal"
3738 " the frequency dimension of *cloudbox_field*.");
3744 cloudbox_field.
nrows(),
3745 cloudbox_field.
ncols());
3747 for (
Index f_index = 0; f_index < cloudbox_field.
nlibraries(); f_index++) {
3748 cloudbox_field_mono =
3758 cloudbox_field_values(f_index,
joker),
3762 cloudbox_field_mono;
3774 const Index& atmosphere_dim,
3775 const Index& stokes_dim,
3777 const Vector& cloudbox_field_values,
3782 out2 <<
" Set initial field to constant values: " << cloudbox_field_values
3792 "The dimension of stokes vector must be"
3796 "Length of *cloudbox_field_values* has to be equal"
3799 "*cloudbox_limits* is a vector which contains the"
3800 "upper and lower limit of the cloud for all "
3801 "atmospheric dimensions. So its dimension must"
3802 "be 2 x *atmosphere_dim*.");
3804 for (
Index i = 0; i < stokes_dim; i++) {
3806 cloudbox_field_values[i];
Declarations for agendas.
This file contains the definition of Array.
The global header file for ARTS.
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.
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.
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 bool isnan(double d) noexcept
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
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.