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