ARTS 2.5.10 (git: 2f1c442c)
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 "sun.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& suns_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 ArrayOfSun& suns,
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" && suns_do,
247 "If suns are present only iy_unit=\"1\" can be used.");
248
249 ARTS_USER_ERROR_IF(jacobian_quantities.nelem() && (suns_do || gas_scattering_do) , R"--(
250Jacobian calculation are not supported when gas scattering or suns 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 (sun) 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_sunlight(nf, ns);
456
457 if (gas_scattering_do) {
458
459 ArrayOfIndex suns_visible(suns.nelem());
460
461 Numeric minP = min(ppvar_p);
462
463 if (suns_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 (sun-)ppath. The influence of the
468 // uppermost level in view of scattering
469 // is negligible due to the low density.
470 ArrayOfPpath sun_ppaths(suns.nelem());
471 ArrayOfVector sun_rte_los(suns.nelem(), Vector(2));
472
473 get_sun_ppaths(wss,
474 sun_ppaths,
475 suns_visible,
476 sun_rte_los,
477 ppath.pos(ip, joker),
478 suns,
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_sunlight;
493 ArrayOfArrayOfTensor3 dtransmitted_sunlight_dummy(suns.nelem(),ArrayOfTensor3(jacobian_quantities.nelem()));
494
495 get_direct_radiation(wss,
496 transmitted_sunlight,
497 dtransmitted_sunlight_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 sun_ppaths,
519 suns,
520 suns_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 suns to get the total scattered starlight
535 RadiationVector scattered_sunlight_isun(nf, ns);
536
537 for (Index i_sun = 0; i_sun < suns.nelem(); i_sun++) {
538 if (suns_visible[i_sun]) {
539
540 // here we calculate how much incoming sun radiation is scattered
541 //into the direction of the ppath
542 get_scattered_sunsource(wss,
543 scattered_sunlight_isun,
544 f_grid,
545 ppvar_p[ip],
546 ppvar_t[ip],
547 ppvar_vmr(joker, ip),
548 transmitted_sunlight[i_sun],
549 sun_rte_los[i_sun],
550 ppath.los(ip, joker),
551 gas_scattering_agenda);
552
553 scattered_sunlight += scattered_sunlight_isun;
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_sunlight is changed within
588 // stepwise source.
589 stepwise_source(src_rad[ip],
590 dsrc_rad[ip],
591 scattered_sunlight,
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 (suns_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 sun radiation at top of the atmosphere. if sun is not visible
708 // in los, iy_direct_toa will be zero
709 get_sun_background(iy_direct_toa,
710 stars_visible,
711 suns,
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 sun 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(wss)
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,
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.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,
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,
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,
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,
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}
3237
3238/* Workspace method: Doxygen documentation will be auto-generated */
3239void y_geo_seriesFromY_geo(Matrix& y_geo_series,
3240 const Matrix& y_geo,
3241 const Vector& sensor_response_f_grid,
3242 const Verbosity&)
3243{
3244 // Sizes
3245 const Index ly = y_geo.nrows();
3246 const Index nchannel = sensor_response_f_grid.nelem();
3247 const Index lseries = ly / nchannel;
3248
3249 ARTS_USER_ERROR_IF (nchannel * lseries != ly,
3250 "Row size of *y_geo* not an even multiple of length of *sensor_response_f_grid*.")
3251
3252 y_geo_series.resize(lseries, y_geo.ncols());
3253
3254 Index i = 0;
3255 for (Index s=0; s<lseries; ++s) {
3256 y_geo_series(s, joker) = y_geo(i, joker);
3257 i += nchannel;
3258 }
3259}
3260
3261/* Workspace method: Doxygen documentation will be auto-generated */
3262void y_geo_swathFromY_geo(Tensor3& y_geo_swath,
3263 const Matrix& y_geo,
3264 const Vector& sensor_response_f_grid,
3265 const Index& npixel,
3266 const Verbosity&)
3267{
3268 // Sizes
3269 const Index ly = y_geo.nrows();
3270 const Index nchannel = sensor_response_f_grid.nelem();
3271 const Index nswath = ly / (nchannel * npixel);
3272
3273 ARTS_USER_ERROR_IF (nchannel * npixel * nswath != ly,
3274 "Row size of *y_geo* does not match given *npixel* and *sensor_response_f_grid*.")
3275
3276 y_geo_swath.resize(nswath, npixel, y_geo.ncols());
3277
3278 Index i = 0;
3279 for (Index s=0; s<nswath; ++s) {
3280 for (Index p=0; p<npixel; ++p) {
3281 y_geo_swath(s, p, joker) = y_geo(i, joker);
3282 i += nchannel;
3283 }
3284 }
3285}
3286
3287/* Workspace method: Doxygen documentation will be auto-generated */
3288void y_seriesFromY(Matrix& y_series,
3289 const Vector& y,
3290 const Vector& y_f,
3291 const Vector& sensor_response_f_grid,
3292 const Index& safe,
3293 const Verbosity&)
3294{
3295 // Sizes
3296 const Index ly = y.nelem();
3297 const Index nchannel = sensor_response_f_grid.nelem();
3298 const Index lseries = ly / nchannel;
3299
3300 ARTS_USER_ERROR_IF (nchannel * lseries != ly,
3301 "Length of *y* not an even multiple of length of *sensor_response_f_grid*.")
3302
3303 y_series.resize(lseries, nchannel);
3304
3305 Index i = 0;
3306 for (Index s=0; s<lseries; ++s) {
3307 for (Index c=0; c<nchannel; ++c) {
3308 if (safe && s > 0) {
3309 ARTS_USER_ERROR_IF (fabs(y_f[i] - y_f[i-nchannel]) > 1,
3310 "At least one channel varies in frequency.")
3311 }
3312 y_series(s, c) = y[i++];
3313 }
3314 }
3315}
3316
3317/* Workspace method: Doxygen documentation will be auto-generated */
3318void y_swathFromY(Tensor3& y_swath,
3319 const Vector& y,
3320 const Vector& y_f,
3321 const Vector& sensor_response_f_grid,
3322 const Index& npixel,
3323 const Index& safe,
3324 const Verbosity&)
3325{
3326 // Sizes
3327 const Index ly = y.nelem();
3328 const Index nchannel = sensor_response_f_grid.nelem();
3329 const Index nswath = ly / (nchannel * npixel);
3330
3331 ARTS_USER_ERROR_IF (nchannel * npixel * nswath != ly,
3332 "Length of *y* does not match given *npixel* and *sensor_response_f_grid*.")
3333
3334 y_swath.resize(nswath, npixel, nchannel);
3335
3336 Index i = 0;
3337 for (Index s=0; s<nswath; ++s) {
3338 for (Index p=0; p<npixel; ++p) {
3339 for (Index c=0; c<nchannel; ++c) {
3340 if (safe && (p > 0 || s > 0)) {
3341 ARTS_USER_ERROR_IF (fabs(y_f[i] - y_f[i-nchannel]) > 1,
3342 "At least one channel varies in frequency.")
3343 }
3344 y_swath(s, p, c) = y[i++];
3345 }
3346 }
3347 }
3348}
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:136
max
base max(const Array< base > &x)
Max function.
Definition: array.h:145
min
base min(const Array< base > &x)
Min function.
Definition: array.h:161
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:25058
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:25318
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:25266
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:25721
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:25180
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:1079
ConstMatrixView::ncols
Index ncols() const noexcept
Definition: matpackI.h:1080
ConstTensor3View::empty
bool empty() const noexcept
Definition: matpackIII.h:157
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
ConstTensor4View::empty
bool empty() const noexcept
Definition: matpackIV.h:152
ConstTensor4View::nbooks
Index nbooks() const noexcept
Definition: matpackIV.h:143
ConstTensor4View::npages
Index npages() const noexcept
Definition: matpackIV.h:144
ConstTensor7View::ncols
Index ncols() const noexcept
Definition: matpackVII.h:164
ConstTensor7View::npages
Index npages() const noexcept
Definition: matpackVII.h:162
ConstTensor7View::nrows
Index nrows() const noexcept
Definition: matpackVII.h:163
ConstTensor7View::nlibraries
Index nlibraries() const noexcept
Definition: matpackVII.h:158
ConstTensor7View::nvitrines
Index nvitrines() const noexcept
Definition: matpackVII.h:159
ConstTensor7View::nshelves
Index nshelves() const noexcept
Definition: matpackVII.h:160
ConstTensor7View::nbooks
Index nbooks() const noexcept
Definition: matpackVII.h:161
ConstVectorView::nelem
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
ConstVectorView::empty
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:536
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:1188
Matrix
The Matrix class.
Definition: matpackI.h:1285
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
Range
The range class.
Definition: matpackI.h:160
Range::get_start
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:341
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1061
Tensor3
The Tensor3 class.
Definition: matpackIII.h:352
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
Tensor4
The Tensor4 class.
Definition: matpackIV.h:435
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:2407
Vector
The Vector class.
Definition: matpackI.h:910
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:102
ARTS_USER_ERROR
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
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:555
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:547
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
yCalcAppend
void yCalcAppend(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, ArrayOfRetrievalQuantity &jacobian_quantities, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const ArrayOfRetrievalQuantity &jacobian_quantities_copy, const Index &append_instrument_wfs, const Verbosity &verbosity)
WORKSPACE METHOD: yCalcAppend.
Definition: m_rte.cc:2810
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
iyClearsky
void iyClearsky(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Matrix &z_surface, const Vector &refellipsoid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &cloudbox_on, const Index &gas_scattering_do, const Index &suns_do, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Agenda &gas_scattering_agenda, const Agenda &ppath_step_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyClearsky.
Definition: m_rte.cc:165
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
y_geo_swathFromY_geo
void y_geo_swathFromY_geo(Tensor3 &y_geo_swath, const Matrix &y_geo, const Vector &sensor_response_f_grid, const Index &npixel, const Verbosity &)
WORKSPACE METHOD: y_geo_swathFromY_geo.
Definition: m_rte.cc:3262
y_geo_seriesFromY_geo
void y_geo_seriesFromY_geo(Matrix &y_geo_series, const Matrix &y_geo, const Vector &sensor_response_f_grid, const Verbosity &)
WORKSPACE METHOD: y_geo_seriesFromY_geo.
Definition: m_rte.cc:3239
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
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
PI
constexpr Numeric PI
Definition: m_rte.cc:53
y_seriesFromY
void y_seriesFromY(Matrix &y_series, const Vector &y, const Vector &y_f, const Vector &sensor_response_f_grid, const Index &safe, const Verbosity &)
WORKSPACE METHOD: y_seriesFromY.
Definition: m_rte.cc:3288
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
y_swathFromY
void y_swathFromY(Tensor3 &y_swath, const Vector &y, const Vector &y_f, const Vector &sensor_response_f_grid, const Index &npixel, const Index &safe, const Verbosity &)
WORKSPACE METHOD: y_swathFromY.
Definition: m_rte.cc:3318
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, 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:479
abs
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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:1433
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:1510
swap
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
Definition: matpack_complex.cc:761
joker
const Joker joker
messages.h
Declarations having to do with the four output streams.
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:206
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
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
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:974
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:1207
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:1135
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:732
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:1324
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:847
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:181
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:1774
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:108
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:1070
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:2105
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:1967
yCalc_mblock_loop_body
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &y_geo, Matrix &jacobian, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity, const Index &mblock_index, const Index &n1y, const Index &j_analytical_do)
Performs calculations for one measurement block, on y-level.
Definition: rte.cc:2147
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:1114
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:57
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:1120
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:966
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.
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
get_direct_radiation
void get_direct_radiation(Workspace &ws, ArrayOfMatrix &direct_radiation, ArrayOfArrayOfTensor3 &ddirect_radiation_dx, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &gas_scattering_do, const Index &irradiance_flag, const ArrayOfPpath &sun_ppaths, const ArrayOfSun &suns, const ArrayOfIndex &suns_visible, const Vector &refellipsoid, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &gas_scattering_agenda, const Numeric &rte_alonglos_v, const Verbosity &verbosity)
Calculates the transmitted sun radiation at the end position of the ppath.
Definition: sun.cc:215
get_sun_ppaths
void get_sun_ppaths(Workspace &ws, ArrayOfPpath &sun_ppaths, ArrayOfIndex &suns_visible, ArrayOfVector &sun_rte_los, const Vector &rte_pos, const ArrayOfSun &suns, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Matrix &z_surface, const Vector &refellipsoid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &ppath_step_agenda, const Verbosity &verbosity)
Calculates the ppath towards the suns from a given position and indicates if sun is visible or not.
Definition: sun.cc:391
get_scattered_sunsource
void get_scattered_sunsource(Workspace &ws, RadiationVector &scattered_sunlight, const Vector &f_grid, const Numeric &p, const Numeric &T, const Vector &vmr, const Matrix &transmitted_sunlight, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Agenda &gas_scattering_agenda)
Calculates the radiance spectrum of sun which is scattered by the atmospheric gases.
Definition: sun.cc:65
get_sun_background
void get_sun_background(Matrix &iy, Index &suns_visible, const ArrayOfSun &suns, const Ppath &ppath, const Vector &f_grid, const Index &stokes_dim, const Index &atmosphere_dim, const Vector &refellipsoid)
Gets the sun background for a given ppath.
Definition: sun.cc:116
sun.h
Declaration of functions in star.cc.
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:2001
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:1803
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:1943
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:1785
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
c
#define c