ARTS 2.5.4 (git: 31ce4f0e)
m_cloudradar.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 Stefan Buehler <sbuehler(at)ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
28/*===========================================================================
29 === External declarations
30 ===========================================================================*/
31
32#include <cmath>
33#include <stdexcept>
34
35
36#include "arts.h"
37#include "arts_constants.h"
38#include "arts_omp.h"
39#include "auto_md.h"
40#include "logic.h"
41#include "matpack_eigen.h"
42#include "messages.h"
43#include "montecarlo.h"
44#include "propagationmatrix.h"
45#include "rte.h"
46#include "sensor.h"
47
48inline constexpr Numeric PI=Constant::pi;
51
52// Index of grids inside radar inversion tables
55
56
57/* Workspace method: Doxygen documentation will be auto-generated */
59 Matrix& iy,
60 ArrayOfMatrix& iy_aux,
61 ArrayOfTensor3& diy_dx,
62 Vector& ppvar_p,
63 Vector& ppvar_t,
64 Matrix& ppvar_vmr,
65 Matrix& ppvar_wind,
66 Matrix& ppvar_mag,
67 Matrix& ppvar_pnd,
68 Matrix& ppvar_f,
69 const Index& stokes_dim,
70 const Vector& f_grid,
71 const Index& atmosphere_dim,
72 const Vector& p_grid,
73 const Tensor3& t_field,
74 const EnergyLevelMap& nlte_field,
75 const Tensor4& vmr_field,
76 const ArrayOfArrayOfSpeciesTag& abs_species,
77 const Tensor3& wind_u_field,
78 const Tensor3& wind_v_field,
79 const Tensor3& wind_w_field,
80 const Tensor3& mag_u_field,
81 const Tensor3& mag_v_field,
82 const Tensor3& mag_w_field,
83 const Index& cloudbox_on,
84 const ArrayOfIndex& cloudbox_limits,
85 const Tensor4& pnd_field,
86 const ArrayOfTensor4& dpnd_field_dx,
87 const ArrayOfString& scat_species,
89 const Index& scat_data_checked,
90 const ArrayOfString& iy_aux_vars,
91 const Index& jacobian_do,
92 const ArrayOfRetrievalQuantity& jacobian_quantities,
93 const Ppath& ppath,
94 const Matrix& iy_transmitter,
95 const Agenda& propmat_clearsky_agenda,
96 const Agenda& water_p_eq_agenda,
97 const Numeric& rte_alonglos_v,
98 const Index& trans_in_jacobian,
99 const Numeric& pext_scaling,
100 const Index& t_interp_order,
101 const Verbosity& verbosity [[maybe_unused]]) {
102 // Init Jacobian quantities?
103 const Index j_analytical_do = jacobian_do ?
104 do_analytical_jacobian<1>(jacobian_quantities) : 0;
105
106 // Some basic sizes
107 const Index nf = f_grid.nelem();
108 const Index ns = stokes_dim;
109 const Index np = ppath.np;
110 const Index ne = pnd_field.nbooks();
111 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
112 const Index naux = iy_aux_vars.nelem();
113
114 // Radiative background index
115 const Index rbi = ppath_what_background(ppath);
116
117 // Checks of input
118 // Throw error if unsupported features are requested
119 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
120 "ppath.background is invalid. Check your "
121 "calculation of *ppath*?");
122 ARTS_USER_ERROR_IF (rbi == 3 || rbi == 4,
123 "The propagation path ends inside or at boundary of "
124 "the cloudbox.\nFor this method, *ppath* must be "
125 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
126 if (cloudbox_on) {
127 ARTS_USER_ERROR_IF (scat_data_checked != 1,
128 "The scattering data must be flagged to have\n"
129 "passed a consistency check (scat_data_checked=1).");
131 "*pnd_field* and *scat_data* inconsistent regarding total number of"
132 " scattering elements.");
133 }
134 if (jacobian_do) {
135 // FIXME: These needs to be tested properly
136 ARTS_USER_ERROR ( "Jacobian calculations *iyActiveSingleScat* need "
137 "revision before safe to use.");
138
139 ARTS_USER_ERROR_IF (dpnd_field_dx.nelem() != jacobian_quantities.nelem(),
140 "*dpnd_field_dx* not properly initialized:\n"
141 "Number of elements in dpnd_field_dx must be equal number of jacobian"
142 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
143 " is calculated/set.");
144 }
145 // iy_aux_vars checked below
146 chk_if_in_range("pext_scaling", pext_scaling, 0, 2);
147
148 //Check transmitter input
149 ARTS_USER_ERROR_IF (iy_transmitter.ncols() != ns || iy_transmitter.nrows() != nf,
150 "The size of *iy* returned from *iy_transmitter_agenda* is\n"
151 "not correct:\n"
152 " expected size = [", nf, ",", stokes_dim, "]\n"
153 " size of iy = [", iy_transmitter.nrows(), ",", iy_transmitter.ncols(), "]\n")
154 for (Index iv = 0; iv < nf; iv++)
155 ARTS_USER_ERROR_IF (iy_transmitter(iv, 0) != 1,
156 "The *iy* returned from *iy_transmitter_agenda* "
157 "must have the value 1 in the first column.");
158
159 // Set diy_dpath if we are doing are doing jacobian calculations
160 ArrayOfTensor3 diy_dpath = j_analytical_do ?
161 get_standard_diy_dpath(jacobian_quantities, np, nf, ns, true) :
163
164 // Set the species pointers if we are doing jacobian
165 const ArrayOfIndex jac_species_i = j_analytical_do ?
166 get_pointers_for_analytical_species(jacobian_quantities, abs_species) :
167 ArrayOfIndex(0);
168
169 // Start diy_dx out if we are doing the first run and are doing jacobian
170 // calculations
171 diy_dx = j_analytical_do ?
172 get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, true) :
174
175 // Checks that the scattering species are treated correctly if their
176 // derivatives are needed
177 const ArrayOfIndex jac_scat_i = j_analytical_do ?
178 get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on) :
179 ArrayOfIndex(0);
180
181 // Init iy_aux
182 iy_aux.resize(naux);
183 //
184 Index auxBackScat = -1;
185 Index auxAbSpAtte = -1;
186 Index auxPartAtte = -1;
187 //
188 for (Index i = 0; i < naux; i++) {
189 iy_aux[i].resize(nf * np, ns);
190
191 if (iy_aux_vars[i] == "Radiative background") {
192 iy_aux[i] = 0;
193 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
194 } else if( iy_aux_vars[i] == "Backscattering" ) {
195 iy_aux[i] = 0;
196 auxBackScat = i;
197 } else if( iy_aux_vars[i] == "Abs species extinction" ) {
198 iy_aux[i] = 0;
199 auxAbSpAtte = i;
200 } else if( iy_aux_vars[i] == "Particle extinction" ) {
201 iy_aux[i] = 0;
202 auxPartAtte = i;
203 } else {
205 "The only allowed strings in *iy_aux_vars* are:\n"
206 " \"Radiative background\"\n"
207 " \"Backscattering\"\n"
208 // " \"Optical depth\"\n"
209 // " \"Particle extinction\"\n"
210 "but you have selected: \"", iy_aux_vars[i], "\"")
211 }
212 }
213
214 // Get atmospheric and radiative variables along the propagation path
215 //
216 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
220 //
222 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(np,
224 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(np,
226
227 // Consider to make Pe an ArrayOfArrayOfPropagationMatrix!
228 Tensor5 Pe(ne, np, nf, ns, ns, 0);
229
230 // Various
231 ArrayOfMatrix ppvar_dpnd_dx(0);
232 ArrayOfIndex clear2cloudy;
233 // Matrix scalar_ext(np,nf,0); // Only used for iy_aux
234 const Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
235 const Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
236 Tensor6 pha_mat_1se(nf_ssd, 1, 1, 1, ns, ns);
237 Vector t_ok(1), t_array(1);
238 Matrix mlos_sca(1, 2), mlos_inc(1, 2);
239 Index ptype;
240 EnergyLevelMap ppvar_nlte;
241
242 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
243 ppvar_p.resize(0);
244 ppvar_t.resize(0);
245 ppvar_vmr.resize(0, 0);
246 ppvar_wind.resize(0, 0);
247 ppvar_mag.resize(0, 0);
248 ppvar_pnd.resize(0, 0);
249 ppvar_f.resize(0, 0);
250 } else {
251 // Basic atmospheric variables
252 get_ppath_atmvars(ppvar_p,
253 ppvar_t,
254 ppvar_nlte,
255 ppvar_vmr,
256 ppvar_wind,
257 ppvar_mag,
258 ppath,
259 atmosphere_dim,
260 p_grid,
261 t_field,
262 nlte_field,
263 vmr_field,
264 wind_u_field,
265 wind_v_field,
266 wind_w_field,
267 mag_u_field,
268 mag_v_field,
269 mag_w_field);
270
272 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
273
274 // pnd_field
275 if (cloudbox_on)
276 get_ppath_cloudvars(clear2cloudy,
277 ppvar_pnd,
278 ppvar_dpnd_dx,
279 ppath,
280 atmosphere_dim,
281 cloudbox_limits,
282 pnd_field,
283 dpnd_field_dx);
284 else {
285 clear2cloudy.resize(np);
286 for (Index ip = 0; ip < np; ip++)
287 clear2cloudy[ip] = -1;
288 }
289
290 // Size radiative variables always used
291 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
292 StokesVector a(nf, ns), S(nf, ns);
293 ArrayOfIndex lte(np);
294
295 // Init variables only used if transmission part of jacobian
296 Vector dB_dT(0);
297 ArrayOfPropagationMatrix dK_this_dx(0), dK_past_dx(0), dKp_dx(0);
298 ArrayOfStokesVector da_dx(0), dS_dx(0);
299
300 // HSE variables
301 Index temperature_derivative_position = -1;
302 bool do_hse = false;
303
304 if (trans_in_jacobian && j_analytical_do) {
305 dK_this_dx.resize(nq);
306 dK_past_dx.resize(nq);
307 dKp_dx.resize(nq);
308 da_dx.resize(nq);
309 dS_dx.resize(nq);
310 dB_dT.resize(nf);
312 dK_this_dx[iq] = PropagationMatrix(nf, ns);
313 dK_past_dx[iq] = PropagationMatrix(nf, ns);
314 dKp_dx[iq] = PropagationMatrix(nf, ns);
315 da_dx[iq] = StokesVector(nf, ns);
316 dS_dx[iq] = StokesVector(nf, ns);
317 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
318 temperature_derivative_position = iq;
319 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
320 })
321 }
322
323 // Loop ppath points and determine radiative properties
324 for (Index ip = 0; ip < np; ip++) {
326 K_this,
327 S,
328 lte[ip],
329 dK_this_dx,
330 dS_dx,
331 propmat_clearsky_agenda,
332 jacobian_quantities,
333 ppvar_f(joker, ip),
334 ppvar_mag(joker, ip),
335 ppath.los(ip, joker),
336 ppvar_nlte[ip],
337 ppvar_vmr(joker, ip),
338 ppvar_t[ip],
339 ppvar_p[ip],
340 trans_in_jacobian && j_analytical_do);
341
342 if (trans_in_jacobian && j_analytical_do)
344 dK_this_dx,
345 dS_dx,
346 jacobian_quantities,
347 ppvar_f(joker, ip),
348 ppath.los(ip, joker),
349 lte[ip],
350 atmosphere_dim,
351 trans_in_jacobian && j_analytical_do);
352
353 if (clear2cloudy[ip] + 1) {
355 Kp,
356 da_dx,
357 dKp_dx,
358 jacobian_quantities,
359 ppvar_pnd(joker, Range(ip, 1)),
360 ppvar_dpnd_dx,
361 ip,
362 scat_data,
363 ppath.los(ip, joker),
364 ppvar_t[Range(ip, 1)],
365 atmosphere_dim,
366 trans_in_jacobian && jacobian_do);
367
368 if (abs(pext_scaling - 1) > 1e-6) {
369 Kp *= pext_scaling;
370 if (trans_in_jacobian && j_analytical_do) {
371 FOR_ANALYTICAL_JACOBIANS_DO(dKp_dx[iq] *= pext_scaling;)
372 }
373 }
374
375 // Some iy_aux quantities
376 if (auxAbSpAtte >= 0) {
377 for (Index iv = 0; iv < nf; iv++) {
378 iy_aux[auxAbSpAtte](iv*np+ip, joker) = K_this.Kjj()[iv];
379 }
380 }
381 if (auxPartAtte >= 0) {
382 for (Index iv = 0; iv < nf; iv++) {
383 iy_aux[auxPartAtte](iv*np+ip, joker) = Kp.Kjj()[iv];
384 }
385 }
386
387 K_this += Kp;
388
389 if (trans_in_jacobian && j_analytical_do)
390 FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] += dKp_dx[iq];);
391
392 // Get back-scattering per particle, where relevant
393 {
394 // Direction of outgoing scattered radiation (which is reverse to LOS).
395 Vector los_sca;
396 mirror_los(los_sca, ppath.los(ip, joker), atmosphere_dim);
397 mlos_sca(0, joker) = los_sca;
398
399 // Obtain a length-2 vector for incoming direction
400 Vector los_inc;
401 if (atmosphere_dim == 3) {
402 los_inc = ppath.los(ip, joker);
403 } else { // Mirror back to get a correct 3D LOS
404 mirror_los(los_inc, los_sca, 3);
405 }
406 mlos_inc(0, joker) = los_inc;
407
408 Index i_se_flat = 0;
409 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
410 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
411 // determine whether we have some valid pnd for this
412 // scatelem (in pnd or dpnd)
413 Index val_pnd = 0;
414 if (ppvar_pnd(i_se_flat, ip) != 0)
415 val_pnd = 1;
416 else if (j_analytical_do)
417 for (Index iq = 0; iq < nq && !val_pnd; iq++)
418 if (jac_scat_i[iq] >= 0)
419 if (ppvar_dpnd_dx[iq](i_se_flat, ip) != 0) {
420 val_pnd = 1;
421 break;
422 }
423 if (val_pnd) {
424 pha_mat_1ScatElem(pha_mat_1se,
425 ptype,
426 t_ok,
427 scat_data[i_ss][i_se],
428 ppvar_t[Range(ip, 1)],
429 mlos_sca,
430 mlos_inc,
431 0,
432 t_interp_order);
433 if (t_ok[0] not_eq 0)
434 if (duplicate_freqs)
435 for (Index iv = 0; iv < nf; iv++)
436 Pe(i_se_flat, ip, iv, joker, joker) =
437 pha_mat_1se(0, 0, 0, 0, joker, joker);
438 else
439 Pe(i_se_flat, ip, joker, joker, joker) =
440 pha_mat_1se(joker, 0, 0, 0, joker, joker);
441 else {
443 "Interpolation error for (flat-array) scattering"
444 " element #", i_se_flat, "\n"
445 "at location/temperature point #", ip, "\n")
446 }
447 }
448 i_se_flat++;
449 } // for i_ss
450
451 } // local scope
452 } // clear2cloudy
453
454 if (ip not_eq 0) {
455 const Numeric dr_dT_past =
456 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
457 const Numeric dr_dT_this =
458 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
459 stepwise_transmission(lyr_tra[ip],
460 dlyr_tra_above[ip],
461 dlyr_tra_below[ip],
462 K_past,
463 K_this,
464 dK_past_dx,
465 dK_this_dx,
466 ppath.lstep[ip - 1],
467 dr_dT_past,
468 dr_dT_this,
469 temperature_derivative_position);
470 }
471
472 swap(K_past, K_this);
473 swap(dK_past_dx, dK_this_dx);
474 }
475 }
476
477 const ArrayOfTransmissionMatrix tot_tra_forward =
479 const ArrayOfTransmissionMatrix tot_tra_reflect =
481 const ArrayOfTransmissionMatrix reflect_matrix =
482 bulk_backscatter(Pe, ppvar_pnd);
483 const ArrayOfArrayOfTransmissionMatrix dreflect_matrix =
484 bulk_backscatter_derivative(Pe, ppvar_dpnd_dx);
485
486 lvl_rad[0] = iy_transmitter;
487 RadiationVector rad_inc = RadiationVector(nf, ns);
488 rad_inc = iy_transmitter;
490 dlvl_rad,
491 rad_inc,
492 lyr_tra,
493 tot_tra_forward,
494 tot_tra_reflect,
495 reflect_matrix,
496 dlyr_tra_above,
497 dlyr_tra_below,
498 dreflect_matrix,
500
501
502 // Fill iy and diy_dpath
503 //
504 iy.resize(nf * np, ns); // iv*np + ip is the desired output order...
505 //
506 for (Index ip = 0; ip < np; ip++) {
507 for (Index iv = 0; iv < nf; iv++) {
508 for (Index is = 0; is < stokes_dim; is++) {
509 iy(iv*np+ip, is) = lvl_rad[ip](iv, is);
510 if (j_analytical_do) {
511 FOR_ANALYTICAL_JACOBIANS_DO(for (Index ip2 = 0; ip2 < np; ip2++)
512 diy_dpath[iq](ip, iv*np+ip2, is) =
513 dlvl_rad[ip][ip2][iq](iv, is););
514 }
515 }
516 }
517 }
518
519 // Finalize analytical Jacobian
520 if (j_analytical_do) {
521 const Index iy_agenda_call1 = 1;
522 const Tensor3 iy_transmittance(0, 0, 0);
523
525 diy_dx,
526 diy_dpath,
527 ns,
528 nf,
529 np,
530 atmosphere_dim,
531 ppath,
532 ppvar_p,
533 ppvar_t,
534 ppvar_vmr,
535 iy_agenda_call1,
536 iy_transmittance,
537 water_p_eq_agenda,
538 jacobian_quantities,
539 jac_species_i);
540 }
541
542 // Remaining part of iy_aux
543 if (auxBackScat >= 0) {
544 for (Index ip = 0; ip < np; ip++) {
545 for (Index iv = 0; iv < nf; iv++) {
546 VectorView stokesvec = VectorView(iy_aux[auxBackScat](iv*np+ip, joker));
547 matpack::eigen::row_vec(stokesvec).noalias() = reflect_matrix[ip].Mat(iv) * rad_inc.Vec(iv);
548 }
549 }
550 }
551}
552
553
554
555/* Workspace method: Doxygen documentation will be auto-generated */
557 Vector& y,
558 Vector& y_f,
559 ArrayOfIndex& y_pol,
560 Matrix& y_pos,
561 Matrix& y_los,
562 ArrayOfVector& y_aux,
563 Matrix& y_geo,
564 Matrix& jacobian,
565 const Index& atmgeom_checked,
566 const Index& atmfields_checked,
567 const String& iy_unit_radar,
568 const ArrayOfString& iy_aux_vars,
569 const Index& stokes_dim,
570 const Vector& f_grid,
571 const Index& cloudbox_on,
572 const Index& cloudbox_checked,
573 const Matrix& sensor_pos,
574 const Matrix& sensor_los,
575 const Index& sensor_checked,
576 const Index& jacobian_do,
577 const ArrayOfRetrievalQuantity& jacobian_quantities,
578 const Agenda& iy_radar_agenda,
579 const ArrayOfArrayOfIndex& instrument_pol_array,
580 const Vector& range_bins,
581 const Numeric& ze_tref,
582 const Numeric& k2,
583 const Numeric& dbze_min,
584 const Verbosity&) {
585 // Important sizes
586 const Index npos = sensor_pos.nrows();
587 const Index nbins = range_bins.nelem() - 1;
588 const Index nf = f_grid.nelem();
589 const Index naux = iy_aux_vars.nelem();
590
591 //---------------------------------------------------------------------------
592 // Input checks
593 //---------------------------------------------------------------------------
594
595 // Basics
596 //
597 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
598 //
599 ARTS_USER_ERROR_IF (f_grid.empty(), "The frequency grid is empty.");
600 chk_if_increasing("f_grid", f_grid);
601 ARTS_USER_ERROR_IF (f_grid[0] <= 0, "All frequencies in *f_grid* must be > 0.");
602 //
603 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
604 ARTS_USER_ERROR_IF (atmfields_checked != 1,
605 "The atmospheric fields must be flagged to have "
606 "passed a consistency check (atmfields_checked=1).");
607 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
608 "The atmospheric geometry must be flagged to have "
609 "passed a consistency check (atmgeom_checked=1).");
610 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
611 "The cloudbox must be flagged to have "
612 "passed a consistency check (cloudbox_checked=1).");
613 ARTS_USER_ERROR_IF (sensor_checked != 1,
614 "The sensor variables must be flagged to have "
615 "passed a consistency check (sensor_checked=1).");
616
617 // Method specific variables
618 bool is_z = max(range_bins) > 1;
619 ARTS_USER_ERROR_IF (!is_increasing(range_bins),
620 "The vector *range_bins* must contain strictly "
621 "increasing values.");
622 ARTS_USER_ERROR_IF (!is_z && min(range_bins) < 0,
623 "The vector *range_bins* is not allowed to contain "
624 "negative times.");
625 ARTS_USER_ERROR_IF (instrument_pol_array.nelem() != nf,
626 "The main length of *instrument_pol_array* must match "
627 "the number of frequencies.");
628
629 // iy_unit_radar and variables to handle conversion to Ze and dBZe
630 Vector cfac(nf, 1.0);
631 Numeric ze_min = 0;
632 const Numeric jfac = 10 * log10(exp(1.0));
633 if (iy_unit_radar == "1") {
634 } else if (iy_unit_radar == "Ze") {
635 ze_cfac(cfac, f_grid, ze_tref, k2);
636 } else if (iy_unit_radar == "dBZe") {
637 ze_cfac(cfac, f_grid, ze_tref, k2);
638 ze_min = pow(10.0, dbze_min / 10);
639 } else {
641 "For this method, *iy_unit_radar* must be set to \"1\", "
642 "\"Ze\" or \"dBZe\".");
643 }
644
645 //---------------------------------------------------------------------------
646 // Various initializations
647 //---------------------------------------------------------------------------
648
649 // Variables to handle conversion from Stokes to instrument_pol
650 ArrayOfIndex npolcum(nf + 1);
651 npolcum[0] = 0;
653 for (Index i = 0; i < nf; i++) {
654 const Index ni = instrument_pol_array[i].nelem();
655 npolcum[i + 1] = npolcum[i] + ni;
656 W[i].resize(ni);
657 for (Index j = 0; j < ni; j++) {
658 W[i][j].resize(stokes_dim);
659 stokes2pol(W[i][j], stokes_dim, instrument_pol_array[i][j], 0.5);
660 }
661 }
662
663 // Resize and init *y* and *y_XXX*
664 //
665 const Index ntot = npos * npolcum[nf] * nbins;
666 y.resize(ntot);
667 y = NAN;
668 y_f.resize(ntot);
669 y_pol.resize(ntot);
670 y_pos.resize(ntot, sensor_pos.ncols());
671 y_los.resize(ntot, sensor_los.ncols());
672 y_geo.resize(ntot, 5);
673 y_geo = NAN; // Will be replaced if relavant data are provided
674
675 // y_aux
676 y_aux.resize(naux);
677 for (Index i = 0; i < naux; i++) {
678 y_aux[i].resize(ntot);
679 y_aux[i] = NAN;
680 }
681
682 // Jacobian variables
683 //
684 Index j_analytical_do = 0;
685 Index njq = 0;
686 ArrayOfArrayOfIndex jacobian_indices;
687 //
688 if (jacobian_do) {
689 bool any_affine;
690 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
691
692 jacobian.resize(ntot,
693 jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
694 jacobian = 0;
695 njq = jacobian_quantities.nelem();
696 //
697 FOR_ANALYTICAL_JACOBIANS_DO(j_analytical_do = 1;)
698 } else {
699 jacobian.resize(0, 0);
700 }
701
702 //---------------------------------------------------------------------------
703 // The calculations
704 //---------------------------------------------------------------------------
705
706 // Loop positions
707 for (Index p = 0; p < npos; p++) {
708 // RT part
709 ArrayOfTensor3 diy_dx;
710 Vector geo_pos;
711 Matrix iy;
712 Ppath ppath;
713 ArrayOfMatrix iy_aux;
714 const Index iy_id = (Index)1e6 * p;
715 //
717 iy,
718 iy_aux,
719 ppath,
720 diy_dx,
721 geo_pos,
722 iy_aux_vars,
723 iy_id,
724 cloudbox_on,
725 jacobian_do,
726 sensor_pos(p, joker),
727 sensor_los(p, joker),
728 iy_radar_agenda);
729
730 // Check if path and size OK
731 const Index np = ppath.np;
732 ARTS_USER_ERROR_IF (np == 1,
733 "A path consisting of a single point found. "
734 "This is not allowed.");
735 error_if_limb_ppath(ppath);
736 ARTS_USER_ERROR_IF (iy.nrows() != nf * np,
737 "The size of *iy* returned from *iy_radar_agenda* "
738 "is not correct (for this method).");
739
740 // Range of ppath, in altitude or time
741 Vector range(np);
742 if (is_z) {
743 range = ppath.pos(joker, 0);
744 } else { // Calculate round-trip time
745 range[0] = 2 * ppath.end_lstep / SPEED_OF_LIGHT;
746 for (Index i = 1; i < np; i++) {
747 range[i] = range[i - 1] + ppath.lstep[i - 1] *
748 (ppath.ngroup[i - 1] + ppath.ngroup[i]) /
750 }
751 }
752 const Numeric range_end1 = min(range[0], range[np - 1]);
753 const Numeric range_end2 = max(range[0], range[np - 1]);
754
755 // Loop radar bins
756 for (Index b = 0; b < nbins; b++) {
757 if (!(range_bins[b] >= range_end2 || // Otherwise bin totally outside
758 range_bins[b + 1] <= range_end1)) // range of ppath
759 {
760 // Bin limits
761 Numeric blim1 = max(range_bins[b], range_end1);
762 Numeric blim2 = min(range_bins[b + 1], range_end2);
763
764 // Determine weight vector to obtain mean inside bin
765 Vector hbin(np);
766 integration_bin_by_vecmult(hbin, range, blim1, blim2);
767 // The function above handles integration over the bin, while we
768 // want the average, so divide weights with bin width
769 hbin /= (blim2 - blim1);
770
771 for (Index iv = 0; iv < nf; iv++) {
772 // Pick out part of iy for frequency
773 Matrix I = iy(Range(iv * np, np), joker);
774 ArrayOfTensor3 dI(njq);
775 if (j_analytical_do) {
777 dI[iq] = diy_dx[iq](joker, Range(iv * np, np), joker);)
778 }
779 ArrayOfMatrix A(naux);
780 for (Index a = 0; a < naux; a++) {
781 A[a] = iy_aux[a](Range(iv * np, np), joker);
782 }
783
784 // Variables to hold data for one freq and one pol
785 Vector refl(np);
786 ArrayOfMatrix drefl(njq);
787 if (j_analytical_do) {
788 FOR_ANALYTICAL_JACOBIANS_DO(drefl[iq].resize(dI[iq].npages(), np);)
789 }
790 ArrayOfVector auxvar(naux);
791 for (Index a = 0; a < naux; a++) {
792 auxvar[a].resize(np);
793 }
794
795 for (Index ip = 0; ip < instrument_pol_array[iv].nelem(); ip++) {
796 // Apply weights on each Stokes element
797 mult(refl, I, W[iv][ip]);
798 if (j_analytical_do) {
800 k < drefl[iq].nrows();
801 k++) {
802 mult(drefl[iq](k, joker), dI[iq](k, joker, joker), W[iv][ip]);
803 })
804 }
805 for (Index a = 0; a < naux; a++) {
806 if (iy_aux_vars[a] == "Backscattering") {
807 mult(auxvar[a], A[a], W[iv][ip]);
808 } else {
809 for (Index j = 0; j < np; j++) {
810 auxvar[a][j] = A[a](j, 0);
811 }
812 }
813 }
814
815 // Apply bin weight vector to get final values.
816 Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) + b;
817 y[iout] = cfac[iv] * (hbin * refl);
818 //
819 if (j_analytical_do) {
821 for (Index k = 0; k < drefl[iq].nrows(); k++) {
822 jacobian(iout, jacobian_indices[iq][0] + k) =
823 cfac[iv] * (hbin * drefl[iq](k, joker));
824 if (iy_unit_radar == "dBZe") {
825 jacobian(iout, jacobian_indices[iq][0] + k) *=
826 jfac / max(y[iout], ze_min);
827 }
828 })
829 }
830
831 if (iy_unit_radar == "dBZe") {
832 y[iout] = y[iout] <= ze_min ? dbze_min : 10 * log10(y[iout]);
833 }
834
835 // Same for aux variables
836 for (Index a = 0; a < naux; a++) {
837 if (iy_aux_vars[a] == "Backscattering") {
838 y_aux[a][iout] = cfac[iv] * (hbin * auxvar[a]);
839 if (iy_unit_radar == "dBZe") {
840 y_aux[a][iout] = y_aux[a][iout] <= ze_min
841 ? dbze_min
842 : 10 * log10(y_aux[a][iout]);
843 }
844 } else {
845 y_aux[a][iout] = hbin * auxvar[a];
846 }
847 }
848
849 // And geo pos
850 if (geo_pos.nelem()) y_geo(iout, joker) = geo_pos;
851 }
852 } // Frequency
853 }
854 }
855
856 // Other aux variables
857 //
858 for (Index b = 0; b < nbins; b++) {
859 for (Index iv = 0; iv < nf; iv++) {
860 for (Index ip = 0; ip < instrument_pol_array[iv].nelem(); ip++) {
861 const Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) + b;
862 y_f[iout] = f_grid[iv];
863 y_pol[iout] = instrument_pol_array[iv][ip];
864 y_pos(iout, joker) = sensor_pos(p, joker);
865 y_los(iout, joker) = sensor_los(p, joker);
866 }
867 }
868 }
869 }
870}
871
872
873
874/* Workspace method: Doxygen documentation will be auto-generated */
876 Workspace& ws,
877 Tensor4& particle_bulkprop_field,
878 ArrayOfString& particle_bulkprop_names,
879 const Index& atmosphere_dim,
880 const Vector& p_grid,
881 const Vector& lat_grid,
882 const Vector& lon_grid,
883 const Tensor3& t_field,
884 const Tensor3& z_field,
885 const Tensor4& vmr_field,
886 const Matrix& z_surface,
887 const Index& atmfields_checked,
888 const Index& atmgeom_checked,
889 const Vector& f_grid,
890 const Agenda& propmat_clearsky_agenda,
891 const ArrayOfString& scat_species,
892 const ArrayOfGriddedField3& invtable,
893 const Matrix& incangles,
894 const Tensor3& dBZe,
895 const Numeric& dbze_noise,
896 const Matrix& h_clutter,
897 const Index& fill_clutter,
898 const Numeric& t_phase,
899 const Numeric& wc_max,
900 const Numeric& wc_clip,
901 const Index& do_atten_abs,
902 const Index& do_atten_hyd,
903 const Numeric& atten_hyd_scaling,
904 const Numeric& atten_hyd_max,
905 const Verbosity&)
906{
907 const Index np = t_field.npages();
908 const Index nlat = t_field.nrows();
909 const Index nlon = t_field.ncols();
910
911 // Checks of input
912 ARTS_USER_ERROR_IF (atmfields_checked != 1,
913 "The atmospheric fields must be flagged to have\n"
914 "passed a consistency check (atmfields_checked=1).");
915 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
916 "The atmospheric geometry must be flagged to have\n"
917 "passed a consistency check (atmgeom_checked=1).");
918 ARTS_USER_ERROR_IF (scat_species.nelem() != 2, "Length of *scat_data* must be two.");
919 ARTS_USER_ERROR_IF (invtable.nelem() != 2, "Length of *invtable* must be two.");
920 ARTS_USER_ERROR_IF (dbze_noise < invtable[0].get_numeric_grid(GFIELD3_DB_GRID)[0],
921 "*dbze_noise* not covered by invtable[0]." );
922 ARTS_USER_ERROR_IF (dbze_noise < invtable[1].get_numeric_grid(GFIELD3_DB_GRID)[0],
923 "*dbze_noise* not covered by invtable[1]." );
924 chk_atm_field("GIN reflectivities", dBZe, atmosphere_dim, p_grid,
925 lat_grid, lon_grid);
926 chk_atm_surface("GIN incangles", incangles, atmosphere_dim, lat_grid,
927 lon_grid);
928 chk_if_in_range("atten_hyd_scaling", atten_hyd_scaling, 0, 2);
929
930 // h_clutter needs special attention as it is allowed to have size 1x1.
931 // And to make the code simpler below, we always switch to matrix
932 Matrix hclutterm;
933 if (h_clutter.nrows() > 1 || h_clutter.ncols() > 1) {
934 chk_atm_surface("GIN h_clutter", h_clutter, atmosphere_dim, lat_grid, lon_grid);
935 hclutterm = h_clutter;
936 } else {
937 hclutterm.resize(nlat, nlon);
938 hclutterm = h_clutter(0,0);
939 }
940
941 // Init output
942 particle_bulkprop_field.resize(2, np, nlat, nlon);
943 particle_bulkprop_field = 0;
944 particle_bulkprop_names = scat_species;
945
946 // We apply huge extrapolation in dbze, to handle values outside the table
947 // range. Extrapolation should be OK dBZe and log(IWC) have close to linear
948 // relationship.
949 const Numeric extrap_fac = 100;
950
951 ArrayOfString fail_msg;
952
954
955 // Loop all profiles
956#pragma omp parallel for if (!arts_omp_in_parallel() && nlat + nlon > 2) \
957 firstprivate(wss) collapse(2)
958 for (Index ilat = 0; ilat < nlat; ilat++) {
959 for (Index ilon = 0; ilon < nlon; ilon++) {
960 if (fail_msg.nelem() != 0) continue;
961 try {
962 // Check if any significant reflectivity at all in profile
963 if (max(dBZe(joker, ilat, ilon)) > dbze_noise) {
964 Numeric dbze_corr_abs = 0; // Corrections for 2-way attenuation
965 Numeric dbze_corr_hyd = 0;
966 Numeric k_part_above = 0, k_abs_above = 0;
967 // Take abs below if incangle wrongly given as viewing angle
968 const Numeric hfac = abs(1 / Conversion::cosd(incangles(ilat, ilon)));
969
970 for (Index ip = np - 1; ip >= 0; ip--) {
971 // Above clutter zone
972 if (z_field(ip, ilat, ilon) >= z_surface(ilat, ilon) +
973 hclutterm(ilat, ilon)) {
974 // Local dBZe, roughly corrected with attenuation for previous point
975 Numeric dbze =
976 dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
977
978 // Local temperature
979 const Numeric t = t_field(ip, ilat, ilon);
980
981 // Prepare for interpolation of invtable, when needed
982 Index phase = -1, it = -1; // Dummy values
983 if (dbze > dbze_noise) {
984 // Liquid or ice?
985 phase = t >= t_phase ? 0 : 1;
986
987 // Find closest temperature
988 const Index tlast =
989 invtable[phase].get_numeric_grid(GFIELD3_T_GRID).nelem() -
990 1;
991 if (t <= invtable[phase].get_numeric_grid(GFIELD3_T_GRID)[0])
992 it = 0;
993 else if (t >= invtable[phase].get_numeric_grid(
994 GFIELD3_T_GRID)[tlast])
995 it = tlast;
996 else {
997 GridPos gp;
998 gridpos(
999 gp, invtable[phase].get_numeric_grid(GFIELD3_T_GRID), t);
1000 it = gp.fd[0] < 0.5 ? gp.idx : gp.idx + 1;
1001 }
1002 }
1003
1004 // Calculate attenuation from previous point
1005 if (ip < np - 1) {
1006 // Attenuation due particles
1007 if (dBZe(ip, ilat, ilon) > dbze_noise && do_atten_hyd &&
1008 dbze_corr_hyd < atten_hyd_max) {
1009 // Get extinction
1010 ARTS_USER_ERROR_IF ( dbze - 20 >
1011 last(invtable[phase].get_numeric_grid(GFIELD3_DB_GRID)),
1012 "Max 20 dB extrapolation allowed. Found a value\n"
1013 "exceeding this limit.");
1014 GridPos gp;
1015 gridpos(gp,
1016 invtable[phase].get_numeric_grid(GFIELD3_DB_GRID),
1017 dbze,
1018 extrap_fac);
1019 Vector itw(2);
1020 interpweights(itw, gp);
1021 Numeric k_this =
1022 interp(itw, invtable[phase].data(1, joker, it), gp);
1023 // Optical thickness
1024 Numeric tau =
1025 0.5 * hfac * (k_part_above + k_this) *
1026 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1027 // This equals -10*log10(exp(-tau)^2)
1028 dbze_corr_hyd +=
1029 20 * LOG10_EULER_NUMBER * (atten_hyd_scaling * tau);
1030 // Make sure we don't pass max value
1031 dbze_corr_hyd = min(dbze_corr_hyd, atten_hyd_max);
1032 // k_this can now be shifted to be "above value"
1033 k_part_above = k_this;
1034 } else { // To handle case: dBZe(ip,ilat,ilon) <= dbze_noise
1035 k_part_above = 0;
1036 }
1037
1038 // Attenuation due to abs_species
1039 if (do_atten_abs) {
1040 // Calculate local attenuation
1041 Numeric k_this;
1042 PropagationMatrix propmat;
1043 ArrayOfPropagationMatrix partial_dummy;
1044 StokesVector nlte_dummy;
1045 ArrayOfStokesVector partial_nlte_dummy;
1046 EnergyLevelMap rtp_nlte_local_dummy;
1048 wss,
1049 propmat,
1050 nlte_dummy,
1051 partial_dummy,
1052 partial_nlte_dummy,
1054 {},
1055 f_grid,
1056 Vector(3, 0),
1057 Vector(0),
1058 p_grid[ip],
1059 t_field(ip, ilat, ilon),
1060 rtp_nlte_local_dummy,
1061 vmr_field(joker, ip, ilat, ilon),
1062 propmat_clearsky_agenda);
1063 k_this = propmat.Kjj()[0];
1064 // Optical thickness
1065 Numeric tau =
1066 0.5 * hfac * (k_abs_above + k_this) *
1067 (z_field(ip + 1, ilat, ilon) - z_field(ip, ilat, ilon));
1068 // This equals -10*log10(exp(-tau)^2)
1069 dbze_corr_abs += 20 * LOG10_EULER_NUMBER * tau;
1070 // k_this can now be shifted to be "above value"
1071 k_abs_above = k_this;
1072 }
1073 }
1074
1075 // Invert
1076 if (dBZe(ip, ilat, ilon) > dbze_noise) {
1077 // Correct reflectivity with (updated) attenuation
1078 dbze = dBZe(ip, ilat, ilon) + dbze_corr_abs + dbze_corr_hyd;
1079
1080 // Interpolate inversion table (table holds log10 of water content)
1081 ARTS_USER_ERROR_IF ( dbze - 20 >
1082 last(invtable[phase].get_numeric_grid(GFIELD3_DB_GRID)),
1083 "Max 20 dB extrapolation allowed. Found a value\n"
1084 "exceeding this limit.");
1085 GridPos gp;
1086 gridpos(gp,
1087 invtable[phase].get_numeric_grid(GFIELD3_DB_GRID),
1088 dbze,
1089 extrap_fac);
1090 Vector itw(2);
1091 interpweights(itw, gp);
1092 Numeric wc =
1093 pow(10.0,
1094 (interp(itw, invtable[phase].data(0, joker, it), gp)));
1095 // Apply max and clip values
1096 if (wc > wc_max)
1097 wc = 0;
1098 else if (wc > wc_clip)
1099 wc = wc_clip;
1100 particle_bulkprop_field(phase, ip, ilat, ilon) = wc;
1101 }
1102
1103 // In clutter zone or below surface
1104 } else {
1105 if (fill_clutter) {
1106 particle_bulkprop_field(0, ip, ilat, ilon) =
1107 particle_bulkprop_field(0, ip + 1, ilat, ilon);
1108 particle_bulkprop_field(1, ip, ilat, ilon) =
1109 particle_bulkprop_field(1, ip + 1, ilat, ilon);
1110 }
1111 }
1112
1113 } // presuure
1114 } // R > 0
1115 } catch (const std::exception& e) {
1116 ostringstream os;
1117 os << "Run-time error at ilat: " << ilat << " ilon: " << ilon << ":\n"
1118 << e.what();
1119#pragma omp critical(particle_bulkpropRadarOnionPeeling_latlon)
1120 fail_msg.push_back(os.str());
1121 }
1122 } // lon
1123 } // lat
1124
1125 if (fail_msg.nelem()) {
1126 ostringstream os;
1127 for (const auto& msg : fail_msg) os << msg << '\n';
1128 ARTS_USER_ERROR(os.str());
1129 }
1130}
1131
1132/* Workspace method: Doxygen documentation will be auto-generated */
1134 Workspace& ws,
1135 ArrayOfGriddedField3& invtable,
1136 const Vector& f_grid,
1137 const ArrayOfString& scat_species,
1138 const ArrayOfArrayOfSingleScatteringData& scat_data,
1139 const ArrayOfArrayOfScatteringMetaData& scat_meta,
1140 const ArrayOfAgenda& pnd_agenda_array,
1141 const ArrayOfArrayOfString& pnd_agenda_array_input_names,
1142 const Index& i_species,
1143 const Vector& dbze_grid,
1144 const Vector& t_grid,
1145 const Numeric& wc_min,
1146 const Numeric& wc_max,
1147 const Numeric& ze_tref,
1148 const Numeric& k2,
1149 const Verbosity& verbosity)
1150{
1151 // Some index and sizes
1152 const Index nss = scat_data.nelem();
1153 const Index ndb = dbze_grid.nelem();
1154 const Index nt = t_grid.nelem();
1155 const Index iss = i_species;
1156
1157 // Check input
1158 ARTS_USER_ERROR_IF (scat_data[0][0].T_grid[0] < 220,
1159 "First element in T_grid of scattering species is "
1160 "below 220 K.\nThat seems to be too low for liquid.\n"
1161 "Please note that the onion peeling function expects "
1162 "rain\nas the first scattering element.");
1163 ARTS_USER_ERROR_IF (f_grid.nelem() != 1,
1164 "This method requires that *f_grid* has length 1.");
1165 ARTS_USER_ERROR_IF (i_species < 0 || i_species > 1,
1166 "*i_species* must either be 0 or 1.");
1167 ARTS_USER_ERROR_IF (nss != 2,
1168 "*scat_data* must contain data for exactly two "
1169 "scattering species.");
1170 ARTS_USER_ERROR_IF (scat_species.nelem() != nss,
1171 "*scat_data* and *scat_species* are inconsistent in size.");
1172 ARTS_USER_ERROR_IF (scat_meta.nelem() != nss,
1173 "*scat_data* and *scat_meta* are inconsistent in size.");
1174 ARTS_USER_ERROR_IF (scat_data[iss].nelem() != scat_meta[iss].nelem(),
1175 "*scat_data* and *scat_meta* have inconsistent sizes.");
1176 ARTS_USER_ERROR_IF (scat_data[iss][0].f_grid.nelem() != 1,
1177 "*scat_data* should just contain one frequency.");
1178 ARTS_USER_ERROR_IF (pnd_agenda_array_input_names[iss].nelem() != 1,
1179 "The PSD applied must be of 1-moment type.");
1180
1181 // Allocate
1182 if (invtable.empty())
1183 invtable.resize(2);
1184 invtable[iss].set_name("Radar inversion table created by *RadarOnionPeelingTableCalc*");
1185 invtable[iss].resize(2, ndb, nt);
1186 invtable[iss].data = 0;
1187 invtable[iss].set_grid_name(0, "Radiative properties");
1188 invtable[iss].set_grid(0, {"Log of water content","Extinction"});
1189 invtable[iss].set_grid_name(1, "Radar reflectivity");
1190 invtable[iss].set_grid(1, dbze_grid);
1191 invtable[iss].set_grid_name(2, "Temperature");
1192 invtable[iss].set_grid(2, t_grid);
1193
1194 // Determine back-scattering and extinction on t_grid, at each size
1195 const Index nse = scat_data[iss].nelem();
1196 Matrix b(nt,nse), e(nt,nse);
1197 {
1198 Matrix itw(nt,2);
1199 ArrayOfGridPos gp(nt);
1200 for (Index i=0; i<nse; i++) {
1201 ARTS_USER_ERROR_IF (scat_data[iss][i].ptype != PTYPE_TOTAL_RND,
1202 "So far only TRO is handled by this method.");
1203 const Index ib = scat_data[iss][i].za_grid.nelem() - 1;
1204 ARTS_USER_ERROR_IF (scat_data[iss][i].za_grid[ib] < 179.999,
1205 "All za_grid in scat_data must end with 180.\n"
1206 "This is not the case for scat_data[", iss, "][",
1207 i, "] (0-based)")
1208 gridpos(gp, scat_data[iss][i].T_grid, t_grid, 1);
1209 interpweights(itw, gp);
1210 interp(e(joker,i), itw, scat_data[iss][i].ext_mat_data(0,joker,0,0,0), gp);
1211 interp(b(joker,i), itw, scat_data[iss][i].pha_mat_data(0,joker,ib,0,0,0,0), gp);
1212 }
1213 }
1214
1215 // Create test grid for water content
1216 Vector wc_grid;
1217 VectorLogSpace(wc_grid, wc_min, wc_max, 0.04, verbosity);
1218 const Index nwc = wc_grid.nelem();
1219
1220 // Calculate dBZe and extinction for wc_grid
1221 //
1222 Tensor3 D(2, nwc, nt, 0);
1223 //
1224 const Vector& pnd_agenda_input_t = t_grid;
1225 Matrix pnd_agenda_input(nt, 1);
1226 ArrayOfString dpnd_data_dx_names(0);
1227 //
1228 Vector cfac(1);
1229 ze_cfac(cfac, f_grid, ze_tref, k2);
1230 //
1231 for (Index w=0; w<nwc; w++) {
1232 // Get pnd
1233 pnd_agenda_input = wc_grid[w];
1234 Matrix pnd_data;
1235 Tensor3 dpnd_data_dx;
1237 pnd_data,
1238 dpnd_data_dx,
1239 iss,
1240 pnd_agenda_input_t,
1241 pnd_agenda_input,
1242 pnd_agenda_array_input_names[iss],
1243 dpnd_data_dx_names,
1244 pnd_agenda_array);
1245
1246 // Sum up to get bulk back-scattering and extinction
1247 for (Index t=0; t<nt; t++) {
1248 for (Index i=0; i<nse; i++) {
1249 ARTS_USER_ERROR_IF (b(t,i) < 0,
1250 "A negative back-scattering found for scat_species ", iss,
1251 ",\ntemperature ", t_grid[t], "K and element ", i)
1252 ARTS_USER_ERROR_IF (e(t,i) < 0,
1253 "A negative extinction found for scat_species ", iss,
1254 ",\ntemperature ", t_grid[t], "K and element ", i)
1255 ARTS_USER_ERROR_IF (pnd_data(t,i) < 0,
1256 "A negative PSD value found for scat_species ", iss,
1257 ",\ntemperature ", t_grid[t], "K and ", wc_grid[w],
1258 " kg/m3")
1259 D(0,w,t) += pnd_data(t,i) * b(t,i);
1260 D(1,w,t) += pnd_data(t,i) * e(t,i);
1261 }
1262 // Convert to dBZe
1263 D(0,w,t) = 10 * log10(cfac[0] * D(0,w,t));
1264 }
1265 }
1266
1267 // Get water content and extinction as a function of dBZe by interpolation
1268 Matrix itw(ndb,2);
1269 ArrayOfGridPos gp(ndb);
1270 // Water content interpolated in log
1271 Vector wc_log(nwc);
1272 transform(wc_log, log10, wc_grid);
1273 for (Index t=0; t<nt; t++) {
1274 if (!is_increasing(D(0,joker,t))) {
1275 for (Index w=0; w<nwc; w++) {
1276 cout << wc_grid[w] << " " << D(0,w,t) << endl;
1277 }
1279 "A case found of non-increasing dBZe.\n"
1280 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1281 }
1282 if (D(0,0,t) > dbze_grid[0]) {
1283 for (Index w=0; w<nwc; w++) {
1284 cout << wc_grid[w] << " " << D(0,w,t) << endl;
1285 }
1287 "A case found where start of dbze_grid not covered.\n"
1288 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1289 }
1290 if (D(0,nwc-1,t) < dbze_grid[ndb-1]) {
1291 for (Index w=0; w<nwc; w++) {
1292 cout << wc_grid[w] << " " << D(0,w,t) << endl;
1293 }
1295 "A case found where end of dbze_grid not covered.\n"
1296 "Found for scat_species ", iss, " and ", t_grid[t], "K.")
1297 }
1298 //
1299 gridpos(gp, D(0,joker,t), dbze_grid);
1300 interpweights(itw, gp);
1301 interp(invtable[iss].data(0,joker,t), itw, wc_log, gp);
1302 interp(invtable[iss].data(1,joker,t), itw, D(1,joker,t), gp);
1303 }
1304}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:136
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:241
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:25521
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:25655
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:25166
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:86
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:111
The Agenda class.
Definition: agenda_class.h:69
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:525
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:159
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The Tensor5 class.
Definition: matpackV.h:516
The Tensor6 class.
Definition: matpackVI.h:1099
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Workspace class.
Definition: workspace_ng.h:53
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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:565
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:606
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:597
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:619
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:46
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:542
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:539
#define abs(x)
#define ns
#define min(a, b)
#define max(a, b)
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
Header file for logic.cc.
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:54
const Index GFIELD3_DB_GRID
Definition: m_cloudradar.cc:53
constexpr Numeric SPEED_OF_LIGHT
Definition: m_cloudradar.cc:49
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:50
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:58
constexpr Numeric PI
Definition: m_cloudradar.cc:48
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:161
Array< Tensor3 > ArrayOfTensor3
An array of Tensor3.
Definition: matpackIII.h:452
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1431
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
const Joker joker
Declarations having to do with the four output streams.
Index nelem(const Lines &l)
Number of lines.
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 k
Boltzmann constant convenience name [J/K].
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 e
Elementary charge convenience name [C].
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
Temperature
Definition: jacobian.h:59
auto row_vec(matpack::vector auto &&x)
Map the input to a non-owning const-correct Eigen Map representing a row vector.
Definition: matpack_eigen.h:63
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
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:55
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition: ppath.cc:541
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1477
Stuff related to the propagation matrix.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
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:1208
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:2287
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 ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_magnetic_field, const ConstVectorView &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const ConstVectorView &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:1136
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:848
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:1071
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:1968
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:1812
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:58
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:967
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:1480
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition: sensor.cc:1033
Sensor modelling functions.
Structure to store a grid position.
Definition: interpolation.h:73
std::array< Numeric, 2 > fd
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74
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