ARTS 2.5.11 (git: 6827797f)
m_cloudradar.cc
Go to the documentation of this file.
1
9/*===========================================================================
10 === External declarations
11 ===========================================================================*/
12
13#include <cmath>
14#include <stdexcept>
15
16
17#include "arts.h"
18#include "arts_constants.h"
19#include "arts_omp.h"
20#include "auto_md.h"
21#include "logic.h"
22#include "matpack_eigen.h"
23#include "messages.h"
24#include "montecarlo.h"
25#include "propagationmatrix.h"
26#include "rte.h"
27#include "sensor.h"
28
29inline constexpr Numeric PI=Constant::pi;
32
33// Index of grids inside radar inversion tables
34const Index GFIELD3_DB_GRID = 1;
35const Index GFIELD3_T_GRID = 2;
36
37
38/* Workspace method: Doxygen documentation will be auto-generated */
40 Matrix& iy,
41 ArrayOfMatrix& iy_aux,
42 ArrayOfTensor3& diy_dx,
43 Vector& ppvar_p,
44 Vector& ppvar_t,
45 Matrix& ppvar_vmr,
46 Matrix& ppvar_wind,
47 Matrix& ppvar_mag,
48 Matrix& ppvar_pnd,
49 Matrix& ppvar_f,
50 const Index& stokes_dim,
51 const Vector& f_grid,
52 const Index& atmosphere_dim,
53 const Vector& p_grid,
54 const Tensor3& t_field,
55 const EnergyLevelMap& nlte_field,
56 const Tensor4& vmr_field,
57 const ArrayOfArrayOfSpeciesTag& abs_species,
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,
65 const ArrayOfIndex& cloudbox_limits,
66 const Tensor4& pnd_field,
67 const ArrayOfTensor4& dpnd_field_dx,
68 const ArrayOfString& scat_species,
70 const Index& scat_data_checked,
71 const ArrayOfString& iy_aux_vars,
72 const Index& jacobian_do,
73 const ArrayOfRetrievalQuantity& jacobian_quantities,
74 const Ppath& ppath,
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]]) {
83 // Init Jacobian quantities?
84 const Index j_analytical_do = jacobian_do ?
85 do_analytical_jacobian<1>(jacobian_quantities) : 0;
86
87 // Some basic sizes
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();
94
95 // Radiative background index
96 const Index rbi = ppath_what_background(ppath);
97
98 // Checks of input
99 // Throw error if unsupported features are requested
100 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
101 "ppath.background is invalid. Check your "
102 "calculation of *ppath*?");
103 ARTS_USER_ERROR_IF (rbi == 3 || rbi == 4,
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 ).");
107 if (cloudbox_on) {
108 ARTS_USER_ERROR_IF (scat_data_checked != 1,
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.");
114 }
115 if (jacobian_do) {
116 // FIXME: These needs to be tested properly
117 ARTS_USER_ERROR ( "Jacobian calculations *iyActiveSingleScat* need "
118 "revision before safe to use.");
119
120 ARTS_USER_ERROR_IF (dpnd_field_dx.nelem() != jacobian_quantities.nelem(),
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.");
125 }
126 // iy_aux_vars checked below
127 chk_if_in_range("pext_scaling", pext_scaling, 0, 2);
128
129 //Check transmitter input
130 ARTS_USER_ERROR_IF (iy_transmitter.ncols() != ns || iy_transmitter.nrows() != nf,
131 "The size of *iy* returned from *iy_transmitter_agenda* is\n"
132 "not correct:\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++)
136 ARTS_USER_ERROR_IF (iy_transmitter(iv, 0) != 1,
137 "The *iy* returned from *iy_transmitter_agenda* "
138 "must have the value 1 in the first column.");
139
140 // Set diy_dpath if we are doing are doing jacobian calculations
141 ArrayOfTensor3 diy_dpath = j_analytical_do ?
142 get_standard_diy_dpath(jacobian_quantities, np, nf, ns, true) :
143 ArrayOfTensor3(0);
144
145 // Set the species pointers if we are doing jacobian
146 const ArrayOfIndex jac_species_i = j_analytical_do ?
147 get_pointers_for_analytical_species(jacobian_quantities, abs_species) :
148 ArrayOfIndex(0);
149
150 // Start diy_dx out if we are doing the first run and are doing jacobian
151 // calculations
152 diy_dx = j_analytical_do ?
153 get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, true) :
154 ArrayOfTensor3(0);
155
156 // Checks that the scattering species are treated correctly if their
157 // derivatives are needed
158 const ArrayOfIndex jac_scat_i = j_analytical_do ?
159 get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on) :
160 ArrayOfIndex(0);
161
162 // Init iy_aux
163 iy_aux.resize(naux);
164 //
165 Index auxBackScat = -1;
166 Index auxAbSpAtte = -1;
167 Index auxPartAtte = -1;
168 //
169 for (Index i = 0; i < naux; i++) {
170 iy_aux[i].resize(nf * np, ns);
171
172 if (iy_aux_vars[i] == "Radiative background") {
173 iy_aux[i] = 0;
174 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
175 } else if( iy_aux_vars[i] == "Backscattering" ) {
176 iy_aux[i] = 0;
177 auxBackScat = i;
178 } else if( iy_aux_vars[i] == "Abs species extinction" ) {
179 iy_aux[i] = 0;
180 auxAbSpAtte = i;
181 } else if( iy_aux_vars[i] == "Particle extinction" ) {
182 iy_aux[i] = 0;
183 auxPartAtte = i;
184 } else {
186 "The only allowed strings in *iy_aux_vars* are:\n"
187 " \"Radiative background\"\n"
188 " \"Backscattering\"\n"
189 // " \"Optical depth\"\n"
190 // " \"Particle extinction\"\n"
191 "but you have selected: \"", iy_aux_vars[i], "\"")
192 }
193 }
194
195 // Get atmospheric and radiative variables along the propagation path
196 //
197 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
201 //
203 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(np,
205 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(np,
207
208 // Consider to make Pe an ArrayOfArrayOfPropagationMatrix!
209 Tensor5 Pe(ne, np, nf, ns, ns, 0);
210
211 // Various
212 ArrayOfMatrix ppvar_dpnd_dx(0);
213 ArrayOfIndex clear2cloudy;
214 // Matrix scalar_ext(np,nf,0); // Only used for iy_aux
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);
220 Index ptype;
221 EnergyLevelMap ppvar_nlte;
222
223 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
224 ppvar_p.resize(0);
225 ppvar_t.resize(0);
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);
231 } else {
232 // Basic atmospheric variables
233 get_ppath_atmvars(ppvar_p,
234 ppvar_t,
235 ppvar_nlte,
236 ppvar_vmr,
237 ppvar_wind,
238 ppvar_mag,
239 ppath,
240 atmosphere_dim,
241 p_grid,
242 t_field,
243 nlte_field,
244 vmr_field,
245 wind_u_field,
246 wind_v_field,
247 wind_w_field,
248 mag_u_field,
249 mag_v_field,
250 mag_w_field);
251
253 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
254
255 // pnd_field
256 if (cloudbox_on)
257 get_ppath_cloudvars(clear2cloudy,
258 ppvar_pnd,
259 ppvar_dpnd_dx,
260 ppath,
261 atmosphere_dim,
262 cloudbox_limits,
263 pnd_field,
264 dpnd_field_dx);
265 else {
266 clear2cloudy.resize(np);
267 for (Index ip = 0; ip < np; ip++)
268 clear2cloudy[ip] = -1;
269 }
270
271 // Size radiative variables always used
272 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
273 StokesVector a(nf, ns), S(nf, ns);
274 ArrayOfIndex lte(np);
275
276 // Init variables only used if transmission part of jacobian
277 Vector dB_dT(0);
278 ArrayOfPropagationMatrix dK_this_dx(0), dK_past_dx(0), dKp_dx(0);
279 ArrayOfStokesVector da_dx(0), dS_dx(0);
280
281 // HSE variables
282 Index temperature_derivative_position = -1;
283 bool do_hse = false;
284
285 if (trans_in_jacobian && j_analytical_do) {
286 dK_this_dx.resize(nq);
287 dK_past_dx.resize(nq);
288 dKp_dx.resize(nq);
289 da_dx.resize(nq);
290 dS_dx.resize(nq);
291 dB_dT.resize(nf);
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";
301 })
302 }
303
304 // Loop ppath points and determine radiative properties
305 for (Index ip = 0; ip < np; ip++) {
307 K_this,
308 S,
309 lte[ip],
310 dK_this_dx,
311 dS_dx,
312 propmat_clearsky_agenda,
313 jacobian_quantities,
314 Vector{ppvar_f(joker, ip)},
315 Vector{ppvar_mag(joker, ip)},
316 Vector{ppath.los(ip, joker)},
317 ppvar_nlte[ip],
318 Vector{ppvar_vmr(joker, ip)},
319 ppvar_t[ip],
320 ppvar_p[ip],
321 trans_in_jacobian && j_analytical_do);
322
323 if (trans_in_jacobian && j_analytical_do)
325 dK_this_dx,
326 dS_dx,
327 jacobian_quantities,
328 ppvar_f(joker, ip),
329 ppath.los(ip, joker),
330 lte[ip],
331 atmosphere_dim,
332 trans_in_jacobian && j_analytical_do);
333
334 if (clear2cloudy[ip] + 1) {
336 Kp,
337 da_dx,
338 dKp_dx,
339 jacobian_quantities,
340 ppvar_pnd(joker, Range(ip, 1)),
341 ppvar_dpnd_dx,
342 ip,
343 scat_data,
344 ppath.los(ip, joker),
345 ppvar_t[Range(ip, 1)],
346 atmosphere_dim,
347 trans_in_jacobian && jacobian_do);
348
349 if (abs(pext_scaling - 1) > 1e-6) {
350 Kp *= pext_scaling;
351 if (trans_in_jacobian && j_analytical_do) {
352 FOR_ANALYTICAL_JACOBIANS_DO(dKp_dx[iq] *= pext_scaling;)
353 }
354 }
355
356 // Some iy_aux quantities
357 if (auxAbSpAtte >= 0) {
358 for (Index iv = 0; iv < nf; iv++) {
359 iy_aux[auxAbSpAtte](iv*np+ip, joker) = K_this.Kjj()[iv];
360 }
361 }
362 if (auxPartAtte >= 0) {
363 for (Index iv = 0; iv < nf; iv++) {
364 iy_aux[auxPartAtte](iv*np+ip, joker) = Kp.Kjj()[iv];
365 }
366 }
367
368 K_this += Kp;
369
370 if (trans_in_jacobian && j_analytical_do)
371 FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] += dKp_dx[iq];);
372
373 // Get back-scattering per particle, where relevant
374 {
375 // Direction of outgoing scattered radiation (which is reverse to LOS).
376 Vector los_sca;
377 mirror_los(los_sca, ppath.los(ip, joker), atmosphere_dim);
378 mlos_sca(0, joker) = los_sca;
379
380 // Obtain a length-2 vector for incoming direction
381 Vector los_inc;
382 if (atmosphere_dim == 3) {
383 los_inc = ppath.los(ip, joker);
384 } else { // Mirror back to get a correct 3D LOS
385 mirror_los(los_inc, los_sca, 3);
386 }
387 mlos_inc(0, joker) = los_inc;
388
389 Index i_se_flat = 0;
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++) {
392 // determine whether we have some valid pnd for this
393 // scatelem (in pnd or dpnd)
394 Index val_pnd = 0;
395 if (ppvar_pnd(i_se_flat, ip) != 0)
396 val_pnd = 1;
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) {
401 val_pnd = 1;
402 break;
403 }
404 if (val_pnd) {
405 pha_mat_1ScatElem(pha_mat_1se,
406 ptype,
407 t_ok,
408 scat_data[i_ss][i_se],
409 Vector{ppvar_t[Range(ip, 1)]},
410 mlos_sca,
411 mlos_inc,
412 0,
413 t_interp_order);
414 if (t_ok[0] not_eq 0)
415 if (duplicate_freqs)
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);
419 else
420 Pe(i_se_flat, ip, joker, joker, joker) =
421 pha_mat_1se(joker, 0, 0, 0, joker, joker);
422 else {
424 "Interpolation error for (flat-array) scattering"
425 " element #", i_se_flat, "\n"
426 "at location/temperature point #", ip, "\n")
427 }
428 }
429 i_se_flat++;
430 } // for i_ss
431
432 } // local scope
433 } // clear2cloudy
434
435 if (ip not_eq 0) {
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;
440 stepwise_transmission(lyr_tra[ip],
441 dlyr_tra_above[ip],
442 dlyr_tra_below[ip],
443 K_past,
444 K_this,
445 dK_past_dx,
446 dK_this_dx,
447 ppath.lstep[ip - 1],
448 dr_dT_past,
449 dr_dT_this,
450 temperature_derivative_position);
451 }
452
453 swap(K_past, K_this);
454 swap(dK_past_dx, dK_this_dx);
455 }
456 }
457
458 const ArrayOfTransmissionMatrix tot_tra_forward =
460 const ArrayOfTransmissionMatrix tot_tra_reflect =
462 const ArrayOfTransmissionMatrix reflect_matrix =
463 bulk_backscatter(Pe, ppvar_pnd);
464 const ArrayOfArrayOfTransmissionMatrix dreflect_matrix =
465 bulk_backscatter_derivative(Pe, ppvar_dpnd_dx);
466
467 lvl_rad[0] = iy_transmitter;
468 RadiationVector rad_inc = RadiationVector(nf, ns);
469 rad_inc = iy_transmitter;
471 dlvl_rad,
472 rad_inc,
473 lyr_tra,
474 tot_tra_forward,
475 tot_tra_reflect,
476 reflect_matrix,
477 dlyr_tra_above,
478 dlyr_tra_below,
479 dreflect_matrix,
481
482
483 // Fill iy and diy_dpath
484 //
485 iy.resize(nf * np, ns); // iv*np + ip is the desired output order...
486 //
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) {
492 FOR_ANALYTICAL_JACOBIANS_DO(for (Index ip2 = 0; ip2 < np; ip2++)
493 diy_dpath[iq](ip, iv*np+ip2, is) =
494 dlvl_rad[ip][ip2][iq](iv, is););
495 }
496 }
497 }
498 }
499
500 // Finalize analytical Jacobian
501 if (j_analytical_do) {
502 const Index iy_agenda_call1 = 1;
503 const Tensor3 iy_transmittance(0, 0, 0);
504
506 diy_dx,
507 diy_dpath,
508 ns,
509 nf,
510 np,
511 atmosphere_dim,
512 ppath,
513 ppvar_p,
514 ppvar_t,
515 ppvar_vmr,
516 iy_agenda_call1,
517 iy_transmittance,
518 water_p_eq_agenda,
519 jacobian_quantities,
520 jac_species_i);
521 }
522
523 // Remaining part of iy_aux
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);
529 }
530 }
531 }
532}
533
534
535
536/* Workspace method: Doxygen documentation will be auto-generated */
538 Vector& y,
539 Vector& y_f,
540 ArrayOfIndex& y_pol,
541 Matrix& y_pos,
542 Matrix& y_los,
543 ArrayOfVector& y_aux,
544 Matrix& y_geo,
545 Matrix& jacobian,
546 const Index& atmgeom_checked,
547 const Index& atmfields_checked,
548 const String& iy_unit_radar,
549 const ArrayOfString& iy_aux_vars,
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,
558 const ArrayOfRetrievalQuantity& jacobian_quantities,
559 const Agenda& iy_radar_agenda,
560 const ArrayOfArrayOfIndex& instrument_pol_array,
561 const Vector& range_bins,
562 const Numeric& ze_tref,
563 const Numeric& k2,
564 const Numeric& dbze_min,
565 const Verbosity&) {
566 // Important sizes
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();
571
572 //---------------------------------------------------------------------------
573 // Input checks
574 //---------------------------------------------------------------------------
575
576 // Basics
577 //
578 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
579 //
580 ARTS_USER_ERROR_IF (f_grid.empty(), "The frequency grid is empty.");
581 chk_if_increasing("f_grid", f_grid);
582 ARTS_USER_ERROR_IF (f_grid[0] <= 0, "All frequencies in *f_grid* must be > 0.");
583 //
584 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
585 ARTS_USER_ERROR_IF (atmfields_checked != 1,
586 "The atmospheric fields must be flagged to have "
587 "passed a consistency check (atmfields_checked=1).");
588 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
589 "The atmospheric geometry must be flagged to have "
590 "passed a consistency check (atmgeom_checked=1).");
591 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
592 "The cloudbox must be flagged to have "
593 "passed a consistency check (cloudbox_checked=1).");
594 ARTS_USER_ERROR_IF (sensor_checked != 1,
595 "The sensor variables must be flagged to have "
596 "passed a consistency check (sensor_checked=1).");
597
598 // Method specific variables
599 bool is_z = max(range_bins) > 1;
600 ARTS_USER_ERROR_IF (!is_increasing(range_bins),
601 "The vector *range_bins* must contain strictly "
602 "increasing values.");
603 ARTS_USER_ERROR_IF (!is_z && min(range_bins) < 0,
604 "The vector *range_bins* is not allowed to contain "
605 "negative times.");
606 ARTS_USER_ERROR_IF (instrument_pol_array.nelem() != nf,
607 "The main length of *instrument_pol_array* must match "
608 "the number of frequencies.");
609
610 // iy_unit_radar and variables to handle conversion to Ze and dBZe
611 Vector cfac(nf, 1.0);
612 Numeric ze_min = 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);
620 } else {
622 "For this method, *iy_unit_radar* must be set to \"1\", "
623 "\"Ze\" or \"dBZe\".");
624 }
625
626 //---------------------------------------------------------------------------
627 // Various initializations
628 //---------------------------------------------------------------------------
629
630 // Variables to handle conversion from Stokes to instrument_pol
631 ArrayOfIndex npolcum(nf + 1);
632 npolcum[0] = 0;
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;
637 W[i].resize(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);
641 }
642 }
643
644 // Resize and init *y* and *y_XXX*
645 //
646 const Index ntot = npos * npolcum[nf] * nbins;
647 y.resize(ntot);
648 y = NAN;
649 y_f.resize(ntot);
650 y_pol.resize(ntot);
651 y_pos.resize(ntot, sensor_pos.ncols());
652 y_los.resize(ntot, sensor_los.ncols());
653 y_geo.resize(ntot, 5);
654 y_geo = NAN; // Will be replaced if relavant data are provided
655
656 // y_aux
657 y_aux.resize(naux);
658 for (Index i = 0; i < naux; i++) {
659 y_aux[i].resize(ntot);
660 y_aux[i] = NAN;
661 }
662
663 // Jacobian variables
664 //
665 Index j_analytical_do = 0;
666 Index njq = 0;
667 ArrayOfArrayOfIndex jacobian_indices;
668 //
669 if (jacobian_do) {
670 bool any_affine;
671 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
672
673 jacobian.resize(ntot,
674 jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
675 jacobian = 0;
676 njq = jacobian_quantities.nelem();
677 //
678 FOR_ANALYTICAL_JACOBIANS_DO(j_analytical_do = 1;)
679 } else {
680 jacobian.resize(0, 0);
681 }
682
683 //---------------------------------------------------------------------------
684 // The calculations
685 //---------------------------------------------------------------------------
686
687 // Loop positions
688 for (Index p = 0; p < npos; p++) {
689 // RT part
690 ArrayOfTensor3 diy_dx;
691 Vector geo_pos;
692 Matrix iy;
693 Ppath ppath;
694 ArrayOfMatrix iy_aux;
695 const Index iy_id = (Index)1e6 * p;
696 //
698 iy,
699 iy_aux,
700 ppath,
701 diy_dx,
702 geo_pos,
703 iy_aux_vars,
704 iy_id,
705 cloudbox_on,
706 jacobian_do,
707 Vector{sensor_pos(p, joker)},
708 Vector{sensor_los(p, joker)},
709 iy_radar_agenda);
710
711 // Check if path and size OK
712 const Index np = ppath.np;
713 ARTS_USER_ERROR_IF (np == 1,
714 "A path consisting of a single point found. "
715 "This is not allowed.");
716 error_if_limb_ppath(ppath);
717 ARTS_USER_ERROR_IF (iy.nrows() != nf * np,
718 "The size of *iy* returned from *iy_radar_agenda* "
719 "is not correct (for this method).");
720
721 // Range of ppath, in altitude or time
722 Vector range(np);
723 if (is_z) {
724 range = ppath.pos(joker, 0);
725 } else { // Calculate round-trip time
726 range[0] = 2 * ppath.end_lstep / SPEED_OF_LIGHT;
727 for (Index i = 1; i < np; i++) {
728 range[i] = range[i - 1] + ppath.lstep[i - 1] *
729 (ppath.ngroup[i - 1] + ppath.ngroup[i]) /
731 }
732 }
733 const Numeric range_end1 = min(range[0], range[np - 1]);
734 const Numeric range_end2 = max(range[0], range[np - 1]);
735
736 // Loop radar bins
737 for (Index b = 0; b < nbins; b++) {
738 if (!(range_bins[b] >= range_end2 || // Otherwise bin totally outside
739 range_bins[b + 1] <= range_end1)) // range of ppath
740 {
741 // Bin limits
742 Numeric blim1 = max(range_bins[b], range_end1);
743 Numeric blim2 = min(range_bins[b + 1], range_end2);
744
745 // Determine weight vector to obtain mean inside bin
746 Vector hbin(np);
747 integration_bin_by_vecmult(hbin, range, blim1, blim2);
748 // The function above handles integration over the bin, while we
749 // want the average, so divide weights with bin width
750 hbin /= (blim2 - blim1);
751
752 for (Index iv = 0; iv < nf; iv++) {
753 // Pick out part of iy for frequency
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);)
759 }
760 ArrayOfMatrix A(naux);
761 for (Index a = 0; a < naux; a++) {
762 A[a] = iy_aux[a](Range(iv * np, np), joker);
763 }
764
765 // Variables to hold data for one freq and one pol
766 Vector refl(np);
767 ArrayOfMatrix drefl(njq);
768 if (j_analytical_do) {
769 FOR_ANALYTICAL_JACOBIANS_DO(drefl[iq].resize(dI[iq].npages(), np);)
770 }
771 ArrayOfVector auxvar(naux);
772 for (Index a = 0; a < naux; a++) {
773 auxvar[a].resize(np);
774 }
775
776 for (Index ip = 0; ip < instrument_pol_array[iv].nelem(); ip++) {
777 // Apply weights on each Stokes element
778 mult(refl, I, W[iv][ip]);
779 if (j_analytical_do) {
780 FOR_ANALYTICAL_JACOBIANS_DO(for (Index k = 0;
781 k < drefl[iq].nrows();
782 k++) {
783 mult(drefl[iq](k, joker), dI[iq](k, joker, joker), W[iv][ip]);
784 })
785 }
786 for (Index a = 0; a < naux; a++) {
787 if (iy_aux_vars[a] == "Backscattering") {
788 mult(auxvar[a], A[a], W[iv][ip]);
789 } else {
790 for (Index j = 0; j < np; j++) {
791 auxvar[a][j] = A[a](j, 0);
792 }
793 }
794 }
795
796 // Apply bin weight vector to get final values.
797 Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) + b;
798 y[iout] = cfac[iv] * (hbin * refl);
799 //
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);
808 }
809 })
810 }
811
812 if (iy_unit_radar == "dBZe") {
813 y[iout] = y[iout] <= ze_min ? dbze_min : 10 * log10(y[iout]);
814 }
815
816 // Same for aux variables
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
822 ? dbze_min
823 : 10 * log10(y_aux[a][iout]);
824 }
825 } else {
826 y_aux[a][iout] = hbin * auxvar[a];
827 }
828 }
829
830 // And geo pos
831 if (geo_pos.nelem()) y_geo(iout, joker) = geo_pos;
832 }
833 } // Frequency
834 }
835 }
836
837 // Other aux variables
838 //
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);
847 }
848 }
849 }
850 }
851}
852
853
854
855/* Workspace method: Doxygen documentation will be auto-generated */
857 Workspace& ws,
858 Tensor4& particle_bulkprop_field,
859 ArrayOfString& particle_bulkprop_names,
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,
872 const ArrayOfString& scat_species,
873 const ArrayOfGriddedField3& invtable,
874 const Matrix& incangles,
875 const Tensor3& dBZe,
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,
886 const Verbosity&)
887{
888 const Index np = t_field.npages();
889 const Index nlat = t_field.nrows();
890 const Index nlon = t_field.ncols();
891
892 // Checks of input
893 ARTS_USER_ERROR_IF (atmfields_checked != 1,
894 "The atmospheric fields must be flagged to have\n"
895 "passed a consistency check (atmfields_checked=1).");
896 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
897 "The atmospheric geometry must be flagged to have\n"
898 "passed a consistency check (atmgeom_checked=1).");
899 ARTS_USER_ERROR_IF (scat_species.nelem() != 2, "Length of *scat_data* must be two.");
900 ARTS_USER_ERROR_IF (invtable.nelem() != 2, "Length of *invtable* must be two.");
901 ARTS_USER_ERROR_IF (dbze_noise < invtable[0].get_numeric_grid(GFIELD3_DB_GRID)[0],
902 "*dbze_noise* not covered by invtable[0]." );
903 ARTS_USER_ERROR_IF (dbze_noise < invtable[1].get_numeric_grid(GFIELD3_DB_GRID)[0],
904 "*dbze_noise* not covered by invtable[1]." );
905 chk_atm_field("GIN reflectivities", dBZe, atmosphere_dim, p_grid,
906 lat_grid, lon_grid);
907 chk_atm_surface("GIN incangles", incangles, atmosphere_dim, lat_grid,
908 lon_grid);
909 chk_if_in_range("atten_hyd_scaling", atten_hyd_scaling, 0, 2);
910
911 // h_clutter needs special attention as it is allowed to have size 1x1.
912 // And to make the code simpler below, we always switch to matrix
913 Matrix hclutterm;
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;
917 } else {
918 hclutterm.resize(nlat, nlon);
919 hclutterm = h_clutter(0,0);
920 }
921
922 // Init output
923 particle_bulkprop_field.resize(2, np, nlat, nlon);
924 particle_bulkprop_field = 0;
925 particle_bulkprop_names = scat_species;
926
927 // We apply huge extrapolation in dbze, to handle values outside the table
928 // range. Extrapolation should be OK dBZe and log(IWC) have close to linear
929 // relationship.
930 const Numeric extrap_fac = 100;
931
932 ArrayOfString fail_msg;
933
935
936 // Loop all profiles
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;
942 try {
943 // Check if any significant reflectivity at all in profile
944 if (max(dBZe(joker, ilat, ilon)) > dbze_noise) {
945 Numeric dbze_corr_abs = 0; // Corrections for 2-way attenuation
946 Numeric dbze_corr_hyd = 0;
947 Numeric k_part_above = 0, k_abs_above = 0;
948 // Take abs below if incangle wrongly given as viewing angle
949 const Numeric hfac = abs(1 / Conversion::cosd(incangles(ilat, ilon)));
950
951 for (Index ip = np - 1; ip >= 0; ip--) {
952 // Above clutter zone
953 if (z_field(ip, ilat, ilon) >= z_surface(ilat, ilon) +
954 hclutterm(ilat, ilon)) {
955 // Local dBZe, roughly corrected with attenuation for previous point
956 Numeric dbze =
957 dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
958
959 // Local temperature
960 const Numeric t = t_field(ip, ilat, ilon);
961
962 // Prepare for interpolation of invtable, when needed
963 Index phase = -1, it = -1; // Dummy values
964 if (dbze > dbze_noise) {
965 // Liquid or ice?
966 phase = t >= t_phase ? 0 : 1;
967
968 // Find closest temperature
969 const Index tlast =
970 invtable[phase].get_numeric_grid(GFIELD3_T_GRID).nelem() -
971 1;
972 if (t <= invtable[phase].get_numeric_grid(GFIELD3_T_GRID)[0])
973 it = 0;
974 else if (t >= invtable[phase].get_numeric_grid(
975 GFIELD3_T_GRID)[tlast])
976 it = tlast;
977 else {
978 GridPos gp;
979 gridpos(
980 gp, invtable[phase].get_numeric_grid(GFIELD3_T_GRID), t);
981 it = gp.fd[0] < 0.5 ? gp.idx : gp.idx + 1;
982 }
983 }
984
985 // Calculate attenuation from previous point
986 if (ip < np - 1) {
987 // Attenuation due particles
988 if (dBZe(ip, ilat, ilon) > dbze_noise && do_atten_hyd &&
989 dbze_corr_hyd < atten_hyd_max) {
990 // Get extinction
991 ARTS_USER_ERROR_IF ( dbze - 20 >
992 last(invtable[phase].get_numeric_grid(GFIELD3_DB_GRID)),
993 "Max 20 dB extrapolation allowed. Found a value\n"
994 "exceeding this limit.");
995 GridPos gp;
996 gridpos(gp,
997 invtable[phase].get_numeric_grid(GFIELD3_DB_GRID),
998 dbze,
999 extrap_fac);
1000 Vector itw(2);
1001 interpweights(itw, gp);
1002 Numeric k_this =
1003 interp(itw, invtable[phase].data(1, joker, it), gp);
1004 // Optical thickness
1005 Numeric tau =
1006 0.5 * hfac * (k_part_above + k_this) *
1007 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1008 // This equals -10*log10(exp(-tau)^2)
1009 dbze_corr_hyd +=
1010 20 * LOG10_EULER_NUMBER * (atten_hyd_scaling * tau);
1011 // Make sure we don't pass max value
1012 dbze_corr_hyd = min(dbze_corr_hyd, atten_hyd_max);
1013 // k_this can now be shifted to be "above value"
1014 k_part_above = k_this;
1015 } else { // To handle case: dBZe(ip,ilat,ilon) <= dbze_noise
1016 k_part_above = 0;
1017 }
1018
1019 // Attenuation due to abs_species
1020 if (do_atten_abs) {
1021 // Calculate local attenuation
1022 Numeric k_this;
1023 PropagationMatrix propmat;
1024 ArrayOfPropagationMatrix partial_dummy;
1025 StokesVector nlte_dummy;
1026 ArrayOfStokesVector partial_nlte_dummy;
1027 EnergyLevelMap rtp_nlte_local_dummy;
1029 wss,
1030 propmat,
1031 nlte_dummy,
1032 partial_dummy,
1033 partial_nlte_dummy,
1035 {},
1036 f_grid,
1037 Vector(3, 0),
1038 Vector(0),
1039 p_grid[ip],
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];
1045 // Optical thickness
1046 Numeric tau =
1047 0.5 * hfac * (k_abs_above + k_this) *
1048 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1049 // This equals -10*log10(exp(-tau)^2)
1050 dbze_corr_abs += 20 * LOG10_EULER_NUMBER * tau;
1051 // k_this can now be shifted to be "above value"
1052 k_abs_above = k_this;
1053 }
1054 }
1055
1056 // Invert
1057 if (dBZe(ip, ilat, ilon) > dbze_noise) {
1058 // Correct reflectivity with (updated) attenuation
1059 dbze = dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
1060
1061 // Interpolate inversion table (table holds log10 of water content)
1062 ARTS_USER_ERROR_IF ( dbze - 20 >
1063 last(invtable[phase].get_numeric_grid(GFIELD3_DB_GRID)),
1064 "Max 20 dB extrapolation allowed. Found a value\n"
1065 "exceeding this limit.");
1066 GridPos gp;
1067 gridpos(gp,
1068 invtable[phase].get_numeric_grid(GFIELD3_DB_GRID),
1069 dbze,
1070 extrap_fac);
1071 Vector itw(2);
1072 interpweights(itw, gp);
1073 Numeric wc =
1074 pow(10.0,
1075 (interp(itw, invtable[phase].data(0, joker, it), gp)));
1076 // Apply max and clip values
1077 if (wc > wc_max)
1078 wc = 0;
1079 else if (wc > wc_clip)
1080 wc = wc_clip;
1081 particle_bulkprop_field(phase, ip, ilat, ilon) = wc;
1082 }
1083
1084 // In clutter zone or below surface
1085 } else {
1086 if (fill_clutter) {
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);
1091 }
1092 }
1093
1094 } // presuure
1095 } // R > 0
1096 } catch (const std::exception& e) {
1097 ostringstream os;
1098 os << "Run-time error at ilat: " << ilat << " ilon: " << ilon << ":\n"
1099 << e.what();
1100#pragma omp critical(particle_bulkpropRadarOnionPeeling_latlon)
1101 fail_msg.push_back(os.str());
1102 }
1103 } // lon
1104 } // lat
1105
1106 if (fail_msg.nelem()) {
1107 ostringstream os;
1108 for (const auto& msg : fail_msg) os << msg << '\n';
1109 ARTS_USER_ERROR(os.str());
1110 }
1111}
1112
1113/* Workspace method: Doxygen documentation will be auto-generated */
1115 Workspace& ws,
1116 ArrayOfGriddedField3& invtable,
1117 const Vector& f_grid,
1118 const ArrayOfString& scat_species,
1119 const ArrayOfArrayOfSingleScatteringData& scat_data,
1120 const ArrayOfArrayOfScatteringMetaData& scat_meta,
1121 const ArrayOfAgenda& pnd_agenda_array,
1122 const ArrayOfArrayOfString& pnd_agenda_array_input_names,
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,
1129 const Numeric& k2,
1130 const Verbosity& verbosity)
1131{
1132 // Some index and sizes
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;
1137
1138 // Check input
1139 ARTS_USER_ERROR_IF (scat_data[0][0].T_grid[0] < 220,
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.");
1144 ARTS_USER_ERROR_IF (f_grid.nelem() != 1,
1145 "This method requires that *f_grid* has length 1.");
1146 ARTS_USER_ERROR_IF (i_species < 0 || i_species > 1,
1147 "*i_species* must either be 0 or 1.");
1148 ARTS_USER_ERROR_IF (nss != 2,
1149 "*scat_data* must contain data for exactly two "
1150 "scattering species.");
1151 ARTS_USER_ERROR_IF (scat_species.nelem() != nss,
1152 "*scat_data* and *scat_species* are inconsistent in size.");
1153 ARTS_USER_ERROR_IF (scat_meta.nelem() != nss,
1154 "*scat_data* and *scat_meta* are inconsistent in size.");
1155 ARTS_USER_ERROR_IF (scat_data[iss].nelem() != scat_meta[iss].nelem(),
1156 "*scat_data* and *scat_meta* have inconsistent sizes.");
1157 ARTS_USER_ERROR_IF (scat_data[iss][0].f_grid.nelem() != 1,
1158 "*scat_data* should just contain one frequency.");
1159 ARTS_USER_ERROR_IF (pnd_agenda_array_input_names[iss].nelem() != 1,
1160 "The PSD applied must be of 1-moment type.");
1161
1162 // Allocate
1163 if (invtable.empty())
1164 invtable.resize(2);
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);
1174
1175 // Determine back-scattering and extinction on t_grid, at each size
1176 const Index nse = scat_data[iss].nelem();
1177 Matrix b(nt,nse), e(nt,nse);
1178 {
1179 Matrix itw(nt,2);
1180 ArrayOfGridPos gp(nt);
1181 for (Index i=0; i<nse; i++) {
1182 ARTS_USER_ERROR_IF (scat_data[iss][i].ptype != PTYPE_TOTAL_RND,
1183 "So far only TRO is handled by this method.");
1184 const Index ib = scat_data[iss][i].za_grid.nelem() - 1;
1185 ARTS_USER_ERROR_IF (scat_data[iss][i].za_grid[ib] < 179.999,
1186 "All za_grid in scat_data must end with 180.\n"
1187 "This is not the case for scat_data[", iss, "][",
1188 i, "] (0-based)")
1189 gridpos(gp, scat_data[iss][i].T_grid, t_grid, 1);
1190 interpweights(itw, gp);
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);
1193 }
1194 }
1195
1196 // Create test grid for water content
1197 Vector wc_grid;
1198 VectorLogSpace(wc_grid, wc_min, wc_max, 0.04, verbosity);
1199 const Index nwc = wc_grid.nelem();
1200
1201 // Calculate dBZe and extinction for wc_grid
1202 //
1203 Tensor3 D(2, nwc, nt, 0);
1204 //
1205 const Vector& pnd_agenda_input_t = t_grid;
1206 Matrix pnd_agenda_input(nt, 1);
1207 ArrayOfString dpnd_data_dx_names(0);
1208 //
1209 Vector cfac(1);
1210 ze_cfac(cfac, f_grid, ze_tref, k2);
1211 //
1212 for (Index w=0; w<nwc; w++) {
1213 // Get pnd
1214 pnd_agenda_input = wc_grid[w];
1215 Matrix pnd_data;
1216 Tensor3 dpnd_data_dx;
1218 pnd_data,
1219 dpnd_data_dx,
1220 iss,
1221 pnd_agenda_input_t,
1222 pnd_agenda_input,
1223 pnd_agenda_array_input_names[iss],
1224 dpnd_data_dx_names,
1225 pnd_agenda_array);
1226
1227 // Sum up to get bulk back-scattering and extinction
1228 for (Index t=0; t<nt; t++) {
1229 for (Index i=0; i<nse; i++) {
1230 ARTS_USER_ERROR_IF (b(t,i) < 0,
1231 "A negative back-scattering found for scat_species ", iss,
1232 ",\ntemperature ", t_grid[t], "K and element ", i)
1233 ARTS_USER_ERROR_IF (e(t,i) < 0,
1234 "A negative extinction found for scat_species ", iss,
1235 ",\ntemperature ", t_grid[t], "K and element ", i)
1236 ARTS_USER_ERROR_IF (pnd_data(t,i) < 0,
1237 "A negative PSD value found for scat_species ", iss,
1238 ",\ntemperature ", t_grid[t], "K and ", wc_grid[w],
1239 " kg/m3")
1240 D(0,w,t) += pnd_data(t,i) * b(t,i);
1241 D(1,w,t) += pnd_data(t,i) * e(t,i);
1242 }
1243 // Convert to dBZe
1244 D(0,w,t) = 10 * log10(cfac[0] * D(0,w,t));
1245 }
1246 }
1247
1248 // Get water content and extinction as a function of dBZe by interpolation
1249 Matrix itw(ndb,2);
1250 ArrayOfGridPos gp(ndb);
1251 // Water content interpolated in log
1252 Vector wc_log(nwc);
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;
1258 }
1260 "A case found of non-increasing dBZe.\n"
1261 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1262 }
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;
1266 }
1268 "A case found where start of dbze_grid not covered.\n"
1269 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1270 }
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;
1274 }
1276 "A case found where end of dbze_grid not covered.\n"
1277 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1278 }
1279 //
1280 gridpos(gp, D(0,joker,t), dbze_grid);
1281 interpweights(itw, gp);
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);
1284 }
1285}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:119
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:224
base max(const Array< base > &x)
Max function.
Definition: array.h:128
base min(const Array< base > &x)
Min function.
Definition: array.h:144
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)
Definition: auto_md.cc:25660
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)
Definition: auto_md.cc:25794
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)
Definition: auto_md.cc:25371
void chk_atm_surface(const String &x_name, const Matrix &x, const Index &dim, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_surface
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:67
void chk_atm_field(const String &x_name, ConstTensor3View x, const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, const bool &chk_lat90)
chk_atm_field (simple fields)
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:92
The Agenda class.
Definition: agenda_class.h:52
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
Workspace class.
Definition: workspace_ng.h:36
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_USER_ERROR(...)
Definition: debug.h:151
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
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.
Definition: jacobian.cc:548
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:589
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:580
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:602
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.
Definition: jacobian.cc:29
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:532
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:529
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
Definition: m_cloudradar.cc:35
const Index GFIELD3_DB_GRID
Definition: m_cloudradar.cc:34
constexpr Numeric SPEED_OF_LIGHT
Definition: m_cloudradar.cc:30
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
Definition: m_cloudradar.cc:31
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.
Definition: m_cloudradar.cc:39
constexpr Numeric PI
Definition: m_cloudradar.cc:29
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:142
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.
@ PTYPE_TOTAL_RND
Definition: optproperties.h:38
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition: ppath.cc:524
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1460
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.
Definition: rte.cc:1188
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
Definition: rte.cc:2267
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.
Definition: rte.cc:828
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.
Definition: rte.cc:1051
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.
Definition: rte.cc:1948
void mirror_los(Vector &los_mirrored, const ConstVectorView &los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:1792
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, 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.
Definition: rte.cc:1116
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 &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition: rte.cc:38
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.
Definition: rte.cc:947
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
Definition: sensor.cc:1419
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition: sensor.cc:978
Sensor modelling functions.
Structure to store a grid position.
Definition: interpolation.h:56
std::array< Numeric, 2 > fd
Definition: interpolation.h:58
Index idx
Definition: interpolation.h:57
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Matrix los
Line-of-sight at each ppath point.
Definition: ppath_struct.h:35
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath_struct.h:45
Vector lstep
The length between ppath points.
Definition: ppath_struct.h:39
Vector ngroup
The group index of refraction.
Definition: ppath_struct.h:49
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.
Array< ArrayOfRadiationVector > ArrayOfArrayOfRadiationVector
Array< RadiationVector > ArrayOfRadiationVector
#define w
#define a
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
#define b