ARTS 2.5.4 (git: 4c0d3b4d)
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_omp.h"
34#include "auto_md.h"
35#include "check_input.h"
36#include "geodetic.h"
37#include "jacobian.h"
38#include "logic.h"
39#include "math_funcs.h"
40#include "messages.h"
41#include "montecarlo.h"
42#include "physics_funcs.h"
43#include "ppath.h"
44#include "rte.h"
45#include "special_interp.h"
46#include "transmissionmatrix.h"
47#include <cmath>
48#include <stdexcept>
49
50extern const Numeric PI;
51extern const Numeric SPEED_OF_LIGHT;
52extern const Index GFIELD4_FIELD_NAMES;
53extern const Index GFIELD4_P_GRID;
54extern const Index GFIELD4_LAT_GRID;
55extern const Index GFIELD4_LON_GRID;
56
57/*===========================================================================
58 === Workspace methods
59 ===========================================================================*/
60
61/* Workspace method: Doxygen documentation will be auto-generated */
63 ArrayOfMatrix& iy_aux,
64 const Index& stokes_dim,
65 const Vector& f_grid,
66 const ArrayOfString& iy_aux_vars,
67 const String& iy_unit,
68 const Verbosity&) {
69 ARTS_USER_ERROR_IF (iy_unit == "1",
70 "No need to use this method with *iy_unit* = \"1\".");
71
72 ARTS_USER_ERROR_IF (max(iy(joker, 0)) > 1e-3,
73 "The spectrum matrix *iy* is required to have original radiance\n"
74 "unit, but this seems not to be the case. This as a value above\n"
75 "1e-3 is found in *iy*.")
76
77 // Polarisation index variable
78 ArrayOfIndex i_pol(stokes_dim);
79 for (Index is = 0; is < stokes_dim; is++) {
80 i_pol[is] = is + 1;
81 }
82
83 apply_iy_unit(iy, iy_unit, f_grid, 1, i_pol);
84
85 for (Index i = 0; i < iy_aux_vars.nelem(); i++) {
86 if (iy_aux_vars[i] == "iy" || iy_aux_vars[i] == "Error" ||
87 iy_aux_vars[i] == "Error (uncorrelated)") {
88 apply_iy_unit(iy_aux[i], iy_unit, f_grid, 1, i_pol);
89 }
90 }
91}
92
93/* Workspace method: Doxygen documentation will be auto-generated */
95 Matrix& iy,
96 ArrayOfMatrix& iy_aux,
97 Ppath& ppath,
98 const Index& atmfields_checked,
99 const Index& atmgeom_checked,
100 const ArrayOfString& iy_aux_vars,
101 const Index& iy_id,
102 const Index& cloudbox_on,
103 const Index& cloudbox_checked,
104 const Index& scat_data_checked,
105 const Vector& f_grid,
106 const EnergyLevelMap& nlte_field,
107 const Vector& rte_pos,
108 const Vector& rte_los,
109 const Vector& rte_pos2,
110 const String& iy_unit,
111 const Agenda& iy_main_agenda,
112 const Verbosity&) {
113 // Basics
114 //
115 ARTS_USER_ERROR_IF (atmfields_checked != 1,
116 "The atmospheric fields must be flagged to have\n"
117 "passed a consistency check (atmfields_checked=1).");
118 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
119 "The atmospheric geometry must be flagged to have\n"
120 "passed a consistency check (atmgeom_checked=1).");
121 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
122 "The cloudbox must be flagged to have\n"
123 "passed a consistency check (cloudbox_checked=1).");
124 if (cloudbox_on)
125 ARTS_USER_ERROR_IF (scat_data_checked != 1,
126 "The scattering data must be flagged to have\n"
127 "passed a consistency check (scat_data_checked=1).");
128
129 // iy_transmittance is just input and can be left empty for first call
130 Tensor3 iy_transmittance(0, 0, 0);
131 ArrayOfTensor3 diy_dx;
132
134 iy,
135 iy_aux,
136 ppath,
137 diy_dx,
138 1,
139 iy_transmittance,
140 iy_aux_vars,
141 iy_id,
142 iy_unit,
143 cloudbox_on,
144 0,
145 f_grid,
146 nlte_field,
147 rte_pos,
148 rte_los,
149 rte_pos2,
150 iy_main_agenda);
151
152 // Don't allow NaNs (should suffice to check first stokes element)
153 for (Index i = 0; i < iy.nrows(); i++) {
155 "One or several NaNs found in *iy*.");
156 }
157}
158
159/* Workspace method: Doxygen documentation will be auto-generated */
161 Matrix& iy,
162 ArrayOfMatrix& iy_aux,
163 ArrayOfTensor3& diy_dx,
164 Vector& ppvar_p,
165 Vector& ppvar_t,
166 EnergyLevelMap& ppvar_nlte,
167 Matrix& ppvar_vmr,
168 Matrix& ppvar_wind,
169 Matrix& ppvar_mag,
170 Matrix& ppvar_pnd,
171 Matrix& ppvar_f,
172 Tensor3& ppvar_iy,
173 Tensor4& ppvar_trans_cumulat,
174 Tensor4& ppvar_trans_partial,
175 const Index& iy_id,
176 const Index& stokes_dim,
177 const Vector& f_grid,
178 const Index& atmosphere_dim,
179 const Vector& p_grid,
180 const Tensor3& t_field,
181 const EnergyLevelMap& nlte_field,
182 const Tensor4& vmr_field,
183 const ArrayOfArrayOfSpeciesTag& abs_species,
184 const Tensor3& wind_u_field,
185 const Tensor3& wind_v_field,
186 const Tensor3& wind_w_field,
187 const Tensor3& mag_u_field,
188 const Tensor3& mag_v_field,
189 const Tensor3& mag_w_field,
190 const Index& cloudbox_on,
191 const ArrayOfIndex& cloudbox_limits,
192 const Tensor4& pnd_field,
193 const ArrayOfTensor4& dpnd_field_dx,
194 const ArrayOfString& scat_species,
195 const ArrayOfArrayOfSingleScatteringData& scat_data,
196 const String& iy_unit,
197 const ArrayOfString& iy_aux_vars,
198 const Index& jacobian_do,
199 const ArrayOfRetrievalQuantity& jacobian_quantities,
200 const Agenda& propmat_clearsky_agenda,
201 const Agenda& water_p_eq_agenda,
202 const String& rt_integration_option,
203 const Agenda& iy_main_agenda,
204 const Agenda& iy_space_agenda,
205 const Agenda& iy_surface_agenda,
206 const Agenda& iy_cloudbox_agenda,
207 const Index& iy_agenda_call1,
208 const Tensor3& iy_transmittance,
209 const Ppath& ppath,
210 const Vector& rte_pos2,
211 const Numeric& rte_alonglos_v,
212 const Tensor3& surface_props_data,
213 const Tensor7& cloudbox_field,
214 const Vector& za_grid,
215 const Index& Naa,
216 const Index& t_interp_order,
217 const Verbosity& verbosity) {
218 // If cloudbox off, switch to use clearsky method
219 if (!cloudbox_on) {
220 Tensor4 dummy;
222 iy,
223 iy_aux,
224 diy_dx,
225 ppvar_p,
226 ppvar_t,
227 ppvar_nlte,
228 ppvar_vmr,
229 ppvar_wind,
230 ppvar_mag,
231 ppvar_f,
232 ppvar_iy,
233 ppvar_trans_cumulat,
234 ppvar_trans_partial,
235 iy_id,
236 stokes_dim,
237 f_grid,
238 atmosphere_dim,
239 p_grid,
240 t_field,
241 nlte_field,
242 vmr_field,
243 abs_species,
244 wind_u_field,
245 wind_v_field,
246 wind_w_field,
247 mag_u_field,
248 mag_v_field,
249 mag_w_field,
250 cloudbox_on,
251 iy_unit,
252 iy_aux_vars,
253 jacobian_do,
254 jacobian_quantities,
255 ppath,
256 rte_pos2,
257 propmat_clearsky_agenda,
258 water_p_eq_agenda,
259 rt_integration_option,
260 iy_main_agenda,
261 iy_space_agenda,
262 iy_surface_agenda,
263 iy_cloudbox_agenda,
264 iy_agenda_call1,
265 iy_transmittance,
266 rte_alonglos_v,
267 surface_props_data,
268 verbosity);
269 return;
270 }
271 // Init Jacobian quantities?
272 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<1>(jacobian_quantities) : 0;
273
274 // Some basic sizes
275 const Index nf = f_grid.nelem();
276 const Index ns = stokes_dim;
277 const Index np = ppath.np;
278 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
279
280 // Radiative background index
281 const Index rbi = ppath_what_background(ppath);
282
283 // Throw error if unsupported features are requested
284 if (atmosphere_dim != 1)
285 throw runtime_error(
286 "With cloudbox on, this method handles only 1D calculations.");
287 if (Naa < 3) throw runtime_error("Naa must be > 2.");
288 if (jacobian_do)
289 if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
290 throw runtime_error(
291 "*dpnd_field_dx* not properly initialized:\n"
292 "Number of elements in dpnd_field_dx must be equal number of jacobian"
293 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
294 " is calculated/set.");
295 if (rbi < 1 || rbi > 9)
296 throw runtime_error(
297 "ppath.background is invalid. Check your "
298 "calculation of *ppath*?");
299 if (rbi == 3 || rbi == 4)
300 throw runtime_error(
301 "The propagation path ends inside or at boundary of "
302 "the cloudbox.\nFor this method, *ppath* must be "
303 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
304 // iy_aux_vars checked below
305 // Checks of i_field
306 if (cloudbox_field.ncols() != stokes_dim)
307 throw runtime_error(
308 "Obtained *cloudbox_field* number of Stokes elements inconsistent with "
309 "*stokes_dim*.");
310 if (cloudbox_field.nrows() != 1)
311 throw runtime_error(
312 "Obtained *cloudbox_field* has wrong number of azimuth angles.");
313 if (cloudbox_field.npages() != za_grid.nelem())
314 throw runtime_error(
315 "Obtained *cloudbox_field* number of zenith angles inconsistent with "
316 "*za_grid*.");
317 if (cloudbox_field.nbooks() != 1)
318 throw runtime_error(
319 "Obtained *cloudbox_field* has wrong number of longitude points.");
320 if (cloudbox_field.nshelves() != 1)
321 throw runtime_error(
322 "Obtained *cloudbox_field* has wrong number of latitude points.");
323 if (cloudbox_field.nvitrines() != cloudbox_limits[1] - cloudbox_limits[0] + 1)
324 throw runtime_error(
325 "Obtained *cloudbox_field* number of pressure points inconsistent with "
326 "*cloudbox_limits*.");
327 if (cloudbox_field.nlibraries() != nf)
328 throw runtime_error(
329 "Obtained *cloudbox_field* number of frequency points inconsistent with "
330 "*f_grid*.");
331
332 // Set diy_dpath if we are doing are doing jacobian calculations
333 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
334
335 // Set the species pointers if we are doing jacobian
336 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
337
338 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
339 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
340
341 // Checks that the scattering species are treated correctly if their derivatives are needed (we can here discard the Array)
342 if (j_analytical_do and iy_agenda_call1) get_pointers_for_scat_species(jacobian_quantities, scat_species, cloudbox_on);
343
344 // Init iy_aux and fill where possible
345 const Index naux = iy_aux_vars.nelem();
346 iy_aux.resize(naux);
347 //
348 for (Index i = 0; i < naux; i++) {
349 iy_aux[i].resize(nf, ns);
350
351 if (iy_aux_vars[i] == "Optical depth") { /*pass*/
352 } // Filled below
353 else if (iy_aux_vars[i] == "Radiative background")
354 iy_aux[i] = (Numeric)min((Index)2, rbi - 1);
355 else {
356 ostringstream os;
357 os << "The only allowed strings in *iy_aux_vars* are:\n"
358 << " \"Radiative background\"\n"
359 << " \"Optical depth\"\n"
360 << "but you have selected: \"" << iy_aux_vars[i] << "\"";
361 throw runtime_error(os.str());
362 }
363 }
364
365 // Get atmospheric and radiative variables along the propagation path
366 ppvar_trans_cumulat.resize(np, nf, ns, ns);
367 ppvar_trans_partial.resize(np, nf, ns, ns);
368 ppvar_iy.resize(nf, ns, np);
369
371 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
374 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
377
382
383 ArrayOfIndex clear2cloudy;
384 //
385 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
386 ppvar_p.resize(0);
387 ppvar_t.resize(0);
388 ppvar_vmr.resize(0, 0);
389 ppvar_wind.resize(0, 0);
390 ppvar_mag.resize(0, 0);
391 ppvar_f.resize(0, 0);
392 ppvar_trans_cumulat = 0;
393 ppvar_trans_partial = 0;
394 for (Index iv = 0; iv < nf; iv++) {
395 for (Index is = 0; is < ns; is++) {
396 ppvar_trans_cumulat(0,iv,is,is) = 1;
397 ppvar_trans_partial(0,iv,is,is) = 1;
398 }
399 }
400 } else {
401 // Basic atmospheric variables
402 get_ppath_atmvars(ppvar_p,
403 ppvar_t,
404 ppvar_nlte,
405 ppvar_vmr,
406 ppvar_wind,
407 ppvar_mag,
408 ppath,
409 atmosphere_dim,
410 p_grid,
411 t_field,
412 nlte_field,
413 vmr_field,
414 wind_u_field,
415 wind_v_field,
416 wind_w_field,
417 mag_u_field,
418 mag_v_field,
419 mag_w_field);
420
422 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
423
424 // here, the cloudbox is on, ie we don't need to check and branch this here
425 // anymore.
426 ArrayOfMatrix ppvar_dpnd_dx;
427 //
428 get_ppath_cloudvars(clear2cloudy,
429 ppvar_pnd,
430 ppvar_dpnd_dx,
431 ppath,
432 atmosphere_dim,
433 cloudbox_limits,
434 pnd_field,
435 dpnd_field_dx);
436
437 // Size radiative variables always used
438 Vector B(nf);
439 PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
440 StokesVector a(nf, ns), S(nf, ns), Sp(nf, ns);
441 ArrayOfIndex lte(np);
442
443 // Init variables only used if analytical jacobians done
444 Vector dB_dT(0);
445 ArrayOfPropagationMatrix dK_this_dx(nq), dK_past_dx(nq), dKp_dx(nq);
446 ArrayOfStokesVector da_dx(nq), dS_dx(nq), dSp_dx(nq);
447
448 // HSE variables
449 Index temperature_derivative_position = -1;
450 bool do_hse = false;
451
452 if (j_analytical_do) {
453 dB_dT.resize(nf);
455 dK_this_dx[iq] = PropagationMatrix(nf, ns);
456 dK_past_dx[iq] = PropagationMatrix(nf, ns);
457 dKp_dx[iq] = PropagationMatrix(nf, ns);
458 da_dx[iq] = StokesVector(nf, ns);
459 dS_dx[iq] = StokesVector(nf, ns);
460 dSp_dx[iq] = StokesVector(nf, ns);
461 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
462 temperature_derivative_position = iq;
463 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
464 })
465 }
466 const bool temperature_jacobian =
467 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
468
469 // Loop ppath points and determine radiative properties
470 for (Index ip = 0; ip < np; ip++) {
472 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
473
475 K_this,
476 S,
477 lte[ip],
478 dK_this_dx,
479 dS_dx,
480 propmat_clearsky_agenda,
481 jacobian_quantities,
482 ppvar_f(joker, ip),
483 ppvar_mag(joker, ip),
484 ppath.los(ip, joker),
485 ppvar_nlte[ip],
486 ppvar_vmr(joker, ip),
487 ppvar_t[ip],
488 ppvar_p[ip],
489 j_analytical_do);
490
491 if (j_analytical_do)
493 dS_dx,
494 jacobian_quantities,
495 ppvar_f(joker, ip),
496 ppath.los(ip, joker),
497 lte[ip],
498 atmosphere_dim,
499 j_analytical_do);
500
501 if (clear2cloudy[ip] + 1) {
503 Kp,
504 da_dx,
505 dKp_dx,
506 jacobian_quantities,
507 ppvar_pnd(joker, Range(ip, 1)),
508 ppvar_dpnd_dx,
509 ip,
510 scat_data,
511 ppath.los(ip, joker),
512 ppvar_t[Range(ip, 1)],
513 atmosphere_dim,
514 jacobian_do);
515 a += K_this;
516 K_this += Kp;
517
518 if (j_analytical_do)
519 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] += dK_this_dx[iq];
520 dK_this_dx[iq] += dKp_dx[iq];)
521
522 Vector aa_grid;
523 nlinspace(aa_grid, 0, 360, Naa);
524 //
526 dSp_dx,
527 jacobian_quantities,
528 ppvar_pnd(joker, ip),
529 ppvar_dpnd_dx,
530 ip,
531 scat_data,
532 cloudbox_field,
533 za_grid,
534 aa_grid,
535 ppath.los(Range(ip, 1), joker),
536 ppath.gp_p[ip],
537 ppvar_t[Range(ip, 1)],
538 atmosphere_dim,
539 jacobian_do,
540 t_interp_order);
541 S += Sp;
542
543 if (j_analytical_do)
544 FOR_ANALYTICAL_JACOBIANS_DO(dS_dx[iq] += dSp_dx[iq];)
545 } else { // no particles present at this level
546 a = K_this;
547 if (j_analytical_do)
548 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_this_dx[iq];)
549 }
550
551 if (ip not_eq 0) {
552 const Numeric dr_dT_past =
553 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
554 const Numeric dr_dT_this =
555 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
556 stepwise_transmission(lyr_tra[ip],
557 dlyr_tra_above[ip],
558 dlyr_tra_below[ip],
559 K_past,
560 K_this,
561 dK_past_dx,
562 dK_this_dx,
563 ppath.lstep[ip - 1],
564 dr_dT_past,
565 dr_dT_this,
566 temperature_derivative_position);
567 }
568
569 stepwise_source(src_rad[ip],
570 dsrc_rad[ip],
571 K_this,
572 a,
573 S,
574 dK_this_dx,
575 da_dx,
576 dS_dx,
577 B,
578 dB_dT,
579 jacobian_quantities,
580 jacobian_do);
581
582 swap(K_past, K_this);
583 swap(dK_past_dx, dK_this_dx);
584 }
585 }
586
587 const ArrayOfTransmissionMatrix tot_tra =
589
590 // iy_transmittance
591 Tensor3 iy_trans_new;
592 if (iy_agenda_call1)
593 iy_trans_new = tot_tra[np - 1];
594 else
595 iy_transmittance_mult(iy_trans_new, iy_transmittance, tot_tra[np - 1]);
596
597 // Copy transmission to iy_aux
598 for (Index i = 0; i < naux; i++)
599 if (iy_aux_vars[i] == "Optical depth")
600 for (Index iv = 0; iv < nf; iv++)
601 iy_aux[i](iv, joker) = -log(ppvar_trans_cumulat(np - 1, iv, 0, 0));
602
603 // Radiative background
605 iy,
606 diy_dx,
607 iy_trans_new,
608 iy_id,
609 jacobian_do,
610 jacobian_quantities,
611 ppath,
612 rte_pos2,
613 atmosphere_dim,
614 nlte_field,
615 cloudbox_on,
616 stokes_dim,
617 f_grid,
618 iy_unit,
619 surface_props_data,
620 iy_main_agenda,
621 iy_space_agenda,
622 iy_surface_agenda,
623 iy_cloudbox_agenda,
624 iy_agenda_call1,
625 verbosity);
626
627 lvl_rad[np - 1] = iy;
628
629 // Radiative transfer calculations
630 for (Index ip = np - 2; ip >= 0; ip--) {
631 lvl_rad[ip] = lvl_rad[ip + 1];
632 update_radiation_vector(lvl_rad[ip],
633 dlvl_rad[ip],
634 dlvl_rad[ip + 1],
635 src_rad[ip],
636 src_rad[ip + 1],
637 dsrc_rad[ip],
638 dsrc_rad[ip + 1],
639 lyr_tra[ip + 1],
640 tot_tra[ip],
641 dlyr_tra_above[ip + 1],
642 dlyr_tra_below[ip + 1],
647 Numeric(),
648 Vector(),
649 Vector(),
650 0,
651 0,
653 }
654
655 // Copy back to ARTS external style
656 iy = lvl_rad[0];
657 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
658 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
659 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
660 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
661 if (j_analytical_do)
662 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
663 dlvl_rad[ip][iq];);
664 }
665
666 // Finalize analytical Jacobians
667 if (j_analytical_do)
669 diy_dx,
670 diy_dpath,
671 ns,
672 nf,
673 np,
674 atmosphere_dim,
675 ppath,
676 ppvar_p,
677 ppvar_t,
678 ppvar_vmr,
679 iy_agenda_call1,
680 iy_transmittance,
681 water_p_eq_agenda,
682 jacobian_quantities,
683 jac_species_i);
684
685 // Unit conversions
686 if (iy_agenda_call1)
688 diy_dx,
689 ppvar_iy,
690 ns,
691 np,
692 f_grid,
693 ppath,
694 jacobian_quantities,
695 j_analytical_do,
696 iy_unit);
697}
698
699/* Workspace method: Doxygen documentation will be auto-generated */
701 Workspace& ws,
702 Matrix& iy,
703 ArrayOfMatrix& iy_aux,
704 ArrayOfTensor3& diy_dx,
705 Vector& ppvar_p,
706 Vector& ppvar_t,
707 EnergyLevelMap& ppvar_nlte,
708 Matrix& ppvar_vmr,
709 Matrix& ppvar_wind,
710 Matrix& ppvar_mag,
711 Matrix& ppvar_f,
712 Tensor3& ppvar_iy,
713 Tensor4& ppvar_trans_cumulat,
714 Tensor4& ppvar_trans_partial,
715 const Index& iy_id,
716 const Index& stokes_dim,
717 const Vector& f_grid,
718 const Index& atmosphere_dim,
719 const Vector& p_grid,
720 const Tensor3& t_field,
721 const EnergyLevelMap& nlte_field,
722 const Tensor4& vmr_field,
723 const ArrayOfArrayOfSpeciesTag& abs_species,
724 const Tensor3& wind_u_field,
725 const Tensor3& wind_v_field,
726 const Tensor3& wind_w_field,
727 const Tensor3& mag_u_field,
728 const Tensor3& mag_v_field,
729 const Tensor3& mag_w_field,
730 const Index& cloudbox_on,
731 const String& iy_unit,
732 const ArrayOfString& iy_aux_vars,
733 const Index& jacobian_do,
734 const ArrayOfRetrievalQuantity& jacobian_quantities,
735 const Ppath& ppath,
736 const Vector& rte_pos2,
737 const Agenda& propmat_clearsky_agenda,
738 const Agenda& water_p_eq_agenda,
739 const String& rt_integration_option,
740 const Agenda& iy_main_agenda,
741 const Agenda& iy_space_agenda,
742 const Agenda& iy_surface_agenda,
743 const Agenda& iy_cloudbox_agenda,
744 const Index& iy_agenda_call1,
745 const Tensor3& iy_transmittance,
746 const Numeric& rte_alonglos_v,
747 const Tensor3& surface_props_data,
748 const Verbosity& verbosity) {
749 // Init Jacobian quantities?
750 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<2>(jacobian_quantities) : 0;
751
752 // Some basic sizes
753 const Index nf = f_grid.nelem();
754 const Index ns = stokes_dim;
755 const Index np = ppath.np;
756 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
757
758 // Radiative background index
759 const Index rbi = ppath_what_background(ppath);
760
761 // Checks of input
762 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
763 "ppath.background is invalid. Check your "
764 "calculation of *ppath*?");
765 ARTS_USER_ERROR_IF (!iy_agenda_call1 && np == 1 && rbi == 2,
766 "A secondary propagation path starting at the "
767 "surface and is going directly into the surface "
768 "is found. This is not allowed.");
769
770 // Set diy_dpath if we are doing are doing jacobian calculations
771 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
772
773 // Set the species pointers if we are doing jacobian
774 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
775
776 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
777 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
778
779 // Init iy_aux and fill where possible
780 const Index naux = iy_aux_vars.nelem();
781 iy_aux.resize(naux);
782 //
783 Index auxOptDepth = -1;
784 //
785 for (Index i = 0; i < naux; i++) {
786 iy_aux[i].resize(nf, ns);
787 iy_aux[i] = 0;
788
789 if (iy_aux_vars[i] == "Radiative background")
790 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
791 else if (iy_aux_vars[i] == "Optical depth")
792 auxOptDepth = i;
793 else {
795 "The only allowed strings in *iy_aux_vars* are:\n"
796 " \"Radiative background\"\n"
797 " \"Optical depth\"\n"
798 "but you have selected: \"", iy_aux_vars[i], "\"")
799 }
800 }
801
802 // Get atmospheric and radiative variables along the propagation path
803 ppvar_trans_cumulat.resize(np, nf, ns, ns);
804 ppvar_trans_partial.resize(np, nf, ns, ns);
805 ppvar_iy.resize(nf, ns, np);
806
807 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
810
811 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
814
820
823 Vector r(np);
824 ArrayOfVector dr_below(np, Vector(nq, 0));
825 ArrayOfVector dr_above(np, Vector(nq, 0));
826
827 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
828 ppvar_p.resize(0);
829 ppvar_t.resize(0);
830 ppvar_vmr.resize(0, 0);
831 ppvar_wind.resize(0, 0);
832 ppvar_mag.resize(0, 0);
833 ppvar_f.resize(0, 0);
834 ppvar_trans_cumulat = 0;
835 ppvar_trans_partial = 0;
836 for (Index iv = 0; iv < nf; iv++) {
837 for (Index is = 0; is < ns; is++) {
838 ppvar_trans_cumulat(0,iv,is,is) = 1;
839 ppvar_trans_partial(0,iv,is,is) = 1;
840 }
841 }
842
843 } else {
844 // Basic atmospheric variables
845 get_ppath_atmvars(ppvar_p,
846 ppvar_t,
847 ppvar_nlte,
848 ppvar_vmr,
849 ppvar_wind,
850 ppvar_mag,
851 ppath,
852 atmosphere_dim,
853 p_grid,
854 t_field,
855 nlte_field,
856 vmr_field,
857 wind_u_field,
858 wind_v_field,
859 wind_w_field,
860 mag_u_field,
861 mag_v_field,
862 mag_w_field);
863
865 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
866
867 const bool temperature_jacobian =
868 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
869
870 // Size radiative variables always used
871 Vector B(nf);
872 StokesVector a(nf, ns), S(nf, ns);
873
874 // Init variables only used if analytical jacobians done
875 Vector dB_dT(temperature_jacobian ? nf : 0);
876 ArrayOfStokesVector da_dx(nq), dS_dx(nq);
877
878 // HSE variables
879 Index temperature_derivative_position = -1;
880 bool do_hse = false;
881
882 if (j_analytical_do) {
883 for (Index ip = 0; ip < np; ip++) {
884 dK_dx[ip].resize(nq);
886 }
888 da_dx[iq] = StokesVector(nf, ns); dS_dx[iq] = StokesVector(nf, ns);
889 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
890 temperature_derivative_position = iq;
891 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
892 })
893 }
894
895 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
896 Workspace l_ws(ws);
897 ArrayOfString fail_msg;
898 bool do_abort = false;
899
900 // Loop ppath points and determine radiative properties
901#pragma omp parallel for if (!arts_omp_in_parallel()) \
902 firstprivate(l_ws, l_propmat_clearsky_agenda, a, B, dB_dT, S, da_dx, dS_dx)
903 for (Index ip = 0; ip < np; ip++) {
904 if (do_abort) continue;
905 try {
907 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
908
909 Index lte;
911 K[ip],
912 S,
913 lte,
914 dK_dx[ip],
915 dS_dx,
916 l_propmat_clearsky_agenda,
917 jacobian_quantities,
918 ppvar_f(joker, ip),
919 ppvar_mag(joker, ip),
920 ppath.los(ip, joker),
921 ppvar_nlte[ip],
922 ppvar_vmr(joker, ip),
923 ppvar_t[ip],
924 ppvar_p[ip],
925 j_analytical_do);
926
927 if (j_analytical_do)
929 dS_dx,
930 jacobian_quantities,
931 ppvar_f(joker, ip),
932 ppath.los(ip, joker),
933 lte,
934 atmosphere_dim,
935 j_analytical_do);
936
937 // Here absorption equals extinction
938 a = K[ip];
939 if (j_analytical_do)
940 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_dx[ip][iq];);
941
942 stepwise_source(src_rad[ip],
943 dsrc_rad[ip],
944 K[ip],
945 a,
946 S,
947 dK_dx[ip],
948 da_dx,
949 dS_dx,
950 B,
951 dB_dT,
952 jacobian_quantities,
953 jacobian_do);
954 } catch (const std::runtime_error& e) {
955 ostringstream os;
956 os << "Runtime-error in source calculation at index " << ip
957 << ": \n";
958 os << e.what();
959#pragma omp critical(iyEmissionStandard_source)
960 {
961 do_abort = true;
962 fail_msg.push_back(os.str());
963 }
964 }
965 }
966
967#pragma omp parallel for if (!arts_omp_in_parallel())
968 for (Index ip = 1; ip < np; ip++) {
969 if (do_abort) continue;
970 try {
971 const Numeric dr_dT_past =
972 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
973 const Numeric dr_dT_this =
974 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
975 stepwise_transmission(lyr_tra[ip],
976 dlyr_tra_above[ip],
977 dlyr_tra_below[ip],
978 K[ip - 1],
979 K[ip],
980 dK_dx[ip - 1],
981 dK_dx[ip],
982 ppath.lstep[ip - 1],
983 dr_dT_past,
984 dr_dT_this,
985 temperature_derivative_position);
986
987 r[ip - 1] = ppath.lstep[ip - 1];
988 if (temperature_derivative_position >= 0){
989 dr_below[ip][temperature_derivative_position] = dr_dT_past;
990 dr_above[ip][temperature_derivative_position] = dr_dT_this;
991 }
992 } catch (const std::runtime_error& e) {
993 ostringstream os;
994 os << "Runtime-error in transmission calculation at index " << ip
995 << ": \n";
996 os << e.what();
997#pragma omp critical(iyEmissionStandard_transmission)
998 {
999 do_abort = true;
1000 fail_msg.push_back(os.str());
1001 }
1002 }
1003 }
1004
1005 ARTS_USER_ERROR_IF (do_abort,
1006 "Error messages from failed cases:\n", fail_msg)
1007 }
1008
1009 const ArrayOfTransmissionMatrix tot_tra =
1011
1012 // iy_transmittance
1013 Tensor3 iy_trans_new;
1014 if (iy_agenda_call1)
1015 iy_trans_new = tot_tra[np - 1];
1016 else
1017 iy_transmittance_mult(iy_trans_new, iy_transmittance, tot_tra[np - 1]);
1018
1019 // iy_aux: Optical depth
1020 if (auxOptDepth >= 0)
1021 for (Index iv = 0; iv < nf; iv++)
1022 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
1023
1024 // Radiative background
1026 iy,
1027 diy_dx,
1028 iy_trans_new,
1029 iy_id,
1030 jacobian_do,
1031 jacobian_quantities,
1032 ppath,
1033 rte_pos2,
1034 atmosphere_dim,
1035 nlte_field,
1036 cloudbox_on,
1037 stokes_dim,
1038 f_grid,
1039 iy_unit,
1040 surface_props_data,
1041 iy_main_agenda,
1042 iy_space_agenda,
1043 iy_surface_agenda,
1044 iy_cloudbox_agenda,
1045 iy_agenda_call1,
1046 verbosity);
1047
1048 lvl_rad[np - 1] = iy;
1049
1050 // Radiative transfer calculations
1051 if (rt_integration_option == "first order" || rt_integration_option == "default") {
1052 for (Index ip = np - 2; ip >= 0; ip--) {
1053 lvl_rad[ip] = lvl_rad[ip + 1];
1054 update_radiation_vector(lvl_rad[ip],
1055 dlvl_rad[ip],
1056 dlvl_rad[ip + 1],
1057 src_rad[ip],
1058 src_rad[ip + 1],
1059 dsrc_rad[ip],
1060 dsrc_rad[ip + 1],
1061 lyr_tra[ip + 1],
1062 tot_tra[ip],
1063 dlyr_tra_above[ip + 1],
1064 dlyr_tra_below[ip + 1],
1069 Numeric(),
1070 Vector(),
1071 Vector(),
1072 0,
1073 0,
1075 }
1076 } else if (rt_integration_option == "second order") {
1077 for (Index ip = np - 2; ip >= 0; ip--) {
1078 lvl_rad[ip] = lvl_rad[ip + 1];
1079 update_radiation_vector(lvl_rad[ip],
1080 dlvl_rad[ip],
1081 dlvl_rad[ip + 1],
1082 src_rad[ip],
1083 src_rad[ip + 1],
1084 dsrc_rad[ip],
1085 dsrc_rad[ip + 1],
1086 lyr_tra[ip + 1],
1087 tot_tra[ip],
1088 dlyr_tra_above[ip + 1],
1089 dlyr_tra_below[ip + 1],
1090 K[ip],
1091 K[ip + 1],
1092 dK_dx[ip + 1],
1093 dK_dx[ip + 1],
1094 r[ip],
1095 dr_above[ip + 1],
1096 dr_below[ip + 1],
1097 0,
1098 0,
1100 }
1101 } else {
1102 ARTS_USER_ERROR ( "Only allowed choices for *integration order* are "
1103 "1 and 2.");
1104 }
1105
1106 // Copy back to ARTS external style
1107 iy = lvl_rad[0];
1108 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
1109 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
1110 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
1111 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
1112 if (j_analytical_do)
1113 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
1114 dlvl_rad[ip][iq];);
1115 }
1116
1117 // Finalize analytical Jacobians
1118 if (j_analytical_do) {
1120 diy_dx,
1121 diy_dpath,
1122 ns,
1123 nf,
1124 np,
1125 atmosphere_dim,
1126 ppath,
1127 ppvar_p,
1128 ppvar_t,
1129 ppvar_vmr,
1130 iy_agenda_call1,
1131 iy_transmittance,
1132 water_p_eq_agenda,
1133 jacobian_quantities,
1134 jac_species_i);
1135 }
1136
1137 // Radiance unit conversions
1138 if (iy_agenda_call1) {
1140 diy_dx,
1141 ppvar_iy,
1142 ns,
1143 np,
1144 f_grid,
1145 ppath,
1146 jacobian_quantities,
1147 j_analytical_do,
1148 iy_unit);
1149 }
1150}
1151
1152/* Workspace method: Doxygen documentation will be auto-generated */
1154 Matrix& iy,
1155 ArrayOfMatrix& iy_aux,
1156 Ppath& ppath,
1157 ArrayOfTensor3& diy_dx,
1158 GriddedField4& atm_fields_compact,
1159 const Index& iy_id,
1160 const Vector& f_grid,
1161 const Index& atmosphere_dim,
1162 const Vector& p_grid,
1163 const Vector& lat_grid,
1164 const Vector& lon_grid,
1165 const Vector& lat_true,
1166 const Vector& lon_true,
1167 const Tensor3& t_field,
1168 const Tensor3& z_field,
1169 const Tensor4& vmr_field,
1170 const EnergyLevelMap& nlte_field,
1171 const Tensor3& wind_u_field,
1172 const Tensor3& wind_v_field,
1173 const Tensor3& wind_w_field,
1174 const Tensor3& mag_u_field,
1175 const Tensor3& mag_v_field,
1176 const Tensor3& mag_w_field,
1177 const Index& cloudbox_on,
1178 const ArrayOfIndex& cloudbox_limits,
1179 const Tensor4& pnd_field,
1180 const Matrix& particle_masses,
1181 const Agenda& ppath_agenda,
1182 const Numeric& ppath_lmax,
1183 const Numeric& ppath_lraytrace,
1184 const Index& iy_agenda_call1,
1185 const String& iy_unit,
1186 const Tensor3& iy_transmittance,
1187 const Vector& rte_pos,
1188 const Vector& rte_los,
1189 const Vector& rte_pos2,
1190 const Index& jacobian_do,
1191 const ArrayOfString& iy_aux_vars,
1192 const Agenda& iy_independent_beam_approx_agenda,
1193 const Index& return_atm1d,
1194 const Index& skip_vmr,
1195 const Index& skip_pnd,
1196 const Index& return_masses,
1197 const Verbosity&) {
1198 // Throw error if unsupported features are requested
1199 ARTS_USER_ERROR_IF (jacobian_do,
1200 "Jacobians not provided by the method, *jacobian_do* "
1201 "must be 0.");
1202 ARTS_USER_ERROR_IF (!nlte_field.value.empty(),
1203 "This method does not yet support non-empty *nlte_field*.");
1204 ARTS_USER_ERROR_IF (!wind_u_field.empty(),
1205 "This method does not yet support non-empty *wind_u_field*.");
1206 ARTS_USER_ERROR_IF (!wind_v_field.empty(),
1207 "This method does not yet support non-empty *wind_v_field*.");
1208 ARTS_USER_ERROR_IF (!wind_w_field.empty(),
1209 "This method does not yet support non-empty *wind_w_field*.");
1210 ARTS_USER_ERROR_IF (!mag_u_field.empty(),
1211 "This method does not yet support non-empty *mag_u_field*.");
1212 ARTS_USER_ERROR_IF (!mag_v_field.empty(),
1213 "This method does not yet support non-empty *mag_v_field*.");
1214 ARTS_USER_ERROR_IF (!mag_w_field.empty(),
1215 "This method does not yet support non-empty *mag_w_field*.");
1216 //
1217 if (return_masses) {
1218 ARTS_USER_ERROR_IF (pnd_field.nbooks() != particle_masses.nrows(),
1219 "Sizes of *pnd_field* and *particle_masses* "
1220 "are inconsistent.");
1221 }
1222 chk_latlon_true(atmosphere_dim,
1223 lat_grid,
1224 lat_true,
1225 lon_true);
1226
1227 // Note that input 1D atmospheres are handled exactly as 2D and 3D, to make
1228 // the function totally general. And 1D must be handled for iterative calls.
1229
1230 // Determine propagation path (with cloudbox deactivated) and check
1231 // that is OK for ICA
1232 //
1234 ppath,
1235 ppath_lmax,
1236 ppath_lraytrace,
1237 rte_pos,
1238 rte_los,
1239 rte_pos2,
1240 0,
1241 0,
1242 f_grid,
1243 ppath_agenda);
1244 //
1245 error_if_limb_ppath(ppath);
1246
1247 // If scattering and sensor inside atmosphere, we need a pseudo-ppath that
1248 // samples altitudes not covered by main ppath. We make this second path
1249 // strictly vertical.
1250 //
1251 Ppath ppath2;
1252 //
1253 if (cloudbox_on && ppath.end_lstep == 0) {
1254 Vector los_tmp = rte_los;
1255 if (abs(rte_los[0]) < 90) {
1256 los_tmp[0] = 180;
1257 } else {
1258 los_tmp[0] = 0;
1259 }
1260 //
1262 ppath2,
1263 ppath_lmax,
1264 ppath_lraytrace,
1265 rte_pos,
1266 los_tmp,
1267 rte_pos2,
1268 0,
1269 0,
1270 f_grid,
1271 ppath_agenda);
1272 } else {
1273 ppath2.np = 1;
1274 }
1275
1276 // Merge grid positions, and sort from bottom to top of atmosphere
1277 const Index np = ppath.np + ppath2.np - 1;
1278 ArrayOfGridPos gp_p(np), gp_lat(np), gp_lon(np);
1279 if (ppath.np > 1 &&
1280 ppath.pos(0, 0) >
1281 ppath.pos(1, 0)) { // Ppath is sorted in downward direction
1282 // Copy ppath in reversed order
1283 for (Index i = 0; i < ppath.np; i++) {
1284 const Index ip = ppath.np - i - 1;
1285 gp_p[i] = ppath.gp_p[ip];
1286 if (atmosphere_dim > 1) {
1287 gp_lat[i] = ppath.gp_lat[ip];
1288 if (atmosphere_dim == 3) {
1289 gp_lon[i] = ppath.gp_lon[ip];
1290 }
1291 }
1292 }
1293 // Append ppath2, but skipping element [0]
1294 for (Index i = ppath.np; i < np; i++) {
1295 const Index ip = i - ppath.np + 1;
1296 gp_p[i] = ppath2.gp_p[ip];
1297 if (atmosphere_dim > 1) {
1298 gp_lat[i] = ppath2.gp_lat[ip];
1299 if (atmosphere_dim == 3) {
1300 gp_lon[i] = ppath2.gp_lon[ip];
1301 }
1302 }
1303 }
1304 } else {
1305 // Copy ppath2 in reversed order, but skipping element [0]
1306 for (Index i = 0; i < ppath2.np - 1; i++) {
1307 const Index ip = ppath2.np - i - 1;
1308 gp_p[i] = ppath2.gp_p[ip];
1309 if (atmosphere_dim > 1) {
1310 gp_lat[i] = ppath2.gp_lat[ip];
1311 if (atmosphere_dim == 3) {
1312 gp_lon[i] = ppath2.gp_lon[ip];
1313 }
1314 }
1315 }
1316 // Append ppath
1317 for (Index i = ppath2.np - 1; i < np; i++) {
1318 const Index ip = i - ppath2.np + 1;
1319 gp_p[i] = ppath.gp_p[ip];
1320 if (atmosphere_dim > 1) {
1321 gp_lat[i] = ppath.gp_lat[ip];
1322 if (atmosphere_dim == 3) {
1323 gp_lon[i] = ppath.gp_lon[ip];
1324 }
1325 }
1326 }
1327 }
1328
1329 // 1D version of p_grid
1330 Matrix itw;
1331 Vector p1(np);
1332 ArrayOfGridPos gp0(0), gp1(1);
1333 interp_atmfield_gp2itw(itw, 1, gp_p, gp0, gp0);
1334 itw2p(p1, p_grid, gp_p, itw);
1335
1336 // 1D version of lat and lon variables
1337 Vector lat1(0), lon1(0);
1338 Vector lat_true1(1), lon_true1(1);
1339 if (atmosphere_dim == 3) {
1340 gp1[0] = gp_lat[0];
1341 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
1342 interp(lat_true1, itw, lat_grid, gp1);
1343 gp1[0] = gp_lon[0];
1344 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
1345 interp(lon_true1, itw, lon_grid, gp1);
1346 } else if (atmosphere_dim == 2) {
1347 gp1[0] = gp_lat[0];
1348 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
1349 interp(lat_true1, itw, lat_true, gp1);
1350 interp(lon_true1, itw, lon_true, gp1);
1351 } else {
1352 lat_true1[0] = lat_true[0];
1353 lon_true1[0] = lon_true[0];
1354 }
1355
1356 // 2D/3D interpolation weights
1357 interp_atmfield_gp2itw(itw, atmosphere_dim, gp_p, gp_lat, gp_lon);
1358
1359 // 1D temperature field
1360 Tensor3 t1(np, 1, 1);
1362 t1(joker, 0, 0), atmosphere_dim, t_field, gp_p, gp_lat, gp_lon, itw);
1363
1364 // 1D altitude field
1365 Tensor3 z1(np, 1, 1);
1367 z1(joker, 0, 0), atmosphere_dim, z_field, gp_p, gp_lat, gp_lon, itw);
1368
1369 // 1D VMR field
1370 Tensor4 vmr1(vmr_field.nbooks(), np, 1, 1);
1371 for (Index is = 0; is < vmr_field.nbooks(); is++) {
1372 interp_atmfield_by_itw(vmr1(is, joker, 0, 0),
1373 atmosphere_dim,
1374 vmr_field(is, joker, joker, joker),
1375 gp_p,
1376 gp_lat,
1377 gp_lon,
1378 itw);
1379 }
1380
1381 // 1D surface altitude
1382 Matrix zsurf1(1, 1);
1383 zsurf1(0, 0) = z1(0, 0, 0);
1384
1385 // 1D version of rte_pos/los
1386 Vector pos1(1);
1387 pos1[0] = rte_pos[0];
1388 Vector los1(1);
1389 los1[0] = abs(rte_los[0]);
1390 Vector pos2(0);
1391 if (rte_pos2.nelem()) {
1392 pos2 = rte_pos2[Range(0, rte_pos2.nelem())];
1393 }
1394
1395 // Cloudbox variables
1396 //
1397 Index cbox_on1 = cloudbox_on;
1398 ArrayOfIndex cbox_lims1(0);
1399 Tensor4 pnd1(0, 0, 0, 0);
1400 //
1401 if (cloudbox_on) {
1402 // Determine what p1-levels that are inside cloudbox
1403 Index ifirst = np;
1404 Index ilast = -1;
1405 for (Index i = 0; i < np; i++) {
1406 if (is_gp_inside_cloudbox(gp_p[i],
1407 gp_lat[i],
1408 gp_lon[i],
1409 cloudbox_limits,
1410 true,
1411 atmosphere_dim)) {
1412 if (i < ifirst) {
1413 ifirst = i;
1414 }
1415 ilast = i;
1416 }
1417 }
1418
1419 // If no hit, deactive cloudbox
1420 if (ifirst == np) {
1421 cbox_on1 = 0;
1422 }
1423
1424 // Otherwise set 1D cloud variables
1425 else {
1426 // We can enter the cloudbox from the side, and we need to add 1
1427 // level on each side to be safe (if not limits already at edges)
1428 //
1429 const Index extra_bot = ifirst == 0 ? 0 : 1;
1430 const Index extra_top = ilast == np - 1 ? 0 : 1;
1431 //
1432 cbox_lims1.resize(2);
1433 cbox_lims1[0] = ifirst - extra_bot;
1434 cbox_lims1[1] = ilast + extra_top;
1435
1436 // pnd_field
1437 //
1438 pnd1.resize(pnd_field.nbooks(), cbox_lims1[1] - cbox_lims1[0] + 1, 1, 1);
1439 pnd1 = 0; // As lowermost and uppermost level not always filled
1440 //
1441 itw.resize(1, Index(pow(2.0, Numeric(atmosphere_dim))));
1442 //
1443 for (Index i = extra_bot; i < pnd1.npages() - extra_top; i++) {
1444 const Index i0 = cbox_lims1[0] + i;
1445 ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
1447 gpc_p[0],
1448 gpc_lat[0],
1449 gpc_lon[0],
1450 gp_p[i0],
1451 gp_lat[i0],
1452 gp_lon[i0],
1453 atmosphere_dim,
1454 cloudbox_limits);
1455 for (Index p = 0; p < pnd_field.nbooks(); p++) {
1456 interp_atmfield_by_itw(pnd1(p, i, 0, 0),
1457 atmosphere_dim,
1458 pnd_field(p, joker, joker, joker),
1459 gpc_p,
1460 gpc_lat,
1461 gpc_lon,
1462 itw);
1463 }
1464 }
1465 }
1466 }
1467
1468 // Call sub agenda
1469 //
1470 {
1471 const Index adim1 = 1;
1472 const Numeric lmax1 = -1;
1473 Ppath ppath1d;
1474 //
1476 iy,
1477 iy_aux,
1478 ppath1d,
1479 diy_dx,
1480 iy_agenda_call1,
1481 iy_unit,
1482 iy_transmittance,
1483 iy_aux_vars,
1484 iy_id,
1485 adim1,
1486 p1,
1487 lat1,
1488 lon1,
1489 lat_true1,
1490 lon_true1,
1491 t1,
1492 z1,
1493 vmr1,
1494 zsurf1,
1495 lmax1,
1496 ppath_lraytrace,
1497 cbox_on1,
1498 cbox_lims1,
1499 pnd1,
1500 jacobian_do,
1501 f_grid,
1502 pos1,
1503 los1,
1504 pos2,
1505 iy_independent_beam_approx_agenda);
1506 }
1507
1508 // Fill *atm_fields_compact*?
1509 if (return_atm1d) {
1510 // Sizes and allocate memory
1511 const Index nvmr = skip_vmr ? 0 : vmr1.nbooks();
1512 const Index npnd = skip_pnd ? 0 : pnd1.nbooks();
1513 const Index nmass = return_masses ? particle_masses.ncols() : 0;
1514 const Index ntot = 2 + nvmr + npnd + nmass;
1515 ArrayOfString field_names(ntot);
1516 atm_fields_compact.resize(ntot, np, 1, 1);
1517
1518 // Altitudes
1519 field_names[0] = "Geometric altitudes";
1520 atm_fields_compact.data(0, joker, 0, 0) = z1(joker, 0, 0);
1521
1522 // Temperature
1523 field_names[1] = "Temperature";
1524 atm_fields_compact.data(1, joker, 0, 0) = t1(joker, 0, 0);
1525
1526 // VMRs
1527 if (nvmr) {
1528 for (Index i = 0; i < nvmr; i++) {
1529 const Index iout = 2 + i;
1530 ostringstream sstr;
1531 sstr << "VMR species " << i;
1532 field_names[iout] = sstr.str();
1533 atm_fields_compact.data(iout, joker, 0, 0) = vmr1(i, joker, 0, 0);
1534 }
1535 }
1536
1537 // PNDs
1538 if (npnd) {
1539 for (Index i = 0; i < npnd; i++) {
1540 const Index iout = 2 + nvmr + i;
1541 ostringstream sstr;
1542 sstr << "Scattering element " << i;
1543 field_names[iout] = sstr.str();
1544 atm_fields_compact.data(iout, joker, 0, 0) = 0;
1545 atm_fields_compact.data(
1546 iout, Range(cbox_lims1[0], pnd1.npages()), 0, 0) =
1547 pnd1(i, joker, 0, 0);
1548 }
1549 }
1550
1551 // Masses
1552 if (nmass) {
1553 for (Index i = 0; i < nmass; i++) {
1554 const Index iout = 2 + nvmr + npnd + i;
1555 ostringstream sstr;
1556 sstr << "Mass category " << i;
1557 field_names[iout] = sstr.str();
1558 atm_fields_compact.data(iout, joker, 0, 0) = 0;
1559 for (Index ip = cbox_lims1[0]; ip < pnd1.npages(); ip++) {
1560 for (Index is = 0; is < pnd1.nbooks(); is++) {
1561 atm_fields_compact.data(iout, ip, 0, 0) +=
1562 particle_masses(is, i) * pnd1(is, ip, 0, 0);
1563 }
1564 }
1565 }
1566 }
1567
1568 // Finally, set grids and names
1569 //
1570 atm_fields_compact.set_name(
1571 "Data created by *iyIndependentBeamApproximation*");
1572 //
1573 atm_fields_compact.set_grid_name(GFIELD4_FIELD_NAMES,
1574 "Atmospheric quantity");
1575 atm_fields_compact.set_grid(GFIELD4_FIELD_NAMES, field_names);
1576 atm_fields_compact.set_grid_name(GFIELD4_P_GRID, "Pressure");
1577 atm_fields_compact.set_grid(GFIELD4_P_GRID, p1);
1578 atm_fields_compact.set_grid_name(GFIELD4_LAT_GRID, "Latitude");
1579 atm_fields_compact.set_grid(GFIELD4_LAT_GRID, lat_true1);
1580 atm_fields_compact.set_grid_name(GFIELD4_LON_GRID, "Longitude");
1581 atm_fields_compact.set_grid(GFIELD4_LON_GRID, lon_true1);
1582 }
1583}
1584
1585/* Workspace method: Doxygen documentation will be auto-generated */
1587 Matrix& iy,
1588 ArrayOfMatrix& iy_aux,
1589 Ppath& ppath,
1590 ArrayOfTensor3& diy_dx,
1591 const ArrayOfString& iy_aux_vars,
1592 const Index& iy_agenda_call1,
1593 const Tensor3& iy_transmittance,
1594 const Vector& rte_pos,
1595 const Vector& rte_los,
1596 const Vector& rte_pos2,
1597 const Index& stokes_dim,
1598 const Vector& f_grid,
1599 const Agenda& iy_loop_freqs_agenda,
1600 const Verbosity&) {
1601 // Throw error if unsupported features are requested
1602 ARTS_USER_ERROR_IF (!iy_agenda_call1,
1603 "Recursive usage not possible (iy_agenda_call1 must be 1).");
1604 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
1605 "*iy_transmittance* must be empty.");
1606
1607 const Index nf = f_grid.nelem();
1608
1609 for (Index i = 0; i < nf; i++) {
1610 // Variables for 1 frequency
1611 Matrix iy1;
1612 ArrayOfMatrix iy_aux1;
1613 ArrayOfTensor3 diy_dx1;
1614
1616 iy1,
1617 iy_aux1,
1618 ppath,
1619 diy_dx1,
1620 iy_agenda_call1,
1621 iy_transmittance,
1622 iy_aux_vars,
1623 0,
1624 Vector(1, f_grid[i]),
1625 rte_pos,
1626 rte_los,
1627 rte_pos2,
1628 iy_loop_freqs_agenda);
1629
1630 // After first frequency, give output its size
1631 if (i == 0) {
1632 iy.resize(nf, stokes_dim);
1633 //
1634 iy_aux.resize(iy_aux1.nelem());
1635 for (Index q = 0; q < iy_aux1.nelem(); q++) {
1636 iy_aux[q].resize(nf, stokes_dim);
1637 }
1638 //
1639 diy_dx.resize(diy_dx1.nelem());
1640 for (Index q = 0; q < diy_dx1.nelem(); q++) {
1641 diy_dx[q].resize(diy_dx1[q].npages(), nf, stokes_dim);
1642 }
1643 }
1644
1645 // Copy to output variables
1646 iy(i, joker) = iy1(0, joker);
1647 for (Index q = 0; q < iy_aux1.nelem(); q++) {
1648 iy_aux[q](i, joker) = iy_aux1[q](0, joker);
1649 }
1650 for (Index q = 0; q < diy_dx1.nelem(); q++) {
1651 diy_dx[q](joker, i, joker) = diy_dx1[q](joker, 0, joker);
1652 }
1653 }
1654}
1655
1656/* Workspace method: Doxygen documentation will be auto-generated */
1658 Matrix& iy,
1659 ArrayOfMatrix& iy_aux,
1660 ArrayOfTensor3& diy_dx,
1661 const Index& iy_agenda_call1,
1662 const Tensor3& iy_transmittance,
1663 const Vector& rte_pos,
1664 const Vector& rte_los,
1665 const ArrayOfString& iy_aux_vars,
1666 const Index& jacobian_do,
1667 const Index& atmosphere_dim,
1668 const Vector& p_grid,
1669 const Vector& lat_grid,
1670 const Vector& lon_grid,
1671 const Tensor3& z_field,
1672 const Tensor3& t_field,
1673 const Tensor4& vmr_field,
1674 const Vector& refellipsoid,
1675 const Matrix& z_surface,
1676 const Index& cloudbox_on,
1677 const ArrayOfIndex& cloudbox_limits,
1678 const Index& stokes_dim,
1679 const Vector& f_grid,
1680 const ArrayOfArrayOfSingleScatteringData& scat_data,
1681 const Agenda& iy_space_agenda,
1682 const Agenda& surface_rtprop_agenda,
1683 const Agenda& propmat_clearsky_agenda,
1684 const Agenda& ppath_step_agenda,
1685 const Numeric& ppath_lmax,
1686 const Numeric& ppath_lraytrace,
1687 const Tensor4& pnd_field,
1688 const String& iy_unit,
1689 const Numeric& mc_std_err,
1690 const Index& mc_max_time,
1691 const Index& mc_max_iter,
1692 const Index& mc_min_iter,
1693 const Numeric& mc_taustep_limit,
1694 const Index& t_interp_order,
1695 const Verbosity& verbosity) {
1696 // Throw error if unsupported features are requested
1697 ARTS_USER_ERROR_IF (atmosphere_dim != 3,
1698 "Only 3D atmospheres are allowed (atmosphere_dim must be 3)");
1699 ARTS_USER_ERROR_IF (!cloudbox_on,
1700 "The cloudbox must be activated (cloudbox_on must be 1)");
1701 ARTS_USER_ERROR_IF (jacobian_do,
1702 "This method does not provide any jacobians (jacobian_do must be 0)");
1703 ARTS_USER_ERROR_IF (!iy_agenda_call1,
1704 "Recursive usage not possible (iy_agenda_call1 must be 1)");
1705 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
1706 "*iy_transmittance* must be empty");
1707
1708 // Size output variables
1709 //
1710 const Index nf = f_grid.nelem();
1711 //
1712 iy.resize(nf, stokes_dim);
1713 diy_dx.resize(0);
1714 //
1715 //=== iy_aux part ===========================================================
1716 Index auxError = -1;
1717 {
1718 const Index naux = iy_aux_vars.nelem();
1719 iy_aux.resize(naux);
1720 //
1721 for (Index i = 0; i < naux; i++) {
1722 if (iy_aux_vars[i] == "Error (uncorrelated)") {
1723 auxError = i;
1724 iy_aux[i].resize(nf, stokes_dim);
1725 } else {
1727 "In *iy_aux_vars* you have included: \"", iy_aux_vars[i],
1728 "\"\nThis choice is not recognised.")
1729 }
1730 }
1731 }
1732 //===========================================================================
1733
1734 // Some MC variables are only local here
1735 //
1736 MCAntenna mc_antenna;
1737 mc_antenna.set_pencil_beam();
1738
1739 // Pos and los must be matrices
1740 Matrix pos(1, 3), los(1, 2);
1741 //
1742 pos(0, joker) = rte_pos;
1743 los(0, joker) = rte_los;
1744
1745 Workspace l_ws(ws);
1746 Agenda l_ppath_step_agenda(ppath_step_agenda);
1747 Agenda l_iy_space_agenda(iy_space_agenda);
1748 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
1749 Agenda l_surface_rtprop_agenda(surface_rtprop_agenda);
1750
1751 String fail_msg;
1752 bool failed = false;
1753
1754 if (nf)
1755#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) \
1756 firstprivate(l_ws, \
1757 l_ppath_step_agenda, \
1758 l_iy_space_agenda, \
1759 l_propmat_clearsky_agenda, \
1760 l_surface_rtprop_agenda)
1761 for (Index f_index = 0; f_index < nf; f_index++) {
1762 if (failed) continue;
1763
1764 try {
1765 // Seed reset for each loop. If not done, the errors
1766 // appear to be highly correlated.
1767 Index mc_seed;
1768 MCSetSeedFromTime(mc_seed, verbosity);
1769
1770 Vector y, mc_error;
1771 Index mc_iteration_count;
1772 Tensor3 mc_points;
1773 ArrayOfIndex mc_scat_order, mc_source_domain;
1774
1775 MCGeneral(l_ws,
1776 y,
1777 mc_iteration_count,
1778 mc_error,
1779 mc_points,
1780 mc_scat_order,
1781 mc_source_domain,
1782 mc_antenna,
1783 f_grid,
1784 f_index,
1785 pos,
1786 los,
1787 stokes_dim,
1788 atmosphere_dim,
1789 l_ppath_step_agenda,
1790 ppath_lmax,
1791 ppath_lraytrace,
1792 l_iy_space_agenda,
1793 l_surface_rtprop_agenda,
1794 l_propmat_clearsky_agenda,
1795 p_grid,
1796 lat_grid,
1797 lon_grid,
1798 z_field,
1799 refellipsoid,
1800 z_surface,
1801 t_field,
1802 vmr_field,
1803 cloudbox_on,
1804 cloudbox_limits,
1805 pnd_field,
1806 scat_data,
1807 1,
1808 1,
1809 1,
1810 1,
1811 iy_unit,
1812 mc_seed,
1813 mc_std_err,
1814 mc_max_time,
1815 mc_max_iter,
1816 mc_min_iter,
1817 mc_taustep_limit,
1818 1,
1819 t_interp_order,
1820 verbosity);
1821
1822 ARTS_ASSERT(y.nelem() == stokes_dim);
1823
1824 iy(f_index, joker) = y;
1825
1826 if (auxError >= 0) {
1827 iy_aux[auxError](f_index, joker) = mc_error;
1828 }
1829 } catch (const std::exception& e) {
1830 ostringstream os;
1831 os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
1832 << ")" << endl
1833 << e.what();
1834#pragma omp critical(iyMC_fail)
1835 {
1836 failed = true;
1837 fail_msg = os.str();
1838 }
1839 continue;
1840 }
1841 }
1842
1843 ARTS_USER_ERROR_IF (failed, fail_msg);
1844}
1845
1846/* Workspace method: Doxygen documentation will be auto-generated */
1848 const ArrayOfMatrix& iy_aux,
1849 const ArrayOfString& iy_aux_vars,
1850 const Index& jacobian_do,
1851 const String& aux_var,
1852 const Verbosity&) {
1853 ARTS_USER_ERROR_IF (iy_aux.nelem() != iy_aux_vars.nelem(),
1854 "*iy_aux* and *iy_aux_vars* must have the same "
1855 "number of elements.");
1856
1857 ARTS_USER_ERROR_IF (jacobian_do,
1858 "This method can not provide any jacobians and "
1859 "*jacobian_do* must be 0.");
1860
1861 bool ready = false;
1862
1863 for (Index i = 0; i < iy_aux.nelem() && !ready; i++) {
1864 if (iy_aux_vars[i] == aux_var) {
1865 iy = iy_aux[i];
1866 ready = true;
1867 }
1868 }
1869
1870 ARTS_USER_ERROR_IF (!ready,
1871 "The selected auxiliary variable to insert in *iy* "
1872 "is either not defined at all or is not set.");
1873}
1874
1875/* Workspace method: Doxygen documentation will be auto-generated */
1877 Matrix& ppvar_optical_depth,
1878 const Tensor4& ppvar_trans_cumulat,
1879 const Verbosity&) {
1880 ppvar_optical_depth = ppvar_trans_cumulat(joker, joker, 0, 0);
1881 transform(ppvar_optical_depth, log, ppvar_optical_depth);
1882 ppvar_optical_depth *= -1;
1883}
1884
1885/* Workspace method: Doxygen documentation will be auto-generated */
1887 Vector& y,
1888 Vector& y_f,
1889 ArrayOfIndex& y_pol,
1890 Matrix& y_pos,
1891 Matrix& y_los,
1892 ArrayOfVector& y_aux,
1893 Matrix& y_geo,
1894 Matrix& jacobian,
1895 const Index& atmgeom_checked,
1896 const Index& atmfields_checked,
1897 const Index& atmosphere_dim,
1898 const EnergyLevelMap& nlte_field,
1899 const Index& cloudbox_on,
1900 const Index& cloudbox_checked,
1901 const Index& scat_data_checked,
1902 const Index& sensor_checked,
1903 const Index& stokes_dim,
1904 const Vector& f_grid,
1905 const Matrix& sensor_pos,
1906 const Matrix& sensor_los,
1907 const Matrix& transmitter_pos,
1908 const Matrix& mblock_dlos_grid,
1909 const Sparse& sensor_response,
1910 const Vector& sensor_response_f,
1911 const ArrayOfIndex& sensor_response_pol,
1912 const Matrix& sensor_response_dlos,
1913 const String& iy_unit,
1914 const Agenda& iy_main_agenda,
1915 const Agenda& geo_pos_agenda,
1916 const Agenda& jacobian_agenda,
1917 const Index& jacobian_do,
1918 const ArrayOfRetrievalQuantity& jacobian_quantities,
1919 const ArrayOfString& iy_aux_vars,
1920 const Verbosity& verbosity) {
1922
1923 // Basics
1924 //
1925 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1926 //
1927 ARTS_USER_ERROR_IF (f_grid.empty(),
1928 "The frequency grid is empty.");
1929 chk_if_increasing("f_grid", f_grid);
1930 ARTS_USER_ERROR_IF (f_grid[0] <= 0,
1931 "All frequencies in *f_grid* must be > 0.");
1932 //
1933 ARTS_USER_ERROR_IF (atmfields_checked != 1,
1934 "The atmospheric fields must be flagged to have\n"
1935 "passed a consistency check (atmfields_checked=1).");
1936 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
1937 "The atmospheric geometry must be flagged to have\n"
1938 "passed a consistency check (atmgeom_checked=1).");
1939 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
1940 "The cloudbox must be flagged to have\n"
1941 "passed a consistency check (cloudbox_checked=1).");
1942 if (cloudbox_on)
1943 ARTS_USER_ERROR_IF (scat_data_checked != 1,
1944 "The scattering data must be flagged to have\n"
1945 "passed a consistency check (scat_data_checked=1).");
1946 ARTS_USER_ERROR_IF (sensor_checked != 1,
1947 "The sensor variables must be flagged to have\n"
1948 "passed a consistency check (sensor_checked=1).");
1949
1950 // Some sizes
1951 const Index nf = f_grid.nelem();
1952 const Index nlos = mblock_dlos_grid.nrows();
1953 const Index n1y = sensor_response.nrows();
1954 const Index nmblock = sensor_pos.nrows();
1955 const Index niyb = nf * nlos * stokes_dim;
1956
1957 //---------------------------------------------------------------------------
1958 // Allocations and resizing
1959 //---------------------------------------------------------------------------
1960
1961 // Resize *y* and *y_XXX*
1962 //
1963 y.resize(nmblock * n1y);
1964 y_f.resize(nmblock * n1y);
1965 y_pol.resize(nmblock * n1y);
1966 y_pos.resize(nmblock * n1y, sensor_pos.ncols());
1967 y_los.resize(nmblock * n1y, sensor_los.ncols());
1968 y_geo.resize(nmblock * n1y, 5);
1969 y_geo = NAN; // Will be replaced if relavant data are provided (*geo_pos*)
1970
1971 // For y_aux we don't know the number of quantities, and we need to
1972 // store all output
1973 ArrayOfArrayOfVector iyb_aux_array(nmblock);
1974
1975 // Jacobian variables
1976 //
1977 Index j_analytical_do = 0;
1978 ArrayOfArrayOfIndex jacobian_indices;
1979 //
1980 if (jacobian_do) {
1981 bool any_affine;
1982 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
1983 jacobian.resize(nmblock * n1y,
1984 jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
1985 jacobian = 0;
1986 //
1987 FOR_ANALYTICAL_JACOBIANS_DO2(j_analytical_do = 1;)
1988 } else {
1989 jacobian.resize(0, 0);
1990 }
1991
1992 //---------------------------------------------------------------------------
1993 // The calculations
1994 //---------------------------------------------------------------------------
1995
1996 String fail_msg;
1997 bool failed = false;
1998
1999 if (nmblock >= arts_omp_get_max_threads() ||
2000 (nf <= nmblock && nmblock >= nlos)) {
2001 out3 << " Parallelizing mblock loop (" << nmblock << " iterations)\n";
2002
2003 // We have to make a local copy of the Workspace and the agendas because
2004 // only non-reference types can be declared firstprivate in OpenMP
2005 Workspace l_ws(ws);
2006 Agenda l_jacobian_agenda(jacobian_agenda);
2007 Agenda l_iy_main_agenda(iy_main_agenda);
2008 Agenda l_geo_pos_agenda(geo_pos_agenda);
2009
2010#pragma omp parallel for firstprivate( \
2011 l_ws, l_jacobian_agenda, l_iy_main_agenda, l_geo_pos_agenda)
2012 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2013 // Skip remaining iterations if an error occurred
2014 if (failed) continue;
2015
2017 fail_msg,
2018 iyb_aux_array,
2019 l_ws,
2020 y,
2021 y_f,
2022 y_pol,
2023 y_pos,
2024 y_los,
2025 y_geo,
2026 jacobian,
2027 atmosphere_dim,
2028 nlte_field,
2029 cloudbox_on,
2030 stokes_dim,
2031 f_grid,
2032 sensor_pos,
2033 sensor_los,
2034 transmitter_pos,
2035 mblock_dlos_grid,
2036 sensor_response,
2037 sensor_response_f,
2038 sensor_response_pol,
2039 sensor_response_dlos,
2040 iy_unit,
2041 l_iy_main_agenda,
2042 l_geo_pos_agenda,
2043 l_jacobian_agenda,
2044 jacobian_do,
2045 jacobian_quantities,
2046 jacobian_indices,
2047 iy_aux_vars,
2048 verbosity,
2049 mblock_index,
2050 n1y,
2051 j_analytical_do);
2052 } // End mblock loop
2053 } else {
2054 out3 << " Not parallelizing mblock loop (" << nmblock << " iterations)\n";
2055
2056 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2057 // Skip remaining iterations if an error occurred
2058 if (failed) continue;
2059
2061 fail_msg,
2062 iyb_aux_array,
2063 ws,
2064 y,
2065 y_f,
2066 y_pol,
2067 y_pos,
2068 y_los,
2069 y_geo,
2070 jacobian,
2071 atmosphere_dim,
2072 nlte_field,
2073 cloudbox_on,
2074 stokes_dim,
2075 f_grid,
2076 sensor_pos,
2077 sensor_los,
2078 transmitter_pos,
2079 mblock_dlos_grid,
2080 sensor_response,
2081 sensor_response_f,
2082 sensor_response_pol,
2083 sensor_response_dlos,
2084 iy_unit,
2085 iy_main_agenda,
2086 geo_pos_agenda,
2087 jacobian_agenda,
2088 jacobian_do,
2089 jacobian_quantities,
2090 jacobian_indices,
2091 iy_aux_vars,
2092 verbosity,
2093 mblock_index,
2094 n1y,
2095 j_analytical_do);
2096 } // End mblock loop
2097 }
2098
2099 // Rethrow exception if a runtime error occurred in the mblock loop
2100 ARTS_USER_ERROR_IF (failed, fail_msg);
2101
2102 // Compile y_aux
2103 //
2104 const Index nq = iyb_aux_array[0].nelem();
2105 y_aux.resize(nq);
2106 //
2107 for (Index q = 0; q < nq; q++) {
2108 y_aux[q].resize(nmblock * n1y);
2109 //
2110 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2111 const Range rowind =
2112 get_rowindex_for_mblock(sensor_response, mblock_index);
2113 const Index row0 = rowind.get_start();
2114
2115 // The sensor response must be applied in a special way for
2116 // uncorrelated errors. Schematically: sqrt( H.^2 * y.^2 )
2117 if (iy_aux_vars[q] == "Error (uncorrelated)") {
2118 for (Index i = 0; i < n1y; i++) {
2119 const Index row = row0 + i;
2120 y_aux[q][row] = 0;
2121 for (Index j = 0; j < niyb; j++) {
2122 y_aux[q][row] +=
2123 pow(sensor_response(i, j) * iyb_aux_array[mblock_index][q][j],
2124 (Numeric)2.0);
2125 }
2126 y_aux[q][row] = sqrt(y_aux[q][row]);
2127 }
2128 } else {
2129 mult(y_aux[q][rowind], sensor_response, iyb_aux_array[mblock_index][q]);
2130 }
2131 }
2132 }
2133}
2134
2135/* Workspace method: Doxygen documentation will be auto-generated */
2137 Vector& y,
2138 Vector& y_f,
2139 ArrayOfIndex& y_pol,
2140 Matrix& y_pos,
2141 Matrix& y_los,
2142 ArrayOfVector& y_aux,
2143 Matrix& y_geo,
2144 Matrix& jacobian,
2145 ArrayOfRetrievalQuantity& jacobian_quantities,
2146 const Index& atmfields_checked,
2147 const Index& atmgeom_checked,
2148 const Index& atmosphere_dim,
2149 const EnergyLevelMap& nlte_field,
2150 const Index& cloudbox_on,
2151 const Index& cloudbox_checked,
2152 const Index& scat_data_checked,
2153 const Index& sensor_checked,
2154 const Index& stokes_dim,
2155 const Vector& f_grid,
2156 const Matrix& sensor_pos,
2157 const Matrix& sensor_los,
2158 const Matrix& transmitter_pos,
2159 const Matrix& mblock_dlos_grid,
2160 const Sparse& sensor_response,
2161 const Vector& sensor_response_f,
2162 const ArrayOfIndex& sensor_response_pol,
2163 const Matrix& sensor_response_dlos,
2164 const String& iy_unit,
2165 const Agenda& iy_main_agenda,
2166 const Agenda& geo_pos_agenda,
2167 const Agenda& jacobian_agenda,
2168 const Index& jacobian_do,
2169 const ArrayOfString& iy_aux_vars,
2170 const ArrayOfRetrievalQuantity& jacobian_quantities_copy,
2171 const Index& append_instrument_wfs,
2172 const Verbosity& verbosity) {
2173 // The jacobian indices of old and new part (without transformations)
2174 ArrayOfArrayOfIndex jacobian_indices, jacobian_indices_copy;
2175 {
2176 bool any_affine;
2178 jacobian_indices_copy, any_affine, jacobian_quantities_copy, true);
2179 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
2180 }
2181
2182 // Check consistency of data representing first measurement
2183 const Index n1 = y.nelem();
2184 Index nrq1 = 0;
2186 "Input *y* is empty. Use *yCalc*");
2187 ARTS_USER_ERROR_IF (y_f.nelem() != n1,
2188 "Lengths of input *y* and *y_f* are inconsistent.");
2189 ARTS_USER_ERROR_IF (y_pol.nelem() != n1,
2190 "Lengths of input *y* and *y_pol* are inconsistent.");
2191 ARTS_USER_ERROR_IF (y_pos.nrows() != n1,
2192 "Sizes of input *y* and *y_pos* are inconsistent.");
2193 ARTS_USER_ERROR_IF (y_los.nrows() != n1,
2194 "Sizes of input *y* and *y_los* are inconsistent.");
2195 ARTS_USER_ERROR_IF (y_geo.nrows() != n1,
2196 "Sizes of input *y* and *y_geo* are inconsistent.");
2197 if (jacobian_do) {
2198 nrq1 = jacobian_quantities_copy.nelem();
2199 ARTS_USER_ERROR_IF (jacobian.nrows() != n1,
2200 "Sizes of *y* and *jacobian* are inconsistent.");
2201 ARTS_USER_ERROR_IF (jacobian.ncols() != jacobian_indices_copy[nrq1 - 1][1] + 1,
2202 "Size of input *jacobian* and size implied "
2203 "*jacobian_quantities_copy* are inconsistent.");
2204 }
2205
2206 // Calculate new measurement
2207 //
2208 Vector y2, y_f2;
2209 Matrix y_pos2, y_los2, y_geo2, jacobian2;
2210 ArrayOfIndex y_pol2;
2211 ArrayOfVector y_aux2;
2212 //
2213 yCalc(ws,
2214 y2,
2215 y_f2,
2216 y_pol2,
2217 y_pos2,
2218 y_los2,
2219 y_aux2,
2220 y_geo2,
2221 jacobian2,
2222 atmfields_checked,
2223 atmgeom_checked,
2224 atmosphere_dim,
2225 nlte_field,
2226 cloudbox_on,
2227 cloudbox_checked,
2228 scat_data_checked,
2229 sensor_checked,
2230 stokes_dim,
2231 f_grid,
2232 sensor_pos,
2233 sensor_los,
2234 transmitter_pos,
2235 mblock_dlos_grid,
2236 sensor_response,
2237 sensor_response_f,
2238 sensor_response_pol,
2239 sensor_response_dlos,
2240 iy_unit,
2241 iy_main_agenda,
2242 geo_pos_agenda,
2243 jacobian_agenda,
2244 jacobian_do,
2245 jacobian_quantities,
2246 iy_aux_vars,
2247 verbosity);
2248
2249 // Consistency checks
2250 ARTS_USER_ERROR_IF (y_pos.ncols() != y_pos2.ncols(),
2251 "Different number of columns in *y_pos* between the measurements.");
2252 ARTS_USER_ERROR_IF (y_los.ncols() != y_los2.ncols(),
2253 "Different number of columns in *y_los* between the measurements.");
2254
2255 // y and y_XXX
2256 //
2257 const Index n2 = y2.nelem();
2258 //
2259 {
2260 // Make copy of old measurement
2261 const Vector y1 = y, y_f1 = y_f;
2262 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
2263 const ArrayOfIndex y_pol1 = y_pol;
2264 const ArrayOfVector y_aux1 = y_aux;
2265 //
2266 y.resize(n1 + n2);
2267 y[Range(0, n1)] = y1;
2268 y[Range(n1, n2)] = y2;
2269 //
2270 y_f.resize(n1 + n2);
2271 y_f[Range(0, n1)] = y_f1;
2272 y_f[Range(n1, n2)] = y_f2;
2273 //
2274 y_pos.resize(n1 + n2, y_pos1.ncols());
2275 y_pos(Range(0, n1), joker) = y_pos1;
2276 y_pos(Range(n1, n2), joker) = y_pos2;
2277 //
2278 y_los.resize(n1 + n2, y_los1.ncols());
2279 y_los(Range(0, n1), joker) = y_los1;
2280 y_los(Range(n1, n2), joker) = y_los2;
2281 //
2282 y_geo.resize(n1 + n2, y_geo1.ncols());
2283 y_geo(Range(0, n1), joker) = y_geo1;
2284 y_geo(Range(n1, n2), joker) = y_geo2;
2285 //
2286 y_pol.resize(n1 + n2);
2287 for (Index i = 0; i < n1; i++) {
2288 y_pol[i] = y_pol1[i];
2289 }
2290 for (Index i = 0; i < n2; i++) {
2291 y_pol[n1 + i] = y_pol2[i];
2292 }
2293
2294 // y_aux
2295 const Index na1 = y_aux1.nelem();
2296 const Index na2 = y_aux2.nelem();
2297 const Index na = max(na1, na2);
2298 //
2299 y_aux.resize(na);
2300 //
2301 for (Index a = 0; a < na; a++) {
2302 y_aux[a].resize(n1 + n2);
2303 if (a < na1) {
2304 y_aux[a][Range(0, n1)] = y_aux1[a];
2305 } else {
2306 y_aux[a][Range(0, n1)] = 0;
2307 }
2308 if (a < na2) {
2309 y_aux[a][Range(n1, n2)] = y_aux2[a];
2310 } else {
2311 y_aux[a][Range(n1, n2)] = 0;
2312 }
2313 }
2314 }
2315
2316 // Jacobian and friends
2317 if (jacobian_do) {
2318 // Put in old jacobian_quantities and jacobian_indices as first part in
2319 // new version of these variables
2320 ArrayOfRetrievalQuantity jacobian_quantities2 = jacobian_quantities;
2321 ArrayOfArrayOfIndex jacobian_indices2 = jacobian_indices;
2322 //
2323 jacobian_quantities = jacobian_quantities_copy;
2324 jacobian_indices = jacobian_indices_copy;
2325
2326 // Loop new jacobian_quantities to determine how new jacobian data shall
2327 // be inserted
2328 //
2329 const Index nrq2 = jacobian_quantities2.nelem();
2330 ArrayOfIndex map_table(nrq2);
2331 //
2332 for (Index q2 = 0; q2 < nrq2; q2++) {
2333 Index pos = -1;
2334
2335 // Compare to old quantities, to determine if append shall be
2336 // considered. Some special checks performed here, grids checked later
2337 if (jacobian_quantities2[q2].Target().isSpeciesVMR() ||
2338 jacobian_quantities2[q2] == Jacobian::Atm::Temperature ||
2339 jacobian_quantities2[q2] == Jacobian::Special::ScatteringString ||
2340 jacobian_quantities2[q2].Target().isWind() ||
2341 jacobian_quantities2[q2] == Jacobian::Special::SurfaceString ||
2342 append_instrument_wfs) {
2343 for (Index q1 = 0; q1 < nrq1; q1++ && pos < 0) { // FIXME: What is with this "&& pos < 0"
2344 if (jacobian_quantities2[q2].Target().sameTargetType(jacobian_quantities_copy[q1].Target())) {
2345 if (jacobian_quantities2[q2].Target().isSpeciesVMR()) {
2346 if (jacobian_quantities2[q2].Subtag() ==
2347 jacobian_quantities_copy[q1].Subtag()) {
2348 if (jacobian_quantities2[q2].Mode() ==
2349 jacobian_quantities_copy[q1].Mode()) {
2350 pos = q1;
2351 } else {
2353 "Jacobians for ", jacobian_quantities2[q2],
2354 " shall be appended.\nThis requires "
2355 "that the same retrieval unit is used "
2356 "but it seems that this requirement is "
2357 "not met.")
2358 }
2359 }
2360 } else if (jacobian_quantities2[q2] == Jacobian::Atm::Temperature) {
2361 if (jacobian_quantities2[q2].Subtag() ==
2362 jacobian_quantities_copy[q1].Subtag()) {
2363 pos = q1;
2364 } else {
2366 "Jacobians for ", jacobian_quantities2[q2],
2367 " shall be appended.\nThis requires "
2368 "that HSE is either ON or OFF for both "
2369 "parts but it seems that this requirement "
2370 "is not met.")
2371 }
2372 } else if (jacobian_quantities[q2] == Jacobian::Special::ScatteringString) {
2373 if ((jacobian_quantities2[q2].Subtag() ==
2374 jacobian_quantities_copy[q1].Subtag()) &&
2375 (jacobian_quantities2[q2].SubSubtag() ==
2376 jacobian_quantities_copy[q1].SubSubtag())) {
2377 pos = q1;
2378 }
2379 } else if (jacobian_quantities2[q2].Subtag() == jacobian_quantities_copy[q1].Subtag()) {
2380 pos = q1;
2381 }
2382 }
2383 }
2384 }
2385
2386 // New quantity
2387 if (pos < 0) {
2388 map_table[q2] = jacobian_quantities.nelem();
2389 jacobian_quantities.push_back(jacobian_quantities2[q2]);
2390 ArrayOfIndex indices(2);
2391 indices[0] = jacobian_indices[jacobian_indices.nelem() - 1][1] + 1;
2392 indices[1] =
2393 indices[0] + jacobian_indices2[q2][1] - jacobian_indices2[q2][0];
2394 jacobian_indices.push_back(indices);
2395 }
2396 // Existing quantity
2397 else {
2398 map_table[q2] = pos;
2399 // Check if grids are equal
2400 ArrayOfVector grids1 = jacobian_quantities_copy[pos].Grids();
2401 ArrayOfVector grids2 = jacobian_quantities2[q2].Grids();
2402 bool any_wrong = false;
2403 if (grids1.nelem() != grids2.nelem()) {
2404 any_wrong = true;
2405 } else {
2406 for (Index g = 0; g < grids1.nelem(); g++) {
2407 if (grids1[g].nelem() != grids2[g].nelem()) {
2408 any_wrong = true;
2409 } else {
2410 for (Index e = 0; e < grids1[g].nelem(); e++) {
2411 const Numeric v1 = grids1[g][e];
2412 const Numeric v2 = grids2[g][e];
2413 if ((v1 == 0 && abs(v2) > 1e-9) || abs(v1 - v2) / v1 > 1e-6) {
2414 any_wrong = true;
2415 }
2416 }
2417 }
2418 }
2419 }
2420 if (any_wrong) {
2422 "Jacobians for ", jacobian_quantities2[q2],
2423 " shall be appended.\nThis requires that the "
2424 "same grids are used for both measurements,\nbut "
2425 "it seems that this requirement is not met.")
2426 }
2427 }
2428 }
2429
2430 // Create and fill *jacobian*
2431 //
2432 const Index nrq = jacobian_quantities.nelem();
2433 const Matrix jacobian1 = jacobian;
2434 //
2435 jacobian.resize(n1 + n2, jacobian_indices[nrq - 1][1] + 1);
2436 jacobian = 0;
2437 //
2438 // Put in old part in top-left corner
2439 jacobian(Range(0, n1), Range(0, jacobian_indices_copy[nrq1 - 1][1] + 1)) =
2440 jacobian1;
2441 // New parts
2442 for (Index q2 = 0; q2 < nrq2; q2++) {
2443 jacobian(Range(n1, n2),
2444 Range(jacobian_indices[map_table[q2]][0],
2445 jacobian_indices[map_table[q2]][1] -
2446 jacobian_indices[map_table[q2]][0] + 1)) =
2447 jacobian2(
2448 joker,
2449 Range(jacobian_indices2[q2][0],
2450 jacobian_indices2[q2][1] - jacobian_indices2[q2][0] + 1));
2451 }
2452 }
2453}
2454
2455/* Workspace method: Doxygen documentation will be auto-generated */
2457 Matrix& jacobian,
2458 const Vector& y_f,
2459 const ArrayOfIndex& y_pol,
2460 const String& iy_unit,
2461 const Verbosity&) {
2462 ARTS_USER_ERROR_IF (iy_unit == "1",
2463 "No need to use this method with *iy_unit* = \"1\".");
2464
2465 ARTS_USER_ERROR_IF (max(y) > 1e-3,
2466 "The spectrum vector *y* is required to have original radiance\n"
2467 "unit, but this seems not to be the case. This as a value above\n"
2468 "1e-3 is found in *y*.")
2469
2470 // Is jacobian set?
2471 //
2472 const Index ny = y.nelem();
2473 //
2474 const bool do_j = jacobian.nrows() == ny;
2475
2476 // Some jacobian quantities can not be handled
2477 ARTS_USER_ERROR_IF (do_j && max(jacobian) > 1e-3,
2478 "The method can not be used with jacobian quantities that are not\n"
2479 "obtained through radiative transfer calculations. One example on\n"
2480 "quantity that can not be handled is *jacobianAddPolyfit*.\n"
2481 "The maximum value of *jacobian* indicates that one or several\n"
2482 "such jacobian quantities are included.")
2483
2484 // Planck-Tb
2485 //--------------------------------------------------------------------------
2486 if (iy_unit == "PlanckBT") {
2487 // Hard to use telescoping here as the data are sorted differently in y
2488 // and jacobian, than what is expected apply_iy_unit. Copy to temporary
2489 // variables instead.
2490
2491 // Handle the elements in "frequency chunks"
2492
2493 Index i0 = 0; // Index of first element for present chunk
2494 //
2495 while (i0 < ny) {
2496 // Find number of values for this chunk
2497 Index n = 1;
2498 //
2499 while (i0 + n < ny && y_f[i0] == y_f[i0 + n]) {
2500 n++;
2501 }
2502
2503 Matrix yv(1, n);
2504 ArrayOfIndex i_pol(n);
2505 bool any_quv = false;
2506 //
2507 for (Index i = 0; i < n; i++) {
2508 const Index ix = i0 + i;
2509 yv(0, i) = y[ix];
2510 i_pol[i] = y_pol[ix];
2511 if (i_pol[i] > 1 && i_pol[i] < 5) {
2512 any_quv = true;
2513 }
2514 }
2515
2516 // Index of elements to convert
2517 Range ii(i0, n);
2518
2519 if (do_j) {
2520 ARTS_USER_ERROR_IF (any_quv && i_pol[0] != 1,
2521 "The conversion to PlanckBT, of the Jacobian and "
2522 "errors for Q, U and V, requires that I (first Stokes "
2523 "element) is at hand and that the data are sorted in "
2524 "such way that I comes first for each frequency.")
2525
2526 // Jacobian
2527 Tensor3 J(jacobian.ncols(), 1, n);
2528 J(joker, 0, joker) = transpose(jacobian(ii, joker));
2529 apply_iy_unit2(J, yv, iy_unit, y_f[i0], 1, i_pol);
2530 jacobian(ii, joker) = transpose(J(joker, 0, joker));
2531 }
2532
2533 // y (must be done last)
2534 apply_iy_unit(yv, iy_unit, y_f[i0], 1, i_pol);
2535 y[ii] = yv(0, joker);
2536
2537 i0 += n;
2538 }
2539 }
2540
2541 // Other conversions
2542 //--------------------------------------------------------------------------
2543 else {
2544 // Here we take each element of y separately.
2545
2546 Matrix yv(1, 1);
2547 ArrayOfIndex i_pol(1);
2548
2549 for (Index i = 0; i < ny; i++) {
2550 yv(0, 0) = y[i];
2551 i_pol[0] = y_pol[i];
2552
2553 // Jacobian
2554 if (do_j) {
2556 MatrixView(jacobian(i, joker)), yv, iy_unit, y_f[i], 1, i_pol);
2557 }
2558
2559 // y (must be done last)
2560 apply_iy_unit(yv, iy_unit, y_f[i], 1, i_pol);
2561 y[i] = yv(0, 0);
2562 }
2563 }
2564}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
Array< Tensor3 > ArrayOfTensor3
An array of Tensor3.
Definition: array.h:74
The global header file for ARTS.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:46
Header file for helper functions for OpenMP.
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:24210
void iy_main_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 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:24276
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:24890
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:24093
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:86
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:111
The Agenda class.
Definition: agenda_class.h:51
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
Index nrows() const noexcept
Definition: matpackI.h:1061
Index ncols() const noexcept
Definition: matpackI.h:1062
bool empty() const noexcept
Definition: matpackIII.h:152
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
bool empty() const noexcept
Definition: matpackIV.h:148
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
Index ncols() const noexcept
Definition: matpackVII.h:159
Index npages() const noexcept
Definition: matpackVII.h:157
Index nrows() const noexcept
Definition: matpackVII.h:158
Index nlibraries() const noexcept
Definition: matpackVII.h:153
Index nvitrines() const noexcept
Definition: matpackVII.h:154
Index nshelves() const noexcept
Definition: matpackVII.h:155
Index nbooks() const noexcept
Definition: matpackVII.h:156
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:530
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
void set_name(const String &s)
Set name of this gridded field.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
The MatrixView class.
Definition: matpackI.h:1169
The Matrix class.
Definition: matpackI.h:1270
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1063
The range class.
Definition: matpackI.h:166
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:347
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:344
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:427
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The Tensor7 class.
Definition: matpackVII.h:2397
The Vector class.
Definition: matpackI.h:908
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
Workspace class.
Definition: workspace_ng.h:40
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:460
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:564
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:605
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1145
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:596
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:618
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Definition: jacobian.cc:45
Routines for setting up the jacobian.
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:556
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:548
#define q1
#define abs(x)
#define q
#define ns
#define min(a, b)
#define max(a, b)
Header file for logic.cc.
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, ArrayOfIndex &mc_source_domain, ArrayOfIndex &mc_scat_order, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &cloudbox_checked, const String &iy_unit, const Index &mc_seed, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Numeric &taustep_limit, const Index &l_mc_scat_order, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
Definition: m_montecarlo.cc:87
const Index GFIELD4_P_GRID
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:1586
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:62
void iyCalc(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, 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:94
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:700
void yCalcAppend(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, ArrayOfRetrievalQuantity &jacobian_quantities, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_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:2136
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:1153
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:1657
const Index GFIELD4_FIELD_NAMES
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:2456
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:1847
const Numeric PI
const Numeric SPEED_OF_LIGHT
void yCalc(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, const Index &atmgeom_checked, const Index &atmfields_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_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:1886
const Index GFIELD4_LAT_GRID
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:1876
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:160
const Index GFIELD4_LON_GRID
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
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:1484
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:207
Index nelem(const Lines &l)
Number of lines.
Temperature
Definition: jacobian.h:58
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ScatteringString
Definition: jacobian.h:102
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:53
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
This file contains declerations of functions of physical character.
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition: ppath.cc:540
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1476
Propagation path structure and functions.
#define a
Array< PropagationMatrix > ArrayOfPropagationMatrix
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
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:1202
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_magnetic_field, const ConstVectorView &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const ConstVectorView &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1131
void get_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:728
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:1319
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition: rte.cc:843
void 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:179
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &y_geo, Matrix &jacobian, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_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:2176
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:1803
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:106
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
Definition: rte.cc:1066
void rtmethods_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:2134
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:1996
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:1110
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition: rte.cc:55
void get_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:1116
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Definition: rte.cc:962
Declaration of functions in rte.cc.
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Header file for special_interp.cc.
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:50
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Matrix los
Line-of-sight at each ppath point.
Definition: ppath_struct.h:35
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath_struct.h:55
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath_struct.h:53
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath_struct.h:45
Vector lstep
The length between ppath points.
Definition: ppath_struct.h:39
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath_struct.h:51
Radiation Vector for Stokes dimension 1-4.
The Sparse class.
Definition: matpackII.h:67
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, 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.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Stuff related to the transmission matrix.
Array< RadiationVector > ArrayOfRadiationVector
Array< TransmissionMatrix > ArrayOfTransmissionMatrix