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