47 using std::ostringstream;
48 using std::runtime_error;
70 const Index& nstreams,
72 const Index& add_straight_angles,
73 const Index& pnd_ncols) {
77 "Cloudbox is off, no scattering calculations to be"
83 "The atmospheric fields must be flagged to have"
84 " passed a consistency check (atmfields_checked=1).");
87 "The atmospheric geometry must be flagged to have"
88 " passed a consistency check (atmgeom_checked=1).");
91 "The cloudbox must be flagged to have"
92 " passed a consistency check (cloudbox_checked=1).");
95 "The scat_data must be flagged to have "
96 "passed a consistency check (scat_data_checked=1).");
100 "For running RT4, atmospheric dimensionality"
103 if (stokes_dim < 0 || stokes_dim > 2)
105 "For running RT4, the dimension of stokes vector"
106 " must be 1 or 2.\n");
110 os <<
"RT4 calculations currently only possible with"
111 <<
" lower cloudbox limit\n"
112 <<
"at 0th atmospheric level"
113 <<
" (assumes surface there, ignoring z_surface).\n";
114 throw runtime_error(os.str());
119 "*cloudbox_limits* is a vector which contains the"
120 " upper and lower limit of the cloud for all"
121 " atmospheric dimensions. So its dimension must"
122 " be 2 x *atmosphere_dim*");
126 "No single scattering data present.\n"
127 "See documentation of WSV *scat_data* for options to"
128 " make single scattering data available.\n");
132 "*pnd_field* is not 1D! \n"
133 "RT4 can only be used for 1D!\n");
135 if (quad_type.length() > 1) {
137 os <<
"Input parameter *quad_type* not allowed to contain more than a"
138 <<
" single character.\n"
139 <<
"Yours has " << quad_type.length() <<
".\n";
140 throw runtime_error(os.str());
143 if (quad_type ==
"D" || quad_type ==
"G") {
144 if (add_straight_angles)
148 }
else if (quad_type ==
"L") {
152 os <<
"Unknown quadrature type: " << quad_type
153 <<
".\nOnly D, G, and L allowed.\n";
154 throw runtime_error(os.str());
161 if (nstreams / 2 * 2 != nstreams) {
163 os <<
"RT4 requires an even number of streams, but yours is " << nstreams
165 throw runtime_error(os.str());
167 nhstreams = nstreams / 2;
170 nummu = nhstreams + nhza;
173 bool no_arb_ori =
true;
181 os <<
"RT4 can only handle scattering elements of type " <<
PTYPE_TOTAL_RND
185 <<
"but at least one element of other type (" <<
PTYPE_GENERAL <<
"="
187 throw runtime_error(os.str());
193 void get_quad_angles(
200 const Index& nhstreams,
202 const Index& nummu) {
203 if (quad_type ==
"D") {
204 double_gauss_quadrature_(
206 }
else if (quad_type ==
"G") {
207 gauss_legendre_quadrature_(
216 if (nhza > 0) mu_values[nhstreams] = 1.;
223 for (
Index imu = 0; imu < nummu; imu++) {
231 void get_rt4surf_props(
237 const String& ground_type,
243 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
245 os <<
"Surface temperature is set to " <<
surface_skin_t <<
" K,\n"
246 <<
"which is not considered a meaningful value.\n";
247 throw runtime_error(os.str());
252 if (ground_type ==
"L")
257 "All values in *surface_scalar_reflectivity* must be inside [0,1].");
267 os <<
"For Lambertian surface reflection, the number of elements in\n"
268 <<
"*surface_scalar_reflectivity* needs to match the length of\n"
269 <<
"*f_grid* or be 1."
270 <<
"\n length of *f_grid* : " <<
f_grid.nelem()
271 <<
"\n length of *surface_scalar_reflectivity* : "
273 throw runtime_error(os.str());
276 }
else if (ground_type ==
"S")
283 os <<
"The number of rows and columnss in *surface_reflectivity*\n"
284 <<
"must match each other.";
285 throw runtime_error(os.str());
291 "All r11 values in *surface_reflectivity* must be inside [0,1].");
313 os <<
"For specular surface reflection, the number of elements in\n"
314 <<
"*surface_reflectivity* needs to match the length of\n"
315 <<
"*f_grid* or be 1."
316 <<
"\n length of *f_grid* : " <<
f_grid.nelem()
317 <<
"\n length of *surface_reflectivity* : "
319 throw runtime_error(os.str());
322 }
else if (ground_type ==
"F")
325 Matrix n_real(nf, 1), n_imag(nf, 1);
329 "surface_complex_refr_index",
338 os <<
"Unknown surface type.\n";
339 throw runtime_error(os.str());
360 const String& ground_type,
372 const Index& auto_inc_nstreams,
374 const Index& za_interp_order,
375 const Index& cos_za_interp,
376 const String& pfct_method,
377 const Index& pfct_aa_grid_size,
417 Vector height(num_layers + 1);
418 Vector temperatures(num_layers + 1);
419 for (
Index i = 0; i < height.nelem(); i++) {
420 height[i] = z[num_layers - i];
421 temperatures[i] = t[num_layers - i];
431 Vector scatlayers(num_layers, 0.);
432 Vector gas_extinct(num_layers, 0.);
443 for (
Index clev = 0; clev < pnd.
ncols(); clev++)
444 pnd_per_level[clev] = pnd(
joker, clev).sum();
445 Numeric pndtot = pnd_per_level.sum();
447 for (
Index i = 0; i < cboxlims[1] - cboxlims[0]; i++) {
448 scatlayers[num_layers - 1 - cboxlims[0] - i] = float(i + 1);
457 if (!auto_inc_nstreams) {
458 extinct_matrix_allf.
resize(
462 par_optpropCalc(emis_vector_allf,
472 if (emis_vector_allf.
nshelves() == 1)
485 if (auto_inc_nstreams) {
510 if (vmr.
ncols() > 0) {
514 t[
Range(0, num_layers + 1)],
516 p[
Range(0, num_layers + 1)],
520 Index pfct_failed = 0;
522 if (nummu_new < nummu) {
523 if (!auto_inc_nstreams)
526 if (emis_vector_allf.
nshelves() != 1) {
529 extinct_matrix = extinct_matrix_allf(
533 par_optpropCalc(emis_vector,
540 t[
Range(0, num_layers + 1)],
544 sca_optpropCalc(scatter_matrix,
565 #pragma omp critical(fortran_rt4)
583 height.get_c_array(),
584 temperatures.get_c_array(),
585 gas_extinct.get_c_array(),
587 scatlayers.get_c_array(),
588 extinct_matrix.get_c_array(),
589 emis_vector.get_c_array(),
590 scatter_matrix.get_c_array(),
594 up_rad.get_c_array(),
595 down_rad.get_c_array());
600 if (nummu_new < nummu) nummu_new = nummu + 1;
603 Vector mu_values_new, quad_weights_new, aa_grid_new;
610 while (pfct_failed && (2 * nummu_new) <= auto_inc_nstreams) {
613 nhstreams_new = nummu_new - nhza;
614 mu_values_new.
resize(nummu_new);
616 quad_weights_new.
resize(nummu_new);
617 quad_weights_new = 0.;
618 get_quad_angles(mu_values_new,
628 extinct_matrix_new.
resize(
630 extinct_matrix_new = 0.;
632 emis_vector_new = 0.;
641 par_optpropCalc(emis_vector_new,
648 t[
Range(0, num_layers + 1)],
653 scatter_matrix_new.
resize(
655 scatter_matrix_new = 0.;
674 if (pfct_failed) nummu_new = nummu_new + 1;
678 nummu_new = nummu_new - 1;
680 os <<
"Could not increase nstreams sufficiently (current: "
681 << 2 * nummu_new <<
")\n"
682 <<
"to satisfy scattering matrix norm at f[" <<
f_index
688 os <<
"Try higher maximum number of allowed streams (ie. higher"
689 <<
" auto_inc_nstreams than " << auto_inc_nstreams <<
").";
690 throw runtime_error(os.str());
693 os <<
"Continuing with nstreams=" << 2 * nummu_new
694 <<
". Output for this frequency might be erroneous.";
718 if (ground_type ==
"A")
740 #pragma omp critical(fortran_rt4)
758 height.get_c_array(),
759 temperatures.get_c_array(),
760 gas_extinct.get_c_array(),
762 scatlayers.get_c_array(),
770 up_rad_new.get_c_array(),
771 down_rad_new.get_c_array());
784 for (
Index j = 0; j < nummu; j++) {
788 gp_za, mu_values_new, mu_values[j], za_interp_order, 0.5);
799 for (
Index k = 0; k < num_layers + 1; k++)
801 up_rad(k, j, ist) =
interp(itw, up_rad_new(k,
joker, ist), gp_za);
802 down_rad(k, j, ist) =
831 for (
Index j = 0; j < nummu; j++) {
833 for (
Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
835 up_rad(num_layers - k, j, ist) * rad_l2f;
837 f_index, k + ncboxremoved, 0, 0, nummu - 1 - j, 0, ist) =
838 down_rad(num_layers - k, j, ist) * rad_l2f;
842 for (
Index k = ncboxremoved - 1; k >= 0; k--) {
857 const Index& nummu) {
858 for (
Index j = 0; j < nummu; j++) {
879 assert(gas_extinct.
nelem() == Np - 1);
890 for (
Index i = 0; i < Np - 1; i++) {
892 rtp_temperature_local = 0.5 * (t_profile[i] + t_profile[i + 1]);
895 for (
Index j = 0; j < vmr_profiles.
nrows(); j++)
896 rtp_vmr_local[j] = 0.5 * (vmr_profiles(j, i) + vmr_profiles(j, i + 1));
898 const Vector rtp_mag_dummy(3, 0);
899 const Vector ppath_los_dummy;
907 propmat_clearsky_local,
910 partial_source_dummy,
917 rtp_temperature_local,
918 rtp_nlte_local_dummy,
923 if (propmat_clearsky_local.
nelem()) {
924 gas_extinct[Np - 2 - i] = propmat_clearsky_local[0].Kjj()[0];
925 for (
Index j = 1; j < propmat_clearsky_local.
nelem(); j++) {
926 gas_extinct[Np - 2 - i] += propmat_clearsky_local[j].Kjj()[0];
949 assert(emis_vector.
nbooks() == Np_cloud - 1);
950 assert(extinct_matrix.
nshelves() == Np_cloud - 1);
995 for (
Index ipc = 0; ipc < Np_cloud - 1; ipc++) {
996 for (
Index fi = 0; fi < abs_vec_bulk.
nbooks(); fi++) {
997 for (
Index imu = 0; imu < nummu; imu++) {
1000 extinct_matrix(fi, ipc, 0, imu, ist2, ist1) =
1001 .5 * (ext_mat_bulk(fi, ipc, imu, ist1, ist2) +
1002 ext_mat_bulk(fi, ipc + 1, imu, ist1, ist2));
1003 extinct_matrix(fi, ipc, 1, imu, ist2, ist1) =
1004 .5 * (ext_mat_bulk(fi, ipc, nummu + imu, ist1, ist2) +
1005 ext_mat_bulk(fi, ipc + 1, nummu + imu, ist1, ist2));
1007 emis_vector(fi, ipc, 0, imu, ist1) =
1008 .5 * (abs_vec_bulk(fi, ipc, imu, ist1) +
1009 abs_vec_bulk(fi, ipc + 1, imu, ist1));
1010 emis_vector(fi, ipc, 1, imu, ist1) =
1011 .5 * (abs_vec_bulk(fi, ipc, nummu + imu, ist1) +
1012 abs_vec_bulk(fi, ipc + 1, nummu + imu, ist1));
1019 void sca_optpropCalc(
1031 const String& pfct_method,
1032 const Index& pfct_aa_grid_size,
1033 const Numeric& pfct_threshold,
1034 const Index& auto_inc_nstreams,
1048 scatter_matrix = 0.;
1054 assert(scatter_matrix.
nvitrines() == Np_cloud - 1);
1062 os <<
"Total number of scattering elements in scat_data ("
1065 throw runtime_error(os.str());
1068 if (pfct_aa_grid_size < 2) {
1070 os <<
"Azimuth grid size for scatt matrix extraction"
1071 <<
" (*pfct_aa_grid_size*) must be >1.\n"
1072 <<
"Yours is " << pfct_aa_grid_size <<
".\n";
1073 throw runtime_error(os.str());
1078 Index i_se_flat = 0;
1080 Matrix ext_fixT_spt(N_se, nza_rt, 0.), abs_fixT_spt(N_se, nza_rt, 0.);
1085 1. / float(pfct_aa_grid_size - 1);
1098 if (pfct_method ==
"low")
1100 else if (pfct_method ==
"high")
1107 for (
Index iza = 0; iza < nza_rt; iza++) {
1108 for (
Index sza = 0; sza < nza_rt; sza++) {
1110 for (
Index saa = 0; saa < pfct_aa_grid_size; saa++) {
1126 if (saa == 0 || saa == pfct_aa_grid_size - 1)
1127 pha_mat *= (daa_totrand / 2.);
1132 sca_mat(i_se_flat, iza, sza,
joker,
joker) = pha_mat_int;
1134 ext_fixT_spt(i_se_flat, iza) =
1136 abs_fixT_spt(i_se_flat, iza) =
1145 assert(aa_datagrid[0] == 0.);
1146 assert(aa_datagrid[naa_se - 1] == 180.);
1155 daa[0] = (aa_datagrid[1] - aa_datagrid[0]) / 360.;
1156 for (
Index saa = 1; saa < naa_se - 1; saa++)
1157 daa[saa] = (aa_datagrid[saa + 1] - aa_datagrid[saa - 1]) / 360.;
1159 (aa_datagrid[naa_se - 1] - aa_datagrid[naa_se - 2]) / 360.;
1164 for (
Index iza = 0; iza < nza_se; iza++)
1165 for (
Index sza = 0; sza < nza_se; sza++) {
1166 for (
Index saa = 0; saa < naa_se; saa++) {
1169 pha_mat_int(sza, iza, ist1, ist2) +=
1182 for (
Index iza = 0; iza < nza_rt; iza++) {
1183 for (
Index sza = 0; sza < nza_rt; sza++) {
1191 gridpos(za_inc_gp, za_datagrid, za_inc);
1192 gridpos(za_sca_gp, za_datagrid, za_sca);
1197 pha_mat_lab(ist1, ist2) =
1206 sca_mat(i_se_flat, iza, sza,
joker,
joker) = pha_mat_lab;
1208 ext_fixT_spt(i_se_flat, iza) =
1210 abs_fixT_spt(i_se_flat, iza) =
1215 os <<
"Unsuitable particle type encountered.";
1216 throw runtime_error(os.str());
1222 assert(i_se_flat == N_se);
1230 Index nummu = nza_rt / 2;
1231 for (
Index scat_p_index_local = 0; scat_p_index_local < Np_cloud - 1;
1232 scat_p_index_local++) {
1233 Vector ext_fixT(nza_rt, 0.), abs_fixT(nza_rt, 0.);
1235 for (
Index i_se = 0; i_se < N_se; i_se++) {
1236 Numeric pnd_mean = 0.5 * (pnd_profiles(i_se, scat_p_index_local + 1) +
1237 pnd_profiles(i_se, scat_p_index_local));
1239 for (
Index iza = 0; iza < nummu; iza++) {
1245 for (
Index sza = 0; sza < nummu; sza++) {
1252 scatter_matrix(scat_p_index_local, 0, iza, ist2, sza, ist1) +=
1253 pnd_mean * sca_mat(i_se, iza, sza, ist1, ist2);
1254 scatter_matrix(scat_p_index_local, 1, iza, ist2, sza, ist1) +=
1255 pnd_mean * sca_mat(i_se, nummu + iza, sza, ist1, ist2);
1256 scatter_matrix(scat_p_index_local, 2, iza, ist2, sza, ist1) +=
1257 pnd_mean * sca_mat(i_se, iza, nummu + sza, ist1, ist2);
1258 scatter_matrix(scat_p_index_local, 3, iza, ist2, sza, ist1) +=
1260 sca_mat(i_se, nummu + iza, nummu + sza, ist1, ist2);
1264 ext_fixT[iza] += pnd_mean * ext_fixT_spt(i_se, iza);
1265 abs_fixT[iza] += pnd_mean * abs_fixT_spt(i_se, iza);
1269 for (
Index iza = 0; iza < nummu; iza++) {
1270 for (
Index ih = 0; ih < 2; ih++) {
1271 if (extinct_matrix(scat_p_index_local, ih, iza, 0, 0) > 0.) {
1280 Numeric ext_nom = ext_fixT[iza];
1281 Numeric sca_nom = ext_nom - abs_fixT[iza];
1282 Numeric w0_nom = sca_nom / ext_nom;
1283 assert(w0_nom >= 0.);
1285 for (
Index sza = 0; sza < nummu; sza++) {
1291 (scatter_matrix(scat_p_index_local, ih, iza, 0, sza, 0) +
1292 scatter_matrix(scat_p_index_local, ih + 2, iza, 0, sza, 0));
1300 Numeric w0_act = 2. *
PI * sca_mat_integ / ext_nom;
1301 Numeric pfct_norm = 2. *
PI * sca_mat_integ / sca_nom;
1307 extinct_matrix(scat_p_index_local, ih, iza, 0, 0) -
1308 emis_vector(scat_p_index_local, ih, iza, 0);
1321 if (
abs(w0_act - w0_nom) > pfct_threshold) {
1322 if (pfct_failed >= 0) {
1323 if (auto_inc_nstreams) {
1328 os <<
"Bulk scattering matrix normalization deviates significantly\n"
1329 <<
"from expected value (" << 1e2 *
abs(1. - pfct_norm)
1331 <<
" resulting in albedo deviation of "
1332 <<
abs(w0_act - w0_nom) <<
").\n"
1333 <<
"Something seems wrong with your scattering data "
1334 <<
" (did you run *scat_dataCheck*?)\n"
1335 <<
"or your RT4 setup (try increasing *nstreams* and in case"
1336 <<
" of randomly oriented particles possibly also"
1337 <<
" pfct_aa_grid_size).";
1338 throw runtime_error(os.str());
1341 }
else if (
abs(w0_act - w0_nom) > pfct_threshold * 0.1 ||
1342 abs(1. - pfct_norm) > 1e-2) {
1345 <<
"Warning: The bulk scattering matrix is not well normalized\n"
1346 <<
"Deviating from expected value by "
1347 << 1e2 *
abs(1. - pfct_norm) <<
"% (and "
1348 <<
abs(w0_act - w0_nom) <<
" in terms of scattering albedo).\n";
1359 pfct_norm = 2. *
PI * sca_mat_integ / sca_nom_paropt;
1430 const String B_unit =
"R";
1435 for (
Index rmu = 0; rmu < nummu; rmu++) {
1477 planck_function_(
surface_skin_t, B_unit.c_str(), wave, B_lambda);
1478 Numeric B_ratio = B_lambda / B_freq;
1504 Vector surf_int_grid(nsl + 1);
1507 surf_int_grid[nsl] =
1510 for (
Index imu = 1; imu < nsl; imu++)
1511 surf_int_grid[imu] =
1514 for (
Index imu = 0; imu < nsl; imu++) {
1518 Numeric w = 0.5 * (cos(2. * surf_int_grid[imu]) -
1519 cos(2. * surf_int_grid[imu + 1]));
1528 for (
Index imu = 0; imu < nummu; imu++) {
1548 (quad_weights[imu] * mu_values[imu]);
1549 }
catch (
const runtime_error& e) {
1569 surf_refl_mat(
f_index, rmu, sto2, rmu, sto1) =
1582 os <<
"Something went wrong.\n"
1583 <<
"At reflected stream #" << rmu
1584 <<
", power reflection coefficient for RT4\n"
1585 <<
"became 0, although the one from surface_rtprop_agenda is "
1587 throw runtime_error(os.str());
1590 R_scale = R_arts[
f_index] / R_rt4;
1597 void rt4_test(
Tensor4& out_rad,
1606 Numeric max_delta_tau = 1.0E-6;
1609 String ground_type =
"L";
1618 Vector height, temperatures, gas_extinct;
1624 ReadXML(gas_extinct,
"gas_extinct", datapath +
"abs_gas.xml",
"",
verbosity);
1629 Index num_scatlayers = 3;
1630 Vector scatlayers(num_layers, 0.);
1640 Tensor6 scatter_matrix(num_scatlayers, 4, nummu, nstokes, nummu, nstokes);
1641 for (
Index ii = 0; ii < 4; ii++)
1642 for (
Index ij = 0; ij < nummu; ij++)
1643 for (
Index ik = 0; ik < nstokes; ik++)
1644 for (
Index il = 0; il < nummu; il++)
1645 for (
Index im = 0; im < nstokes; im++)
1646 for (
Index in = 0; in < num_scatlayers; in++)
1647 scatter_matrix(in, ii, ij, ik, il, im) =
1648 sca_data(im, il, ik, ij, ii);
1649 Tensor5 extinct_matrix(num_scatlayers, 2, nummu, nstokes, nstokes);
1650 for (
Index ii = 0; ii < 2; ii++)
1651 for (
Index ij = 0; ij < nummu; ij++)
1652 for (
Index ik = 0; ik < nstokes; ik++)
1653 for (
Index il = 0; il < nstokes; il++)
1654 for (
Index im = 0; im < num_scatlayers; im++)
1655 extinct_matrix(im, ii, ij, ik, il) = ext_data(il, ik, ij, ii);
1656 Tensor4 emis_vector(num_scatlayers, 2, nummu, nstokes);
1657 for (
Index ii = 0; ii < 2; ii++)
1658 for (
Index ij = 0; ij < nummu; ij++)
1659 for (
Index ik = 0; ik < nstokes; ik++)
1660 for (
Index il = 0; il < num_scatlayers; il++)
1661 emis_vector(il, ii, ij, ik) = abs_data(ik, ij, ii);
1664 Tensor4 surf_refl_mat(nummu, nstokes, nummu, nstokes, 0.);
1665 Matrix surf_emis_vec(nummu, nstokes, 0.);
1666 Matrix ground_reflec(nstokes, nstokes, 0.);
1670 Tensor3 up_rad(num_layers + 1, nummu, nstokes, 0.);
1671 Tensor3 down_rad(num_layers + 1, nummu, nstokes, 0.);
1679 ground_type.c_str(),
1682 ground_reflec.get_c_array(),
1692 scatlayers.get_c_array(),
1693 extinct_matrix.get_c_array(),
1694 emis_vector.get_c_array(),
1698 mu_values.get_c_array(),
1699 up_rad.get_c_array(),
1700 down_rad.get_c_array());
1714 out_rad.
resize(num_layers + 1, 2, nummu, nstokes);
1715 for (
Index ii = 0; ii < nummu; ii++)