22#include "matpack_eigen.h"
25#include "propagationmatrix.h"
41 ArrayOfMatrix& iy_aux,
42 ArrayOfTensor3& diy_dx,
50 const Index& stokes_dim,
52 const Index& atmosphere_dim,
54 const Tensor3& t_field,
56 const Tensor4& vmr_field,
58 const Tensor3& wind_u_field,
59 const Tensor3& wind_v_field,
60 const Tensor3& wind_w_field,
61 const Tensor3& mag_u_field,
62 const Tensor3& mag_v_field,
63 const Tensor3& mag_w_field,
64 const Index& cloudbox_on,
66 const Tensor4& pnd_field,
67 const ArrayOfTensor4& dpnd_field_dx,
70 const Index& scat_data_checked,
72 const Index& jacobian_do,
75 const Matrix& iy_transmitter,
76 const Agenda& propmat_clearsky_agenda,
77 const Agenda& water_p_eq_agenda,
78 const Numeric& rte_alonglos_v,
79 const Index& trans_in_jacobian,
80 const Numeric& pext_scaling,
81 const Index& t_interp_order,
82 const Verbosity& verbosity [[maybe_unused]]) {
84 const Index j_analytical_do = jacobian_do ?
85 do_analytical_jacobian<1>(jacobian_quantities) : 0;
88 const Index nf = f_grid.nelem();
89 const Index ns = stokes_dim;
90 const Index np = ppath.
np;
91 const Index ne = pnd_field.nbooks();
92 const Index nq = j_analytical_do ? jacobian_quantities.
nelem() : 0;
93 const Index naux = iy_aux_vars.
nelem();
101 "ppath.background is invalid. Check your "
102 "calculation of *ppath*?");
104 "The propagation path ends inside or at boundary of "
105 "the cloudbox.\nFor this method, *ppath* must be "
106 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
109 "The scattering data must be flagged to have\n"
110 "passed a consistency check (scat_data_checked=1).");
112 "*pnd_field* and *scat_data* inconsistent regarding total number of"
113 " scattering elements.");
118 "revision before safe to use.");
121 "*dpnd_field_dx* not properly initialized:\n"
122 "Number of elements in dpnd_field_dx must be equal number of jacobian"
123 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
124 " is calculated/set.");
131 "The size of *iy* returned from *iy_transmitter_agenda* is\n"
133 " expected size = [", nf,
",", stokes_dim,
"]\n"
134 " size of iy = [", iy_transmitter.nrows(),
",", iy_transmitter.ncols(),
"]\n")
135 for (Index iv = 0; iv < nf; iv++)
137 "The *iy* returned from *iy_transmitter_agenda* "
138 "must have the value 1 in the first column.");
141 ArrayOfTensor3 diy_dpath = j_analytical_do ?
152 diy_dx = j_analytical_do ?
165 Index auxBackScat = -1;
166 Index auxAbSpAtte = -1;
167 Index auxPartAtte = -1;
169 for (Index i = 0; i < naux; i++) {
170 iy_aux[i].resize(nf * np, ns);
172 if (iy_aux_vars[i] ==
"Radiative background") {
174 iy_aux[i](joker, 0) = (Numeric)
min((Index)2, rbi - 1);
175 }
else if( iy_aux_vars[i] ==
"Backscattering" ) {
178 }
else if( iy_aux_vars[i] ==
"Abs species extinction" ) {
181 }
else if( iy_aux_vars[i] ==
"Particle extinction" ) {
186 "The only allowed strings in *iy_aux_vars* are:\n"
187 " \"Radiative background\"\n"
188 " \"Backscattering\"\n"
191 "but you have selected: \"", iy_aux_vars[i],
"\"")
209 Tensor5 Pe(ne, np, nf, ns, ns, 0);
212 ArrayOfMatrix ppvar_dpnd_dx(0);
215 const Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
216 const Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
217 Tensor6 pha_mat_1se(nf_ssd, 1, 1, 1, ns, ns);
218 Vector t_ok(1), t_array(1);
219 Matrix mlos_sca(1, 2), mlos_inc(1, 2);
223 if (np == 1 && rbi == 1) {
226 ppvar_vmr.resize(0, 0);
227 ppvar_wind.resize(0, 0);
228 ppvar_mag.resize(0, 0);
229 ppvar_pnd.resize(0, 0);
230 ppvar_f.resize(0, 0);
253 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
266 clear2cloudy.resize(np);
267 for (Index ip = 0; ip < np; ip++)
268 clear2cloudy[ip] = -1;
272 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
273 StokesVector
a(nf, ns), S(nf, ns);
278 ArrayOfPropagationMatrix dK_this_dx(0), dK_past_dx(0), dKp_dx(0);
279 ArrayOfStokesVector da_dx(0), dS_dx(0);
282 Index temperature_derivative_position = -1;
285 if (trans_in_jacobian && j_analytical_do) {
286 dK_this_dx.resize(nq);
287 dK_past_dx.resize(nq);
293 dK_this_dx[iq] = PropagationMatrix(nf, ns);
294 dK_past_dx[iq] = PropagationMatrix(nf, ns);
295 dKp_dx[iq] = PropagationMatrix(nf, ns);
296 da_dx[iq] = StokesVector(nf, ns);
297 dS_dx[iq] = StokesVector(nf, ns);
298 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
299 temperature_derivative_position = iq;
300 do_hse = jacobian_quantities[iq].Subtag() ==
"HSE on";
305 for (Index ip = 0; ip < np; ip++) {
312 propmat_clearsky_agenda,
314 Vector{ppvar_f(joker, ip)},
315 Vector{ppvar_mag(joker, ip)},
316 Vector{ppath.
los(ip, joker)},
318 Vector{ppvar_vmr(joker, ip)},
321 trans_in_jacobian && j_analytical_do);
323 if (trans_in_jacobian && j_analytical_do)
329 ppath.
los(ip, joker),
332 trans_in_jacobian && j_analytical_do);
334 if (clear2cloudy[ip] + 1) {
340 ppvar_pnd(joker, Range(ip, 1)),
344 ppath.
los(ip, joker),
345 ppvar_t[Range(ip, 1)],
347 trans_in_jacobian && jacobian_do);
349 if (abs(pext_scaling - 1) > 1e-6) {
351 if (trans_in_jacobian && j_analytical_do) {
357 if (auxAbSpAtte >= 0) {
358 for (Index iv = 0; iv < nf; iv++) {
359 iy_aux[auxAbSpAtte](iv*np+ip, joker) = K_this.Kjj()[iv];
362 if (auxPartAtte >= 0) {
363 for (Index iv = 0; iv < nf; iv++) {
364 iy_aux[auxPartAtte](iv*np+ip, joker) = Kp.Kjj()[iv];
370 if (trans_in_jacobian && j_analytical_do)
378 mlos_sca(0, joker) = los_sca;
382 if (atmosphere_dim == 3) {
383 los_inc = ppath.
los(ip, joker);
387 mlos_inc(0, joker) = los_inc;
390 for (Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
391 for (Index i_se = 0; i_se < scat_data[i_ss].
nelem(); i_se++) {
395 if (ppvar_pnd(i_se_flat, ip) != 0)
397 else if (j_analytical_do)
398 for (Index iq = 0; iq < nq && !val_pnd; iq++)
399 if (jac_scat_i[iq] >= 0)
400 if (ppvar_dpnd_dx[iq](i_se_flat, ip) != 0) {
408 scat_data[i_ss][i_se],
409 Vector{ppvar_t[Range(ip, 1)]},
414 if (t_ok[0] not_eq 0)
416 for (Index iv = 0; iv < nf; iv++)
417 Pe(i_se_flat, ip, iv, joker, joker) =
418 pha_mat_1se(0, 0, 0, 0, joker, joker);
420 Pe(i_se_flat, ip, joker, joker, joker) =
421 pha_mat_1se(joker, 0, 0, 0, joker, joker);
424 "Interpolation error for (flat-array) scattering"
425 " element #", i_se_flat,
"\n"
426 "at location/temperature point #", ip,
"\n")
436 const Numeric dr_dT_past =
437 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
438 const Numeric dr_dT_this =
439 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
450 temperature_derivative_position);
453 swap(K_past, K_this);
454 swap(dK_past_dx, dK_this_dx);
467 lvl_rad[0] = iy_transmitter;
469 rad_inc = iy_transmitter;
485 iy.resize(nf * np, ns);
487 for (Index ip = 0; ip < np; ip++) {
488 for (Index iv = 0; iv < nf; iv++) {
489 for (Index is = 0; is < stokes_dim; is++) {
490 iy(iv*np+ip, is) = lvl_rad[ip](iv, is);
491 if (j_analytical_do) {
493 diy_dpath[iq](ip, iv*np+ip2, is) =
494 dlvl_rad[ip][ip2][iq](iv, is););
501 if (j_analytical_do) {
502 const Index iy_agenda_call1 = 1;
503 const Tensor3 iy_transmittance(0, 0, 0);
524 if (auxBackScat >= 0) {
525 for (Index ip = 0; ip < np; ip++) {
526 for (Index iv = 0; iv < nf; iv++) {
527 VectorView stokesvec = VectorView(iy_aux[auxBackScat](iv*np+ip, joker));
528 matpack::eigen::row_vec(stokesvec).noalias() = reflect_matrix[ip].Mat(iv) * rad_inc.
Vec(iv);
543 ArrayOfVector& y_aux,
546 const Index& atmgeom_checked,
547 const Index& atmfields_checked,
548 const String& iy_unit_radar,
550 const Index& stokes_dim,
551 const Vector& f_grid,
552 const Index& cloudbox_on,
553 const Index& cloudbox_checked,
554 const Matrix& sensor_pos,
555 const Matrix& sensor_los,
556 const Index& sensor_checked,
557 const Index& jacobian_do,
559 const Agenda& iy_radar_agenda,
561 const Vector& range_bins,
562 const Numeric& ze_tref,
564 const Numeric& dbze_min,
567 const Index npos = sensor_pos.nrows();
568 const Index nbins = range_bins.nelem() - 1;
569 const Index nf = f_grid.nelem();
570 const Index naux = iy_aux_vars.
nelem();
586 "The atmospheric fields must be flagged to have "
587 "passed a consistency check (atmfields_checked=1).");
589 "The atmospheric geometry must be flagged to have "
590 "passed a consistency check (atmgeom_checked=1).");
592 "The cloudbox must be flagged to have "
593 "passed a consistency check (cloudbox_checked=1).");
595 "The sensor variables must be flagged to have "
596 "passed a consistency check (sensor_checked=1).");
599 bool is_z =
max(range_bins) > 1;
601 "The vector *range_bins* must contain strictly "
602 "increasing values.");
604 "The vector *range_bins* is not allowed to contain "
607 "The main length of *instrument_pol_array* must match "
608 "the number of frequencies.");
611 Vector cfac(nf, 1.0);
613 const Numeric jfac = 10 * log10(exp(1.0));
614 if (iy_unit_radar ==
"1") {
615 }
else if (iy_unit_radar ==
"Ze") {
616 ze_cfac(cfac, f_grid, ze_tref, k2);
617 }
else if (iy_unit_radar ==
"dBZe") {
618 ze_cfac(cfac, f_grid, ze_tref, k2);
619 ze_min = pow(10.0, dbze_min / 10);
622 "For this method, *iy_unit_radar* must be set to \"1\", "
623 "\"Ze\" or \"dBZe\".");
633 ArrayOfArrayOfVector W(nf);
634 for (Index i = 0; i < nf; i++) {
635 const Index ni = instrument_pol_array[i].
nelem();
636 npolcum[i + 1] = npolcum[i] + ni;
638 for (Index j = 0; j < ni; j++) {
639 W[i][j].resize(stokes_dim);
640 stokes2pol(W[i][j], stokes_dim, instrument_pol_array[i][j], 0.5);
646 const Index ntot = npos * npolcum[nf] * nbins;
651 y_pos.resize(ntot, sensor_pos.ncols());
652 y_los.resize(ntot, sensor_los.ncols());
653 y_geo.resize(ntot, 5);
658 for (Index i = 0; i < naux; i++) {
659 y_aux[i].resize(ntot);
665 Index j_analytical_do = 0;
673 jacobian.resize(ntot,
674 jacobian_indices[jacobian_indices.
nelem() - 1][1] + 1);
676 njq = jacobian_quantities.
nelem();
680 jacobian.resize(0, 0);
688 for (Index p = 0; p < npos; p++) {
690 ArrayOfTensor3 diy_dx;
694 ArrayOfMatrix iy_aux;
695 const Index iy_id = (Index)1e6 * p;
707 Vector{sensor_pos(p, joker)},
708 Vector{sensor_los(p, joker)},
712 const Index np = ppath.
np;
714 "A path consisting of a single point found. "
715 "This is not allowed.");
718 "The size of *iy* returned from *iy_radar_agenda* "
719 "is not correct (for this method).");
724 range = ppath.
pos(joker, 0);
727 for (Index i = 1; i < np; i++) {
728 range[i] = range[i - 1] + ppath.
lstep[i - 1] *
733 const Numeric range_end1 =
min(range[0], range[np - 1]);
734 const Numeric range_end2 =
max(range[0], range[np - 1]);
737 for (Index
b = 0;
b < nbins;
b++) {
738 if (!(range_bins[
b] >= range_end2 ||
739 range_bins[
b + 1] <= range_end1))
742 Numeric blim1 =
max(range_bins[
b], range_end1);
743 Numeric blim2 =
min(range_bins[
b + 1], range_end2);
750 hbin /= (blim2 - blim1);
752 for (Index iv = 0; iv < nf; iv++) {
754 Matrix I{iy(Range(iv * np, np), joker)};
755 ArrayOfTensor3 dI(njq);
756 if (j_analytical_do) {
758 dI[iq] = diy_dx[iq](joker, Range(iv * np, np), joker);)
760 ArrayOfMatrix A(naux);
761 for (Index
a = 0;
a < naux;
a++) {
762 A[
a] = iy_aux[
a](Range(iv * np, np), joker);
767 ArrayOfMatrix drefl(njq);
768 if (j_analytical_do) {
771 ArrayOfVector auxvar(naux);
772 for (Index
a = 0;
a < naux;
a++) {
773 auxvar[
a].resize(np);
776 for (Index ip = 0; ip < instrument_pol_array[iv].
nelem(); ip++) {
778 mult(refl, I, W[iv][ip]);
779 if (j_analytical_do) {
781 k < drefl[iq].nrows();
783 mult(drefl[iq](k, joker), dI[iq](k, joker, joker), W[iv][ip]);
786 for (Index
a = 0;
a < naux;
a++) {
787 if (iy_aux_vars[
a] ==
"Backscattering") {
788 mult(auxvar[
a], A[
a], W[iv][ip]);
790 for (Index j = 0; j < np; j++) {
791 auxvar[
a][j] = A[
a](j, 0);
797 Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) +
b;
798 y[iout] = cfac[iv] * (hbin * refl);
800 if (j_analytical_do) {
802 for (Index k = 0; k < drefl[iq].nrows(); k++) {
803 jacobian(iout, jacobian_indices[iq][0] + k) =
804 cfac[iv] * (hbin * drefl[iq](k, joker));
805 if (iy_unit_radar ==
"dBZe") {
806 jacobian(iout, jacobian_indices[iq][0] + k) *=
807 jfac /
max(y[iout], ze_min);
812 if (iy_unit_radar ==
"dBZe") {
813 y[iout] = y[iout] <= ze_min ? dbze_min : 10 * log10(y[iout]);
817 for (Index
a = 0;
a < naux;
a++) {
818 if (iy_aux_vars[
a] ==
"Backscattering") {
819 y_aux[
a][iout] = cfac[iv] * (hbin * auxvar[
a]);
820 if (iy_unit_radar ==
"dBZe") {
821 y_aux[
a][iout] = y_aux[
a][iout] <= ze_min
823 : 10 * log10(y_aux[
a][iout]);
826 y_aux[
a][iout] = hbin * auxvar[
a];
831 if (geo_pos.nelem()) y_geo(iout, joker) = geo_pos;
839 for (Index
b = 0;
b < nbins;
b++) {
840 for (Index iv = 0; iv < nf; iv++) {
841 for (Index ip = 0; ip < instrument_pol_array[iv].
nelem(); ip++) {
842 const Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) +
b;
843 y_f[iout] = f_grid[iv];
844 y_pol[iout] = instrument_pol_array[iv][ip];
845 y_pos(iout, joker) = sensor_pos(p, joker);
846 y_los(iout, joker) = sensor_los(p, joker);
858 Tensor4& particle_bulkprop_field,
860 const Index& atmosphere_dim,
861 const Vector& p_grid,
862 const Vector& lat_grid,
863 const Vector& lon_grid,
864 const Tensor3& t_field,
865 const Tensor3& z_field,
866 const Tensor4& vmr_field,
867 const Matrix& z_surface,
868 const Index& atmfields_checked,
869 const Index& atmgeom_checked,
870 const Vector& f_grid,
871 const Agenda& propmat_clearsky_agenda,
874 const Matrix& incangles,
876 const Numeric& dbze_noise,
877 const Matrix& h_clutter,
878 const Index& fill_clutter,
879 const Numeric& t_phase,
880 const Numeric& wc_max,
881 const Numeric& wc_clip,
882 const Index& do_atten_abs,
883 const Index& do_atten_hyd,
884 const Numeric& atten_hyd_scaling,
885 const Numeric& atten_hyd_max,
888 const Index np = t_field.npages();
889 const Index nlat = t_field.nrows();
890 const Index nlon = t_field.ncols();
894 "The atmospheric fields must be flagged to have\n"
895 "passed a consistency check (atmfields_checked=1).");
897 "The atmospheric geometry must be flagged to have\n"
898 "passed a consistency check (atmgeom_checked=1).");
902 "*dbze_noise* not covered by invtable[0]." );
904 "*dbze_noise* not covered by invtable[1]." );
905 chk_atm_field(
"GIN reflectivities", dBZe, atmosphere_dim, p_grid,
914 if (h_clutter.nrows() > 1 || h_clutter.ncols() > 1) {
915 chk_atm_surface(
"GIN h_clutter", h_clutter, atmosphere_dim, lat_grid, lon_grid);
916 hclutterm = h_clutter;
918 hclutterm.resize(nlat, nlon);
919 hclutterm = h_clutter(0,0);
923 particle_bulkprop_field.resize(2, np, nlat, nlon);
924 particle_bulkprop_field = 0;
925 particle_bulkprop_names = scat_species;
930 const Numeric extrap_fac = 100;
937#pragma omp parallel for if (!arts_omp_in_parallel() && nlat + nlon > 2) \
938 firstprivate(wss) collapse(2)
939 for (Index ilat = 0; ilat < nlat; ilat++) {
940 for (Index ilon = 0; ilon < nlon; ilon++) {
941 if (fail_msg.
nelem() != 0)
continue;
944 if (
max(dBZe(joker, ilat, ilon)) > dbze_noise) {
945 Numeric dbze_corr_abs = 0;
946 Numeric dbze_corr_hyd = 0;
947 Numeric k_part_above = 0, k_abs_above = 0;
951 for (Index ip = np - 1; ip >= 0; ip--) {
953 if (z_field(ip, ilat, ilon) >= z_surface(ilat, ilon) +
954 hclutterm(ilat, ilon)) {
957 dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
960 const Numeric t = t_field(ip, ilat, ilon);
963 Index phase = -1, it = -1;
964 if (dbze > dbze_noise) {
966 phase = t >= t_phase ? 0 : 1;
974 else if (t >= invtable[phase].get_numeric_grid(
981 it = gp.
fd[0] < 0.5 ? gp.
idx : gp.
idx + 1;
988 if (dBZe(ip, ilat, ilon) > dbze_noise && do_atten_hyd &&
989 dbze_corr_hyd < atten_hyd_max) {
993 "Max 20 dB extrapolation allowed. Found a value\n"
994 "exceeding this limit.");
1003 interp(itw, invtable[phase].data(1, joker, it), gp);
1006 0.5 * hfac * (k_part_above + k_this) *
1007 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1012 dbze_corr_hyd =
min(dbze_corr_hyd, atten_hyd_max);
1014 k_part_above = k_this;
1023 PropagationMatrix propmat;
1024 ArrayOfPropagationMatrix partial_dummy;
1025 StokesVector nlte_dummy;
1026 ArrayOfStokesVector partial_nlte_dummy;
1040 t_field(ip, ilat, ilon),
1041 rtp_nlte_local_dummy,
1042 Vector{vmr_field(joker, ip, ilat, ilon)},
1043 propmat_clearsky_agenda);
1044 k_this = propmat.Kjj()[0];
1047 0.5 * hfac * (k_abs_above + k_this) *
1048 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1052 k_abs_above = k_this;
1057 if (dBZe(ip, ilat, ilon) > dbze_noise) {
1059 dbze = dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
1064 "Max 20 dB extrapolation allowed. Found a value\n"
1065 "exceeding this limit.");
1075 (
interp(itw, invtable[phase].data(0, joker, it), gp)));
1079 else if (wc > wc_clip)
1081 particle_bulkprop_field(phase, ip, ilat, ilon) = wc;
1087 particle_bulkprop_field(0, ip, ilat, ilon) =
1088 particle_bulkprop_field(0, ip + 1, ilat, ilon);
1089 particle_bulkprop_field(1, ip, ilat, ilon) =
1090 particle_bulkprop_field(1, ip + 1, ilat, ilon);
1096 }
catch (
const std::exception& e) {
1098 os <<
"Run-time error at ilat: " << ilat <<
" ilon: " << ilon <<
":\n"
1100#pragma omp critical(particle_bulkpropRadarOnionPeeling_latlon)
1101 fail_msg.push_back(os.str());
1106 if (fail_msg.
nelem()) {
1108 for (
const auto& msg : fail_msg) os << msg <<
'\n';
1117 const Vector& f_grid,
1123 const Index& i_species,
1124 const Vector& dbze_grid,
1125 const Vector& t_grid,
1126 const Numeric& wc_min,
1127 const Numeric& wc_max,
1128 const Numeric& ze_tref,
1133 const Index nss = scat_data.
nelem();
1134 const Index ndb = dbze_grid.nelem();
1135 const Index nt = t_grid.nelem();
1136 const Index iss = i_species;
1140 "First element in T_grid of scattering species is "
1141 "below 220 K.\nThat seems to be too low for liquid.\n"
1142 "Please note that the onion peeling function expects "
1143 "rain\nas the first scattering element.");
1145 "This method requires that *f_grid* has length 1.");
1147 "*i_species* must either be 0 or 1.");
1149 "*scat_data* must contain data for exactly two "
1150 "scattering species.");
1152 "*scat_data* and *scat_species* are inconsistent in size.");
1154 "*scat_data* and *scat_meta* are inconsistent in size.");
1156 "*scat_data* and *scat_meta* have inconsistent sizes.");
1158 "*scat_data* should just contain one frequency.");
1160 "The PSD applied must be of 1-moment type.");
1163 if (invtable.empty())
1165 invtable[iss].set_name(
"Radar inversion table created by *RadarOnionPeelingTableCalc*");
1166 invtable[iss].resize(2, ndb, nt);
1167 invtable[iss].data = 0;
1168 invtable[iss].set_grid_name(0,
"Radiative properties");
1169 invtable[iss].set_grid(0, {
"Log of water content",
"Extinction"});
1170 invtable[iss].set_grid_name(1,
"Radar reflectivity");
1171 invtable[iss].set_grid(1, dbze_grid);
1172 invtable[iss].set_grid_name(2,
"Temperature");
1173 invtable[iss].set_grid(2, t_grid);
1176 const Index nse = scat_data[iss].
nelem();
1177 Matrix
b(nt,nse), e(nt,nse);
1181 for (Index i=0; i<nse; i++) {
1183 "So far only TRO is handled by this method.");
1184 const Index ib = scat_data[iss][i].za_grid.
nelem() - 1;
1186 "All za_grid in scat_data must end with 180.\n"
1187 "This is not the case for scat_data[", iss,
"][",
1189 gridpos(gp, scat_data[iss][i].T_grid, t_grid, 1);
1191 interp(e(joker,i), itw, scat_data[iss][i].ext_mat_data(0,joker,0,0,0), gp);
1192 interp(
b(joker,i), itw, scat_data[iss][i].pha_mat_data(0,joker,ib,0,0,0,0), gp);
1199 const Index nwc = wc_grid.nelem();
1203 Tensor3 D(2, nwc, nt, 0);
1205 const Vector& pnd_agenda_input_t = t_grid;
1206 Matrix pnd_agenda_input(nt, 1);
1210 ze_cfac(cfac, f_grid, ze_tref, k2);
1212 for (Index
w=0;
w<nwc;
w++) {
1214 pnd_agenda_input = wc_grid[
w];
1216 Tensor3 dpnd_data_dx;
1223 pnd_agenda_array_input_names[iss],
1228 for (Index t=0; t<nt; t++) {
1229 for (Index i=0; i<nse; i++) {
1231 "A negative back-scattering found for scat_species ", iss,
1232 ",\ntemperature ", t_grid[t],
"K and element ", i)
1234 "A negative extinction found for scat_species ", iss,
1235 ",\ntemperature ", t_grid[t],
"K and element ", i)
1237 "A negative PSD value found for scat_species ", iss,
1238 ",\ntemperature ", t_grid[t],
"K and ", wc_grid[
w],
1240 D(0,
w,t) += pnd_data(t,i) *
b(t,i);
1241 D(1,
w,t) += pnd_data(t,i) * e(t,i);
1244 D(0,
w,t) = 10 * log10(cfac[0] * D(0,
w,t));
1253 transform(wc_log, log10, wc_grid);
1254 for (Index t=0; t<nt; t++) {
1255 if (!is_increasing(D(0,joker,t))) {
1256 for (Index
w=0;
w<nwc;
w++) {
1257 cout << wc_grid[
w] <<
" " << D(0,
w,t) << endl;
1260 "A case found of non-increasing dBZe.\n"
1261 "Found for scat_species ", iss,
" and ", t_grid[t],
"K.")
1263 if (D(0,0,t) > dbze_grid[0]) {
1264 for (Index
w=0;
w<nwc;
w++) {
1265 cout << wc_grid[
w] <<
" " << D(0,
w,t) << endl;
1268 "A case found where start of dbze_grid not covered.\n"
1269 "Found for scat_species ", iss,
" and ", t_grid[t],
"K.")
1271 if (D(0,nwc-1,t) < dbze_grid[ndb-1]) {
1272 for (Index
w=0;
w<nwc;
w++) {
1273 cout << wc_grid[
w] <<
" " << D(0,
w,t) << endl;
1276 "A case found where end of dbze_grid not covered.\n"
1277 "Found for scat_species ", iss,
" and ", t_grid[t],
"K.")
1280 gridpos(gp, D(0,joker,t), dbze_grid);
1282 interp(invtable[iss].data(0,joker,t), itw, wc_log, gp);
1283 interp(invtable[iss].data(1,joker,t), itw, D(1,joker,t), gp);
Array< Index > ArrayOfIndex
An array of Index.
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
The global header file for ARTS.
Constants of physical expressions as constexpr.
Header file for helper functions for OpenMP.
void pnd_agenda_arrayExecute(Workspace &ws, Matrix &pnd_data, Tensor3 &dpnd_data_dx, const Index agenda_array_index, const Vector &pnd_agenda_input_t, const Matrix &pnd_agenda_input, const ArrayOfString &pnd_agenda_input_names, const ArrayOfString &dpnd_data_dx_names, const ArrayOfAgenda &input_agenda_array)
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 iy_radar_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, Vector &geo_pos, const ArrayOfString &iy_aux_vars, const Index iy_id, const Index cloudbox_on, const Index jacobian_do, const Vector &rte_pos, const Vector &rte_los, const Agenda &input_agenda)
Index nelem() const ARTS_NOEXCEPT
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
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.
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
void VectorLogSpace(Vector &x, const Numeric &start, const Numeric &stop, const Numeric &step, const Verbosity &verbosity)
WORKSPACE METHOD: VectorLogSpace.
void yRadar(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, const Index &atmgeom_checked, const Index &atmfields_checked, const String &iy_unit_radar, const ArrayOfString &iy_aux_vars, const Index &stokes_dim, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &sensor_checked, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &iy_radar_agenda, const ArrayOfArrayOfIndex &instrument_pol_array, const Vector &range_bins, const Numeric &ze_tref, const Numeric &k2, const Numeric &dbze_min, const Verbosity &)
WORKSPACE METHOD: yRadar.
const Index GFIELD3_T_GRID
const Index GFIELD3_DB_GRID
constexpr Numeric SPEED_OF_LIGHT
void RadarOnionPeelingTableCalc(Workspace &ws, ArrayOfGriddedField3 &invtable, const Vector &f_grid, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfArrayOfScatteringMetaData &scat_meta, const ArrayOfAgenda &pnd_agenda_array, const ArrayOfArrayOfString &pnd_agenda_array_input_names, const Index &i_species, const Vector &dbze_grid, const Vector &t_grid, const Numeric &wc_min, const Numeric &wc_max, const Numeric &ze_tref, const Numeric &k2, const Verbosity &verbosity)
WORKSPACE METHOD: RadarOnionPeelingTableCalc.
void particle_bulkpropRadarOnionPeeling(Workspace &ws, Tensor4 &particle_bulkprop_field, ArrayOfString &particle_bulkprop_names, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Matrix &z_surface, const Index &atmfields_checked, const Index &atmgeom_checked, const Vector &f_grid, const Agenda &propmat_clearsky_agenda, const ArrayOfString &scat_species, const ArrayOfGriddedField3 &invtable, const Matrix &incangles, const Tensor3 &dBZe, const Numeric &dbze_noise, const Matrix &h_clutter, const Index &fill_clutter, const Numeric &t_phase, const Numeric &wc_max, const Numeric &wc_clip, const Index &do_atten_abs, const Index &do_atten_hyd, const Numeric &atten_hyd_scaling, const Numeric &atten_hyd_max, const Verbosity &)
WORKSPACE METHOD: particle_bulkpropRadarOnionPeeling.
constexpr Numeric LOG10_EULER_NUMBER
void iyRadarSingleScat(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_pnd, Matrix &ppvar_f, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &scat_data_checked, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Matrix &iy_transmitter, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Numeric &rte_alonglos_v, const Index &trans_in_jacobian, const Numeric &pext_scaling, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyRadarSingleScat.
Numeric last(ConstVectorView x)
last
Declarations having to do with the four output streams.
constexpr Numeric log10_euler
Ten's logarithm of Euler's number.
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...
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
void pha_mat_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstMatrixView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstVectorView &ppath_line_of_sight, const ConstVectorView &ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
void mirror_los(Vector &los_mirrored, const ConstVectorView &los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index <e, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &ppath_f_grid, const Vector &ppath_magnetic_field, const Vector &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const Vector &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index <e, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Declaration of functions in rte.cc.
void integration_bin_by_vecmult(VectorView h, ConstVectorView x_g_in, const Numeric &limit1, const Numeric &limit2)
integration_bin_by_vecmult
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Sensor modelling functions.
Structure to store a grid position.
std::array< Numeric, 2 > fd
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
Numeric end_lstep
The distance between end pos and the first position in pos.
Vector lstep
The length between ppath points.
Vector ngroup
The group index of refraction.
Radiation Vector for Stokes dimension 1-4.
Eigen::VectorXd Vec(size_t i) const
Return Vector at position by copy.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
ArrayOfTransmissionMatrix bulk_backscatter(const ConstTensor5View &Pe, const ConstMatrixView &pnd)
Bulk back-scattering.
ArrayOfArrayOfTransmissionMatrix bulk_backscatter_derivative(const ConstTensor5View &Pe, const ArrayOfMatrix &dpnd_dx)
Derivatives of bulk back-scattering
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void set_backscatter_radiation_vector(ArrayOfRadiationVector &I, ArrayOfArrayOfArrayOfRadiationVector &dI, const RadiationVector &I_incoming, const ArrayOfTransmissionMatrix &T, const ArrayOfTransmissionMatrix &PiTf, const ArrayOfTransmissionMatrix &PiTr, const ArrayOfTransmissionMatrix &Z, const ArrayOfArrayOfTransmissionMatrix &dT1, const ArrayOfArrayOfTransmissionMatrix &dT2, const ArrayOfArrayOfTransmissionMatrix &dZ, const BackscatterSolver solver)
Set the backscatter radiation vector.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
@ CommutativeTransmission
Array< ArrayOfRadiationVector > ArrayOfArrayOfRadiationVector
Array< RadiationVector > ArrayOfRadiationVector
Array< TransmissionMatrix > ArrayOfTransmissionMatrix