ARTS 2.5.11 (git: 725533f0)
m_rte.cc
Go to the documentation of this file.
1
9/*===========================================================================
10 === External declarations
11 ===========================================================================*/
12
13#include "arts.h"
14#include "arts_constants.h"
15#include "arts_omp.h"
16#include "auto_md.h"
17#include "check_input.h"
18#include "geodetic.h"
19#include "gridded_fields.h"
20#include "jacobian.h"
21#include "logic.h"
22#include "math_funcs.h"
23#include "messages.h"
24#include "montecarlo.h"
25#include "physics_funcs.h"
26#include "ppath.h"
27#include "rte.h"
28#include "special_interp.h"
29#include "sun.h"
30#include "transmissionmatrix.h"
31#include <cmath>
32#include <stdexcept>
33
34inline constexpr Numeric PI=Constant::pi;
40
41/*===========================================================================
42 === Workspace methods
43 ===========================================================================*/
44
45/* Workspace method: Doxygen documentation will be auto-generated */
46void iyApplyUnit(Matrix& iy,
47 ArrayOfMatrix& iy_aux,
48 const Index& stokes_dim,
49 const Vector& f_grid,
50 const ArrayOfString& iy_aux_vars,
51 const String& iy_unit,
52 const Verbosity&) {
53 ARTS_USER_ERROR_IF (iy_unit == "1",
54 "No need to use this method with *iy_unit* = \"1\".");
55
56 ARTS_USER_ERROR_IF (max(iy(joker, 0)) > 1e-3,
57 "The spectrum matrix *iy* is required to have original radiance\n"
58 "unit, but this seems not to be the case. This as a value above\n"
59 "1e-3 is found in *iy*.")
60
61 // Polarisation index variable
62 ArrayOfIndex i_pol(stokes_dim);
63 for (Index is = 0; is < stokes_dim; is++) {
64 i_pol[is] = is + 1;
65 }
66
67 apply_iy_unit(iy, iy_unit, f_grid, 1, i_pol);
68
69 for (Index i = 0; i < iy_aux_vars.nelem(); i++) {
70 if (iy_aux_vars[i] == "iy" || iy_aux_vars[i] == "Error" ||
71 iy_aux_vars[i] == "Error (uncorrelated)") {
72 apply_iy_unit(iy_aux[i], iy_unit, f_grid, 1, i_pol);
73 }
74 }
75}
76
77/* Workspace method: Doxygen documentation will be auto-generated */
79 Matrix& iy,
80 ArrayOfMatrix& iy_aux,
81 Ppath& ppath,
82 Vector& geo_pos,
83 const Index& atmfields_checked,
84 const Index& atmgeom_checked,
85 const ArrayOfString& iy_aux_vars,
86 const Index& iy_id,
87 const Index& cloudbox_on,
88 const Index& cloudbox_checked,
89 const Index& scat_data_checked,
90 const Vector& f_grid,
91 const EnergyLevelMap& nlte_field,
92 const Vector& rte_pos,
93 const Vector& rte_los,
94 const Vector& rte_pos2,
95 const String& iy_unit,
96 const Agenda& iy_main_agenda,
97 const Verbosity&) {
98 // Basics
99 //
100 ARTS_USER_ERROR_IF (atmfields_checked != 1,
101 "The atmospheric fields must be flagged to have\n"
102 "passed a consistency check (atmfields_checked=1).");
103 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
104 "The atmospheric geometry must be flagged to have\n"
105 "passed a consistency check (atmgeom_checked=1).");
106 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
107 "The cloudbox must be flagged to have\n"
108 "passed a consistency check (cloudbox_checked=1).");
109 if (cloudbox_on)
110 ARTS_USER_ERROR_IF (scat_data_checked != 1,
111 "The scattering data must be flagged to have\n"
112 "passed a consistency check (scat_data_checked=1).");
113
114 // iy_transmittance is just input and can be left empty for first call
115 Tensor3 iy_transmittance(0, 0, 0);
116 ArrayOfTensor3 diy_dx;
117
119 iy,
120 iy_aux,
121 ppath,
122 diy_dx,
123 geo_pos,
124 1,
125 iy_transmittance,
126 iy_aux_vars,
127 iy_id,
128 iy_unit,
129 cloudbox_on,
130 0,
131 f_grid,
132 nlte_field,
133 rte_pos,
134 rte_los,
135 rte_pos2,
136 iy_main_agenda);
137
138 // Don't allow NaNs (should suffice to check first stokes element)
139 for (Index i = 0; i < iy.nrows(); i++) {
140 ARTS_USER_ERROR_IF (std::isnan(iy(i, 0)),
141 "One or several NaNs found in *iy*.");
142 }
143}
144
145/* Workspace method: Doxygen documentation will be auto-generated */
147 Workspace& ws,
148 Matrix& iy,
149 ArrayOfMatrix& iy_aux,
150 ArrayOfTensor3& diy_dx,
151 Vector& ppvar_p,
152 Vector& ppvar_t,
153 EnergyLevelMap& ppvar_nlte,
154 Matrix& ppvar_vmr,
155 Matrix& ppvar_wind,
156 Matrix& ppvar_mag,
157 Matrix& ppvar_f,
158 Tensor3& ppvar_iy,
159 Tensor4& ppvar_trans_cumulat,
160 Tensor4& ppvar_trans_partial,
161 const Index& iy_id,
162 const Index& stokes_dim,
163 const Vector& f_grid,
164 const Index& atmosphere_dim,
165 const Vector& p_grid,
166 const Vector& lat_grid,
167 const Vector& lon_grid,
168 const Tensor3& z_field,
169 const Tensor3& t_field,
170 const EnergyLevelMap& nlte_field,
171 const Tensor4& vmr_field,
172 const ArrayOfArrayOfSpeciesTag& abs_species,
173 const Tensor3& wind_u_field,
174 const Tensor3& wind_v_field,
175 const Tensor3& wind_w_field,
176 const Tensor3& mag_u_field,
177 const Tensor3& mag_v_field,
178 const Tensor3& mag_w_field,
179 const Matrix& z_surface,
180 const Vector& refellipsoid,
181 const Numeric& ppath_lmax,
182 const Numeric& ppath_lraytrace,
183 const Index& cloudbox_on,
184 const Index& gas_scattering_do,
185 const Index& suns_do,
186 const String& iy_unit,
187 const ArrayOfString& iy_aux_vars,
188 const Index& jacobian_do,
189 const ArrayOfRetrievalQuantity& jacobian_quantities,
190 const Ppath& ppath,
191 const Vector& rte_pos2,
192 const ArrayOfSun& suns,
193 const Agenda& propmat_clearsky_agenda,
194 const Agenda& water_p_eq_agenda,
195 const String& rt_integration_option,
196 const Agenda& iy_main_agenda,
197 const Agenda& iy_space_agenda,
198 const Agenda& iy_surface_agenda,
199 const Agenda& iy_cloudbox_agenda,
200 const Agenda& gas_scattering_agenda,
201 const Agenda& ppath_step_agenda,
202 const Index& iy_agenda_call1,
203 const Tensor3& iy_transmittance,
204 const Numeric& rte_alonglos_v,
205 const Tensor3& surface_props_data,
206 const Verbosity& verbosity) {
207 // Init Jacobian quantities?
208 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<2>(jacobian_quantities) : 0;
209
210 // Some basic sizes
211 const Index nf = f_grid.nelem();
212 const Index ns = stokes_dim;
213 const Index np = ppath.np;
214 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
215
216 // Radiative background index
217 const Index rbi = ppath_what_background(ppath);
218
219 // Checks of input
220 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
221 "ppath.background is invalid. Check your "
222 "calculation of *ppath*?");
223 ARTS_USER_ERROR_IF (!iy_agenda_call1 && np == 1 && rbi == 2,
224 "A secondary propagation path starting at the "
225 "surface and is going directly into the surface "
226 "is found. This is not allowed.");
227 ARTS_USER_ERROR_IF (iy_unit != "1" && suns_do,
228 "If suns are present only iy_unit=\"1\" can be used.");
229
230 ARTS_USER_ERROR_IF(jacobian_quantities.nelem() && (suns_do || gas_scattering_do) , R"--(
231Jacobian calculation are not supported when gas scattering or suns are included.
232This feature will be added in a future version.
233)--");
234
235
236 // Set diy_dpath if we are doing are doing jacobian calculations
237 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
238
239 // Set the species pointers if we are doing jacobian
240 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
241
242 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
243 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
244
245 // Init iy_aux and fill where possible
246 const Index naux = iy_aux_vars.nelem();
247// iy_aux.resize(naux);
248 //
249 Index auxOptDepth = -1;
250 Index auxDirectRad = -1;
251 Index auxRadBackGrnd = -1;
252
253 // Allocate
254 iy_aux.resize(iy_aux_vars.nelem());
255
256 // Allocate and set (if possible here) iy_aux
257 Index cnt=-1;
258 for (Index i = 0; i < naux; i++) {
259
260
261 if (iy_aux_vars[i] == "Radiative background"){
262 cnt+=1;
263 iy_aux[cnt].resize(nf, ns);
264 iy_aux[cnt](joker, 0) = (Numeric)min((Index)2, rbi - 1);
265 }
266 else if (iy_aux_vars[i] == "Optical depth"){
267 // we set it further below
268 cnt+=1;
269 auxOptDepth = cnt;
270 iy_aux[cnt].resize(nf, ns);
271 iy_aux[cnt] = 0;
272 }
273 else if (iy_aux_vars[i] == "Direct radiation"){
274 cnt+=1;
275 auxDirectRad = cnt;
276 iy_aux[cnt].resize(nf, ns);
277 iy_aux[cnt] = 0;
278 }
279 else if (iy_aux_vars[i] == "Radiation Background"){
280 cnt+=1;
281 auxRadBackGrnd = cnt;
282 iy_aux[cnt].resize(nf, ns);
283 iy_aux[cnt] = 0;
284 }
285 else {
287 "The only allowed strings in *iy_aux_vars* are:\n"
288 " \"Radiative background\"\n"
289 " \"Optical depth\"\n"
290 " \"Direct radiation\"\n"
291 " \"Radiation Background\"\n"
292 "but you have selected: \"", iy_aux_vars[i], "\"")
293 }
294 }
295
296 // Get atmospheric and radiative variables along the propagation path
297 ppvar_trans_cumulat.resize(np, nf, ns, ns);
298 ppvar_trans_partial.resize(np, nf, ns, ns);
299 ppvar_iy.resize(nf, ns, np);
300
301 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
303 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
304
305 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
307 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
308
314
315 ArrayOfPropagationMatrix K(np, PropagationMatrix(nf, ns));
316 ArrayOfArrayOfPropagationMatrix dK_dx(np);
317 Vector r(np);
318 ArrayOfVector dr_below(np, Vector(nq, 0));
319 ArrayOfVector dr_above(np, Vector(nq, 0));
320
321 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
322 ppvar_p.resize(0);
323 ppvar_t.resize(0);
324 ppvar_vmr.resize(0, 0);
325 ppvar_wind.resize(0, 0);
326 ppvar_mag.resize(0, 0);
327 ppvar_f.resize(0, 0);
328 ppvar_trans_cumulat = 1;
329 } else {
330 // Basic atmospheric variables
331 get_ppath_atmvars(ppvar_p,
332 ppvar_t,
333 ppvar_nlte,
334 ppvar_vmr,
335 ppvar_wind,
336 ppvar_mag,
337 ppath,
338 atmosphere_dim,
339 p_grid,
340 t_field,
341 nlte_field,
342 vmr_field,
343 wind_u_field,
344 wind_v_field,
345 wind_w_field,
346 mag_u_field,
347 mag_v_field,
348 mag_w_field);
349
351 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
352
353 const bool temperature_jacobian =
354 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
355
356 // Size radiative variables always used
357 Vector B(nf);
358 StokesVector a(nf, ns), S(nf, ns);
359
360 // Init variables only used if analytical jacobians done
361 Vector dB_dT(temperature_jacobian ? nf : 0);
362 ArrayOfStokesVector da_dx(nq), dS_dx(nq);
363
364 // HSE variables
365 Index temperature_derivative_position = -1;
366 bool do_hse = false;
367
368 if (j_analytical_do) {
369 for (Index ip = 0; ip < np; ip++) {
370 dK_dx[ip].resize(nq);
371 FOR_ANALYTICAL_JACOBIANS_DO(dK_dx[ip][iq] = PropagationMatrix(nf, ns);)
372 }
374 da_dx[iq] = StokesVector(nf, ns); dS_dx[iq] = StokesVector(nf, ns);
375 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
376 temperature_derivative_position = iq;
377 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
378 })
379 }
380
381
382 //allocate Varibale for direct (sun) source, that is needed outside ppath loop.
383
384 //dummy variables needed for the output and input of
385 // gas_scattering_agenda
386 const Vector in_los_dummy;
387 const Vector out_los_dummy;
388 const ArrayOfIndex cloudbox_limits_dummy;
389 const Tensor4 pnd_field_dummy;
390 const ArrayOfTensor4 dpnd_field_dx_dummy(jacobian_quantities.nelem());
391 const ArrayOfString scat_species_dummy;
392 const ArrayOfArrayOfSingleScatteringData scat_data_dummy;
393
395 ArrayOfString fail_msg;
396 bool do_abort = false;
397
398 // Loop ppath points and determine radiative properties
399#pragma omp parallel for if (!arts_omp_in_parallel()) \
400 firstprivate(wss, a, B, dB_dT, S, da_dx, dS_dx)
401 for (Index ip = 0; ip < np; ip++) {
402 if (do_abort) continue;
403 try {
405 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
406
407 Index lte;
409 K[ip],
410 S,
411 lte,
412 dK_dx[ip],
413 dS_dx,
414 propmat_clearsky_agenda,
415 jacobian_quantities,
416 Vector{ppvar_f(joker, ip)},
417 Vector{ppvar_mag(joker, ip)},
418 Vector{ppath.los(ip, joker)},
419 ppvar_nlte[ip],
420 Vector{ppvar_vmr(joker, ip)},
421 ppvar_t[ip],
422 ppvar_p[ip],
423 j_analytical_do);
424
425 if (j_analytical_do)
427 dS_dx,
428 jacobian_quantities,
429 ppvar_f(joker, ip),
430 ppath.los(ip, joker),
431 lte,
432 atmosphere_dim,
433 j_analytical_do);
434
435
436 RadiationVector scattered_sunlight(nf, ns);
437
438 if (gas_scattering_do) {
439
440 ArrayOfIndex suns_visible(suns.nelem());
441
442 Numeric minP = min(ppvar_p);
443
444 if (suns_do && ppvar_p[ip] > minP) {
445 // We skip the uppermost altitude
446 // level as there can be sometimes issue due
447 // to the finite precision when calculating
448 // the (sun-)ppath. The influence of the
449 // uppermost level in view of scattering
450 // is negligible due to the low density.
451 ArrayOfPpath sun_ppaths(suns.nelem());
452 ArrayOfVector sun_rte_los(suns.nelem(), Vector(2));
453
454 get_sun_ppaths(wss,
455 sun_ppaths,
456 suns_visible,
457 sun_rte_los,
458 Vector{ppath.pos(ip, joker)},
459 suns,
460 f_grid,
461 atmosphere_dim,
462 p_grid,
463 lat_grid,
464 lon_grid,
465 z_field,
466 z_surface,
467 refellipsoid,
468 ppath_lmax,
469 ppath_lraytrace,
470 ppath_step_agenda,
471 verbosity);
472
473 ArrayOfMatrix transmitted_sunlight;
474 ArrayOfArrayOfTensor3 dtransmitted_sunlight_dummy(suns.nelem(),ArrayOfTensor3(jacobian_quantities.nelem()));
475
477 transmitted_sunlight,
478 dtransmitted_sunlight_dummy,
479 stokes_dim,
480 f_grid,
481 atmosphere_dim,
482 p_grid,
483 lat_grid,
484 lon_grid,
485 t_field,
486 nlte_field,
487 vmr_field,
488 abs_species,
489 wind_u_field,
490 wind_v_field,
491 wind_w_field,
492 mag_u_field,
493 mag_v_field,
494 mag_w_field,
495 cloudbox_on,
496 cloudbox_limits_dummy,
497 gas_scattering_do,
498 1,
499 sun_ppaths,
500 suns,
501 suns_visible,
502 refellipsoid,
503 pnd_field_dummy,
504 dpnd_field_dx_dummy,
505 scat_species_dummy,
506 scat_data_dummy,
507 0,
508 jacobian_quantities,
509 propmat_clearsky_agenda,
510 water_p_eq_agenda,
511 gas_scattering_agenda,
512 rte_alonglos_v,
513 verbosity);
514
515 //Loop over the different suns to get the total scattered starlight
516 RadiationVector scattered_sunlight_isun(nf, ns);
517
518 for (Index i_sun = 0; i_sun < suns.nelem(); i_sun++) {
519 if (suns_visible[i_sun]) {
520
521 // here we calculate how much incoming sun radiation is scattered
522 //into the direction of the ppath
524 scattered_sunlight_isun,
525 f_grid,
526 ppvar_p[ip],
527 ppvar_t[ip],
528 Vector{ppvar_vmr(joker, ip)},
529 transmitted_sunlight[i_sun],
530 sun_rte_los[i_sun],
531 Vector{ppath.los(ip, joker)},
532 gas_scattering_agenda);
533
534 scattered_sunlight += scattered_sunlight_isun;
535 }
536 }
537 }
538
539 // Calculate gas scattering extiction
540 PropagationMatrix K_sca;
541 TransmissionMatrix sca_mat_dummy;
542 Vector sca_fct_dummy;
543
545 K_sca,
546 sca_mat_dummy,
547 sca_fct_dummy,
548 f_grid,
549 ppvar_p[ip],
550 ppvar_t[ip],
551 Vector{ppvar_vmr(joker, ip)},
552 in_los_dummy,
553 out_los_dummy,
554 0,
555 gas_scattering_agenda);
556
557 // absorption equals extinction only for the gas absorption part.
558 a = K[ip];
559 K[ip] += K_sca;
560
561 } else {
562 // Here absorption equals extinction
563 a = K[ip];
564 if (j_analytical_do)
565 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_dx[ip][iq];);
566 }
567
568 // scattered_sunlight is changed within
569 // stepwise source.
570 stepwise_source(src_rad[ip],
571 dsrc_rad[ip],
572 scattered_sunlight,
573 K[ip],
574 a,
575 S,
576 dK_dx[ip],
577 da_dx,
578 dS_dx,
579 B,
580 dB_dT,
581 jacobian_quantities,
582 jacobian_do);
583
584
585 } catch (const std::runtime_error& e) {
586 ostringstream os;
587 os << "Runtime-error in source calculation at index " << ip
588 << ": \n";
589 os << e.what();
590#pragma omp critical(iyEmissionStandard_source)
591 {
592 do_abort = true;
593 fail_msg.push_back(os.str());
594 }
595 }
596 }
597
598#pragma omp parallel for if (!arts_omp_in_parallel())
599 for (Index ip = 1; ip < np; ip++) {
600 if (do_abort) continue;
601 try {
602 const Numeric dr_dT_past =
603 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
604 const Numeric dr_dT_this =
605 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
606 stepwise_transmission(lyr_tra[ip],
607 dlyr_tra_above[ip],
608 dlyr_tra_below[ip],
609 K[ip - 1],
610 K[ip],
611 dK_dx[ip - 1],
612 dK_dx[ip],
613 ppath.lstep[ip - 1],
614 dr_dT_past,
615 dr_dT_this,
616 temperature_derivative_position);
617
618 r[ip - 1] = ppath.lstep[ip - 1];
619 if (temperature_derivative_position >= 0){
620 dr_below[ip][temperature_derivative_position] = dr_dT_past;
621 dr_above[ip][temperature_derivative_position] = dr_dT_this;
622 }
623 } catch (const std::runtime_error& e) {
624 ostringstream os;
625 os << "Runtime-error in transmission calculation at index " << ip
626 << ": \n";
627 os << e.what();
628#pragma omp critical(iyEmissionStandard_transmission)
629 {
630 do_abort = true;
631 fail_msg.push_back(os.str());
632 }
633 }
634 }
635
636 ARTS_USER_ERROR_IF (do_abort,
637 "Error messages from failed cases:\n", fail_msg)
638 }
639
640 const ArrayOfTransmissionMatrix tot_tra =
642
643 // iy_transmittance
644 Tensor3 iy_trans_new;
645 if (iy_agenda_call1)
646 iy_trans_new = tot_tra[np - 1];
647 else
648 iy_transmittance_mult(iy_trans_new, iy_transmittance, Tensor3{tot_tra[np - 1]});
649
650 // iy_aux: Optical depth
651 if (auxOptDepth >= 0)
652 for (Index iv = 0; iv < nf; iv++)
653 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
654
655 // Diffuse radiative background
657 iy,
658 diy_dx,
659 iy_trans_new,
660 iy_id,
661 jacobian_do,
662 jacobian_quantities,
663 ppath,
664 rte_pos2,
665 atmosphere_dim,
666 nlte_field,
667 cloudbox_on,
668 stokes_dim,
669 f_grid,
670 iy_unit,
671 surface_props_data,
672 iy_main_agenda,
673 iy_space_agenda,
674 iy_surface_agenda,
675 iy_cloudbox_agenda,
676 iy_agenda_call1,
677 verbosity);
678
679 // Direct radiative background
680 Matrix iy_direct(nf, ns, 0.);
681
682 if (suns_do) {
683 Matrix iy_direct_toa;
684 Tensor3 total_transmission;
685 total_transmission = tot_tra[np - 1];
686 Index stars_visible;
687
688 // Get incoming sun radiation at top of the atmosphere. if sun is not visible
689 // in los, iy_direct_toa will be zero
690 get_sun_background(iy_direct_toa,
691 stars_visible,
692 suns,
693 ppath,
694 f_grid,
695 stokes_dim,
696 atmosphere_dim,
697 refellipsoid);
698
699 if (stars_visible) {
700 for (Index iv = 0; iv < nf; iv++) {
701 mult(iy_direct(iv, joker),
702 total_transmission(iv, joker, joker),
703 iy_direct_toa(iv, joker));
704 }
705
706 //Add sun background
707 iy += iy_direct_toa;
708 }
709 }
710
711 // iy_aux: Direct radiation
712 if (auxDirectRad >= 0)
713 iy_aux[auxDirectRad] = iy_direct;
714
715 // iy_aux: Background Radiation
716 if (auxRadBackGrnd >= 0)
717 iy_aux[auxRadBackGrnd] = iy;
718
719 // set the radiation at the start of the ppath
720 lvl_rad[np - 1] = iy;
721
722 // Radiative transfer calculations
723 if (rt_integration_option == "first order" || rt_integration_option == "default") {
724 for (Index ip = np - 2; ip >= 0; ip--) {
725 lvl_rad[ip] = lvl_rad[ip + 1];
726 update_radiation_vector(lvl_rad[ip],
727 dlvl_rad[ip],
728 dlvl_rad[ip + 1],
729 src_rad[ip],
730 src_rad[ip + 1],
731 dsrc_rad[ip],
732 dsrc_rad[ip + 1],
733 lyr_tra[ip + 1],
734 tot_tra[ip],
735 dlyr_tra_above[ip + 1],
736 dlyr_tra_below[ip + 1],
737 PropagationMatrix(),
738 PropagationMatrix(),
739 ArrayOfPropagationMatrix(),
740 ArrayOfPropagationMatrix(),
741 Numeric(),
742 Vector(),
743 Vector(),
744 0,
745 0,
747
748 }
749 } else if (rt_integration_option == "second order") {
750 for (Index ip = np - 2; ip >= 0; ip--) {
751 lvl_rad[ip] = lvl_rad[ip + 1];
752 update_radiation_vector(lvl_rad[ip],
753 dlvl_rad[ip],
754 dlvl_rad[ip + 1],
755 src_rad[ip],
756 src_rad[ip + 1],
757 dsrc_rad[ip],
758 dsrc_rad[ip + 1],
759 lyr_tra[ip + 1],
760 tot_tra[ip],
761 dlyr_tra_above[ip + 1],
762 dlyr_tra_below[ip + 1],
763 K[ip],
764 K[ip + 1],
765 dK_dx[ip + 1],
766 dK_dx[ip + 1],
767 r[ip],
768 dr_above[ip + 1],
769 dr_below[ip + 1],
770 0,
771 0,
773 }
774 } else {
775 ARTS_USER_ERROR ( "Only allowed choices for *integration order* are "
776 "1 and 2.");
777 }
778
779 // Copy back to ARTS external style
780 iy = lvl_rad[0];
781 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
782 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
783 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
784 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
785 if (j_analytical_do)
786 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
787 dlvl_rad[ip][iq];);
788 }
789
790 // Finalize analytical Jacobians
791 if (j_analytical_do) {
793 diy_dx,
794 diy_dpath,
795 ns,
796 nf,
797 np,
798 atmosphere_dim,
799 ppath,
800 ppvar_p,
801 ppvar_t,
802 ppvar_vmr,
803 iy_agenda_call1,
804 iy_transmittance,
805 water_p_eq_agenda,
806 jacobian_quantities,
807 jac_species_i);
808 }
809
810 // Radiance unit conversions
811 if (iy_agenda_call1) {
813 diy_dx,
814 ppvar_iy,
815 ns,
816 np,
817 f_grid,
818 ppath,
819 jacobian_quantities,
820 j_analytical_do,
821 iy_unit);
822 }
823}
824
825/* Workspace method: Doxygen documentation will be auto-generated */
827 Matrix& iy,
828 ArrayOfMatrix& iy_aux,
829 ArrayOfTensor3& diy_dx,
830 Vector& ppvar_p,
831 Vector& ppvar_t,
832 EnergyLevelMap& ppvar_nlte,
833 Matrix& ppvar_vmr,
834 Matrix& ppvar_wind,
835 Matrix& ppvar_mag,
836 Matrix& ppvar_pnd,
837 Matrix& ppvar_f,
838 Tensor3& ppvar_iy,
839 Tensor4& ppvar_trans_cumulat,
840 Tensor4& ppvar_trans_partial,
841 const Index& iy_id,
842 const Index& stokes_dim,
843 const Vector& f_grid,
844 const Index& atmosphere_dim,
845 const Vector& p_grid,
846 const Tensor3& t_field,
847 const EnergyLevelMap& nlte_field,
848 const Tensor4& vmr_field,
849 const ArrayOfArrayOfSpeciesTag& abs_species,
850 const Tensor3& wind_u_field,
851 const Tensor3& wind_v_field,
852 const Tensor3& wind_w_field,
853 const Tensor3& mag_u_field,
854 const Tensor3& mag_v_field,
855 const Tensor3& mag_w_field,
856 const Index& cloudbox_on,
857 const ArrayOfIndex& cloudbox_limits,
858 const Tensor4& pnd_field,
859 const ArrayOfTensor4& dpnd_field_dx,
860 const ArrayOfString& scat_species,
861 const ArrayOfArrayOfSingleScatteringData& scat_data,
862 const String& iy_unit,
863 const ArrayOfString& iy_aux_vars,
864 const Index& jacobian_do,
865 const ArrayOfRetrievalQuantity& jacobian_quantities,
866 const Agenda& propmat_clearsky_agenda,
867 const Agenda& water_p_eq_agenda,
868 const String& rt_integration_option,
869 const Agenda& iy_main_agenda,
870 const Agenda& iy_space_agenda,
871 const Agenda& iy_surface_agenda,
872 const Agenda& iy_cloudbox_agenda,
873 const Index& iy_agenda_call1,
874 const Tensor3& iy_transmittance,
875 const Ppath& ppath,
876 const Vector& rte_pos2,
877 const Numeric& rte_alonglos_v,
878 const Tensor3& surface_props_data,
879 const Tensor7& cloudbox_field,
880 const Vector& za_grid,
881 const Index& Naa,
882 const Index& t_interp_order,
883 const Verbosity& verbosity) {
884 // If cloudbox off, switch to use clearsky method
885 if (!cloudbox_on) {
886 Tensor4 dummy;
888 iy,
889 iy_aux,
890 diy_dx,
891 ppvar_p,
892 ppvar_t,
893 ppvar_nlte,
894 ppvar_vmr,
895 ppvar_wind,
896 ppvar_mag,
897 ppvar_f,
898 ppvar_iy,
899 ppvar_trans_cumulat,
900 ppvar_trans_partial,
901 iy_id,
902 stokes_dim,
903 f_grid,
904 atmosphere_dim,
905 p_grid,
906 t_field,
907 nlte_field,
908 vmr_field,
909 abs_species,
910 wind_u_field,
911 wind_v_field,
912 wind_w_field,
913 mag_u_field,
914 mag_v_field,
915 mag_w_field,
916 cloudbox_on,
917 iy_unit,
918 iy_aux_vars,
919 jacobian_do,
920 jacobian_quantities,
921 ppath,
922 rte_pos2,
923 propmat_clearsky_agenda,
924 water_p_eq_agenda,
925 rt_integration_option,
926 iy_main_agenda,
927 iy_space_agenda,
928 iy_surface_agenda,
929 iy_cloudbox_agenda,
930 iy_agenda_call1,
931 iy_transmittance,
932 rte_alonglos_v,
933 surface_props_data,
934 verbosity);
935 return;
936 }
937 // Init Jacobian quantities?
938 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<1>(jacobian_quantities) : 0;
939
940 // Some basic sizes
941 const Index nf = f_grid.nelem();
942 const Index ns = stokes_dim;
943 const Index np = ppath.np;
944 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
945
946 // Radiative background index
947 const Index rbi = ppath_what_background(ppath);
948
949 // Throw error if unsupported features are requested
950 if (atmosphere_dim != 1)
951 throw runtime_error(
952 "With cloudbox on, this method handles only 1D calculations.");
953 if (Naa < 3) throw runtime_error("Naa must be > 2.");
954 if (jacobian_do)
955 if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
956 throw runtime_error(
957 "*dpnd_field_dx* not properly initialized:\n"
958 "Number of elements in dpnd_field_dx must be equal number of jacobian"
959 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
960 " is calculated/set.");
961 if (rbi < 1 || rbi > 9)
962 throw runtime_error(
963 "ppath.background is invalid. Check your "
964 "calculation of *ppath*?");
965 if (rbi == 3 || rbi == 4)
966 throw runtime_error(
967 "The propagation path ends inside or at boundary of "
968 "the cloudbox.\nFor this method, *ppath* must be "
969 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
970 // iy_aux_vars checked below
971 // Checks of i_field
972 if (cloudbox_field.ncols() != stokes_dim)
973 throw runtime_error(
974 "Obtained *cloudbox_field* number of Stokes elements inconsistent with "
975 "*stokes_dim*.");
976 if (cloudbox_field.nrows() != 1)
977 throw runtime_error(
978 "Obtained *cloudbox_field* has wrong number of azimuth angles.");
979 if (cloudbox_field.npages() != za_grid.nelem())
980 throw runtime_error(
981 "Obtained *cloudbox_field* number of zenith angles inconsistent with "
982 "*za_grid*.");
983 if (cloudbox_field.nbooks() != 1)
984 throw runtime_error(
985 "Obtained *cloudbox_field* has wrong number of longitude points.");
986 if (cloudbox_field.nshelves() != 1)
987 throw runtime_error(
988 "Obtained *cloudbox_field* has wrong number of latitude points.");
989 if (cloudbox_field.nvitrines() != cloudbox_limits[1] - cloudbox_limits[0] + 1)
990 throw runtime_error(
991 "Obtained *cloudbox_field* number of pressure points inconsistent with "
992 "*cloudbox_limits*.");
993 if (cloudbox_field.nlibraries() != nf)
994 throw runtime_error(
995 "Obtained *cloudbox_field* number of frequency points inconsistent with "
996 "*f_grid*.");
997
998 // Set diy_dpath if we are doing are doing jacobian calculations
999 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
1000
1001 // Set the species pointers if we are doing jacobian
1002 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
1003
1004 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
1005 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
1006
1007 // Checks that the scattering species are treated correctly if their derivatives are needed (we can here discard the Array)
1008 if (j_analytical_do and iy_agenda_call1) get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on);
1009
1010 // Init iy_aux and fill where possible
1011 const Index naux = iy_aux_vars.nelem();
1012 iy_aux.resize(naux);
1013 //
1014 for (Index i = 0; i < naux; i++) {
1015 iy_aux[i].resize(nf, ns);
1016
1017 if (iy_aux_vars[i] == "Optical depth") { /*pass*/
1018 } // Filled below
1019 else if (iy_aux_vars[i] == "Radiative background")
1020 iy_aux[i] = (Numeric)min((Index)2, rbi - 1);
1021 else {
1022 ostringstream os;
1023 os << "The only allowed strings in *iy_aux_vars* are:\n"
1024 << " \"Radiative background\"\n"
1025 << " \"Optical depth\"\n"
1026 << "but you have selected: \"" << iy_aux_vars[i] << "\"";
1027 throw runtime_error(os.str());
1028 }
1029 }
1030
1031 // Get atmospheric and radiative variables along the propagation path
1032 ppvar_trans_cumulat.resize(np, nf, ns, ns);
1033 ppvar_trans_partial.resize(np, nf, ns, ns);
1034 ppvar_iy.resize(nf, ns, np);
1035
1036 ArrayOfTransmissionMatrix lyr_tra(np, TransmissionMatrix(nf, ns));
1037 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
1039 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1040 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
1042 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1043
1044 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
1046 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
1048
1049 ArrayOfIndex clear2cloudy;
1050 //
1051 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
1052 ppvar_p.resize(0);
1053 ppvar_t.resize(0);
1054 ppvar_vmr.resize(0, 0);
1055 ppvar_wind.resize(0, 0);
1056 ppvar_mag.resize(0, 0);
1057 ppvar_f.resize(0, 0);
1058 ppvar_trans_cumulat = 0;
1059 ppvar_trans_partial = 0;
1060 for (Index iv = 0; iv < nf; iv++) {
1061 for (Index is = 0; is < ns; is++) {
1062 ppvar_trans_cumulat(0,iv,is,is) = 1;
1063 ppvar_trans_partial(0,iv,is,is) = 1;
1064 }
1065 }
1066 } else {
1067 // Basic atmospheric variables
1068 get_ppath_atmvars(ppvar_p,
1069 ppvar_t,
1070 ppvar_nlte,
1071 ppvar_vmr,
1072 ppvar_wind,
1073 ppvar_mag,
1074 ppath,
1075 atmosphere_dim,
1076 p_grid,
1077 t_field,
1078 nlte_field,
1079 vmr_field,
1080 wind_u_field,
1081 wind_v_field,
1082 wind_w_field,
1083 mag_u_field,
1084 mag_v_field,
1085 mag_w_field);
1086
1088 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
1089
1090 // here, the cloudbox is on, ie we don't need to check and branch this here
1091 // anymore.
1092 ArrayOfMatrix ppvar_dpnd_dx;
1093 //
1094 get_ppath_cloudvars(clear2cloudy,
1095 ppvar_pnd,
1096 ppvar_dpnd_dx,
1097 ppath,
1098 atmosphere_dim,
1099 cloudbox_limits,
1100 pnd_field,
1101 dpnd_field_dx);
1102
1103 // Size radiative variables always used
1104 Vector B(nf);
1105 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
1106 StokesVector a(nf, ns), S(nf, ns), Sp(nf, ns);
1107 ArrayOfIndex lte(np);
1108 RadiationVector J_add_dummy;
1109 ArrayOfRadiationVector dJ_add_dummy;
1110
1111 // Init variables only used if analytical jacobians done
1112 Vector dB_dT(0);
1113 ArrayOfPropagationMatrix dK_this_dx(nq), dK_past_dx(nq), dKp_dx(nq);
1114 ArrayOfStokesVector da_dx(nq), dS_dx(nq), dSp_dx(nq);
1115
1116 // HSE variables
1117 Index temperature_derivative_position = -1;
1118 bool do_hse = false;
1119
1120 if (j_analytical_do) {
1121 dB_dT.resize(nf);
1123 dK_this_dx[iq] = PropagationMatrix(nf, ns);
1124 dK_past_dx[iq] = PropagationMatrix(nf, ns);
1125 dKp_dx[iq] = PropagationMatrix(nf, ns);
1126 da_dx[iq] = StokesVector(nf, ns);
1127 dS_dx[iq] = StokesVector(nf, ns);
1128 dSp_dx[iq] = StokesVector(nf, ns);
1129 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
1130 temperature_derivative_position = iq;
1131 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
1132 })
1133 }
1134 const bool temperature_jacobian =
1135 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
1136
1137 // Loop ppath points and determine radiative properties
1138 for (Index ip = 0; ip < np; ip++) {
1140 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
1141
1143 K_this,
1144 S,
1145 lte[ip],
1146 dK_this_dx,
1147 dS_dx,
1148 propmat_clearsky_agenda,
1149 jacobian_quantities,
1150 Vector{ppvar_f(joker, ip)},
1151 Vector{ppvar_mag(joker, ip)},
1152 Vector{ppath.los(ip, joker)},
1153 ppvar_nlte[ip],
1154 Vector{ppvar_vmr(joker, ip)},
1155 ppvar_t[ip],
1156 ppvar_p[ip],
1157 j_analytical_do);
1158
1159 if (j_analytical_do)
1161 dS_dx,
1162 jacobian_quantities,
1163 ppvar_f(joker, ip),
1164 ppath.los(ip, joker),
1165 lte[ip],
1166 atmosphere_dim,
1167 j_analytical_do);
1168
1169 if (clear2cloudy[ip] + 1) {
1171 Kp,
1172 da_dx,
1173 dKp_dx,
1174 jacobian_quantities,
1175 ppvar_pnd(joker, Range(ip, 1)),
1176 ppvar_dpnd_dx,
1177 ip,
1178 scat_data,
1179 ppath.los(ip, joker),
1180 ppvar_t[Range(ip, 1)],
1181 atmosphere_dim,
1182 jacobian_do);
1183 a += K_this;
1184 K_this += Kp;
1185
1186 if (j_analytical_do)
1187 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] += dK_this_dx[iq];
1188 dK_this_dx[iq] += dKp_dx[iq];)
1189
1190 Vector aa_grid;
1191 nlinspace(aa_grid, 0, 360, Naa);
1192 //
1194 dSp_dx,
1195 jacobian_quantities,
1196 ppvar_pnd(joker, ip),
1197 ppvar_dpnd_dx,
1198 ip,
1199 scat_data,
1200 cloudbox_field,
1201 za_grid,
1202 aa_grid,
1203 ppath.los(Range(ip, 1), joker),
1204 ppath.gp_p[ip],
1205 Vector{ppvar_t[Range(ip, 1)]},
1206 atmosphere_dim,
1207 jacobian_do,
1208 t_interp_order);
1209 S += Sp;
1210
1211 if (j_analytical_do)
1212 FOR_ANALYTICAL_JACOBIANS_DO(dS_dx[iq] += dSp_dx[iq];)
1213 } else { // no particles present at this level
1214 a = K_this;
1215 if (j_analytical_do)
1216 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_this_dx[iq];)
1217 }
1218
1219 if (ip not_eq 0) {
1220 const Numeric dr_dT_past =
1221 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
1222 const Numeric dr_dT_this =
1223 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
1224 stepwise_transmission(lyr_tra[ip],
1225 dlyr_tra_above[ip],
1226 dlyr_tra_below[ip],
1227 K_past,
1228 K_this,
1229 dK_past_dx,
1230 dK_this_dx,
1231 ppath.lstep[ip - 1],
1232 dr_dT_past,
1233 dr_dT_this,
1234 temperature_derivative_position);
1235 }
1236
1237 stepwise_source(src_rad[ip],
1238 dsrc_rad[ip],
1239 J_add_dummy,
1240 K_this,
1241 a,
1242 S,
1243 dK_this_dx,
1244 da_dx,
1245 dS_dx,
1246 B,
1247 dB_dT,
1248 jacobian_quantities,
1249 jacobian_do);
1250
1251 swap(K_past, K_this);
1252 swap(dK_past_dx, dK_this_dx);
1253 }
1254 }
1255
1256 const ArrayOfTransmissionMatrix tot_tra =
1258
1259 // iy_transmittance
1260 Tensor3 iy_trans_new;
1261 if (iy_agenda_call1)
1262 iy_trans_new = tot_tra[np - 1];
1263 else
1264 iy_transmittance_mult(iy_trans_new, iy_transmittance, Tensor3{tot_tra[np - 1]});
1265
1266 // Copy transmission to iy_aux
1267 for (Index i = 0; i < naux; i++)
1268 if (iy_aux_vars[i] == "Optical depth")
1269 for (Index iv = 0; iv < nf; iv++)
1270 iy_aux[i](iv, joker) = -log(ppvar_trans_cumulat(np - 1, iv, 0, 0));
1271
1272 // Radiative background
1274 iy,
1275 diy_dx,
1276 iy_trans_new,
1277 iy_id,
1278 jacobian_do,
1279 jacobian_quantities,
1280 ppath,
1281 rte_pos2,
1282 atmosphere_dim,
1283 nlte_field,
1284 cloudbox_on,
1285 stokes_dim,
1286 f_grid,
1287 iy_unit,
1288 surface_props_data,
1289 iy_main_agenda,
1290 iy_space_agenda,
1291 iy_surface_agenda,
1292 iy_cloudbox_agenda,
1293 iy_agenda_call1,
1294 verbosity);
1295
1296 lvl_rad[np - 1] = iy;
1297
1298 // Radiative transfer calculations
1299 for (Index ip = np - 2; ip >= 0; ip--) {
1300 lvl_rad[ip] = lvl_rad[ip + 1];
1301 update_radiation_vector(lvl_rad[ip],
1302 dlvl_rad[ip],
1303 dlvl_rad[ip + 1],
1304 src_rad[ip],
1305 src_rad[ip + 1],
1306 dsrc_rad[ip],
1307 dsrc_rad[ip + 1],
1308 lyr_tra[ip + 1],
1309 tot_tra[ip],
1310 dlyr_tra_above[ip + 1],
1311 dlyr_tra_below[ip + 1],
1312 PropagationMatrix(),
1313 PropagationMatrix(),
1314 ArrayOfPropagationMatrix(),
1315 ArrayOfPropagationMatrix(),
1316 Numeric(),
1317 Vector(),
1318 Vector(),
1319 0,
1320 0,
1322 }
1323
1324 // Copy back to ARTS external style
1325 iy = lvl_rad[0];
1326 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
1327 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
1328 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
1329 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
1330 if (j_analytical_do)
1331 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
1332 dlvl_rad[ip][iq];);
1333 }
1334
1335 // Finalize analytical Jacobians
1336 if (j_analytical_do)
1338 diy_dx,
1339 diy_dpath,
1340 ns,
1341 nf,
1342 np,
1343 atmosphere_dim,
1344 ppath,
1345 ppvar_p,
1346 ppvar_t,
1347 ppvar_vmr,
1348 iy_agenda_call1,
1349 iy_transmittance,
1350 water_p_eq_agenda,
1351 jacobian_quantities,
1352 jac_species_i);
1353
1354 // Unit conversions
1355 if (iy_agenda_call1)
1357 diy_dx,
1358 ppvar_iy,
1359 ns,
1360 np,
1361 f_grid,
1362 ppath,
1363 jacobian_quantities,
1364 j_analytical_do,
1365 iy_unit);
1366}
1367
1368/* Workspace method: Doxygen documentation will be auto-generated */
1370 Workspace& ws,
1371 Matrix& iy,
1372 ArrayOfMatrix& iy_aux,
1373 ArrayOfTensor3& diy_dx,
1374 Vector& ppvar_p,
1375 Vector& ppvar_t,
1376 EnergyLevelMap& ppvar_nlte,
1377 Matrix& ppvar_vmr,
1378 Matrix& ppvar_wind,
1379 Matrix& ppvar_mag,
1380 Matrix& ppvar_f,
1381 Tensor3& ppvar_iy,
1382 Tensor4& ppvar_trans_cumulat,
1383 Tensor4& ppvar_trans_partial,
1384 const Index& iy_id,
1385 const Index& stokes_dim,
1386 const Vector& f_grid,
1387 const Index& atmosphere_dim,
1388 const Vector& p_grid,
1389 const Tensor3& t_field,
1390 const EnergyLevelMap& nlte_field,
1391 const Tensor4& vmr_field,
1392 const ArrayOfArrayOfSpeciesTag& abs_species,
1393 const Tensor3& wind_u_field,
1394 const Tensor3& wind_v_field,
1395 const Tensor3& wind_w_field,
1396 const Tensor3& mag_u_field,
1397 const Tensor3& mag_v_field,
1398 const Tensor3& mag_w_field,
1399 const Index& cloudbox_on,
1400 const String& iy_unit,
1401 const ArrayOfString& iy_aux_vars,
1402 const Index& jacobian_do,
1403 const ArrayOfRetrievalQuantity& jacobian_quantities,
1404 const Ppath& ppath,
1405 const Vector& rte_pos2,
1406 const Agenda& propmat_clearsky_agenda,
1407 const Agenda& water_p_eq_agenda,
1408 const String& rt_integration_option,
1409 const Agenda& iy_main_agenda,
1410 const Agenda& iy_space_agenda,
1411 const Agenda& iy_surface_agenda,
1412 const Agenda& iy_cloudbox_agenda,
1413 const Index& iy_agenda_call1,
1414 const Tensor3& iy_transmittance,
1415 const Numeric& rte_alonglos_v,
1416 const Tensor3& surface_props_data,
1417 const Verbosity& verbosity) {
1418 // Init Jacobian quantities?
1419 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<2>(jacobian_quantities) : 0;
1420
1421 // Some basic sizes
1422 const Index nf = f_grid.nelem();
1423 const Index ns = stokes_dim;
1424 const Index np = ppath.np;
1425 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
1426
1427 // Radiative background index
1428 const Index rbi = ppath_what_background(ppath);
1429
1430 // Checks of input
1431 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
1432 "ppath.background is invalid. Check your "
1433 "calculation of *ppath*?");
1434 ARTS_USER_ERROR_IF (!iy_agenda_call1 && np == 1 && rbi == 2,
1435 "A secondary propagation path starting at the "
1436 "surface and is going directly into the surface "
1437 "is found. This is not allowed.");
1438
1439 // Set diy_dpath if we are doing are doing jacobian calculations
1440 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
1441
1442 // Set the species pointers if we are doing jacobian
1443 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
1444
1445 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
1446 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
1447
1448 // Init iy_aux and fill where possible
1449 const Index naux = iy_aux_vars.nelem();
1450 iy_aux.resize(naux);
1451 //
1452 Index auxOptDepth = -1;
1453 //
1454 for (Index i = 0; i < naux; i++) {
1455 iy_aux[i].resize(nf, ns);
1456 iy_aux[i] = 0;
1457
1458 if (iy_aux_vars[i] == "Radiative background")
1459 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
1460 else if (iy_aux_vars[i] == "Optical depth")
1461 auxOptDepth = i;
1462 else {
1464 "The only allowed strings in *iy_aux_vars* are:\n"
1465 " \"Radiative background\"\n"
1466 " \"Optical depth\"\n"
1467 "but you have selected: \"", iy_aux_vars[i], "\"")
1468 }
1469 }
1470
1471 // Get atmospheric and radiative variables along the propagation path
1472 ppvar_trans_cumulat.resize(np, nf, ns, ns);
1473 ppvar_trans_partial.resize(np, nf, ns, ns);
1474 ppvar_iy.resize(nf, ns, np);
1475
1476 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
1478 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1479
1480 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
1482 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1483
1484 ArrayOfTransmissionMatrix lyr_tra(np, TransmissionMatrix(nf, ns));
1485 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
1487 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
1489
1490 ArrayOfPropagationMatrix K(np, PropagationMatrix(nf, ns));
1491 ArrayOfArrayOfPropagationMatrix dK_dx(np);
1492 Vector r(np);
1493 ArrayOfVector dr_below(np, Vector(nq, 0));
1494 ArrayOfVector dr_above(np, Vector(nq, 0));
1495
1496 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
1497 ppvar_p.resize(0);
1498 ppvar_t.resize(0);
1499 ppvar_vmr.resize(0, 0);
1500 ppvar_wind.resize(0, 0);
1501 ppvar_mag.resize(0, 0);
1502 ppvar_f.resize(0, 0);
1503 ppvar_trans_cumulat = 0;
1504 ppvar_trans_partial = 0;
1505 for (Index iv = 0; iv < nf; iv++) {
1506 for (Index is = 0; is < ns; is++) {
1507 ppvar_trans_cumulat(0,iv,is,is) = 1;
1508 ppvar_trans_partial(0,iv,is,is) = 1;
1509 }
1510 }
1511
1512 } else {
1513 // Basic atmospheric variables
1514 get_ppath_atmvars(ppvar_p,
1515 ppvar_t,
1516 ppvar_nlte,
1517 ppvar_vmr,
1518 ppvar_wind,
1519 ppvar_mag,
1520 ppath,
1521 atmosphere_dim,
1522 p_grid,
1523 t_field,
1524 nlte_field,
1525 vmr_field,
1526 wind_u_field,
1527 wind_v_field,
1528 wind_w_field,
1529 mag_u_field,
1530 mag_v_field,
1531 mag_w_field);
1532
1534 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
1535
1536 const bool temperature_jacobian =
1537 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
1538
1539 // Size radiative variables always used
1540 Vector B(nf);
1541 StokesVector a(nf, ns), S(nf, ns);
1542 RadiationVector J_add_dummy;
1543 ArrayOfRadiationVector dJ_add_dummy;
1544
1545 // Init variables only used if analytical jacobians done
1546 Vector dB_dT(temperature_jacobian ? nf : 0);
1547 ArrayOfStokesVector da_dx(nq), dS_dx(nq);
1548
1549 // HSE variables
1550 Index temperature_derivative_position = -1;
1551 bool do_hse = false;
1552
1553 if (j_analytical_do) {
1554 for (Index ip = 0; ip < np; ip++) {
1555 dK_dx[ip].resize(nq);
1556 FOR_ANALYTICAL_JACOBIANS_DO(dK_dx[ip][iq] = PropagationMatrix(nf, ns);)
1557 }
1559 da_dx[iq] = StokesVector(nf, ns); dS_dx[iq] = StokesVector(nf, ns);
1560 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
1561 temperature_derivative_position = iq;
1562 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
1563 })
1564 }
1565
1566 ArrayOfString fail_msg;
1567 bool do_abort = false;
1568
1570
1571 // Loop ppath points and determine radiative properties
1572#pragma omp parallel for if (!arts_omp_in_parallel()) \
1573 firstprivate(wss, a, B, dB_dT, S, da_dx, dS_dx)
1574 for (Index ip = 0; ip < np; ip++) {
1575 if (do_abort) continue;
1576 try {
1578 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
1579
1580 Index lte;
1582 K[ip],
1583 S,
1584 lte,
1585 dK_dx[ip],
1586 dS_dx,
1587 propmat_clearsky_agenda,
1588 jacobian_quantities,
1589 Vector{ppvar_f(joker, ip)},
1590 Vector{ppvar_mag(joker, ip)},
1591 Vector{ppath.los(ip, joker)},
1592 ppvar_nlte[ip],
1593 Vector{ppvar_vmr(joker, ip)},
1594 ppvar_t[ip],
1595 ppvar_p[ip],
1596 j_analytical_do);
1597
1598 if (j_analytical_do)
1600 dS_dx,
1601 jacobian_quantities,
1602 ppvar_f(joker, ip),
1603 ppath.los(ip, joker),
1604 lte,
1605 atmosphere_dim,
1606 j_analytical_do);
1607
1608 // Here absorption equals extinction
1609 a = K[ip];
1610 if (j_analytical_do)
1611 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_dx[ip][iq];);
1612
1613 stepwise_source(src_rad[ip],
1614 dsrc_rad[ip],
1615 J_add_dummy,
1616 K[ip],
1617 a,
1618 S,
1619 dK_dx[ip],
1620 da_dx,
1621 dS_dx,
1622 B,
1623 dB_dT,
1624 jacobian_quantities,
1625 jacobian_do);
1626 } catch (const std::runtime_error& e) {
1627 ostringstream os;
1628 os << "Runtime-error in source calculation at index " << ip
1629 << ": \n";
1630 os << e.what();
1631#pragma omp critical(iyEmissionStandard_source)
1632 {
1633 do_abort = true;
1634 fail_msg.push_back(os.str());
1635 }
1636 }
1637 }
1638
1639#pragma omp parallel for if (!arts_omp_in_parallel())
1640 for (Index ip = 1; ip < np; ip++) {
1641 if (do_abort) continue;
1642 try {
1643 const Numeric dr_dT_past =
1644 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
1645 const Numeric dr_dT_this =
1646 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
1647 stepwise_transmission(lyr_tra[ip],
1648 dlyr_tra_above[ip],
1649 dlyr_tra_below[ip],
1650 K[ip - 1],
1651 K[ip],
1652 dK_dx[ip - 1],
1653 dK_dx[ip],
1654 ppath.lstep[ip - 1],
1655 dr_dT_past,
1656 dr_dT_this,
1657 temperature_derivative_position);
1658
1659 r[ip - 1] = ppath.lstep[ip - 1];
1660 if (temperature_derivative_position >= 0){
1661 dr_below[ip][temperature_derivative_position] = dr_dT_past;
1662 dr_above[ip][temperature_derivative_position] = dr_dT_this;
1663 }
1664 } catch (const std::runtime_error& e) {
1665 ostringstream os;
1666 os << "Runtime-error in transmission calculation at index " << ip
1667 << ": \n";
1668 os << e.what();
1669#pragma omp critical(iyEmissionStandard_transmission)
1670 {
1671 do_abort = true;
1672 fail_msg.push_back(os.str());
1673 }
1674 }
1675 }
1676
1677 ARTS_USER_ERROR_IF (do_abort,
1678 "Error messages from failed cases:\n", fail_msg)
1679 }
1680
1681 const ArrayOfTransmissionMatrix tot_tra =
1683
1684 // iy_transmittance
1685 Tensor3 iy_trans_new;
1686 if (iy_agenda_call1)
1687 iy_trans_new = tot_tra[np - 1];
1688 else
1689 iy_transmittance_mult(iy_trans_new, iy_transmittance, Tensor3{tot_tra[np - 1]});
1690
1691 // iy_aux: Optical depth
1692 if (auxOptDepth >= 0)
1693 for (Index iv = 0; iv < nf; iv++)
1694 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
1695
1696 // Radiative background
1698 iy,
1699 diy_dx,
1700 iy_trans_new,
1701 iy_id,
1702 jacobian_do,
1703 jacobian_quantities,
1704 ppath,
1705 rte_pos2,
1706 atmosphere_dim,
1707 nlte_field,
1708 cloudbox_on,
1709 stokes_dim,
1710 f_grid,
1711 iy_unit,
1712 surface_props_data,
1713 iy_main_agenda,
1714 iy_space_agenda,
1715 iy_surface_agenda,
1716 iy_cloudbox_agenda,
1717 iy_agenda_call1,
1718 verbosity);
1719
1720 lvl_rad[np - 1] = iy;
1721
1722 // Radiative transfer calculations
1723 if (rt_integration_option == "first order" || rt_integration_option == "default") {
1724 for (Index ip = np - 2; ip >= 0; ip--) {
1725 lvl_rad[ip] = lvl_rad[ip + 1];
1726 update_radiation_vector(lvl_rad[ip],
1727 dlvl_rad[ip],
1728 dlvl_rad[ip + 1],
1729 src_rad[ip],
1730 src_rad[ip + 1],
1731 dsrc_rad[ip],
1732 dsrc_rad[ip + 1],
1733 lyr_tra[ip + 1],
1734 tot_tra[ip],
1735 dlyr_tra_above[ip + 1],
1736 dlyr_tra_below[ip + 1],
1737 PropagationMatrix(),
1738 PropagationMatrix(),
1739 ArrayOfPropagationMatrix(),
1740 ArrayOfPropagationMatrix(),
1741 Numeric(),
1742 Vector(),
1743 Vector(),
1744 0,
1745 0,
1747 }
1748 } else if (rt_integration_option == "second order") {
1749 for (Index ip = np - 2; ip >= 0; ip--) {
1750 lvl_rad[ip] = lvl_rad[ip + 1];
1751 update_radiation_vector(lvl_rad[ip],
1752 dlvl_rad[ip],
1753 dlvl_rad[ip + 1],
1754 src_rad[ip],
1755 src_rad[ip + 1],
1756 dsrc_rad[ip],
1757 dsrc_rad[ip + 1],
1758 lyr_tra[ip + 1],
1759 tot_tra[ip],
1760 dlyr_tra_above[ip + 1],
1761 dlyr_tra_below[ip + 1],
1762 K[ip],
1763 K[ip + 1],
1764 dK_dx[ip + 1],
1765 dK_dx[ip + 1],
1766 r[ip],
1767 dr_above[ip + 1],
1768 dr_below[ip + 1],
1769 0,
1770 0,
1772 }
1773 } else {
1774 ARTS_USER_ERROR ( "Only allowed choices for *integration order* are "
1775 "1 and 2.");
1776 }
1777
1778 // Copy back to ARTS external style
1779 iy = lvl_rad[0];
1780 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
1781 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
1782 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
1783 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
1784 if (j_analytical_do)
1785 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
1786 dlvl_rad[ip][iq];);
1787 }
1788
1789 // Finalize analytical Jacobians
1790 if (j_analytical_do) {
1792 diy_dx,
1793 diy_dpath,
1794 ns,
1795 nf,
1796 np,
1797 atmosphere_dim,
1798 ppath,
1799 ppvar_p,
1800 ppvar_t,
1801 ppvar_vmr,
1802 iy_agenda_call1,
1803 iy_transmittance,
1804 water_p_eq_agenda,
1805 jacobian_quantities,
1806 jac_species_i);
1807 }
1808
1809 // Radiance unit conversions
1810 if (iy_agenda_call1) {
1812 diy_dx,
1813 ppvar_iy,
1814 ns,
1815 np,
1816 f_grid,
1817 ppath,
1818 jacobian_quantities,
1819 j_analytical_do,
1820 iy_unit);
1821 }
1822}
1823
1824/* Workspace method: Doxygen documentation will be auto-generated */
1826 Matrix& iy,
1827 ArrayOfMatrix& iy_aux,
1828 Ppath& ppath,
1829 ArrayOfTensor3& diy_dx,
1830 GriddedField4& atm_fields_compact,
1831 const Index& iy_id,
1832 const Vector& f_grid,
1833 const Index& atmosphere_dim,
1834 const Vector& p_grid,
1835 const Vector& lat_grid,
1836 const Vector& lon_grid,
1837 const Vector& lat_true,
1838 const Vector& lon_true,
1839 const Tensor3& t_field,
1840 const Tensor3& z_field,
1841 const Tensor4& vmr_field,
1842 const EnergyLevelMap& nlte_field,
1843 const Tensor3& wind_u_field,
1844 const Tensor3& wind_v_field,
1845 const Tensor3& wind_w_field,
1846 const Tensor3& mag_u_field,
1847 const Tensor3& mag_v_field,
1848 const Tensor3& mag_w_field,
1849 const Index& cloudbox_on,
1850 const ArrayOfIndex& cloudbox_limits,
1851 const Tensor4& pnd_field,
1852 const Matrix& particle_masses,
1853 const Agenda& ppath_agenda,
1854 const Numeric& ppath_lmax,
1855 const Numeric& ppath_lraytrace,
1856 const Index& iy_agenda_call1,
1857 const String& iy_unit,
1858 const Tensor3& iy_transmittance,
1859 const Vector& rte_pos,
1860 const Vector& rte_los,
1861 const Vector& rte_pos2,
1862 const Index& jacobian_do,
1863 const ArrayOfString& iy_aux_vars,
1864 const Agenda& iy_independent_beam_approx_agenda,
1865 const Index& return_atm1d,
1866 const Index& skip_vmr,
1867 const Index& skip_pnd,
1868 const Index& return_masses,
1869 const Verbosity&) {
1870 // Throw error if unsupported features are requested
1871 ARTS_USER_ERROR_IF (jacobian_do,
1872 "Jacobians not provided by the method, *jacobian_do* "
1873 "must be 0.");
1874 ARTS_USER_ERROR_IF (!nlte_field.value.empty(),
1875 "This method does not yet support non-empty *nlte_field*.");
1876 ARTS_USER_ERROR_IF (!wind_u_field.empty(),
1877 "This method does not yet support non-empty *wind_u_field*.");
1878 ARTS_USER_ERROR_IF (!wind_v_field.empty(),
1879 "This method does not yet support non-empty *wind_v_field*.");
1880 ARTS_USER_ERROR_IF (!wind_w_field.empty(),
1881 "This method does not yet support non-empty *wind_w_field*.");
1882 ARTS_USER_ERROR_IF (!mag_u_field.empty(),
1883 "This method does not yet support non-empty *mag_u_field*.");
1884 ARTS_USER_ERROR_IF (!mag_v_field.empty(),
1885 "This method does not yet support non-empty *mag_v_field*.");
1886 ARTS_USER_ERROR_IF (!mag_w_field.empty(),
1887 "This method does not yet support non-empty *mag_w_field*.");
1888 //
1889 if (return_masses) {
1890 ARTS_USER_ERROR_IF (pnd_field.nbooks() != particle_masses.nrows(),
1891 "Sizes of *pnd_field* and *particle_masses* "
1892 "are inconsistent.");
1893 }
1894 chk_latlon_true(atmosphere_dim,
1895 lat_grid,
1896 lat_true,
1897 lon_true);
1898
1899 // Note that input 1D atmospheres are handled exactly as 2D and 3D, to make
1900 // the function totally general. And 1D must be handled for iterative calls.
1901
1902 // Determine propagation path (with cloudbox deactivated) and check
1903 // that is OK for ICA
1904 //
1906 ppath,
1907 ppath_lmax,
1908 ppath_lraytrace,
1909 rte_pos,
1910 rte_los,
1911 rte_pos2,
1912 0,
1913 0,
1914 f_grid,
1915 ppath_agenda);
1916 //
1917 error_if_limb_ppath(ppath);
1918
1919 // If scattering and sensor inside atmosphere, we need a pseudo-ppath that
1920 // samples altitudes not covered by main ppath. We make this second path
1921 // strictly vertical.
1922 //
1923 Ppath ppath2;
1924 //
1925 if (cloudbox_on && ppath.end_lstep == 0) {
1926 Vector los_tmp = rte_los;
1927 if (abs(rte_los[0]) < 90) {
1928 los_tmp[0] = 180;
1929 } else {
1930 los_tmp[0] = 0;
1931 }
1932 //
1934 ppath2,
1935 ppath_lmax,
1936 ppath_lraytrace,
1937 rte_pos,
1938 los_tmp,
1939 rte_pos2,
1940 0,
1941 0,
1942 f_grid,
1943 ppath_agenda);
1944 } else {
1945 ppath2.np = 1;
1946 }
1947
1948 // Merge grid positions, and sort from bottom to top of atmosphere
1949 const Index np = ppath.np + ppath2.np - 1;
1950 ArrayOfGridPos gp_p(np), gp_lat(np), gp_lon(np);
1951 if (ppath.np > 1 &&
1952 ppath.pos(0, 0) >
1953 ppath.pos(1, 0)) { // Ppath is sorted in downward direction
1954 // Copy ppath in reversed order
1955 for (Index i = 0; i < ppath.np; i++) {
1956 const Index ip = ppath.np - i - 1;
1957 gp_p[i] = ppath.gp_p[ip];
1958 if (atmosphere_dim > 1) {
1959 gp_lat[i] = ppath.gp_lat[ip];
1960 if (atmosphere_dim == 3) {
1961 gp_lon[i] = ppath.gp_lon[ip];
1962 }
1963 }
1964 }
1965 // Append ppath2, but skipping element [0]
1966 for (Index i = ppath.np; i < np; i++) {
1967 const Index ip = i - ppath.np + 1;
1968 gp_p[i] = ppath2.gp_p[ip];
1969 if (atmosphere_dim > 1) {
1970 gp_lat[i] = ppath2.gp_lat[ip];
1971 if (atmosphere_dim == 3) {
1972 gp_lon[i] = ppath2.gp_lon[ip];
1973 }
1974 }
1975 }
1976 } else {
1977 // Copy ppath2 in reversed order, but skipping element [0]
1978 for (Index i = 0; i < ppath2.np - 1; i++) {
1979 const Index ip = ppath2.np - i - 1;
1980 gp_p[i] = ppath2.gp_p[ip];
1981 if (atmosphere_dim > 1) {
1982 gp_lat[i] = ppath2.gp_lat[ip];
1983 if (atmosphere_dim == 3) {
1984 gp_lon[i] = ppath2.gp_lon[ip];
1985 }
1986 }
1987 }
1988 // Append ppath
1989 for (Index i = ppath2.np - 1; i < np; i++) {
1990 const Index ip = i - ppath2.np + 1;
1991 gp_p[i] = ppath.gp_p[ip];
1992 if (atmosphere_dim > 1) {
1993 gp_lat[i] = ppath.gp_lat[ip];
1994 if (atmosphere_dim == 3) {
1995 gp_lon[i] = ppath.gp_lon[ip];
1996 }
1997 }
1998 }
1999 }
2000
2001 // 1D version of p_grid
2002 Matrix itw;
2003 Vector p1(np);
2004 ArrayOfGridPos gp0(0), gp1(1);
2005 interp_atmfield_gp2itw(itw, 1, gp_p, gp0, gp0);
2006 itw2p(p1, p_grid, gp_p, itw);
2007
2008 // 1D version of lat and lon variables
2009 Vector lat1(0), lon1(0);
2010 Vector lat_true1(1), lon_true1(1);
2011 if (atmosphere_dim == 3) {
2012 gp1[0] = gp_lat[0];
2013 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
2014 interp(lat_true1, itw, lat_grid, gp1);
2015 gp1[0] = gp_lon[0];
2016 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
2017 interp(lon_true1, itw, lon_grid, gp1);
2018 } else if (atmosphere_dim == 2) {
2019 gp1[0] = gp_lat[0];
2020 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
2021 interp(lat_true1, itw, lat_true, gp1);
2022 interp(lon_true1, itw, lon_true, gp1);
2023 } else {
2024 lat_true1[0] = lat_true[0];
2025 lon_true1[0] = lon_true[0];
2026 }
2027
2028 // 2D/3D interpolation weights
2029 interp_atmfield_gp2itw(itw, atmosphere_dim, gp_p, gp_lat, gp_lon);
2030
2031 // 1D temperature field
2032 Tensor3 t1(np, 1, 1);
2034 t1(joker, 0, 0), atmosphere_dim, t_field, gp_p, gp_lat, gp_lon, itw);
2035
2036 // 1D altitude field
2037 Tensor3 z1(np, 1, 1);
2039 z1(joker, 0, 0), atmosphere_dim, z_field, gp_p, gp_lat, gp_lon, itw);
2040
2041 // 1D VMR field
2042 Tensor4 vmr1(vmr_field.nbooks(), np, 1, 1);
2043 for (Index is = 0; is < vmr_field.nbooks(); is++) {
2044 interp_atmfield_by_itw(vmr1(is, joker, 0, 0),
2045 atmosphere_dim,
2046 vmr_field(is, joker, joker, joker),
2047 gp_p,
2048 gp_lat,
2049 gp_lon,
2050 itw);
2051 }
2052
2053 // 1D surface altitude
2054 Matrix zsurf1(1, 1);
2055 zsurf1(0, 0) = z1(0, 0, 0);
2056
2057 // 1D version of rte_pos/los
2058 Vector pos1(1);
2059 pos1[0] = rte_pos[0];
2060 Vector los1(1);
2061 los1[0] = abs(rte_los[0]);
2062 Vector pos2(0);
2063 if (rte_pos2.nelem()) {
2064 pos2 = rte_pos2[Range(0, rte_pos2.nelem())];
2065 }
2066
2067 // Cloudbox variables
2068 //
2069 Index cbox_on1 = cloudbox_on;
2070 ArrayOfIndex cbox_lims1(0);
2071 Tensor4 pnd1(0, 0, 0, 0);
2072 //
2073 if (cloudbox_on) {
2074 // Determine what p1-levels that are inside cloudbox
2075 Index ifirst = np;
2076 Index ilast = -1;
2077 for (Index i = 0; i < np; i++) {
2078 if (is_gp_inside_cloudbox(gp_p[i],
2079 gp_lat[i],
2080 gp_lon[i],
2081 cloudbox_limits,
2082 true,
2083 atmosphere_dim)) {
2084 if (i < ifirst) {
2085 ifirst = i;
2086 }
2087 ilast = i;
2088 }
2089 }
2090
2091 // If no hit, deactive cloudbox
2092 if (ifirst == np) {
2093 cbox_on1 = 0;
2094 }
2095
2096 // Otherwise set 1D cloud variables
2097 else {
2098 // We can enter the cloudbox from the side, and we need to add 1
2099 // level on each side to be safe (if not limits already at edges)
2100 //
2101 const Index extra_bot = ifirst == 0 ? 0 : 1;
2102 const Index extra_top = ilast == np - 1 ? 0 : 1;
2103 //
2104 cbox_lims1.resize(2);
2105 cbox_lims1[0] = ifirst - extra_bot;
2106 cbox_lims1[1] = ilast + extra_top;
2107
2108 // pnd_field
2109 //
2110 pnd1.resize(pnd_field.nbooks(), cbox_lims1[1] - cbox_lims1[0] + 1, 1, 1);
2111 pnd1 = 0; // As lowermost and uppermost level not always filled
2112 //
2113 itw.resize(1, Index(pow(2.0, Numeric(atmosphere_dim))));
2114 //
2115 for (Index i = extra_bot; i < pnd1.npages() - extra_top; i++) {
2116 const Index i0 = cbox_lims1[0] + i;
2117 ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
2118 interp_cloudfield_gp2itw(itw(0, joker),
2119 gpc_p[0],
2120 gpc_lat[0],
2121 gpc_lon[0],
2122 gp_p[i0],
2123 gp_lat[i0],
2124 gp_lon[i0],
2125 atmosphere_dim,
2126 cloudbox_limits);
2127 for (Index p = 0; p < pnd_field.nbooks(); p++) {
2128 interp_atmfield_by_itw(ExhaustiveVectorView{pnd1(p, i, 0, 0)},
2129 atmosphere_dim,
2130 pnd_field(p, joker, joker, joker),
2131 gpc_p,
2132 gpc_lat,
2133 gpc_lon,
2134 itw);
2135 }
2136 }
2137 }
2138 }
2139
2140 // Call sub agenda
2141 //
2142 {
2143 const Index adim1 = 1;
2144 const Numeric lmax1 = -1;
2145 Ppath ppath1d;
2146 //
2148 iy,
2149 iy_aux,
2150 ppath1d,
2151 diy_dx,
2152 iy_agenda_call1,
2153 iy_unit,
2154 iy_transmittance,
2155 iy_aux_vars,
2156 iy_id,
2157 adim1,
2158 p1,
2159 lat1,
2160 lon1,
2161 lat_true1,
2162 lon_true1,
2163 t1,
2164 z1,
2165 vmr1,
2166 zsurf1,
2167 lmax1,
2168 ppath_lraytrace,
2169 cbox_on1,
2170 cbox_lims1,
2171 pnd1,
2172 jacobian_do,
2173 f_grid,
2174 pos1,
2175 los1,
2176 pos2,
2177 iy_independent_beam_approx_agenda);
2178 }
2179
2180 // Fill *atm_fields_compact*?
2181 if (return_atm1d) {
2182 // Sizes and allocate memory
2183 const Index nvmr = skip_vmr ? 0 : vmr1.nbooks();
2184 const Index npnd = skip_pnd ? 0 : pnd1.nbooks();
2185 const Index nmass = return_masses ? particle_masses.ncols() : 0;
2186 const Index ntot = 2 + nvmr + npnd + nmass;
2187 ArrayOfString field_names(ntot);
2188 atm_fields_compact.resize(ntot, np, 1, 1);
2189
2190 // Altitudes
2191 field_names[0] = "Geometric altitudes";
2192 atm_fields_compact.data(0, joker, 0, 0) = z1(joker, 0, 0);
2193
2194 // Temperature
2195 field_names[1] = "Temperature";
2196 atm_fields_compact.data(1, joker, 0, 0) = t1(joker, 0, 0);
2197
2198 // VMRs
2199 if (nvmr) {
2200 for (Index i = 0; i < nvmr; i++) {
2201 const Index iout = 2 + i;
2202 ostringstream sstr;
2203 sstr << "VMR species " << i;
2204 field_names[iout] = sstr.str();
2205 atm_fields_compact.data(iout, joker, 0, 0) = vmr1(i, joker, 0, 0);
2206 }
2207 }
2208
2209 // PNDs
2210 if (npnd) {
2211 for (Index i = 0; i < npnd; i++) {
2212 const Index iout = 2 + nvmr + i;
2213 ostringstream sstr;
2214 sstr << "Scattering element " << i;
2215 field_names[iout] = sstr.str();
2216 atm_fields_compact.data(iout, joker, 0, 0) = 0;
2217 atm_fields_compact.data(
2218 iout, Range(cbox_lims1[0], pnd1.npages()), 0, 0) =
2219 pnd1(i, joker, 0, 0);
2220 }
2221 }
2222
2223 // Masses
2224 if (nmass) {
2225 for (Index i = 0; i < nmass; i++) {
2226 const Index iout = 2 + nvmr + npnd + i;
2227 ostringstream sstr;
2228 sstr << "Mass category " << i;
2229 field_names[iout] = sstr.str();
2230 atm_fields_compact.data(iout, joker, 0, 0) = 0;
2231 for (Index ip = cbox_lims1[0]; ip < pnd1.npages(); ip++) {
2232 for (Index is = 0; is < pnd1.nbooks(); is++) {
2233 atm_fields_compact.data(iout, ip, 0, 0) +=
2234 particle_masses(is, i) * pnd1(is, ip, 0, 0);
2235 }
2236 }
2237 }
2238 }
2239
2240 // Finally, set grids and names
2241 //
2242 atm_fields_compact.set_name(
2243 "Data created by *iyIndependentBeamApproximation*");
2244 //
2245 atm_fields_compact.set_grid_name(GFIELD4_FIELD_NAMES,
2246 "Atmospheric quantity");
2247 atm_fields_compact.set_grid(GFIELD4_FIELD_NAMES, field_names);
2248 atm_fields_compact.set_grid_name(GFIELD4_P_GRID, "Pressure");
2249 atm_fields_compact.set_grid(GFIELD4_P_GRID, p1);
2250 atm_fields_compact.set_grid_name(GFIELD4_LAT_GRID, "Latitude");
2251 atm_fields_compact.set_grid(GFIELD4_LAT_GRID, lat_true1);
2252 atm_fields_compact.set_grid_name(GFIELD4_LON_GRID, "Longitude");
2253 atm_fields_compact.set_grid(GFIELD4_LON_GRID, lon_true1);
2254 }
2255}
2256
2257/* Workspace method: Doxygen documentation will be auto-generated */
2259 Matrix& iy,
2260 ArrayOfMatrix& iy_aux,
2261 Ppath& ppath,
2262 ArrayOfTensor3& diy_dx,
2263 const ArrayOfString& iy_aux_vars,
2264 const Index& iy_agenda_call1,
2265 const Tensor3& iy_transmittance,
2266 const Vector& rte_pos,
2267 const Vector& rte_los,
2268 const Vector& rte_pos2,
2269 const Index& stokes_dim,
2270 const Vector& f_grid,
2271 const Agenda& iy_loop_freqs_agenda,
2272 const Verbosity&) {
2273 // Throw error if unsupported features are requested
2274 ARTS_USER_ERROR_IF (!iy_agenda_call1,
2275 "Recursive usage not possible (iy_agenda_call1 must be 1).");
2276 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
2277 "*iy_transmittance* must be empty.");
2278
2279 const Index nf = f_grid.nelem();
2280
2281 for (Index i = 0; i < nf; i++) {
2282 // Variables for 1 frequency
2283 Matrix iy1;
2284 ArrayOfMatrix iy_aux1;
2285 ArrayOfTensor3 diy_dx1;
2286
2288 iy1,
2289 iy_aux1,
2290 ppath,
2291 diy_dx1,
2292 iy_agenda_call1,
2293 iy_transmittance,
2294 iy_aux_vars,
2295 0,
2296 Vector(1, f_grid[i]),
2297 rte_pos,
2298 rte_los,
2299 rte_pos2,
2300 iy_loop_freqs_agenda);
2301
2302 // After first frequency, give output its size
2303 if (i == 0) {
2304 iy.resize(nf, stokes_dim);
2305 //
2306 iy_aux.resize(iy_aux1.nelem());
2307 for (Index q = 0; q < iy_aux1.nelem(); q++) {
2308 iy_aux[q].resize(nf, stokes_dim);
2309 }
2310 //
2311 diy_dx.resize(diy_dx1.nelem());
2312 for (Index q = 0; q < diy_dx1.nelem(); q++) {
2313 diy_dx[q].resize(diy_dx1[q].npages(), nf, stokes_dim);
2314 }
2315 }
2316
2317 // Copy to output variables
2318 iy(i, joker) = iy1(0, joker);
2319 for (Index q = 0; q < iy_aux1.nelem(); q++) {
2320 iy_aux[q](i, joker) = iy_aux1[q](0, joker);
2321 }
2322 for (Index q = 0; q < diy_dx1.nelem(); q++) {
2323 diy_dx[q](joker, i, joker) = diy_dx1[q](joker, 0, joker);
2324 }
2325 }
2326}
2327
2328/* Workspace method: Doxygen documentation will be auto-generated */
2330 Matrix& iy,
2331 ArrayOfMatrix& iy_aux,
2332 ArrayOfTensor3& diy_dx,
2333 const Index& iy_agenda_call1,
2334 const Tensor3& iy_transmittance,
2335 const Vector& rte_pos,
2336 const Vector& rte_los,
2337 const ArrayOfString& iy_aux_vars,
2338 const Index& jacobian_do,
2339 const Index& atmosphere_dim,
2340 const Vector& p_grid,
2341 const Vector& lat_grid,
2342 const Vector& lon_grid,
2343 const Tensor3& z_field,
2344 const Tensor3& t_field,
2345 const Tensor4& vmr_field,
2346 const Vector& refellipsoid,
2347 const Matrix& z_surface,
2348 const Index& cloudbox_on,
2349 const ArrayOfIndex& cloudbox_limits,
2350 const Index& stokes_dim,
2351 const Vector& f_grid,
2352 const ArrayOfArrayOfSingleScatteringData& scat_data,
2353 const Agenda& iy_space_agenda,
2354 const Agenda& surface_rtprop_agenda,
2355 const Agenda& propmat_clearsky_agenda,
2356 const Agenda& ppath_step_agenda,
2357 const Numeric& ppath_lmax,
2358 const Numeric& ppath_lraytrace,
2359 const Tensor4& pnd_field,
2360 const String& iy_unit,
2361 const Numeric& mc_std_err,
2362 const Index& mc_max_time,
2363 const Index& mc_max_iter,
2364 const Index& mc_min_iter,
2365 const Numeric& mc_taustep_limit,
2366 const Index& t_interp_order,
2367 const Verbosity& verbosity) {
2368 // Throw error if unsupported features are requested
2369 ARTS_USER_ERROR_IF (atmosphere_dim != 3,
2370 "Only 3D atmospheres are allowed (atmosphere_dim must be 3)");
2371 ARTS_USER_ERROR_IF (!cloudbox_on,
2372 "The cloudbox must be activated (cloudbox_on must be 1)");
2373 ARTS_USER_ERROR_IF (jacobian_do,
2374 "This method does not provide any jacobians (jacobian_do must be 0)");
2375 ARTS_USER_ERROR_IF (!iy_agenda_call1,
2376 "Recursive usage not possible (iy_agenda_call1 must be 1)");
2377 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
2378 "*iy_transmittance* must be empty");
2379
2380 // Size output variables
2381 //
2382 const Index nf = f_grid.nelem();
2383 //
2384 iy.resize(nf, stokes_dim);
2385 diy_dx.resize(0);
2386 //
2387 //=== iy_aux part ===========================================================
2388 Index auxError = -1;
2389 {
2390 const Index naux = iy_aux_vars.nelem();
2391 iy_aux.resize(naux);
2392 //
2393 for (Index i = 0; i < naux; i++) {
2394 if (iy_aux_vars[i] == "Error (uncorrelated)") {
2395 auxError = i;
2396 iy_aux[i].resize(nf, stokes_dim);
2397 } else {
2399 "In *iy_aux_vars* you have included: \"", iy_aux_vars[i],
2400 "\"\nThis choice is not recognised.")
2401 }
2402 }
2403 }
2404 //===========================================================================
2405
2406 // Some MC variables are only local here
2407 //
2408 MCAntenna mc_antenna;
2409 mc_antenna.set_pencil_beam();
2410
2411 // Pos and los must be matrices
2412 Matrix pos(1, 3), los(1, 2);
2413 //
2414 pos(0, joker) = rte_pos;
2415 los(0, joker) = rte_los;
2416
2417 String fail_msg;
2418 bool failed = false;
2419
2420 if (nf) {
2422
2423#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) firstprivate(wss)
2424 for (Index f_index = 0; f_index < nf; f_index++) {
2425 if (failed) continue;
2426
2427 try {
2428 // Seed reset for each loop. If not done, the errors
2429 // appear to be highly correlated.
2430 Index mc_seed;
2431 MCSetSeedFromTime(mc_seed, verbosity);
2432
2433 Vector y, mc_error;
2434 Index mc_iteration_count;
2435 Tensor3 mc_points;
2436 ArrayOfIndex mc_scat_order, mc_source_domain;
2437
2438 MCGeneral(wss,
2439 y,
2440 mc_iteration_count,
2441 mc_error,
2442 mc_points,
2443 mc_scat_order,
2444 mc_source_domain,
2445 mc_antenna,
2446 f_grid,
2447 f_index,
2448 pos,
2449 los,
2450 stokes_dim,
2451 atmosphere_dim,
2452 ppath_step_agenda,
2453 ppath_lmax,
2454 ppath_lraytrace,
2455 iy_space_agenda,
2456 surface_rtprop_agenda,
2457 propmat_clearsky_agenda,
2458 p_grid,
2459 lat_grid,
2460 lon_grid,
2461 z_field,
2462 refellipsoid,
2463 z_surface,
2464 t_field,
2465 vmr_field,
2466 cloudbox_on,
2467 cloudbox_limits,
2468 pnd_field,
2469 scat_data,
2470 1,
2471 1,
2472 1,
2473 1,
2474 iy_unit,
2475 mc_seed,
2476 mc_std_err,
2477 mc_max_time,
2478 mc_max_iter,
2479 mc_min_iter,
2480 mc_taustep_limit,
2481 1,
2482 t_interp_order,
2483 verbosity);
2484
2485 ARTS_ASSERT(y.nelem() == stokes_dim);
2486
2487 iy(f_index, joker) = y;
2488
2489 if (auxError >= 0) {
2490 iy_aux[auxError](f_index, joker) = mc_error;
2491 }
2492 } catch (const std::exception& e) {
2493 ostringstream os;
2494 os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
2495 << ")" << endl
2496 << e.what();
2497#pragma omp critical(iyMC_fail)
2498 {
2499 failed = true;
2500 fail_msg = os.str();
2501 }
2502 continue;
2503 }
2504 }
2505 }
2506
2507 ARTS_USER_ERROR_IF (failed, fail_msg);
2508}
2509
2510/* Workspace method: Doxygen documentation will be auto-generated */
2511void iyReplaceFromAux(Matrix& iy,
2512 const ArrayOfMatrix& iy_aux,
2513 const ArrayOfString& iy_aux_vars,
2514 const Index& jacobian_do,
2515 const String& aux_var,
2516 const Verbosity&) {
2517 ARTS_USER_ERROR_IF (iy_aux.nelem() != iy_aux_vars.nelem(),
2518 "*iy_aux* and *iy_aux_vars* must have the same "
2519 "number of elements.");
2520
2521 ARTS_USER_ERROR_IF (jacobian_do,
2522 "This method can not provide any jacobians and "
2523 "*jacobian_do* must be 0.");
2524
2525 bool ready = false;
2526
2527 for (Index i = 0; i < iy_aux.nelem() && !ready; i++) {
2528 if (iy_aux_vars[i] == aux_var) {
2529 iy = iy_aux[i];
2530 ready = true;
2531 }
2532 }
2533
2534 ARTS_USER_ERROR_IF (!ready,
2535 "The selected auxiliary variable to insert in *iy* "
2536 "is either not defined at all or is not set.");
2537}
2538
2539/* Workspace method: Doxygen documentation will be auto-generated */
2541 Matrix& ppvar_optical_depth,
2542 const Tensor4& ppvar_trans_cumulat,
2543 const Verbosity&) {
2544 ppvar_optical_depth = ppvar_trans_cumulat(joker, joker, 0, 0);
2545 transform(ppvar_optical_depth, log, ppvar_optical_depth);
2546 ppvar_optical_depth *= -1;
2547}
2548
2549/* Workspace method: Doxygen documentation will be auto-generated */
2551 Vector& y,
2552 Vector& y_f,
2553 ArrayOfIndex& y_pol,
2554 Matrix& y_pos,
2555 Matrix& y_los,
2556 ArrayOfVector& y_aux,
2557 Matrix& y_geo,
2558 Matrix& jacobian,
2559 const Index& atmgeom_checked,
2560 const Index& atmfields_checked,
2561 const Index& atmosphere_dim,
2562 const EnergyLevelMap& nlte_field,
2563 const Index& cloudbox_on,
2564 const Index& cloudbox_checked,
2565 const Index& scat_data_checked,
2566 const Index& sensor_checked,
2567 const Index& stokes_dim,
2568 const Vector& f_grid,
2569 const Matrix& sensor_pos,
2570 const Matrix& sensor_los,
2571 const Matrix& transmitter_pos,
2572 const Matrix& mblock_dlos,
2573 const Sparse& sensor_response,
2574 const Vector& sensor_response_f,
2575 const ArrayOfIndex& sensor_response_pol,
2576 const Matrix& sensor_response_dlos,
2577 const String& iy_unit,
2578 const Agenda& iy_main_agenda,
2579 const Agenda& jacobian_agenda,
2580 const Index& jacobian_do,
2581 const ArrayOfRetrievalQuantity& jacobian_quantities,
2582 const ArrayOfString& iy_aux_vars,
2583 const Verbosity& verbosity) {
2585
2586 // Basics
2587 //
2588 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2589 //
2590 ARTS_USER_ERROR_IF (f_grid.empty(),
2591 "The frequency grid is empty.");
2592 chk_if_increasing("f_grid", f_grid);
2593 ARTS_USER_ERROR_IF (f_grid[0] <= 0,
2594 "All frequencies in *f_grid* must be > 0.");
2595 //
2596 ARTS_USER_ERROR_IF (atmfields_checked != 1,
2597 "The atmospheric fields must be flagged to have\n"
2598 "passed a consistency check (atmfields_checked=1).");
2599 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
2600 "The atmospheric geometry must be flagged to have\n"
2601 "passed a consistency check (atmgeom_checked=1).");
2602 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
2603 "The cloudbox must be flagged to have\n"
2604 "passed a consistency check (cloudbox_checked=1).");
2605 if (cloudbox_on)
2606 ARTS_USER_ERROR_IF (scat_data_checked != 1,
2607 "The scattering data must be flagged to have\n"
2608 "passed a consistency check (scat_data_checked=1).");
2609 ARTS_USER_ERROR_IF (sensor_checked != 1,
2610 "The sensor variables must be flagged to have\n"
2611 "passed a consistency check (sensor_checked=1).");
2612
2613 // Some sizes
2614 const Index nf = f_grid.nelem();
2615 const Index nlos = mblock_dlos.nrows();
2616 const Index n1y = sensor_response.nrows();
2617 const Index nmblock = sensor_pos.nrows();
2618 const Index niyb = nf * nlos * stokes_dim;
2619
2620 //---------------------------------------------------------------------------
2621 // Allocations and resizing
2622 //---------------------------------------------------------------------------
2623
2624 // Resize *y* and *y_XXX*
2625 //
2626 y.resize(nmblock * n1y);
2627 y_f.resize(nmblock * n1y);
2628 y_pol.resize(nmblock * n1y);
2629 y_pos.resize(nmblock * n1y, sensor_pos.ncols());
2630 y_los.resize(nmblock * n1y, sensor_los.ncols());
2631 y_geo.resize(nmblock * n1y, 5);
2632 y_geo = NAN; // Will be replaced if relavant data are provided (*geo_pos*)
2633
2634 // For y_aux we don't know the number of quantities, and we need to
2635 // store all output
2636 ArrayOfArrayOfVector iyb_aux_array(nmblock);
2637
2638 // Jacobian variables
2639 //
2640 Index j_analytical_do = 0;
2641 ArrayOfArrayOfIndex jacobian_indices;
2642 //
2643 if (jacobian_do) {
2644 bool any_affine;
2645 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
2646 jacobian.resize(nmblock * n1y,
2647 jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
2648 jacobian = 0;
2649 //
2650 FOR_ANALYTICAL_JACOBIANS_DO2(j_analytical_do = 1;)
2651 } else {
2652 jacobian.resize(0, 0);
2653 }
2654
2655 //---------------------------------------------------------------------------
2656 // The calculations
2657 //---------------------------------------------------------------------------
2658
2659 String fail_msg;
2660 bool failed = false;
2661
2662 if (nmblock >= arts_omp_get_max_threads() ||
2663 (nf <= nmblock && nmblock >= nlos)) {
2664 out3 << " Parallelizing mblock loop (" << nmblock << " iterations)\n";
2665
2667
2668#pragma omp parallel for firstprivate(wss)
2669 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2670 // Skip remaining iterations if an error occurred
2671 if (failed) continue;
2672
2674 fail_msg,
2675 iyb_aux_array,
2676 wss,
2677 y,
2678 y_f,
2679 y_pol,
2680 y_pos,
2681 y_los,
2682 y_geo,
2683 jacobian,
2684 atmosphere_dim,
2685 nlte_field,
2686 cloudbox_on,
2687 stokes_dim,
2688 f_grid,
2689 sensor_pos,
2690 sensor_los,
2691 transmitter_pos,
2692 mblock_dlos,
2693 sensor_response,
2694 sensor_response_f,
2695 sensor_response_pol,
2696 sensor_response_dlos,
2697 iy_unit,
2698 iy_main_agenda,
2699 jacobian_agenda,
2700 jacobian_do,
2701 jacobian_quantities,
2702 jacobian_indices,
2703 iy_aux_vars,
2704 verbosity,
2705 mblock_index,
2706 n1y,
2707 j_analytical_do);
2708 } // End mblock loop
2709 } else {
2710 out3 << " Not parallelizing mblock loop (" << nmblock << " iterations)\n";
2711
2712 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2713 // Skip remaining iterations if an error occurred
2714 if (failed) continue;
2715
2717 fail_msg,
2718 iyb_aux_array,
2719 ws,
2720 y,
2721 y_f,
2722 y_pol,
2723 y_pos,
2724 y_los,
2725 y_geo,
2726 jacobian,
2727 atmosphere_dim,
2728 nlte_field,
2729 cloudbox_on,
2730 stokes_dim,
2731 f_grid,
2732 sensor_pos,
2733 sensor_los,
2734 transmitter_pos,
2735 mblock_dlos,
2736 sensor_response,
2737 sensor_response_f,
2738 sensor_response_pol,
2739 sensor_response_dlos,
2740 iy_unit,
2741 iy_main_agenda,
2742 jacobian_agenda,
2743 jacobian_do,
2744 jacobian_quantities,
2745 jacobian_indices,
2746 iy_aux_vars,
2747 verbosity,
2748 mblock_index,
2749 n1y,
2750 j_analytical_do);
2751 } // End mblock loop
2752 }
2753
2754 // Rethrow exception if a runtime error occurred in the mblock loop
2755 ARTS_USER_ERROR_IF (failed, fail_msg);
2756
2757 // Compile y_aux
2758 //
2759 const Index nq = iyb_aux_array[0].nelem();
2760 y_aux.resize(nq);
2761 //
2762 for (Index q = 0; q < nq; q++) {
2763 y_aux[q].resize(nmblock * n1y);
2764 //
2765 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2766 const Range rowind =
2767 get_rowindex_for_mblock(sensor_response, mblock_index);
2768 const Index row0 = rowind.offset;
2769
2770 // The sensor response must be applied in a special way for
2771 // uncorrelated errors. Schematically: sqrt( H.^2 * y.^2 )
2772 if (iy_aux_vars[q] == "Error (uncorrelated)") {
2773 for (Index i = 0; i < n1y; i++) {
2774 const Index row = row0 + i;
2775 y_aux[q][row] = 0;
2776 for (Index j = 0; j < niyb; j++) {
2777 y_aux[q][row] +=
2778 pow(sensor_response(i, j) * iyb_aux_array[mblock_index][q][j],
2779 (Numeric)2.0);
2780 }
2781 y_aux[q][row] = sqrt(y_aux[q][row]);
2782 }
2783 } else {
2784 mult(y_aux[q][rowind], sensor_response, iyb_aux_array[mblock_index][q]);
2785 }
2786 }
2787 }
2788}
2789
2790/* Workspace method: Doxygen documentation will be auto-generated */
2792 Vector& y,
2793 Vector& y_f,
2794 ArrayOfIndex& y_pol,
2795 Matrix& y_pos,
2796 Matrix& y_los,
2797 ArrayOfVector& y_aux,
2798 Matrix& y_geo,
2799 Matrix& jacobian,
2800 ArrayOfRetrievalQuantity& jacobian_quantities,
2801 const Index& atmfields_checked,
2802 const Index& atmgeom_checked,
2803 const Index& atmosphere_dim,
2804 const EnergyLevelMap& nlte_field,
2805 const Index& cloudbox_on,
2806 const Index& cloudbox_checked,
2807 const Index& scat_data_checked,
2808 const Index& sensor_checked,
2809 const Index& stokes_dim,
2810 const Vector& f_grid,
2811 const Matrix& sensor_pos,
2812 const Matrix& sensor_los,
2813 const Matrix& transmitter_pos,
2814 const Matrix& mblock_dlos,
2815 const Sparse& sensor_response,
2816 const Vector& sensor_response_f,
2817 const ArrayOfIndex& sensor_response_pol,
2818 const Matrix& sensor_response_dlos,
2819 const String& iy_unit,
2820 const Agenda& iy_main_agenda,
2821 const Agenda& jacobian_agenda,
2822 const Index& jacobian_do,
2823 const ArrayOfString& iy_aux_vars,
2824 const ArrayOfRetrievalQuantity& jacobian_quantities_copy,
2825 const Index& append_instrument_wfs,
2826 const Verbosity& verbosity) {
2827 // The jacobian indices of old and new part (without transformations)
2828 ArrayOfArrayOfIndex jacobian_indices, jacobian_indices_copy;
2829 {
2830 bool any_affine;
2832 jacobian_indices_copy, any_affine, jacobian_quantities_copy, true);
2833 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
2834 }
2835
2836 // Check consistency of data representing first measurement
2837 const Index n1 = y.nelem();
2838 Index nrq1 = 0;
2839 ARTS_USER_ERROR_IF (y.empty(),
2840 "Input *y* is empty. Use *yCalc*");
2841 ARTS_USER_ERROR_IF (y_f.nelem() != n1,
2842 "Lengths of input *y* and *y_f* are inconsistent.");
2843 ARTS_USER_ERROR_IF (y_pol.nelem() != n1,
2844 "Lengths of input *y* and *y_pol* are inconsistent.");
2845 ARTS_USER_ERROR_IF (y_pos.nrows() != n1,
2846 "Sizes of input *y* and *y_pos* are inconsistent.");
2847 ARTS_USER_ERROR_IF (y_los.nrows() != n1,
2848 "Sizes of input *y* and *y_los* are inconsistent.");
2849 ARTS_USER_ERROR_IF (y_geo.nrows() != n1,
2850 "Sizes of input *y* and *y_geo* are inconsistent.");
2851 if (jacobian_do) {
2852 nrq1 = jacobian_quantities_copy.nelem();
2853 ARTS_USER_ERROR_IF (jacobian.nrows() != n1,
2854 "Sizes of *y* and *jacobian* are inconsistent.");
2855 ARTS_USER_ERROR_IF (jacobian.ncols() != jacobian_indices_copy[nrq1 - 1][1] + 1,
2856 "Size of input *jacobian* and size implied "
2857 "*jacobian_quantities_copy* are inconsistent.");
2858 }
2859
2860 // Calculate new measurement
2861 //
2862 Vector y2, y_f2;
2863 Matrix y_pos2, y_los2, y_geo2, jacobian2;
2864 ArrayOfIndex y_pol2;
2865 ArrayOfVector y_aux2;
2866 //
2867 yCalc(ws,
2868 y2,
2869 y_f2,
2870 y_pol2,
2871 y_pos2,
2872 y_los2,
2873 y_aux2,
2874 y_geo2,
2875 jacobian2,
2876 atmfields_checked,
2877 atmgeom_checked,
2878 atmosphere_dim,
2879 nlte_field,
2880 cloudbox_on,
2881 cloudbox_checked,
2882 scat_data_checked,
2883 sensor_checked,
2884 stokes_dim,
2885 f_grid,
2886 sensor_pos,
2887 sensor_los,
2888 transmitter_pos,
2889 mblock_dlos,
2890 sensor_response,
2891 sensor_response_f,
2892 sensor_response_pol,
2893 sensor_response_dlos,
2894 iy_unit,
2895 iy_main_agenda,
2896 jacobian_agenda,
2897 jacobian_do,
2898 jacobian_quantities,
2899 iy_aux_vars,
2900 verbosity);
2901
2902 // Consistency checks
2903 ARTS_USER_ERROR_IF (y_pos.ncols() != y_pos2.ncols(),
2904 "Different number of columns in *y_pos* between the measurements.");
2905 ARTS_USER_ERROR_IF (y_los.ncols() != y_los2.ncols(),
2906 "Different number of columns in *y_los* between the measurements.");
2907
2908 // y and y_XXX
2909 //
2910 const Index n2 = y2.nelem();
2911 //
2912 {
2913 // Make copy of old measurement
2914 const Vector y1 = y, y_f1 = y_f;
2915 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
2916 const ArrayOfIndex y_pol1 = y_pol;
2917 const ArrayOfVector y_aux1 = y_aux;
2918 //
2919 y.resize(n1 + n2);
2920 y[Range(0, n1)] = y1;
2921 y[Range(n1, n2)] = y2;
2922 //
2923 y_f.resize(n1 + n2);
2924 y_f[Range(0, n1)] = y_f1;
2925 y_f[Range(n1, n2)] = y_f2;
2926 //
2927 y_pos.resize(n1 + n2, y_pos1.ncols());
2928 y_pos(Range(0, n1), joker) = y_pos1;
2929 y_pos(Range(n1, n2), joker) = y_pos2;
2930 //
2931 y_los.resize(n1 + n2, y_los1.ncols());
2932 y_los(Range(0, n1), joker) = y_los1;
2933 y_los(Range(n1, n2), joker) = y_los2;
2934 //
2935 y_geo.resize(n1 + n2, y_geo1.ncols());
2936 y_geo(Range(0, n1), joker) = y_geo1;
2937 y_geo(Range(n1, n2), joker) = y_geo2;
2938 //
2939 y_pol.resize(n1 + n2);
2940 for (Index i = 0; i < n1; i++) {
2941 y_pol[i] = y_pol1[i];
2942 }
2943 for (Index i = 0; i < n2; i++) {
2944 y_pol[n1 + i] = y_pol2[i];
2945 }
2946
2947 // y_aux
2948 const Index na1 = y_aux1.nelem();
2949 const Index na2 = y_aux2.nelem();
2950 const Index na = max(na1, na2);
2951 //
2952 y_aux.resize(na);
2953 //
2954 for (Index a = 0; a < na; a++) {
2955 y_aux[a].resize(n1 + n2);
2956 if (a < na1) {
2957 y_aux[a][Range(0, n1)] = y_aux1[a];
2958 } else {
2959 y_aux[a][Range(0, n1)] = 0;
2960 }
2961 if (a < na2) {
2962 y_aux[a][Range(n1, n2)] = y_aux2[a];
2963 } else {
2964 y_aux[a][Range(n1, n2)] = 0;
2965 }
2966 }
2967 }
2968
2969 // Jacobian and friends
2970 if (jacobian_do) {
2971 // Put in old jacobian_quantities and jacobian_indices as first part in
2972 // new version of these variables
2973 ArrayOfRetrievalQuantity jacobian_quantities2 = jacobian_quantities;
2974 ArrayOfArrayOfIndex jacobian_indices2 = jacobian_indices;
2975 //
2976 jacobian_quantities = jacobian_quantities_copy;
2977 jacobian_indices = jacobian_indices_copy;
2978
2979 // Loop new jacobian_quantities to determine how new jacobian data shall
2980 // be inserted
2981 //
2982 const Index nrq2 = jacobian_quantities2.nelem();
2983 ArrayOfIndex map_table(nrq2);
2984 //
2985 for (Index q2 = 0; q2 < nrq2; q2++) {
2986 Index pos = -1;
2987
2988 // Compare to old quantities, to determine if append shall be
2989 // considered. Some special checks performed here, grids checked later
2990 if (jacobian_quantities2[q2].Target().isSpeciesVMR() ||
2991 jacobian_quantities2[q2] == Jacobian::Atm::Temperature ||
2992 jacobian_quantities2[q2] == Jacobian::Special::ScatteringString ||
2993 jacobian_quantities2[q2].Target().isWind() ||
2994 jacobian_quantities2[q2] == Jacobian::Special::SurfaceString ||
2995 append_instrument_wfs) {
2996 for (Index q1 = 0; q1 < nrq1; q1++ && pos < 0) { // FIXME: What is with this "&& pos < 0"
2997 if (jacobian_quantities2[q2].Target().sameTargetType(jacobian_quantities_copy[q1].Target())) {
2998 if (jacobian_quantities2[q2].Target().isSpeciesVMR()) {
2999 if (jacobian_quantities2[q2].Subtag() ==
3000 jacobian_quantities_copy[q1].Subtag()) {
3001 if (jacobian_quantities2[q2].Mode() ==
3002 jacobian_quantities_copy[q1].Mode()) {
3003 pos = q1;
3004 } else {
3006 "Jacobians for ", jacobian_quantities2[q2],
3007 " shall be appended.\nThis requires "
3008 "that the same retrieval unit is used "
3009 "but it seems that this requirement is "
3010 "not met.")
3011 }
3012 }
3013 } else if (jacobian_quantities2[q2] == Jacobian::Atm::Temperature) {
3014 if (jacobian_quantities2[q2].Subtag() ==
3015 jacobian_quantities_copy[q1].Subtag()) {
3016 pos = q1;
3017 } else {
3019 "Jacobians for ", jacobian_quantities2[q2],
3020 " shall be appended.\nThis requires "
3021 "that HSE is either ON or OFF for both "
3022 "parts but it seems that this requirement "
3023 "is not met.")
3024 }
3025 } else if (jacobian_quantities[q2] == Jacobian::Special::ScatteringString) {
3026 if ((jacobian_quantities2[q2].Subtag() ==
3027 jacobian_quantities_copy[q1].Subtag()) &&
3028 (jacobian_quantities2[q2].SubSubtag() ==
3029 jacobian_quantities_copy[q1].SubSubtag())) {
3030 pos = q1;
3031 }
3032 } else if (jacobian_quantities2[q2].Subtag() == jacobian_quantities_copy[q1].Subtag()) {
3033 pos = q1;
3034 }
3035 }
3036 }
3037 }
3038
3039 // New quantity
3040 if (pos < 0) {
3041 map_table[q2] = jacobian_quantities.nelem();
3042 jacobian_quantities.push_back(jacobian_quantities2[q2]);
3043 ArrayOfIndex indices(2);
3044 indices[0] = jacobian_indices[jacobian_indices.nelem() - 1][1] + 1;
3045 indices[1] =
3046 indices[0] + jacobian_indices2[q2][1] - jacobian_indices2[q2][0];
3047 jacobian_indices.push_back(indices);
3048 }
3049 // Existing quantity
3050 else {
3051 map_table[q2] = pos;
3052 // Check if grids are equal
3053 ArrayOfVector grids1 = jacobian_quantities_copy[pos].Grids();
3054 ArrayOfVector grids2 = jacobian_quantities2[q2].Grids();
3055 bool any_wrong = false;
3056 if (grids1.nelem() != grids2.nelem()) {
3057 any_wrong = true;
3058 } else {
3059 for (Index g = 0; g < grids1.nelem(); g++) {
3060 if (grids1[g].nelem() != grids2[g].nelem()) {
3061 any_wrong = true;
3062 } else {
3063 for (Index e = 0; e < grids1[g].nelem(); e++) {
3064 const Numeric v1 = grids1[g][e];
3065 const Numeric v2 = grids2[g][e];
3066 if ((v1 == 0 && abs(v2) > 1e-9) || abs(v1 - v2) / v1 > 1e-6) {
3067 any_wrong = true;
3068 }
3069 }
3070 }
3071 }
3072 }
3073 if (any_wrong) {
3075 "Jacobians for ", jacobian_quantities2[q2],
3076 " shall be appended.\nThis requires that the "
3077 "same grids are used for both measurements,\nbut "
3078 "it seems that this requirement is not met.")
3079 }
3080 }
3081 }
3082
3083 // Create and fill *jacobian*
3084 //
3085 const Index nrq = jacobian_quantities.nelem();
3086 const Matrix jacobian1 = jacobian;
3087 //
3088 jacobian.resize(n1 + n2, jacobian_indices[nrq - 1][1] + 1);
3089 jacobian = 0;
3090 //
3091 // Put in old part in top-left corner
3092 jacobian(Range(0, n1), Range(0, jacobian_indices_copy[nrq1 - 1][1] + 1)) =
3093 jacobian1;
3094 // New parts
3095 for (Index q2 = 0; q2 < nrq2; q2++) {
3096 jacobian(Range(n1, n2),
3097 Range(jacobian_indices[map_table[q2]][0],
3098 jacobian_indices[map_table[q2]][1] -
3099 jacobian_indices[map_table[q2]][0] + 1)) =
3100 jacobian2(
3101 joker,
3102 Range(jacobian_indices2[q2][0],
3103 jacobian_indices2[q2][1] - jacobian_indices2[q2][0] + 1));
3104 }
3105 }
3106}
3107
3108/* Workspace method: Doxygen documentation will be auto-generated */
3109void yApplyUnit(Vector& y,
3110 Matrix& jacobian,
3111 const Vector& y_f,
3112 const ArrayOfIndex& y_pol,
3113 const String& iy_unit,
3114 const Verbosity&) {
3115 ARTS_USER_ERROR_IF (iy_unit == "1",
3116 "No need to use this method with *iy_unit* = \"1\".");
3117
3118 ARTS_USER_ERROR_IF (max(y) > 1e-3,
3119 "The spectrum vector *y* is required to have original radiance\n"
3120 "unit, but this seems not to be the case. This as a value above\n"
3121 "1e-3 is found in *y*.")
3122
3123 // Is jacobian set?
3124 //
3125 const Index ny = y.nelem();
3126 //
3127 const bool do_j = jacobian.nrows() == ny;
3128
3129 // Some jacobian quantities can not be handled
3130 ARTS_USER_ERROR_IF (do_j && max(jacobian) > 1e-3,
3131 "The method can not be used with jacobian quantities that are not\n"
3132 "obtained through radiative transfer calculations. One example on\n"
3133 "quantity that can not be handled is *jacobianAddPolyfit*.\n"
3134 "The maximum value of *jacobian* indicates that one or several\n"
3135 "such jacobian quantities are included.")
3136
3137 // Planck-Tb
3138 //--------------------------------------------------------------------------
3139 if (iy_unit == "PlanckBT") {
3140 // Hard to use telescoping here as the data are sorted differently in y
3141 // and jacobian, than what is expected apply_iy_unit. Copy to temporary
3142 // variables instead.
3143
3144 // Handle the elements in "frequency chunks"
3145
3146 Index i0 = 0; // Index of first element for present chunk
3147 //
3148 while (i0 < ny) {
3149 // Find number of values for this chunk
3150 Index n = 1;
3151 //
3152 while (i0 + n < ny && y_f[i0] == y_f[i0 + n]) {
3153 n++;
3154 }
3155
3156 Matrix yv(1, n);
3157 ArrayOfIndex i_pol(n);
3158 bool any_quv = false;
3159 //
3160 for (Index i = 0; i < n; i++) {
3161 const Index ix = i0 + i;
3162 yv(0, i) = y[ix];
3163 i_pol[i] = y_pol[ix];
3164 if (i_pol[i] > 1 && i_pol[i] < 5) {
3165 any_quv = true;
3166 }
3167 }
3168
3169 // Index of elements to convert
3170 Range ii(i0, n);
3171
3172 if (do_j) {
3173 ARTS_USER_ERROR_IF (any_quv && i_pol[0] != 1,
3174 "The conversion to PlanckBT, of the Jacobian and "
3175 "errors for Q, U and V, requires that I (first Stokes "
3176 "element) is at hand and that the data are sorted in "
3177 "such way that I comes first for each frequency.")
3178
3179 // Jacobian
3180 Tensor3 J(jacobian.ncols(), 1, n);
3181 J(joker, 0, joker) = transpose(jacobian(ii, joker));
3182 apply_iy_unit2(J, yv, iy_unit, ExhaustiveConstVectorView{y_f[i0]}, 1, i_pol);
3183 jacobian(ii, joker) = transpose(J(joker, 0, joker));
3184 }
3185
3186 // y (must be done last)
3187 apply_iy_unit(yv, iy_unit, ExhaustiveConstVectorView{y_f[i0]}, 1, i_pol);
3188 y[ii] = yv(0, joker);
3189
3190 i0 += n;
3191 }
3192 }
3193
3194 // Other conversions
3195 //--------------------------------------------------------------------------
3196 else {
3197 // Here we take each element of y separately.
3198
3199 Matrix yv(1, 1);
3200 ArrayOfIndex i_pol(1);
3201
3202 for (Index i = 0; i < ny; i++) {
3203 yv(0, 0) = y[i];
3204 i_pol[0] = y_pol[i];
3205
3206 // Jacobian
3207 if (do_j) {
3209 ExhaustiveTensor3View{ExhaustiveMatrixView{jacobian(i, joker)}}, yv, iy_unit, ExhaustiveConstVectorView{y_f[i]}, 1, i_pol);
3210 }
3211
3212 // y (must be done last)
3213 apply_iy_unit(yv, iy_unit, ExhaustiveConstVectorView{y_f[i]}, 1, i_pol);
3214 y[i] = yv(0, 0);
3215 }
3216 }
3217}
3218
3219/* Workspace method: Doxygen documentation will be auto-generated */
3220void y_geo_seriesFromY_geo(Matrix& y_geo_series,
3221 const Matrix& y_geo,
3222 const Vector& sensor_response_f_grid,
3223 const Verbosity&)
3224{
3225 // Sizes
3226 const Index ly = y_geo.nrows();
3227 const Index nchannel = sensor_response_f_grid.nelem();
3228 const Index lseries = ly / nchannel;
3229
3230 ARTS_USER_ERROR_IF (nchannel * lseries != ly,
3231 "Row size of *y_geo* not an even multiple of length of *sensor_response_f_grid*.")
3232
3233 y_geo_series.resize(lseries, y_geo.ncols());
3234
3235 Index i = 0;
3236 for (Index s=0; s<lseries; ++s) {
3237 y_geo_series(s, joker) = y_geo(i, joker);
3238 i += nchannel;
3239 }
3240}
3241
3242/* Workspace method: Doxygen documentation will be auto-generated */
3243void y_geo_swathFromY_geo(Tensor3& y_geo_swath,
3244 const Matrix& y_geo,
3245 const Vector& sensor_response_f_grid,
3246 const Index& npixel,
3247 const Verbosity&)
3248{
3249 // Sizes
3250 const Index ly = y_geo.nrows();
3251 const Index nchannel = sensor_response_f_grid.nelem();
3252 const Index nswath = ly / (nchannel * npixel);
3253
3254 ARTS_USER_ERROR_IF (nchannel * npixel * nswath != ly,
3255 "Row size of *y_geo* does not match given *npixel* and *sensor_response_f_grid*.")
3256
3257 y_geo_swath.resize(nswath, npixel, y_geo.ncols());
3258
3259 Index i = 0;
3260 for (Index s=0; s<nswath; ++s) {
3261 for (Index p=0; p<npixel; ++p) {
3262 y_geo_swath(s, p, joker) = y_geo(i, joker);
3263 i += nchannel;
3264 }
3265 }
3266}
3267
3268/* Workspace method: Doxygen documentation will be auto-generated */
3269void y_seriesFromY(Matrix& y_series,
3270 const Vector& y,
3271 const Vector& y_f,
3272 const Vector& sensor_response_f_grid,
3273 const Index& safe,
3274 const Verbosity&)
3275{
3276 // Sizes
3277 const Index ly = y.nelem();
3278 const Index nchannel = sensor_response_f_grid.nelem();
3279 const Index lseries = ly / nchannel;
3280
3281 ARTS_USER_ERROR_IF (nchannel * lseries != ly,
3282 "Length of *y* not an even multiple of length of *sensor_response_f_grid*.")
3283
3284 y_series.resize(lseries, nchannel);
3285
3286 Index i = 0;
3287 for (Index s=0; s<lseries; ++s) {
3288 for (Index c=0; c<nchannel; ++c) {
3289 if (safe && s > 0) {
3290 ARTS_USER_ERROR_IF (fabs(y_f[i] - y_f[i-nchannel]) > 1,
3291 "At least one channel varies in frequency.")
3292 }
3293 y_series(s, c) = y[i++];
3294 }
3295 }
3296}
3297
3298/* Workspace method: Doxygen documentation will be auto-generated */
3299void y_swathFromY(Tensor3& y_swath,
3300 const Vector& y,
3301 const Vector& y_f,
3302 const Vector& sensor_response_f_grid,
3303 const Index& npixel,
3304 const Index& safe,
3305 const Verbosity&)
3306{
3307 // Sizes
3308 const Index ly = y.nelem();
3309 const Index nchannel = sensor_response_f_grid.nelem();
3310 const Index nswath = ly / (nchannel * npixel);
3311
3312 ARTS_USER_ERROR_IF (nchannel * npixel * nswath != ly,
3313 "Length of *y* does not match given *npixel* and *sensor_response_f_grid*.")
3314
3315 y_swath.resize(nswath, npixel, nchannel);
3316
3317 Index i = 0;
3318 for (Index s=0; s<nswath; ++s) {
3319 for (Index p=0; p<npixel; ++p) {
3320 for (Index c=0; c<nchannel; ++c) {
3321 if (safe && (p > 0 || s > 0)) {
3322 ARTS_USER_ERROR_IF (fabs(y_f[i] - y_f[i-nchannel]) > 1,
3323 "At least one channel varies in frequency.")
3324 }
3325 y_swath(s, p, c) = y[i++];
3326 }
3327 }
3328 }
3329}
Array< Index > ArrayOfIndex
An array of Index.
Definition array.h:119
base max(const Array< base > &x)
Max function.
Definition array.h:128
base min(const Array< base > &x)
Min function.
Definition array.h:144
The global header file for ARTS.
Constants of physical expressions as constexpr.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition arts_omp.cc:29
Header file for helper functions for OpenMP.
void gas_scattering_agendaExecute(Workspace &ws, PropagationMatrix &gas_scattering_coef, TransmissionMatrix &gas_scattering_mat, Vector &gas_scattering_fct_legendre, const Vector &f_grid, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Index gas_scattering_output_type, const Agenda &input_agenda)
Definition auto_md.cc:26834
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, Vector &geo_pos, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const String &iy_unit, const Index cloudbox_on, const Index jacobian_do, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition auto_md.cc:27094
void iy_loop_freqs_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition auto_md.cc:27042
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Vector &f_grid, const Agenda &input_agenda)
Definition auto_md.cc:27497
void iy_independent_beam_approx_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const String &iy_unit, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const Index atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Matrix &z_surface, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Index cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Index jacobian_do, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
Definition auto_md.cc:26956
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
The Agenda class.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
void set_name(const String &s)
Set name of this gridded field.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
Workspace class.
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
Definition cloudbox.cc:443
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
#define ARTS_USER_ERROR(...)
Definition debug.h:153
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
Implementation of gridded fields.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition jacobian.cc:548
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition jacobian.cc:589
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition jacobian.cc:1133
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition jacobian.cc:580
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition jacobian.cc:602
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Definition jacobian.cc:29
Routines for setting up the jacobian.
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition jacobian.h:540
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition jacobian.h:532
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, ArrayOfIndex &mc_source_domain, ArrayOfIndex &mc_scat_order, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &cloudbox_checked, const String &iy_unit, const Index &mc_seed, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Numeric &taustep_limit, const Index &l_mc_scat_order, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
void yCalcAppend(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, ArrayOfRetrievalQuantity &jacobian_quantities, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const ArrayOfRetrievalQuantity &jacobian_quantities_copy, const Index &append_instrument_wfs, const Verbosity &verbosity)
WORKSPACE METHOD: yCalcAppend.
Definition m_rte.cc:2791
void iyLoopFrequencies(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const ArrayOfString &iy_aux_vars, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index &stokes_dim, const Vector &f_grid, const Agenda &iy_loop_freqs_agenda, const Verbosity &)
WORKSPACE METHOD: iyLoopFrequencies.
Definition m_rte.cc:2258
void iyApplyUnit(Matrix &iy, ArrayOfMatrix &iy_aux, const Index &stokes_dim, const Vector &f_grid, const ArrayOfString &iy_aux_vars, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: iyApplyUnit.
Definition m_rte.cc:46
void iyEmissionStandard(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, 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 String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandard.
Definition m_rte.cc:1369
void iyIndependentBeamApproximation(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, GriddedField4 &atm_fields_compact, const Index &iy_id, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const EnergyLevelMap &nlte_field, 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 Matrix &particle_masses, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &iy_agenda_call1, const String &iy_unit, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const Agenda &iy_independent_beam_approx_agenda, const Index &return_atm1d, const Index &skip_vmr, const Index &skip_pnd, const Index &return_masses, const Verbosity &)
WORKSPACE METHOD: iyIndependentBeamApproximation.
Definition m_rte.cc:1825
void iyMC(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const ArrayOfArrayOfSingleScatteringData &scat_data, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Tensor4 &pnd_field, const String &iy_unit, const Numeric &mc_std_err, const Index &mc_max_time, const Index &mc_max_iter, const Index &mc_min_iter, const Numeric &mc_taustep_limit, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyMC.
Definition m_rte.cc:2329
constexpr Numeric SPEED_OF_LIGHT
Definition m_rte.cc:35
void iyClearsky(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, 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 Matrix &z_surface, const Vector &refellipsoid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &cloudbox_on, const Index &gas_scattering_do, const Index &suns_do, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyClearsky.
Definition m_rte.cc:146
void iyCalc(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, Vector &geo_pos, const Index &atmfields_checked, const Index &atmgeom_checked, const ArrayOfString &iy_aux_vars, const Index &iy_id, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda, const Verbosity &)
WORKSPACE METHOD: iyCalc.
Definition m_rte.cc:78
void y_geo_swathFromY_geo(Tensor3 &y_geo_swath, const Matrix &y_geo, const Vector &sensor_response_f_grid, const Index &npixel, const Verbosity &)
WORKSPACE METHOD: y_geo_swathFromY_geo.
Definition m_rte.cc:3243
void y_geo_seriesFromY_geo(Matrix &y_geo_series, const Matrix &y_geo, const Vector &sensor_response_f_grid, const Verbosity &)
WORKSPACE METHOD: y_geo_seriesFromY_geo.
Definition m_rte.cc:3220
void yApplyUnit(Vector &y, Matrix &jacobian, const Vector &y_f, const ArrayOfIndex &y_pol, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: yApplyUnit.
Definition m_rte.cc:3109
void iyReplaceFromAux(Matrix &iy, const ArrayOfMatrix &iy_aux, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const String &aux_var, const Verbosity &)
WORKSPACE METHOD: iyReplaceFromAux.
Definition m_rte.cc:2511
constexpr Numeric PI
Definition m_rte.cc:34
void y_seriesFromY(Matrix &y_series, const Vector &y, const Vector &y_f, const Vector &sensor_response_f_grid, const Index &safe, const Verbosity &)
WORKSPACE METHOD: y_seriesFromY.
Definition m_rte.cc:3269
void ppvar_optical_depthFromPpvar_trans_cumulat(Matrix &ppvar_optical_depth, const Tensor4 &ppvar_trans_cumulat, const Verbosity &)
WORKSPACE METHOD: ppvar_optical_depthFromPpvar_trans_cumulat.
Definition m_rte.cc:2540
void y_swathFromY(Tensor3 &y_swath, const Vector &y, const Vector &y_f, const Vector &sensor_response_f_grid, const Index &npixel, const Index &safe, const Verbosity &)
WORKSPACE METHOD: y_swathFromY.
Definition m_rte.cc:3299
void iyEmissionHybrid(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_pnd, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, 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 String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Ppath &ppath, const Vector &rte_pos2, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Tensor7 &cloudbox_field, const Vector &za_grid, const Index &Naa, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionHybrid.
Definition m_rte.cc:826
void yCalc(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 Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
WORKSPACE METHOD: yCalc.
Definition m_rte.cc:2550
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition messages.h:189
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Index GFIELD4_LAT_GRID
Global constant, Index of the latitude grid in GriddedField4.
constexpr Index GFIELD4_P_GRID
Global constant, Index of the pressure grid in GriddedField4.
constexpr Index GFIELD4_FIELD_NAMES
Global constant, Index of the field names in GriddedField4.
constexpr Index GFIELD4_LON_GRID
Global constant, Index of the longitude grid in GriddedField4.
This file contains declerations of functions of physical character.
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition ppath.cc:524
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition ppath.cc:1460
Propagation path structure and functions.
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstMatrixView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstVectorView &ppath_line_of_sight, const ConstVectorView &ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
Definition rte.cc:1188
void get_stepwise_scattersky_source(StokesVector &Sp, ArrayOfStokesVector &dSp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstTensor7View &cloudbox_field, const ConstVectorView &za_grid, const ConstVectorView &aa_grid, const ConstMatrixView &ppath_line_of_sight, const GridPos &ppath_pressure, const Vector &temperature, const Index &atmosphere_dim, const bool &jacobian_do, const Index &t_interp_order)
Calculates the stepwise scattering source terms.
Definition rte.cc:1305
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition rte.cc:828
void apply_iy_unit2(Tensor3View J, const ConstMatrixView &iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Largely as apply_iy_unit but operates on jacobian data.
Definition rte.cc:162
void iy_transmittance_mult(Tensor3 &iy_trans_total, const ConstTensor3View &iy_trans_old, const ConstTensor3View &iy_trans_new)
Multiplicates iy_transmittance with transmissions.
Definition rte.cc:1755
void apply_iy_unit(MatrixView iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Performs conversion from radiance to other units, as well as applies refractive index to fulfill the ...
Definition rte.cc:89
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const Tensor3 &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const String &iy_unit, const Tensor3 &surface_props_data, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Verbosity &verbosity)
Determines iy of the "background" of a propgation path.
Definition rte.cc:713
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
Definition rte.cc:1051
void rtmethods_unit_conversion(Matrix &iy, ArrayOfTensor3 &diy_dx, Tensor3 &ppvar_iy, const Index &ns, const Index &np, const Vector &f_grid, const Ppath &ppath, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &j_analytical_do, const String &iy_unit)
This function handles the unit conversion to be done at the end of some radiative transfer WSMs.
Definition rte.cc:2086
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Definition rte.cc:1948
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &y_geo, Matrix &jacobian, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity, const Index &mblock_index, const Index &n1y, const Index &j_analytical_do)
Performs calculations for one measurement block, on y-level.
Definition rte.cc:2128
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &ppath_f_grid, const Vector &ppath_magnetic_field, const Vector &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const Vector &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition rte.cc:1116
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
Definition rte.cc:1095
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition rte.cc:38
void get_stepwise_blackbody_radiation(VectorView B, VectorView dB_dT, const ConstVectorView &ppath_f_grid, const Numeric &ppath_temperature, const bool &do_temperature_derivative)
Get the blackbody radiation at propagation path point.
Definition rte.cc:1101
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Definition rte.cc:947
Declaration of functions in rte.cc.
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Header file for special_interp.cc.
An Antenna object used by MCGeneral.
Definition mc_antenna.h:33
void set_pencil_beam()
set_pencil_beam
Definition mc_antenna.cc:80
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Numeric end_lstep
The distance between end pos and the first position in pos.
Vector lstep
The length between ppath points.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Radiation Vector for Stokes dimension 1-4.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
void get_direct_radiation(Workspace &ws, ArrayOfMatrix &direct_radiation, ArrayOfArrayOfTensor3 &ddirect_radiation_dx, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_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 Index &gas_scattering_do, const Index &irradiance_flag, const ArrayOfPpath &sun_ppaths, const ArrayOfSun &suns, const ArrayOfIndex &suns_visible, const Vector &refellipsoid, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Numeric &rte_alonglos_v, const Verbosity &verbosity)
Calculates the transmitted sun radiation at the end position of the ppath.
Definition sun.cc:195
void get_sun_ppaths(Workspace &ws, ArrayOfPpath &sun_ppaths, ArrayOfIndex &suns_visible, ArrayOfVector &sun_rte_los, const Vector &rte_pos, const ArrayOfSun &suns, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Matrix &z_surface, const Vector &refellipsoid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
Calculates the ppath towards the suns from a given position and indicates if sun is visible or not.
Definition sun.cc:371
void get_scattered_sunsource(Workspace &ws, RadiationVector &scattered_sunlight, const Vector &f_grid, const Numeric &p, const Numeric &T, const Vector &vmr, const Matrix &transmitted_sunlight, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Agenda &gas_scattering_agenda)
Calculates the radiance spectrum of sun which is scattered by the atmospheric gases.
Definition sun.cc:45
void get_sun_background(Matrix &iy, Index &suns_visible, const ArrayOfSun &suns, const Ppath &ppath, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &refellipsoid)
Gets the sun background for a given ppath.
Definition sun.cc:96
Declaration of functions in star.cc.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, RadiationVector &J_add, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK, const ArrayOfStokesVector &da, const ArrayOfStokesVector &dS, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Set the stepwise source.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the 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.
Stuff related to the transmission matrix.
Array< RadiationVector > ArrayOfRadiationVector
#define a
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
#define c