75 for (
Index i = 0; i < np1; i++) {
76 for (
Index j = 0; j < nr1 ; j++) {
77 for (
Index k = 0; k < nc1; k++)
80 (sca1(i, j) * pfct1(i, j, k) + sca2(i, j) * pfct2(j, k)) /
81 (sca1(i, j) + sca2(i, j));
87 const Index& cloudbox_on,
88 const Index& atmfields_checked,
89 const Index& atmgeom_checked,
90 const Index& cloudbox_checked,
91 const Index& scat_data_checked,
92 const Index& atmosphere_dim,
93 const Index& stokes_dim,
97 const Index& nstreams) {
100 "Cloudbox is off, no scattering calculations to be"
104 if (atmfields_checked != 1)
106 "The atmospheric fields must be flagged to have "
107 "passed a consistency check (atmfields_checked=1).");
108 if (atmgeom_checked != 1)
110 "The atmospheric geometry must be flagged to have "
111 "passed a consistency check (atmgeom_checked=1).");
112 if (cloudbox_checked != 1)
114 "The cloudbox must be flagged to have "
115 "passed a consistency check (cloudbox_checked=1).");
116 if (scat_data_checked != 1)
118 "The scat_data must be flagged to have "
119 "passed a consistency check (scat_data_checked=1).");
121 if (atmosphere_dim != 1)
123 "For running DISORT, atmospheric dimensionality "
126 if (stokes_dim < 0 || stokes_dim > 1)
128 "For running DISORT, the dimension of stokes vector "
131 if (cloudbox_limits.
nelem() != 2 * atmosphere_dim)
133 "*cloudbox_limits* is a vector which contains the"
134 "upper and lower limit of the cloud for all "
135 "atmospheric dimensions. So its dimension must"
136 "be 2 x *atmosphere_dim*");
138 if (cloudbox_limits[0] != 0) {
140 os <<
"DISORT calculations currently only possible with "
141 <<
"lower cloudbox limit\n"
142 <<
"at 0th atmospheric level "
143 <<
"(assumes surface there, ignoring z_surface).\n";
144 throw runtime_error(os.str());
147 if (scat_data.empty())
149 "No single scattering data present.\n"
150 "See documentation of WSV *scat_data* for options to "
151 "make single scattering data available.\n");
158 if (nstreams / 2 * 2 != nstreams) {
160 os <<
"DISORT requires an even number of streams, but yours is " << nstreams
162 throw runtime_error(os.str());
176 os <<
"We require size of za_grid to be >= 20, to ensure a\n"
177 <<
"reasonable interpolation of the calculated cloudbox field.\n"
178 <<
"Note that for DISORT additional computation costs for\n"
179 <<
"larger numbers of angles are negligible.";
180 throw runtime_error(os.str());
183 if (za_grid[0] != 0. || za_grid[nza - 1] != 180.)
184 throw runtime_error(
"The range of *za_grid* must [0 180].");
187 throw runtime_error(
"*za_grid* must be increasing.");
190 while (za_grid[i] <= 90) {
191 if (za_grid[i] == 90)
192 throw runtime_error(
"*za_grid* is not allowed to contain the value 90");
197 bool all_totrand =
true;
198 for (
Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
199 for (
Index i_se = 0; i_se < scat_data[i_ss].
nelem(); i_se++)
200 if (scat_data[i_ss][i_se].ptype !=
PTYPE_TOTAL_RND) all_totrand =
false;
203 os <<
"DISORT can only handle scattering elements of type "
208 throw runtime_error(os.str());
213 const Index& atmfields_checked,
214 const Index& atmgeom_checked,
215 const Index& scat_data_checked,
216 const Index& atmosphere_dim,
217 const Index& stokes_dim,
219 const Index& nstreams) {
220 if (atmfields_checked != 1)
222 "The atmospheric fields must be flagged to have "
223 "passed a consistency check (atmfields_checked=1).");
225 if (atmgeom_checked != 1)
227 "The atmospheric geometry must be flagged to have "
228 "passed a consistency check (atmgeom_checked=1).");
230 if (scat_data_checked != 1)
232 "The scat_data must be flagged to have "
233 "passed a consistency check (scat_data_checked=1).");
235 if (atmosphere_dim != 1)
237 "For running DISORT, atmospheric dimensionality "
240 if (stokes_dim < 0 || stokes_dim > 1)
242 "For running DISORT, the dimension of stokes vector "
245 if (scat_data.empty())
247 "No single scattering data present.\n"
248 "See documentation of WSV *scat_data* for options to "
249 "make single scattering data available.\n");
256 if (nstreams / 2 * 2 != nstreams) {
258 os <<
"DISORT requires an even number of streams, but yours is " << nstreams
260 throw runtime_error(os.str());
264 bool all_totrand =
true;
265 for (
Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
266 for (
Index i_se = 0; i_se < scat_data[i_ss].
nelem(); i_se++)
267 if (scat_data[i_ss][i_se].ptype !=
PTYPE_TOTAL_RND) all_totrand =
false;
270 os <<
"DISORT can only handle scattering elements of type "
275 throw runtime_error(os.str());
286 const Index& stokes_dim) {
288 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
292 cloudbox_field.
resize(Nf, Np_cloud, 1, 1, n_za, n_aa, stokes_dim);
293 cloudbox_field = NAN;
304 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
306 os <<
"Surface temperature has been set or derived as " << btemp <<
" K,\n"
307 <<
"which is not considered a meaningful value.\n"
308 <<
"For surface method 'L', *surface_skin_t* needs to\n"
309 <<
"be set and passed explicitly. Maybe you didn't do this?";
310 throw runtime_error(os.str());
312 btemp = surface_skin_t;
315 if (surface_scalar_reflectivity.
nelem() != f_grid.
nelem() &&
316 surface_scalar_reflectivity.
nelem() != 1) {
318 os <<
"The number of elements in *surface_scalar_reflectivity*\n"
319 <<
"should match length of *f_grid* or be 1."
320 <<
"\n length of *f_grid* : " << f_grid.
nelem()
321 <<
"\n length of *surface_scalar_reflectivity* : "
322 << surface_scalar_reflectivity.
nelem() <<
"\n";
323 throw runtime_error(os.str());
326 if (
min(surface_scalar_reflectivity) < 0 ||
327 max(surface_scalar_reflectivity) > 1) {
329 "All values in *surface_scalar_reflectivity*"
330 " must be inside [0,1].");
333 if (surface_scalar_reflectivity.
nelem() > 1)
334 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++)
335 albedo[f_index] = surface_scalar_reflectivity[f_index];
337 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++)
338 albedo[f_index] = surface_scalar_reflectivity[0];
343 const Agenda& propmat_clearsky_agenda,
358 const Vector rtp_mag_dummy(3, 0);
359 const Vector ppath_los_dummy;
365 for (
Index ip = 0; ip < Np; ip++) {
367 propmat_clearsky_local,
379 vmr_profiles(
joker, ip),
380 propmat_clearsky_agenda);
381 ext_bulk_gas(
joker, ip) += propmat_clearsky_local.
Kjj();
393 const Agenda& gas_scattering_agenda) {
402 Matrix pmom_gas_level( Np, Nl, 0.);
406 for (
Index ip = 0; ip < Np; ip++) {
418 gas_scattering_agenda);
421 sca_coeff_gas_level(
joker, ip) = K_sca_gas_temp.
Kjj(0, 0);
424 N_polys =
min(Nl, sca_fct_temp.
nelem());
425 for (
Index k = 0; k < N_polys; k++) {
426 pmom_gas_level( ip, k) = sca_fct_temp[k];
431 for (
Index ip = 0; ip < Np - 1; ip++) {
432 for (
Index f_index = 0; f_index < Nf; f_index++) {
434 sca_coeff_gas(f_index, Np - 2 - ip) =
436 (sca_coeff_gas_level(f_index, ip) +
437 sca_coeff_gas_level(f_index, ip + 1));
440 for (
Index l_index = 0; l_index < Nl; l_index++) {
441 pfct_gas( Np - 2 - ip, l_index) =
442 0.5 * (pmom_gas_level( ip, l_index) +
443 pmom_gas_level( ip + 1, l_index));
470 Vector T_array = t_profile[
Range(cloudbox_limits[0], Np_cloud)];
471 Matrix dir_array(1, 2, 0.);
512 bool pf = (abs_vec_bulk.
nbooks() != 1);
513 for (
Index ip = 0; ip < Np_cloud; ip++)
514 for (
Index f_index = 0; f_index < nf; f_index++) {
515 if (pf) f_this = f_index;
516 ext_bulk_par(f_index, ip + cloudbox_limits[0]) =
517 ext_mat_bulk(f_this, ip, 0, 0, 0);
518 abs_bulk_par(f_index, ip + cloudbox_limits[0]) =
519 abs_vec_bulk(f_this, ip, 0, 0);
541 for (
Index ip = 0; ip < Np - 1; ip++)
543 for (
Index f_index = 0; f_index < nf; f_index++) {
545 0.5 * (ext_bulk_gas(f_index, ip) + ext_bulk_par(f_index, ip) +
546 ext_bulk_gas(f_index, ip + 1) + ext_bulk_par(f_index, ip + 1));
550 (ext_bulk_gas(f_index, ip) + abs_bulk_par(f_index, ip) +
551 ext_bulk_gas(f_index, ip + 1) + abs_bulk_par(f_index, ip + 1));
552 ssalb(f_index, Np - 2 - ip) = (ext -
abs) / ext;
556 "pressure level = ", ip,
"\n",
557 "frequency number = ", f_index,
"\n");
561 dtauc(f_index, Np - 2 - ip) = ext * (z_profile[ip + 1] - z_profile[ip]);
567 const Index& Npfct) {
568 const Index min_nang = 3;
572 Index this_ss = 0, this_se = 0;
574 for (
Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
575 for (
Index i_se = scat_data[i_ss].nelem() - 1; i_se >= 0; i_se--)
580 if (nang < scat_data[i_ss][i_se].za_grid.
nelem()) {
581 nang = scat_data[i_ss][i_se].za_grid.
nelem();
585 pfct_angs = scat_data[this_ss][this_se].za_grid;
586 }
else if (Npfct < min_nang) {
588 os <<
"Number of requested angular grid points (Npfct=" << Npfct
589 <<
") is insufficient.\n"
590 <<
"At least " << min_nang <<
" points required.\n";
591 throw runtime_error(os.str());
610 Vector T_array = t_profile[
Range(cloudbox_limits[0], Np_cloud)];
611 Matrix idir_array(1, 2, 0.);
613 Matrix pdir_array(nang, 2, 0.);
614 pdir_array(
joker, 0) = pfct_angs;
643 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
654 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
666 for (
Index ip = cloudbox_limits[0]; ip < Np_cloud - 1; ip++)
667 for (
Index f_index = 0; f_index < nf; f_index++) {
671 (ext_bulk_par(f_index, ip) + ext_bulk_par(f_index, ip + 1)) -
672 (abs_bulk_par(f_index, ip) + abs_bulk_par(f_index, ip + 1));
676 for (
Index ia = 0; ia < nang; ia++)
677 pfct_bulk_par(f_index, Np - 2 - ip, ia) +=
678 pha_bulk_par(f_index, ip, ia) + pha_bulk_par(f_index, ip + 1, ia);
679 pfct_bulk_par(f_index, Np - 2 - ip,
joker) *= 4 *
PI / sca;
687 const Index& Nlegendre) {
704 Vector u(nang), adu(nang - 1), ad_angs(nang - 1);
705 Tensor3 px(nang - 1, Nlegendre, 2, 0.);
706 u[0] = cos(pfct_angs[0] *
PI / 180.);
708 for (
Index ia = 1; ia < nang; ia++) {
709 u[ia] = cos(pfct_angs[ia] *
PI / 180.);
710 adu[ia - 1] =
abs(
u[ia] -
u[ia - 1]);
712 abs(pfct_angs[ia] *
PI / 180. - pfct_angs[ia - 1] *
PI / 180.);
713 px(ia - 1, 1, 0) =
u[ia - 1];
714 px(ia - 1, 1, 1) =
u[ia];
715 for (
Index l = 2; l < Nlegendre; l++) {
717 px(ia - 1, l, 0) = (2 * dl - 1) / dl *
u[ia - 1] * px(ia - 1, l - 1, 0) -
718 (dl - 1) / dl * px(ia - 1, l - 2, 0);
719 px(ia - 1, l, 1) = (2 * dl - 1) / dl *
u[ia] * px(ia - 1, l - 1, 1) -
720 (dl - 1) / dl * px(ia - 1, l - 2, 1);
724 for (
Index il = 0; il < nlyr; il++)
725 if (pfct_bulk_par(
joker, il, 0).sum() != 0.)
726 for (
Index f_index = 0; f_index < nf; f_index++) {
727 if (pfct_bulk_par(f_index, il, 0) != 0) {
736 for (
Index ia = 0; ia < nang - 1; ia++) {
737 pint += 0.5 * ad_angs[ia] *
738 (pfct[ia] * sin(pfct_angs[ia] *
PI / 180.) +
739 pfct[ia + 1] * sin(pfct_angs[ia + 1] *
PI / 180.));
742 if (
abs(pint / 2. - 1.) > pfct_threshold) {
744 os <<
"Phase function normalization deviates from expected value by\n"
745 << 1e2 * pint / 2. - 1e2 <<
"(allowed: " << pfct_threshold * 1e2
747 <<
"Occurs at layer #" << il <<
" and frequency #" << f_index
749 <<
"Something is wrong with your scattering data. Check!\n";
750 throw runtime_error(os.str());
756 pmom(f_index, il, 0) = 1.;
757 for (
Index ia = 0; ia < nang - 1; ia++) {
758 for (
Index l = 1; l < Nlegendre; l++) {
759 pmom(f_index, il, l) +=
761 (px(ia, l, 0) * pfct[ia] * sin(pfct_angs[ia] *
PI / 180.) +
762 px(ia, l, 1) * pfct[ia + 1] *
763 sin(pfct_angs[ia + 1] *
PI / 180.));
784 for (
Index ip = 0; ip < Np - 1; ip++)
786 for (
Index f_index = 0; f_index < nf; f_index++) {
788 0.5 * (ext_bulk(f_index, ip) - abs_bulk(f_index, ip) +
789 ext_bulk(f_index, ip + 1) - abs_bulk(f_index, ip + 1));
791 sca_bulk_layer(f_index, Np - 2 - ip) = sca;
802#define MAX_WARNINGS 100
807 static int warning_limit = FALSE, num_warnings = 0;
811 if (warning_limit)
return;
815 out1 <<
"DISORT WARNING >>> " << messag <<
"\n";
818 out1 <<
"DISORT TOO MANY WARNING MESSAGES -- They will no longer be "
819 "printed <<<<<<<\n\n";
820 warning_limit = TRUE;
830 const int maxmsg = 50;
831 static int nummsg = 0;
834 if (quiet != QUIET) {
837 out1 <<
" **** Input variable " << varnam <<
" in error ****\n";
838 if (nummsg == maxmsg) {
839 c_errmsg(
"Too many input errors. Aborting...", DS_ERROR);
848 if (quiet != QUIET) {
851 out1 <<
" **** Symbolic dimension " << dimnam
852 <<
" should be increased to at least " << minval <<
" ****\n";
873 if (
abs(z_surface - z_profile[0]) < 1e-3) {
879 cboxlims = cloudbox_limits;
887 for (; z_surface >= z_profile[ifirst + 1]; ifirst++) {
891 Range ind(ifirst, np);
895 vmr = vmr_profiles(
joker, ind);
900 gridpos(gp[0], z_profile, z_surface);
904 t[0] =
interp(itw, t, gp[0]);
905 for (
int i = 0; i < vmr.
nrows(); i++) {
912 itw2p(p[0], p, gp, itw2);
914 cboxlims = cloudbox_limits;
915 if (ifirst < cloudbox_limits[0]) {
916 cboxlims[0] -= ifirst;
917 cboxlims[1] -= ifirst;
921 ncboxremoved = ifirst - cboxlims[0];
923 cboxlims[1] = cloudbox_limits[1] - cloudbox_limits[0] - ncboxremoved;
924 ind =
Range(ncboxremoved, cboxlims[1] + 1);
925 pnd = pnd_profiles(
joker, ind);
926 gp[0].idx -= cloudbox_limits[0] + ncboxremoved;
927 for (
int i = 0; i < pnd.
nrows(); i++) {
946 const Agenda& propmat_clearsky_agenda,
947 const Agenda& gas_scattering_agenda,
950 const Vector& surface_scalar_reflectivity,
954 const Index& gas_scattering_do,
955 const Index& suns_do,
958 const Index& nstreams,
960 const Index& only_tro,
962 const Index& emission,
963 const Index& intensity_correction,
1007 nphi = aa_grid.
nelem();
1009 phi0 = sun_rte_los[1];
1016 ds.flag.prnt[0] = FALSE;
1017 ds.flag.prnt[1] = FALSE;
1018 ds.flag.prnt[2] = FALSE;
1019 ds.flag.prnt[3] = FALSE;
1020 ds.flag.prnt[4] = TRUE;
1022 ds.flag.usrtau = FALSE;
1023 ds.flag.usrang = TRUE;
1024 ds.flag.spher = FALSE;
1025 ds.flag.general_source = FALSE;
1026 ds.flag.output_uum = FALSE;
1028 ds.nlyr =
static_cast<int>(p.
nelem() - 1);
1030 ds.flag.brdf_type = BRDF_NONE;
1032 ds.flag.ibcnd = GENERAL_BC;
1033 ds.flag.usrang = TRUE;
1036 ds.flag.planck = TRUE;
1038 ds.flag.planck = FALSE;
1040 ds.flag.onlyfl = FALSE;
1041 ds.flag.lamber = TRUE;
1042 ds.flag.quiet = FALSE;
1043 if (intensity_correction) {
1044 ds.flag.intensity_correction = TRUE;
1045 ds.flag.old_intensity_correction = FALSE;
1047 ds.flag.intensity_correction = FALSE;
1048 ds.flag.old_intensity_correction = FALSE;
1051 ds.nstr =
static_cast<int>(nstreams);
1052 ds.nphase = ds.nstr;
1055 ds.numu =
static_cast<int>(za_grid.
nelem());
1056 ds.nphi =
static_cast<int>(nphi);
1057 Index Nlegendre = nstreams + 1;
1060 c_disort_state_alloc(&ds);
1061 c_disort_out_alloc(&ds, &out);
1071 for (
Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i];
1073 if (ds.flag.planck==TRUE){
1074 for (
Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1078 for (
Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] *
PI / 180);
1081 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1082 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1087 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1088 Index nf_ssd = scat_data[0][0].f_grid.
nelem();
1091 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1095 pha_bulk_par.
resize(nf_ssd, ds.nlyr + 1, nang);
1103 for (
Index iss = 0; iss < scat_data.
nelem(); iss++) {
1117 get_angs(pfct_angs, scat_data, Npfct);
1118 nang = pfct_angs.
nelem();
1120 pha_bulk_par.
resize(nf_ssd, ds.nlyr + 1, nang);
1123 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1124 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1127 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1128 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1131 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1132 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1134 if (gas_scattering_do) {
1138 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1142 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1143 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1144 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1147 sca_coeff_gas_layer,
1148 sca_coeff_gas_level,
1154 gas_scattering_agenda);
1158 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1161 ext_bulk_par += sca_coeff_gas_level;
1165 Matrix dtauc(nf, ds.nlyr);
1167 Matrix ssalb(nf, ds.nlyr);
1168 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1187 ds.bc.btemp = surface_skin_t;
1190 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++) {
1191 snprintf(ds.header, 128,
"ARTS Calc f_index = %" PRId64, f_index);
1193 std::memcpy(ds.dtauc,
1194 dtauc(f_index,
joker).get_c_array(),
1196 std::memcpy(ds.ssalb,
1197 ssalb(f_index,
joker).get_c_array(),
1201 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. *
SPEED_OF_LIGHT);
1202 ds.wvnmhi += ds.wvnmhi * 1e-7;
1203 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1206 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1210 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1213 ds.bc.fbeam = fbeam;
1215 std::memcpy(ds.pmom,
1219 enum class Status { FIRST_TRY, RETRY, SUCCESS };
1220 Status tries = Status::FIRST_TRY;
1224 c_disort(&ds, &out);
1225 tries = Status::SUCCESS;
1226 }
catch (
const std::runtime_error& e) {
1228 if (tries == Status::FIRST_TRY) {
1230 if (umu0 < 1 - eps) {
1232 }
else if (umu0 > 1 - eps) {
1240 <<
"Solar zenith angle coincided with one of the quadrature angles\n"
1241 <<
"We needed to shift the solar sun angle by " << shift
1245 tries = Status::RETRY;
1249 }
while (tries != Status::SUCCESS);
1251 for (
Index i = 0; i < ds.nphi; i++) {
1252 for (
Index j = 0; j < ds.numu; j++) {
1253 for (
Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1254 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, i, 0) =
1255 out.uu[j + ((ds.nlyr - k - cboxlims[0]) + i * (ds.nlyr + 1)) *
1261 for (
Index k = ncboxremoved - 1; k >= 0; k--) {
1262 cloudbox_field(f_index, k, 0, 0, j, i, 0) =
1263 cloudbox_field(f_index, k + 1, 0, 0, j, i, 0);
1270 disort_aux.resize(disort_aux_vars.
nelem());
1273 for (
Index i = 0; i < disort_aux_vars.
nelem(); i++) {
1276 if (disort_aux_vars[i] ==
"Layer optical thickness"){
1278 Matrix deltatau(nf, ds.nlyr, 0);
1279 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++) {
1280 for (
Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1281 deltatau(f_index, k - 1 + ncboxremoved) =
1282 dtauc(f_index, ds.nlyr - k + ncboxremoved);
1285 disort_aux[cnt] = deltatau;
1287 else if (disort_aux_vars[i] ==
"Single scattering albedo"){
1289 Matrix snglsctalbedo(nf, ds.nlyr,0);
1290 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++) {
1291 for (
Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1292 snglsctalbedo(f_index, k - 1 + ncboxremoved) =
1293 ssalb(f_index, ds.nlyr - k + ncboxremoved);
1296 disort_aux[cnt]=snglsctalbedo;
1298 else if (disort_aux_vars[i] ==
"Direct beam") {
1300 Matrix directbeam(nf, ds.nlyr + 1, 0);
1302 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++) {
1303 directbeam(f_index, cboxlims[1] - cboxlims[0] + ncboxremoved) =
1304 suns[0].spectrum(f_index, 0)/
PI;
1306 for (
Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1307 directbeam(f_index, k - 1 + ncboxremoved) =
1308 directbeam(f_index, k + ncboxremoved) *
1309 exp(-dtauc(f_index, ds.nlyr - k + ncboxremoved)/umu0);
1313 disort_aux[cnt]=directbeam;
1316 "The only allowed strings in *disort_aux_vars* are:\n"
1317 " \"Layer optical thickness\"\n"
1318 " \"Single scattering albedo\"\n"
1319 " \"Direct beam\"\n"
1320 "but you have selected: \"", disort_aux_vars[i],
"\"\n");
1325 c_disort_out_free(&ds, &out);
1326 c_disort_state_free(&ds);
1330 Tensor5& spectral_irradiance_field,
1341 const Agenda& propmat_clearsky_agenda,
1342 const Agenda& gas_scattering_agenda,
1344 const Numeric& surface_skin_t,
1345 const Vector& surface_scalar_reflectivity,
1347 const Index& gas_scattering_do,
1348 const Index& suns_do,
1351 const Index& nstreams,
1353 const Index& only_tro,
1355 const Index& emission,
1356 const Index& intensity_correction,
1399 phi0 = sun_rte_los[1];
1406 ds.flag.prnt[0] = FALSE;
1407 ds.flag.prnt[1] = FALSE;
1408 ds.flag.prnt[2] = FALSE;
1409 ds.flag.prnt[3] = FALSE;
1410 ds.flag.prnt[4] = TRUE;
1412 ds.flag.usrtau = FALSE;
1413 ds.flag.usrang = TRUE;
1414 ds.flag.spher = FALSE;
1415 ds.flag.general_source = FALSE;
1416 ds.flag.output_uum = FALSE;
1418 ds.nlyr =
static_cast<int>(p.
nelem() - 1);
1420 ds.flag.brdf_type = BRDF_NONE;
1422 ds.flag.ibcnd = GENERAL_BC;
1423 ds.flag.usrang = FALSE;
1426 ds.flag.planck = TRUE;
1428 ds.flag.planck = FALSE;
1430 ds.flag.onlyfl = TRUE;
1431 ds.flag.lamber = TRUE;
1432 ds.flag.quiet = FALSE;
1433 if (intensity_correction) {
1434 ds.flag.intensity_correction = TRUE;
1435 ds.flag.old_intensity_correction = FALSE;
1437 ds.flag.intensity_correction = FALSE;
1438 ds.flag.old_intensity_correction = FALSE;
1441 ds.nstr =
static_cast<int>(nstreams);
1442 ds.nphase = ds.nstr;
1450 Index Nlegendre = nstreams + 1;
1453 c_disort_state_alloc(&ds);
1454 c_disort_out_alloc(&ds, &out);
1464 if (ds.flag.planck==TRUE){
1465 for (
Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1469 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1470 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1475 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1476 Index nf_ssd = scat_data[0][0].f_grid.
nelem();
1479 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1483 pha_bulk_par.
resize(nf_ssd, ds.nlyr + 1, nang);
1491 for (
Index iss = 0; iss < scat_data.
nelem(); iss++) {
1505 get_angs(pfct_angs, scat_data, Npfct);
1506 nang = pfct_angs.
nelem();
1508 pha_bulk_par.
resize(nf_ssd, ds.nlyr + 1, nang);
1511 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1512 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1515 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1516 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1519 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1520 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1522 if (gas_scattering_do) {
1526 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1530 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1531 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1532 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1535 sca_coeff_gas_layer,
1536 sca_coeff_gas_level,
1542 gas_scattering_agenda);
1546 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1549 ext_bulk_par += sca_coeff_gas_level;
1553 Matrix dtauc(nf, ds.nlyr);
1555 Matrix ssalb(nf, ds.nlyr);
1556 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1575 ds.bc.btemp = surface_skin_t;
1578 Matrix spectral_direct_irradiance_field;
1579 Matrix dFdtau(nf, ds.nlyr+1,0);
1580 Matrix deltatau(nf, ds.nlyr,0);
1581 Matrix snglsctalbedo(nf, ds.nlyr,0);
1585 spectral_direct_irradiance_field.
resize(nf, ds.nlyr+1);
1586 spectral_direct_irradiance_field = 0;
1589 for (
Index f_index = 0; f_index < f_grid.
nelem(); f_index++) {
1590 snprintf(ds.header, 128,
"ARTS Calc f_index = %" PRId64, f_index);
1592 std::memcpy(ds.dtauc,
1593 dtauc(f_index,
joker).get_c_array(),
1595 std::memcpy(ds.ssalb,
1596 ssalb(f_index,
joker).get_c_array(),
1600 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. *
SPEED_OF_LIGHT);
1601 ds.wvnmhi += ds.wvnmhi * 1e-7;
1602 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1605 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1609 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1612 ds.bc.fbeam = fbeam;
1614 std::memcpy(ds.pmom,
1618 c_disort(&ds, &out);
1623 for (
Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1626 spectral_direct_irradiance_field(f_index, k + ncboxremoved) =
1627 -out.rad[ds.nlyr - k - cboxlims[0]].rfldir/conv_fac;
1630 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1631 -(out.rad[ds.nlyr - k - cboxlims[0]].rfldir +
1632 out.rad[ds.nlyr - k - cboxlims[0]].rfldn)/conv_fac;
1636 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1637 -out.rad[ds.nlyr - k - cboxlims[0]].rfldn/conv_fac;
1641 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 1) =
1642 out.rad[ds.nlyr - k - cboxlims[0]].flup/conv_fac;
1645 dFdtau(f_index, k + ncboxremoved) =
1646 -out.rad[ds.nlyr - k - cboxlims[0]].dfdt;
1651 deltatau(f_index, k - 1 + ncboxremoved) = ds.dtauc[ds.nlyr - k - 1 - cboxlims[0]];
1652 snglsctalbedo(f_index, k - 1 + ncboxremoved) = ds.ssalb[ds.nlyr - k - 1 - cboxlims[0]];
1658 for (
Index k = ncboxremoved - 1; k >= 0; k--) {
1659 spectral_irradiance_field(f_index, k, 0, 0,
joker) =
1660 spectral_irradiance_field(f_index, k + 1, 0, 0,
joker);
1663 spectral_direct_irradiance_field(f_index, k) =
1664 spectral_direct_irradiance_field(f_index, k + 1);
1666 dFdtau(f_index, k) = dFdtau(f_index, k + 1);
1671 disort_aux.resize(disort_aux_vars.
nelem());
1674 for (
Index i = 0; i < disort_aux_vars.
nelem(); i++) {
1677 if (disort_aux_vars[i] ==
"Layer optical thickness"){
1679 disort_aux[cnt]=deltatau;
1681 else if (disort_aux_vars[i] ==
"Single scattering albedo"){
1683 disort_aux[cnt]=snglsctalbedo;
1685 else if (disort_aux_vars[i] ==
"Direct downward spectral irradiance") {
1687 disort_aux[cnt] = spectral_direct_irradiance_field;
1689 else if (disort_aux_vars[i] ==
"dFdtau") {
1691 disort_aux[cnt] = dFdtau;
1695 "The only allowed strings in *disort_aux_vars* are:\n"
1696 " \"Layer optical thickness\"\n"
1697 " \"Single scattering albedo\"\n"
1698 " \"Direct downward spectral irradiance\"\n"
1700 "but you have selected: \"", disort_aux_vars[i],
"\"\n");
1707 c_disort_out_free(&ds, &out);
1708 c_disort_state_free(&ds);
1716 const Agenda& surface_rtprop_agenda,
1747 chk_not_empty(
"surface_rtprop_agenda", surface_rtprop_agenda);
1751 while (frza < scat_za_grid.
nelem() && scat_za_grid[frza] < 90.) frza++;
1752 if (frza == scat_za_grid.
nelem()) {
1754 os <<
"No upwelling direction found in scat_za_grid.\n";
1755 throw runtime_error(os.str());
1757 const Index nrza = scat_za_grid.
nelem() - frza;
1760 Matrix dir_refl_coeff(nrza, nf, 0.);
1763 Vector rtp_pos(1, surf_alt);
1767 for (
Index rza = 0; rza < nrza; rza++) {
1774 Vector rtp_los(1, scat_za_grid[rza + frza]);
1775 out2 <<
"Doing reflected dir #" << rza <<
" at " << rtp_los[0] <<
" degs\n";
1785 surface_rtprop_agenda);
1786 if (surface_los.
nrows() > 1) {
1788 os <<
"For this method, *surface_rtprop_agenda* must be set up to\n"
1789 <<
"return zero or one direction in *surface_los*.";
1790 throw runtime_error(os.str());
1794 btemp = surface_skin_t;
1795 else if (surface_skin_t != btemp) {
1797 os <<
"Something went wrong.\n"
1798 <<
" *surface_rtprop_agenda* returned different surface_skin_t\n"
1799 <<
" for different LOS.\n";
1800 throw runtime_error(os.str());
1804 if (!surface_los.
empty()) {
1805 for (
Index f_index = 0; f_index < nf; f_index++)
1806 dir_refl_coeff(rza, f_index) =
1807 surface_rmatrix(
joker, f_index, 0, 0).sum();
1809 out2 <<
" directional albedos[f_grid] = " << dir_refl_coeff(rza,
joker)
1813 if (btemp < 0. || btemp > 1000.) {
1815 os <<
"Surface temperature has been derived as " << btemp <<
" K,\n"
1816 <<
"which is not considered a meaningful value.\n";
1817 throw runtime_error(os.str());
1824 Vector surf_int_grid(nrza + 1);
1831 os <<
"Looks like scat_za_grid contains the 90deg direction,\n"
1832 <<
"which it shouldn't for running Disort.\n";
1833 throw runtime_error(os.str());
1835 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1836 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1837 const Numeric ok_extrapol = 0.5;
1838 if ((za_extrapol - ok_extrapol) > 1e-6) {
1840 os <<
"Extrapolation range from shallowest scat_za_grid point\n"
1841 <<
"to horizon is too big.\n"
1842 <<
" Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1843 <<
", which is " << za_extrapol - 0.5 <<
" too big.\n";
1844 throw runtime_error(os.str());
1847 scat_za_grid[scat_za_grid.
nelem() - 1], 180., 1e-6)) {
1849 os <<
"Looks like last point in scat_za_grid is not 180deg.\n";
1850 throw runtime_error(os.str());
1853 surf_int_grid[0] = 90.;
1854 surf_int_grid[nrza] = 180.;
1855 for (
Index rza = 1; rza < nrza; rza++)
1856 surf_int_grid[rza] =
1857 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1861 for (
Index rza = 0; rza < nrza; rza++) {
1866 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1872 dir_refl_coeff(rza,
joker) *=
w;
1877 for (
Index f_index = 0; f_index < nf; f_index++) {
1878 albedo[f_index] = dir_refl_coeff(
joker, f_index).sum();
1879 out2 <<
"at f=" << f_grid[f_index] * 1e-9
1880 <<
" GHz, ending up with albedo=" << albedo[f_index] <<
"\n";
1881 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1883 os <<
"Something went wrong: Albedo must be inside [0,1],\n"
1884 <<
" but is not at freq #" << f_index <<
" , where it is "
1885 << albedo[f_index] <<
".\n";
1886 throw runtime_error(os.str());
1896 const Agenda& surface_rtprop_agenda,
1901 Vector rtp_pos(1, surf_alt);
1902 Vector rtp_los(1, 180-inc_angle);
1917 surface_rtprop_agenda);
1919 if (surface_los.
nrows() > 1) {
1921 os <<
"For this method, *surface_rtprop_agenda* must be set up to\n"
1922 <<
"return zero or one direction in *surface_los*.";
1923 throw runtime_error(os.str());
1925 if (surface_los.
empty()) {
1928 albedo[
joker] = surface_rmatrix(0,
joker,0,0);
1931 btemp = surface_skin_t;
1932 if (btemp < 0. || btemp > 1000.) {
1934 os <<
"Surface temperature has been derived as " << btemp <<
" K,\n"
1935 <<
"which is not considered a meaningful value.\n";
1936 throw runtime_error(os.str());
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
Constants of physical expressions as constexpr.
void gas_scattering_agendaExecute(Workspace &ws, PropagationMatrix &gas_scattering_coef, TransmissionMatrix &gas_scattering_mat, Vector &gas_scattering_fct_legendre, const Vector &f_grid, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Index gas_scattering_output_type, const Agenda &input_agenda)
void propmat_clearsky_agendaExecute(Workspace &ws, PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfSpeciesTag &select_abs_species, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
void surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Index nelem() const ARTS_NOEXCEPT
A constant view of a Matrix.
bool empty() const noexcept
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Tensor3.
Index npages() const
Returns the number of pages.
Index nrows() const
Returns the number of rows.
Index ncols() const
Returns the number of columns.
Index nbooks() const noexcept
A constant view of a Vector.
Index nelem() const noexcept
Returns the number of elements.
void resize(Index r, Index c)
Resize function.
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 l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
void surf_albedoCalc(Workspace &ws, VectorView albedo, Numeric &btemp, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, ConstVectorView scat_za_grid, const Numeric &surf_alt, const Verbosity &verbosity)
surf_albedoCalc
constexpr Numeric PLANCK_CONST
void c_errmsg(const char *messag, int type)
Verbosity enabled replacement for the original cdisort function.
void get_gas_scattering_properties(Workspace &ws, MatrixView &sca_coeff_gas, MatrixView &sca_coeff_gas_level, MatrixView &pfct_gas, const ConstVectorView &f_grid, const VectorView &p, const VectorView &t, const MatrixView &vmr, const Agenda &gas_scattering_agenda)
get_gas_scattering_properties
void get_pfct(Tensor3 &pfct_bulk_par, ConstTensor3View &pha_bulk_par, ConstMatrixView ext_bulk_par, ConstMatrixView abs_bulk_par, const ArrayOfIndex &cloudbox_limits)
get_pfct.
void reduced_1datm(Vector &p, Vector &z, Vector &t, Matrix &vmr, Matrix &pnd, ArrayOfIndex &cboxlims, Index &ncboxremoved, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfIndex &cloudbox_limits)
reduced_1datm
constexpr Numeric SPEED_OF_LIGHT
void run_cdisort(Workspace &ws, Tensor7 &cloudbox_field, ArrayOfMatrix &disort_aux, ConstVectorView f_grid, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &gas_scattering_agenda, const ArrayOfIndex &cloudbox_limits, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView sun_rte_los, const Index &gas_scattering_do, const Index &suns_do, const ArrayOfString &disort_aux_vars, const Numeric &scale_factor, const Index &nstreams, const Index &Npfct, const Index &only_tro, const Index &quiet, const Index &emission, const Index &intensity_correction, const Verbosity &verbosity)
Calculate doit_i_field with Disort including a sun source.
constexpr Numeric DEG2RAD
void get_parZ(Tensor3 &pha_bulk_par, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstMatrixView pnd_profiles, ConstVectorView t_profile, ConstVectorView pfct_angs, const ArrayOfIndex &cloudbox_limits)
get_parZ.
void run_cdisort_flux(Workspace &ws, Tensor5 &spectral_irradiance_field, ArrayOfMatrix &disort_aux, ConstVectorView f_grid, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &gas_scattering_agenda, const ArrayOfIndex &cloudbox_limits, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, ConstVectorView sun_rte_los, const Index &gas_scattering_do, const Index &suns_do, const ArrayOfString &disort_aux_vars, const Numeric &scale_factor, const Index &nstreams, const Index &Npfct, const Index &only_tro, const Index &quiet, const Index &emission, const Index &intensity_correction, const Verbosity &verbosity)
Calculate spectral_irradiance_field with Disort including a sun source.
void get_disortsurf_props(Vector &albedo, Numeric &btemp, ConstVectorView f_grid, const Numeric &surface_skin_t, ConstVectorView surface_scalar_reflectivity)
get_disortsurf_props.
void surf_albedoCalcSingleAngle(Workspace &ws, VectorView albedo, Numeric &btemp, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, const Numeric &surf_alt, const Numeric &inc_angle)
surf_albedoCalcSingleAngle
void get_scat_bulk_layer(MatrixView &sca_bulk_layer, const MatrixView &ext_bulk, const MatrixView &abs_bulk)
get_scat_bulk_layer
void check_disort_irradiance_input(const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &nstreams)
check_disort_input.
int c_write_bad_var(int quiet, const char *varnam)
Verbosity enabled replacement for the original cdisort function.
void get_dtauc_ssalb(MatrixView dtauc, MatrixView ssalb, ConstMatrixView ext_bulk_gas, ConstMatrixView ext_bulk_par, ConstMatrixView abs_bulk_par, ConstVectorView z_profile)
get_dtauc_ssalb
void check_disort_input(const Index &cloudbox_on, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfIndex &cloudbox_limits, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstVectorView za_grid, const Index &nstreams)
check_disort_input.
void init_ifield(Tensor7 &cloudbox_field, const Vector &f_grid, const ArrayOfIndex &cloudbox_limits, const Index &n_za, const Index &n_aa, const Index &stokes_dim)
init_ifield.
thread_local Verbosity disort_verbosity
void get_angs(Vector &pfct_angs, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &Npfct)
get_angs.
void get_pmom(Tensor3View pmom, ConstTensor3View pfct_bulk_par, ConstVectorView pfct_angs, const Index &Nlegendre)
get_pmom
void get_gasoptprop(Workspace &ws, MatrixView ext_bulk_gas, const Agenda &propmat_clearsky_agenda, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstVectorView p_grid, ConstVectorView f_grid)
get_gasoptprop.
void add_normed_phase_functions(Tensor3View &pfct1, const MatrixView &sca1, const MatrixView &pfct2, const MatrixView &sca2)
add_normed_phase_functions
constexpr Numeric COSMIC_BG_TEMP
void get_paroptprop(MatrixView ext_bulk_par, MatrixView abs_bulk_par, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstMatrixView pnd_profiles, ConstVectorView t_profile, ConstVectorView DEBUG_ONLY(p_grid), const ArrayOfIndex &cloudbox_limits, ConstVectorView f_grid)
int c_write_too_small_dim(int quiet, const char *dimnam, int minval)
Verbosity enabled replacement for the original cdisort function.
Functions for disort interface.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Header file for logic.cc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
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 cosmic_microwave_background_temperature
Global constant, Planck temperature for cosmic background radiation [K].
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric planck_constant
Planck constant [J s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 20...
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto acosd(auto x) noexcept
Returns rad2deg of the arc-cosine of the input.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void pha_mat_ScatSpecBulk(ArrayOfTensor6 &pha_mat, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor6 &pha_mat_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk phase matrices.
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
void pha_mat_NScatElems(ArrayOfArrayOfTensor6 &pha_mat, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_index, const Index &t_interp_order)
Phase matrices from all scattering elements.
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
void ext_abs_pfun_from_tro(MatrixView ext_data, MatrixView abs_data, Tensor3View pfun_data, const ArrayOfSingleScatteringData &scat_data, const Index &iss, ConstMatrixView pnd_data, ArrayOfIndex &cloudbox_limits, ConstVectorView T_grid, ConstVectorView sa_grid)
Extinction, absorption and phase function for one scattering species, 1D and TRO.
Declaration of functions in rte.cc.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
This file contains basic functions to handle XML data files.