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;
35inline constexpr Numeric SPEED_OF_LIGHT=Constant::speed_of_light;
36using GriddedFieldGrids::GFIELD4_FIELD_NAMES;
37using GriddedFieldGrids::GFIELD4_P_GRID;
38using GriddedFieldGrids::GFIELD4_LAT_GRID;
39using GriddedFieldGrids::GFIELD4_LON_GRID;
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 */
78void iyCalc(Workspace& ws,
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
118 iy_main_agendaExecute(ws,
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 */
146void iyClearsky(
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 {
286 ARTS_USER_ERROR (
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));
302 ArrayOfArrayOfRadiationVector dlvl_rad(
303 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
304
305 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
306 ArrayOfArrayOfRadiationVector dsrc_rad(
307 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
308
309 ArrayOfTransmissionMatrix lyr_tra(np, TransmissionMatrix(nf, ns));
310 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
311 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
312 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
313 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
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
350 get_ppath_f(
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 }
373 FOR_ANALYTICAL_JACOBIANS_DO(
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
394 WorkspaceOmpParallelCopyGuard wss{ws};;
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 {
404 get_stepwise_blackbody_radiation(
405 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
406
407 Index lte;
408 get_stepwise_clearsky_propmat(wss,
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)
426 adapt_stepwise_partial_derivatives(dK_dx[ip],
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
476 get_direct_radiation(wss,
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
523 get_scattered_sunsource(wss,
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
544 gas_scattering_agendaExecute(wss,
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 =
641 cumulative_transmission(lyr_tra, CumulativeTransmission::Forward);
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
656 get_iy_of_background(ws,
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,
746 RadiativeTransferSolver::Emission);
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,
772 RadiativeTransferSolver::LinearWeightedEmission);
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) {
792 rtmethods_jacobian_finalisation(ws,
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) {
812 rtmethods_unit_conversion(iy,
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 */
826void iyEmissionHybrid(Workspace& ws,
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;
887 iyEmissionStandard(ws,
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));
1038 ArrayOfArrayOfRadiationVector dlvl_rad(
1039 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1040 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
1041 ArrayOfArrayOfRadiationVector dsrc_rad(
1042 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1043
1044 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
1045 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
1046 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
1047 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
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
1087 get_ppath_f(
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);
1122 FOR_ANALYTICAL_JACOBIANS_DO(
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++) {
1139 get_stepwise_blackbody_radiation(
1140 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
1141
1142 get_stepwise_clearsky_propmat(ws,
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)
1160 adapt_stepwise_partial_derivatives(dK_this_dx,
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) {
1170 get_stepwise_scattersky_propmat(a,
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 //
1193 get_stepwise_scattersky_source(Sp,
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 =
1257 cumulative_transmission(lyr_tra, CumulativeTransmission::Forward);
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
1273 get_iy_of_background(ws,
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,
1321 RadiativeTransferSolver::Emission);
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)
1337 rtmethods_jacobian_finalisation(ws,
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)
1356 rtmethods_unit_conversion(iy,
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 */
1369void iyEmissionStandard(
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 {
1463 ARTS_USER_ERROR (
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));
1477 ArrayOfArrayOfRadiationVector dlvl_rad(
1478 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1479
1480 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
1481 ArrayOfArrayOfRadiationVector dsrc_rad(
1482 np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns)));
1483
1484 ArrayOfTransmissionMatrix lyr_tra(np, TransmissionMatrix(nf, ns));
1485 ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
1486 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
1487 ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
1488 np, ArrayOfTransmissionMatrix(nq, TransmissionMatrix(nf, ns)));
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
1533 get_ppath_f(
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 }
1558 FOR_ANALYTICAL_JACOBIANS_DO(
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
1569 WorkspaceOmpParallelCopyGuard wss{ws};
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 {
1577 get_stepwise_blackbody_radiation(
1578 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
1579
1580 Index lte;
1581 get_stepwise_clearsky_propmat(wss,
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)
1599 adapt_stepwise_partial_derivatives(dK_dx[ip],
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 =
1682 cumulative_transmission(lyr_tra, CumulativeTransmission::Forward);
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
1697 get_iy_of_background(ws,
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,
1746 RadiativeTransferSolver::Emission);
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,
1771 RadiativeTransferSolver::LinearWeightedEmission);
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) {
1791 rtmethods_jacobian_finalisation(ws,
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) {
1811 rtmethods_unit_conversion(iy,
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 */
1825void iyIndependentBeamApproximation(Workspace& ws,
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 //
1905 ppath_agendaExecute(ws,
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 //
1933 ppath_agendaExecute(ws,
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);
2033 interp_atmfield_by_itw(
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);
2038 interp_atmfield_by_itw(
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 //
2147 iy_independent_beam_approx_agendaExecute(ws,
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 */
2258void iyLoopFrequencies(Workspace& ws,
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
2287 iy_loop_freqs_agendaExecute(ws,
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 */
2329void iyMC(Workspace& ws,
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 {
2398 ARTS_USER_ERROR (
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) {
2421 WorkspaceOmpParallelCopyGuard wss{ws};
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 */
2540void ppvar_optical_depthFromPpvar_trans_cumulat(
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 */
2550void yCalc(Workspace& ws,
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) {
2584 CREATE_OUT3;
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
2666 WorkspaceOmpParallelCopyGuard wss{ws};
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
2673 yCalc_mblock_loop_body(failed,
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
2716 yCalc_mblock_loop_body(failed,
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 */
2791void yCalcAppend(Workspace& ws,
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;
2831 jac_ranges_indices(
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 {
3005 ARTS_USER_ERROR (
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 {
3018 ARTS_USER_ERROR (
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) {
3074 ARTS_USER_ERROR (
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) {
3208 apply_iy_unit2(
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}
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition array.h:119
max
base max(const Array< base > &x)
Max function.
Definition array.h:128
min
base min(const Array< base > &x)
Min function.
Definition array.h:144
arts.h
The global header file for ARTS.
arts_constants.h
Constants of physical expressions as constexpr.
arts_omp_get_max_threads
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition arts_omp.cc:29
arts_omp.h
Header file for helper functions for OpenMP.
gas_scattering_agendaExecute
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
iy_main_agendaExecute
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
iy_loop_freqs_agendaExecute
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
ppath_agendaExecute
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
iy_independent_beam_approx_agendaExecute
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
auto_md.h
chk_latlon_true
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
Definition check_input.cc:1484
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition check_input.cc:67
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition check_input.cc:92
check_input.h
Agenda
The Agenda class.
Definition agenda_class.h:52
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
GriddedField4::resize
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
Definition gridded_fields.h:425
GriddedField::set_name
void set_name(const String &s)
Set name of this gridded field.
Definition gridded_fields.h:267
GriddedField::set_grid_name
void set_grid_name(Index i, const String &s)
Set grid name.
Definition gridded_fields.h:260
GriddedField::set_grid
void set_grid(Index i, const Vector &g)
Set a numeric grid.
Definition gridded_fields.cc:139
Workspace
Workspace class.
Definition workspace_ng.h:36
my_basic_string< char >
is_gp_inside_cloudbox
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
mult
void mult(MatrixView C, ConstMatrixView A, const Block &B)
Definition covariance_matrix.cc:22
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR
#define ARTS_USER_ERROR(...)
Definition debug.h:153
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
gridded_fields.h
Implementation of gridded fields.
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition interpolation.cc:955
get_pointers_for_analytical_species
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition jacobian.cc:548
get_standard_starting_diy_dx
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
do_temperature_jacobian
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition jacobian.cc:1133
get_standard_diy_dpath
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
get_pointers_for_scat_species
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
jac_ranges_indices
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
jacobian.h
Routines for setting up the jacobian.
FOR_ANALYTICAL_JACOBIANS_DO2
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition jacobian.h:540
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition jacobian.h:532
MCSetSeedFromTime
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
Definition m_montecarlo.cc:1084
MCGeneral
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.
Definition m_montecarlo.cc:72
yCalcAppend
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
iyLoopFrequencies
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
iyApplyUnit
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
iyEmissionStandard
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
iyIndependentBeamApproximation
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
iyMC
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
SPEED_OF_LIGHT
constexpr Numeric SPEED_OF_LIGHT
Definition m_rte.cc:35
iyClearsky
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
iyCalc
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
y_geo_swathFromY_geo
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
y_geo_seriesFromY_geo
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
yApplyUnit
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
iyReplaceFromAux
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
PI
constexpr Numeric PI
Definition m_rte.cc:34
y_seriesFromY
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
ppvar_optical_depthFromPpvar_trans_cumulat
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
y_swathFromY
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
iyEmissionHybrid
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
yCalc
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
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition math_funcs.cc:208
math_funcs.h
messages.h
Declarations having to do with the four output streams.
CREATE_OUT3
#define CREATE_OUT3
Definition messages.h:189
Constant::pi
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
Definition arts_constants.h:57
Constant::speed_of_light
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
Definition arts_constants.h:116
GriddedFieldGrids::GFIELD4_LAT_GRID
constexpr Index GFIELD4_LAT_GRID
Global constant, Index of the latitude grid in GriddedField4.
Definition gridded_fields.h:79
GriddedFieldGrids::GFIELD4_P_GRID
constexpr Index GFIELD4_P_GRID
Global constant, Index of the pressure grid in GriddedField4.
Definition gridded_fields.h:67
GriddedFieldGrids::GFIELD4_FIELD_NAMES
constexpr Index GFIELD4_FIELD_NAMES
Global constant, Index of the field names in GriddedField4.
Definition gridded_fields.h:55
GriddedFieldGrids::GFIELD4_LON_GRID
constexpr Index GFIELD4_LON_GRID
Global constant, Index of the longitude grid in GriddedField4.
Definition gridded_fields.h:91
physics_funcs.h
This file contains declerations of functions of physical character.
error_if_limb_ppath
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition ppath.cc:524
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition ppath.cc:1460
ppath.h
Propagation path structure and functions.
get_stepwise_scattersky_propmat
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
get_stepwise_scattersky_source
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
get_ppath_atmvars
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
apply_iy_unit2
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
iy_transmittance_mult
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
apply_iy_unit
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
get_iy_of_background
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
get_ppath_f
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
rtmethods_unit_conversion
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
rtmethods_jacobian_finalisation
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
yCalc_mblock_loop_body
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
get_stepwise_clearsky_propmat
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
get_rowindex_for_mblock
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
adapt_stepwise_partial_derivatives
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
get_stepwise_blackbody_radiation
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
get_ppath_cloudvars
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
rte.h
Declaration of functions in rte.cc.
interp_atmfield_gp2itw
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.
Definition special_interp.cc:58
itw2p
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Definition special_interp.cc:701
interp_cloudfield_gp2itw
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 ...
Definition special_interp.cc:148
interp_atmfield_by_itw
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.
Definition special_interp.cc:84
special_interp.h
Header file for special_interp.cc.
MCAntenna
An Antenna object used by MCGeneral.
Definition mc_antenna.h:33
MCAntenna::set_pencil_beam
void set_pencil_beam()
set_pencil_beam
Definition mc_antenna.cc:80
Ppath
The structure to describe a propagation path and releated quantities.
Definition ppath_struct.h:17
Ppath::los
Matrix los
Line-of-sight at each ppath point.
Definition ppath_struct.h:35
Ppath::gp_lon
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition ppath_struct.h:55
Ppath::np
Index np
Number of points describing the ppath.
Definition ppath_struct.h:21
Ppath::pos
Matrix pos
The distance between start pos and the last position in pos.
Definition ppath_struct.h:33
Ppath::gp_lat
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition ppath_struct.h:53
Ppath::end_lstep
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition ppath_struct.h:45
Ppath::lstep
Vector lstep
The length between ppath points.
Definition ppath_struct.h:39
Ppath::gp_p
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition ppath_struct.h:51
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition transmissionmatrix.h:480
TransmissionMatrix
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
Definition transmissionmatrix.h:148
get_direct_radiation
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
get_sun_ppaths
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
get_scattered_sunsource
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
get_sun_background
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
sun.h
Declaration of functions in star.cc.
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition transmissionmatrix.cc:1983
stepwise_source
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.
Definition transmissionmatrix.cc:1785
update_radiation_vector
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.
Definition transmissionmatrix.cc:1925
stepwise_transmission
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.
Definition transmissionmatrix.cc:1767
transmissionmatrix.h
Stuff related to the transmission matrix.
ArrayOfRadiationVector
Array< RadiationVector > ArrayOfRadiationVector
Definition transmissionmatrix.h:773
a
#define a
ArrayOfTransmissionMatrix
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
Definition transmissionmatrix.h:769
c
#define c