ARTS 2.5.0 (git: 9ee3ac6c)
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 Workspace& ws,
162 Matrix& iy,
163 ArrayOfMatrix& iy_aux,
164 ArrayOfTensor3& diy_dx,
165 Vector& ppvar_p,
166 Vector& ppvar_t,
167 EnergyLevelMap& ppvar_nlte,
168 Matrix& ppvar_vmr,
169 Matrix& ppvar_wind,
170 Matrix& ppvar_mag,
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 String& iy_unit,
192 const ArrayOfString& iy_aux_vars,
193 const Index& jacobian_do,
194 const ArrayOfRetrievalQuantity& jacobian_quantities,
195 const Ppath& ppath,
196 const Vector& rte_pos2,
197 const Agenda& propmat_clearsky_agenda,
198 const Agenda& water_p_eq_agenda,
199 const String& rt_integration_option,
200 const Agenda& iy_main_agenda,
201 const Agenda& iy_space_agenda,
202 const Agenda& iy_surface_agenda,
203 const Agenda& iy_cloudbox_agenda,
204 const Index& iy_agenda_call1,
205 const Tensor3& iy_transmittance,
206 const Numeric& rte_alonglos_v,
207 const Tensor3& surface_props_data,
208 const Verbosity& verbosity) {
209 // Init Jacobian quantities?
210 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<2>(jacobian_quantities) : 0;
211
212 // Some basic sizes
213 const Index nf = f_grid.nelem();
214 const Index ns = stokes_dim;
215 const Index np = ppath.np;
216 const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
217
218 // Radiative background index
219 const Index rbi = ppath_what_background(ppath);
220
221 // Checks of input
222 ARTS_USER_ERROR_IF (rbi < 1 || rbi > 9,
223 "ppath.background is invalid. Check your "
224 "calculation of *ppath*?");
225 ARTS_USER_ERROR_IF (!iy_agenda_call1 && np == 1 && rbi == 2,
226 "A secondary propagation path starting at the "
227 "surface and is going directly into the surface "
228 "is found. This is not allowed.");
229
230 // Set diy_dpath if we are doing are doing jacobian calculations
231 ArrayOfTensor3 diy_dpath = j_analytical_do ? get_standard_diy_dpath(jacobian_quantities, np, nf, ns, false) : ArrayOfTensor3(0);
232
233 // Set the species pointers if we are doing jacobian
234 const ArrayOfIndex jac_species_i = j_analytical_do ? get_pointers_for_analytical_species(jacobian_quantities, abs_species) : ArrayOfIndex(0);
235
236 // Start diy_dx out if we are doing the first run and are doing jacobian calculations
237 if (j_analytical_do and iy_agenda_call1) diy_dx = get_standard_starting_diy_dx(jacobian_quantities, np, nf, ns, false);
238
239 // Init iy_aux and fill where possible
240 const Index naux = iy_aux_vars.nelem();
241 iy_aux.resize(naux);
242 //
243 Index auxOptDepth = -1;
244 //
245 for (Index i = 0; i < naux; i++) {
246 iy_aux[i].resize(nf, ns);
247 iy_aux[i] = 0;
248
249 if (iy_aux_vars[i] == "Radiative background")
250 iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
251 else if (iy_aux_vars[i] == "Optical depth")
252 auxOptDepth = i;
253 else {
255 "The only allowed strings in *iy_aux_vars* are:\n"
256 " \"Radiative background\"\n"
257 " \"Optical depth\"\n"
258 "but you have selected: \"", iy_aux_vars[i], "\"")
259 }
260 }
261
262 // Get atmospheric and radiative variables along the propagation path
263 ppvar_trans_cumulat.resize(np, nf, ns, ns);
264 ppvar_trans_partial.resize(np, nf, ns, ns);
265 ppvar_iy.resize(nf, ns, np);
266
267 ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
270
271 ArrayOfRadiationVector src_rad(np, RadiationVector(nf, ns));
274
280
283 Vector r(np);
284 ArrayOfVector dr_below(np, Vector(nq, 0));
285 ArrayOfVector dr_above(np, Vector(nq, 0));
286
287 if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
288 ppvar_p.resize(0);
289 ppvar_t.resize(0);
290 ppvar_vmr.resize(0, 0);
291 ppvar_wind.resize(0, 0);
292 ppvar_mag.resize(0, 0);
293 ppvar_f.resize(0, 0);
294 ppvar_trans_cumulat = 0;
295 ppvar_trans_partial = 0;
296 for (Index iv = 0; iv < nf; iv++) {
297 for (Index is = 0; is < ns; is++) {
298 ppvar_trans_cumulat(0,iv,is,is) = 1;
299 ppvar_trans_partial(0,iv,is,is) = 1;
300 }
301 }
302
303 } else {
304 // Basic atmospheric variables
305 get_ppath_atmvars(ppvar_p,
306 ppvar_t,
307 ppvar_nlte,
308 ppvar_vmr,
309 ppvar_wind,
310 ppvar_mag,
311 ppath,
312 atmosphere_dim,
313 p_grid,
314 t_field,
315 nlte_field,
316 vmr_field,
317 wind_u_field,
318 wind_v_field,
319 wind_w_field,
320 mag_u_field,
321 mag_v_field,
322 mag_w_field);
323
325 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
326
327 const bool temperature_jacobian =
328 j_analytical_do and do_temperature_jacobian(jacobian_quantities);
329
330 // Size radiative variables always used
331 Vector B(nf);
332 StokesVector a(nf, ns), S(nf, ns);
333
334 // Init variables only used if analytical jacobians done
335 Vector dB_dT(temperature_jacobian ? nf : 0);
336 ArrayOfStokesVector da_dx(nq), dS_dx(nq);
337
338 // HSE variables
339 Index temperature_derivative_position = -1;
340 bool do_hse = false;
341
342 if (j_analytical_do) {
343 for (Index ip = 0; ip < np; ip++) {
344 dK_dx[ip].resize(nq);
346 }
348 da_dx[iq] = StokesVector(nf, ns); dS_dx[iq] = StokesVector(nf, ns);
349 if (jacobian_quantities[iq] == Jacobian::Atm::Temperature) {
350 temperature_derivative_position = iq;
351 do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
352 })
353 }
354
355 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
356 Workspace l_ws(ws);
357 ArrayOfString fail_msg;
358 bool do_abort = false;
359
360 // Loop ppath points and determine radiative properties
361#pragma omp parallel for if (!arts_omp_in_parallel()) \
362 firstprivate(l_ws, l_propmat_clearsky_agenda, a, B, dB_dT, S, da_dx, dS_dx)
363 for (Index ip = 0; ip < np; ip++) {
364 if (do_abort) continue;
365 try {
367 B, dB_dT, ppvar_f(joker, ip), ppvar_t[ip], temperature_jacobian);
368
369 Index lte;
371 K[ip],
372 S,
373 lte,
374 dK_dx[ip],
375 dS_dx,
376 l_propmat_clearsky_agenda,
377 jacobian_quantities,
378 ppvar_f(joker, ip),
379 ppvar_mag(joker, ip),
380 ppath.los(ip, joker),
381 ppvar_nlte[ip],
382 ppvar_vmr(joker, ip),
383 ppvar_t[ip],
384 ppvar_p[ip],
385 j_analytical_do);
386
387 if (j_analytical_do)
389 dS_dx,
390 jacobian_quantities,
391 ppvar_f(joker, ip),
392 ppath.los(ip, joker),
393 lte,
394 atmosphere_dim,
395 j_analytical_do);
396
397 // Here absorption equals extinction
398 a = K[ip];
399 if (j_analytical_do)
400 FOR_ANALYTICAL_JACOBIANS_DO(da_dx[iq] = dK_dx[ip][iq];);
401
402 stepwise_source(src_rad[ip],
403 dsrc_rad[ip],
404 K[ip],
405 a,
406 S,
407 dK_dx[ip],
408 da_dx,
409 dS_dx,
410 B,
411 dB_dT,
412 jacobian_quantities,
413 jacobian_do);
414 } catch (const std::runtime_error& e) {
415 ostringstream os;
416 os << "Runtime-error in source calculation at index " << ip
417 << ": \n";
418 os << e.what();
419#pragma omp critical(iyEmissionStandard_source)
420 {
421 do_abort = true;
422 fail_msg.push_back(os.str());
423 }
424 }
425 }
426
427#pragma omp parallel for if (!arts_omp_in_parallel())
428 for (Index ip = 1; ip < np; ip++) {
429 if (do_abort) continue;
430 try {
431 const Numeric dr_dT_past =
432 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
433 const Numeric dr_dT_this =
434 do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
435 stepwise_transmission(lyr_tra[ip],
436 dlyr_tra_above[ip],
437 dlyr_tra_below[ip],
438 K[ip - 1],
439 K[ip],
440 dK_dx[ip - 1],
441 dK_dx[ip],
442 ppath.lstep[ip - 1],
443 dr_dT_past,
444 dr_dT_this,
445 temperature_derivative_position);
446
447 r[ip - 1] = ppath.lstep[ip - 1];
448 if (temperature_derivative_position >= 0){
449 dr_below[ip][temperature_derivative_position] = dr_dT_past;
450 dr_above[ip][temperature_derivative_position] = dr_dT_this;
451 }
452 } catch (const std::runtime_error& e) {
453 ostringstream os;
454 os << "Runtime-error in transmission calculation at index " << ip
455 << ": \n";
456 os << e.what();
457#pragma omp critical(iyEmissionStandard_transmission)
458 {
459 do_abort = true;
460 fail_msg.push_back(os.str());
461 }
462 }
463 }
464
465 ARTS_USER_ERROR_IF (do_abort,
466 "Error messages from failed cases:\n", fail_msg)
467 }
468
469 const ArrayOfTransmissionMatrix tot_tra =
471
472 // iy_transmittance
473 Tensor3 iy_trans_new;
474 if (iy_agenda_call1)
475 iy_trans_new = tot_tra[np - 1];
476 else
477 iy_transmittance_mult(iy_trans_new, iy_transmittance, tot_tra[np - 1]);
478
479 // iy_aux: Optical depth
480 if (auxOptDepth >= 0)
481 for (Index iv = 0; iv < nf; iv++)
482 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
483
484 // Radiative background
486 iy,
487 diy_dx,
488 iy_trans_new,
489 iy_id,
490 jacobian_do,
491 jacobian_quantities,
492 ppath,
493 rte_pos2,
494 atmosphere_dim,
495 nlte_field,
496 cloudbox_on,
497 stokes_dim,
498 f_grid,
499 iy_unit,
500 surface_props_data,
501 iy_main_agenda,
502 iy_space_agenda,
503 iy_surface_agenda,
504 iy_cloudbox_agenda,
505 iy_agenda_call1,
506 verbosity);
507
508 lvl_rad[np - 1] = iy;
509
510 // Radiative transfer calculations
511 if (rt_integration_option == "first order" || rt_integration_option == "default") {
512 for (Index ip = np - 2; ip >= 0; ip--) {
513 lvl_rad[ip] = lvl_rad[ip + 1];
514 update_radiation_vector(lvl_rad[ip],
515 dlvl_rad[ip],
516 dlvl_rad[ip + 1],
517 src_rad[ip],
518 src_rad[ip + 1],
519 dsrc_rad[ip],
520 dsrc_rad[ip + 1],
521 lyr_tra[ip + 1],
522 tot_tra[ip],
523 dlyr_tra_above[ip + 1],
524 dlyr_tra_below[ip + 1],
529 Numeric(),
530 Vector(),
531 Vector(),
532 0,
533 0,
535 }
536 } else if (rt_integration_option == "second order") {
537 for (Index ip = np - 2; ip >= 0; ip--) {
538 lvl_rad[ip] = lvl_rad[ip + 1];
539 update_radiation_vector(lvl_rad[ip],
540 dlvl_rad[ip],
541 dlvl_rad[ip + 1],
542 src_rad[ip],
543 src_rad[ip + 1],
544 dsrc_rad[ip],
545 dsrc_rad[ip + 1],
546 lyr_tra[ip + 1],
547 tot_tra[ip],
548 dlyr_tra_above[ip + 1],
549 dlyr_tra_below[ip + 1],
550 K[ip],
551 K[ip + 1],
552 dK_dx[ip + 1],
553 dK_dx[ip + 1],
554 r[ip],
555 dr_above[ip + 1],
556 dr_below[ip + 1],
557 0,
558 0,
560 }
561 } else {
562 ARTS_USER_ERROR ( "Only allowed choices for *integration order* are "
563 "1 and 2.");
564 }
565
566 // Copy back to ARTS external style
567 iy = lvl_rad[0];
568 for (Index ip = 0; ip < lvl_rad.nelem(); ip++) {
569 ppvar_trans_cumulat(ip, joker, joker, joker) = tot_tra[ip];
570 ppvar_trans_partial(ip, joker, joker, joker) = lyr_tra[ip];
571 ppvar_iy(joker, joker, ip) = lvl_rad[ip];
572 if (j_analytical_do)
573 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq](ip, joker, joker) =
574 dlvl_rad[ip][iq];);
575 }
576
577 // Finalize analytical Jacobians
578 if (j_analytical_do) {
580 diy_dx,
581 diy_dpath,
582 ns,
583 nf,
584 np,
585 atmosphere_dim,
586 ppath,
587 ppvar_p,
588 ppvar_t,
589 ppvar_vmr,
590 iy_agenda_call1,
591 iy_transmittance,
592 water_p_eq_agenda,
593 jacobian_quantities,
594 jac_species_i);
595 }
596
597 // Radiance unit conversions
598 if (iy_agenda_call1) {
600 diy_dx,
601 ppvar_iy,
602 ns,
603 np,
604 f_grid,
605 ppath,
606 jacobian_quantities,
607 j_analytical_do,
608 iy_unit);
609 }
610}
611
612/* Workspace method: Doxygen documentation will be auto-generated */
614 Matrix& iy,
615 ArrayOfMatrix& iy_aux,
616 Ppath& ppath,
617 ArrayOfTensor3& diy_dx,
618 GriddedField4& atm_fields_compact,
619 const Index& iy_id,
620 const Vector& f_grid,
621 const Index& atmosphere_dim,
622 const Vector& p_grid,
623 const Vector& lat_grid,
624 const Vector& lon_grid,
625 const Vector& lat_true,
626 const Vector& lon_true,
627 const Tensor3& t_field,
628 const Tensor3& z_field,
629 const Tensor4& vmr_field,
630 const EnergyLevelMap& nlte_field,
631 const Tensor3& wind_u_field,
632 const Tensor3& wind_v_field,
633 const Tensor3& wind_w_field,
634 const Tensor3& mag_u_field,
635 const Tensor3& mag_v_field,
636 const Tensor3& mag_w_field,
637 const Index& cloudbox_on,
638 const ArrayOfIndex& cloudbox_limits,
639 const Tensor4& pnd_field,
640 const Matrix& particle_masses,
641 const Agenda& ppath_agenda,
642 const Numeric& ppath_lmax,
643 const Numeric& ppath_lraytrace,
644 const Index& iy_agenda_call1,
645 const String& iy_unit,
646 const Tensor3& iy_transmittance,
647 const Vector& rte_pos,
648 const Vector& rte_los,
649 const Vector& rte_pos2,
650 const Index& jacobian_do,
651 const ArrayOfString& iy_aux_vars,
652 const Agenda& iy_independent_beam_approx_agenda,
653 const Index& return_atm1d,
654 const Index& skip_vmr,
655 const Index& skip_pnd,
656 const Index& return_masses,
657 const Verbosity&) {
658 // Throw error if unsupported features are requested
659 ARTS_USER_ERROR_IF (jacobian_do,
660 "Jacobians not provided by the method, *jacobian_do* "
661 "must be 0.");
662 ARTS_USER_ERROR_IF (!nlte_field.Data().empty(),
663 "This method does not yet support non-empty *nlte_field*.");
664 ARTS_USER_ERROR_IF (!wind_u_field.empty(),
665 "This method does not yet support non-empty *wind_u_field*.");
666 ARTS_USER_ERROR_IF (!wind_v_field.empty(),
667 "This method does not yet support non-empty *wind_v_field*.");
668 ARTS_USER_ERROR_IF (!wind_w_field.empty(),
669 "This method does not yet support non-empty *wind_w_field*.");
670 ARTS_USER_ERROR_IF (!mag_u_field.empty(),
671 "This method does not yet support non-empty *mag_u_field*.");
672 ARTS_USER_ERROR_IF (!mag_v_field.empty(),
673 "This method does not yet support non-empty *mag_v_field*.");
674 ARTS_USER_ERROR_IF (!mag_w_field.empty(),
675 "This method does not yet support non-empty *mag_w_field*.");
676 //
677 if (return_masses) {
678 ARTS_USER_ERROR_IF (pnd_field.nbooks() != particle_masses.nrows(),
679 "Sizes of *pnd_field* and *particle_masses* "
680 "are inconsistent.");
681 }
682 chk_latlon_true(atmosphere_dim,
683 lat_grid,
684 lat_true,
685 lon_true);
686
687 // Note that input 1D atmospheres are handled exactly as 2D and 3D, to make
688 // the function totally general. And 1D must be handled for iterative calls.
689
690 // Determine propagation path (with cloudbox deactivated) and check
691 // that is OK for ICA
692 //
694 ppath,
695 ppath_lmax,
696 ppath_lraytrace,
697 rte_pos,
698 rte_los,
699 rte_pos2,
700 0,
701 0,
702 f_grid,
703 ppath_agenda);
704 //
705 error_if_limb_ppath(ppath);
706
707 // If scattering and sensor inside atmosphere, we need a pseudo-ppath that
708 // samples altitudes not covered by main ppath. We make this second path
709 // strictly vertical.
710 //
711 Ppath ppath2;
712 //
713 if (cloudbox_on && ppath.end_lstep == 0) {
714 Vector los_tmp = rte_los;
715 if (abs(rte_los[0]) < 90) {
716 los_tmp[0] = 180;
717 } else {
718 los_tmp[0] = 0;
719 }
720 //
722 ppath2,
723 ppath_lmax,
724 ppath_lraytrace,
725 rte_pos,
726 los_tmp,
727 rte_pos2,
728 0,
729 0,
730 f_grid,
731 ppath_agenda);
732 } else {
733 ppath2.np = 1;
734 }
735
736 // Merge grid positions, and sort from bottom to top of atmosphere
737 const Index np = ppath.np + ppath2.np - 1;
738 ArrayOfGridPos gp_p(np), gp_lat(np), gp_lon(np);
739 if (ppath.np > 1 &&
740 ppath.pos(0, 0) >
741 ppath.pos(1, 0)) { // Ppath is sorted in downward direction
742 // Copy ppath in reversed order
743 for (Index i = 0; i < ppath.np; i++) {
744 const Index ip = ppath.np - i - 1;
745 gp_p[i] = ppath.gp_p[ip];
746 if (atmosphere_dim > 1) {
747 gp_lat[i] = ppath.gp_lat[ip];
748 if (atmosphere_dim == 3) {
749 gp_lon[i] = ppath.gp_lon[ip];
750 }
751 }
752 }
753 // Append ppath2, but skipping element [0]
754 for (Index i = ppath.np; i < np; i++) {
755 const Index ip = i - ppath.np + 1;
756 gp_p[i] = ppath2.gp_p[ip];
757 if (atmosphere_dim > 1) {
758 gp_lat[i] = ppath2.gp_lat[ip];
759 if (atmosphere_dim == 3) {
760 gp_lon[i] = ppath2.gp_lon[ip];
761 }
762 }
763 }
764 } else {
765 // Copy ppath2 in reversed order, but skipping element [0]
766 for (Index i = 0; i < ppath2.np - 1; i++) {
767 const Index ip = ppath2.np - i - 1;
768 gp_p[i] = ppath2.gp_p[ip];
769 if (atmosphere_dim > 1) {
770 gp_lat[i] = ppath2.gp_lat[ip];
771 if (atmosphere_dim == 3) {
772 gp_lon[i] = ppath2.gp_lon[ip];
773 }
774 }
775 }
776 // Append ppath
777 for (Index i = ppath2.np - 1; i < np; i++) {
778 const Index ip = i - ppath2.np + 1;
779 gp_p[i] = ppath.gp_p[ip];
780 if (atmosphere_dim > 1) {
781 gp_lat[i] = ppath.gp_lat[ip];
782 if (atmosphere_dim == 3) {
783 gp_lon[i] = ppath.gp_lon[ip];
784 }
785 }
786 }
787 }
788
789 // 1D version of p_grid
790 Matrix itw;
791 Vector p1(np);
792 ArrayOfGridPos gp0(0), gp1(1);
793 interp_atmfield_gp2itw(itw, 1, gp_p, gp0, gp0);
794 itw2p(p1, p_grid, gp_p, itw);
795
796 // 1D version of lat and lon variables
797 Vector lat1(0), lon1(0);
798 Vector lat_true1(1), lon_true1(1);
799 if (atmosphere_dim == 3) {
800 gp1[0] = gp_lat[0];
801 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
802 interp(lat_true1, itw, lat_grid, gp1);
803 gp1[0] = gp_lon[0];
804 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
805 interp(lon_true1, itw, lon_grid, gp1);
806 } else if (atmosphere_dim == 2) {
807 gp1[0] = gp_lat[0];
808 interp_atmfield_gp2itw(itw, 1, gp1, gp0, gp0);
809 interp(lat_true1, itw, lat_true, gp1);
810 interp(lon_true1, itw, lon_true, gp1);
811 } else {
812 lat_true1[0] = lat_true[0];
813 lon_true1[0] = lon_true[0];
814 }
815
816 // 2D/3D interpolation weights
817 interp_atmfield_gp2itw(itw, atmosphere_dim, gp_p, gp_lat, gp_lon);
818
819 // 1D temperature field
820 Tensor3 t1(np, 1, 1);
822 t1(joker, 0, 0), atmosphere_dim, t_field, gp_p, gp_lat, gp_lon, itw);
823
824 // 1D altitude field
825 Tensor3 z1(np, 1, 1);
827 z1(joker, 0, 0), atmosphere_dim, z_field, gp_p, gp_lat, gp_lon, itw);
828
829 // 1D VMR field
830 Tensor4 vmr1(vmr_field.nbooks(), np, 1, 1);
831 for (Index is = 0; is < vmr_field.nbooks(); is++) {
832 interp_atmfield_by_itw(vmr1(is, joker, 0, 0),
833 atmosphere_dim,
834 vmr_field(is, joker, joker, joker),
835 gp_p,
836 gp_lat,
837 gp_lon,
838 itw);
839 }
840
841 // 1D surface altitude
842 Matrix zsurf1(1, 1);
843 zsurf1(0, 0) = z1(0, 0, 0);
844
845 // 1D version of rte_pos/los
846 Vector pos1(1);
847 pos1[0] = rte_pos[0];
848 Vector los1(1);
849 los1[0] = abs(rte_los[0]);
850 Vector pos2(0);
851 if (rte_pos2.nelem()) {
852 pos2 = rte_pos2[Range(0, rte_pos2.nelem())];
853 }
854
855 // Cloudbox variables
856 //
857 Index cbox_on1 = cloudbox_on;
858 ArrayOfIndex cbox_lims1(0);
859 Tensor4 pnd1(0, 0, 0, 0);
860 //
861 if (cloudbox_on) {
862 // Determine what p1-levels that are inside cloudbox
863 Index ifirst = np;
864 Index ilast = -1;
865 for (Index i = 0; i < np; i++) {
866 if (is_gp_inside_cloudbox(gp_p[i],
867 gp_lat[i],
868 gp_lon[i],
869 cloudbox_limits,
870 true,
871 atmosphere_dim)) {
872 if (i < ifirst) {
873 ifirst = i;
874 }
875 ilast = i;
876 }
877 }
878
879 // If no hit, deactive cloudbox
880 if (ifirst == np) {
881 cbox_on1 = 0;
882 }
883
884 // Otherwise set 1D cloud variables
885 else {
886 // We can enter the cloudbox from the side, and we need to add 1
887 // level on each side to be safe (if not limits already at edges)
888 //
889 const Index extra_bot = ifirst == 0 ? 0 : 1;
890 const Index extra_top = ilast == np - 1 ? 0 : 1;
891 //
892 cbox_lims1.resize(2);
893 cbox_lims1[0] = ifirst - extra_bot;
894 cbox_lims1[1] = ilast + extra_top;
895
896 // pnd_field
897 //
898 pnd1.resize(pnd_field.nbooks(), cbox_lims1[1] - cbox_lims1[0] + 1, 1, 1);
899 pnd1 = 0; // As lowermost and uppermost level not always filled
900 //
901 itw.resize(1, Index(pow(2.0, Numeric(atmosphere_dim))));
902 //
903 for (Index i = extra_bot; i < pnd1.npages() - extra_top; i++) {
904 const Index i0 = cbox_lims1[0] + i;
905 ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
907 gpc_p[0],
908 gpc_lat[0],
909 gpc_lon[0],
910 gp_p[i0],
911 gp_lat[i0],
912 gp_lon[i0],
913 atmosphere_dim,
914 cloudbox_limits);
915 for (Index p = 0; p < pnd_field.nbooks(); p++) {
916 interp_atmfield_by_itw(pnd1(p, i, 0, 0),
917 atmosphere_dim,
918 pnd_field(p, joker, joker, joker),
919 gpc_p,
920 gpc_lat,
921 gpc_lon,
922 itw);
923 }
924 }
925 }
926 }
927
928 // Call sub agenda
929 //
930 {
931 const Index adim1 = 1;
932 const Numeric lmax1 = -1;
933 Ppath ppath1d;
934 //
936 iy,
937 iy_aux,
938 ppath1d,
939 diy_dx,
940 iy_agenda_call1,
941 iy_unit,
942 iy_transmittance,
943 iy_aux_vars,
944 iy_id,
945 adim1,
946 p1,
947 lat1,
948 lon1,
949 lat_true1,
950 lon_true1,
951 t1,
952 z1,
953 vmr1,
954 zsurf1,
955 lmax1,
956 ppath_lraytrace,
957 cbox_on1,
958 cbox_lims1,
959 pnd1,
960 jacobian_do,
961 f_grid,
962 pos1,
963 los1,
964 pos2,
965 iy_independent_beam_approx_agenda);
966 }
967
968 // Fill *atm_fields_compact*?
969 if (return_atm1d) {
970 // Sizes and allocate memory
971 const Index nvmr = skip_vmr ? 0 : vmr1.nbooks();
972 const Index npnd = skip_pnd ? 0 : pnd1.nbooks();
973 const Index nmass = return_masses ? particle_masses.ncols() : 0;
974 const Index ntot = 2 + nvmr + npnd + nmass;
975 ArrayOfString field_names(ntot);
976 atm_fields_compact.resize(ntot, np, 1, 1);
977
978 // Altitudes
979 field_names[0] = "Geometric altitudes";
980 atm_fields_compact.data(0, joker, 0, 0) = z1(joker, 0, 0);
981
982 // Temperature
983 field_names[1] = "Temperature";
984 atm_fields_compact.data(1, joker, 0, 0) = t1(joker, 0, 0);
985
986 // VMRs
987 if (nvmr) {
988 for (Index i = 0; i < nvmr; i++) {
989 const Index iout = 2 + i;
990 ostringstream sstr;
991 sstr << "VMR species " << i;
992 field_names[iout] = sstr.str();
993 atm_fields_compact.data(iout, joker, 0, 0) = vmr1(i, joker, 0, 0);
994 }
995 }
996
997 // PNDs
998 if (npnd) {
999 for (Index i = 0; i < npnd; i++) {
1000 const Index iout = 2 + nvmr + i;
1001 ostringstream sstr;
1002 sstr << "Scattering element " << i;
1003 field_names[iout] = sstr.str();
1004 atm_fields_compact.data(iout, joker, 0, 0) = 0;
1005 atm_fields_compact.data(
1006 iout, Range(cbox_lims1[0], pnd1.npages()), 0, 0) =
1007 pnd1(i, joker, 0, 0);
1008 }
1009 }
1010
1011 // Masses
1012 if (nmass) {
1013 for (Index i = 0; i < nmass; i++) {
1014 const Index iout = 2 + nvmr + npnd + i;
1015 ostringstream sstr;
1016 sstr << "Mass category " << i;
1017 field_names[iout] = sstr.str();
1018 atm_fields_compact.data(iout, joker, 0, 0) = 0;
1019 for (Index ip = cbox_lims1[0]; ip < pnd1.npages(); ip++) {
1020 for (Index is = 0; is < pnd1.nbooks(); is++) {
1021 atm_fields_compact.data(iout, ip, 0, 0) +=
1022 particle_masses(is, i) * pnd1(is, ip, 0, 0);
1023 }
1024 }
1025 }
1026 }
1027
1028 // Finally, set grids and names
1029 //
1030 atm_fields_compact.set_name(
1031 "Data created by *iyIndependentBeamApproximation*");
1032 //
1033 atm_fields_compact.set_grid_name(GFIELD4_FIELD_NAMES,
1034 "Atmospheric quantity");
1035 atm_fields_compact.set_grid(GFIELD4_FIELD_NAMES, field_names);
1036 atm_fields_compact.set_grid_name(GFIELD4_P_GRID, "Pressure");
1037 atm_fields_compact.set_grid(GFIELD4_P_GRID, p1);
1038 atm_fields_compact.set_grid_name(GFIELD4_LAT_GRID, "Latitude");
1039 atm_fields_compact.set_grid(GFIELD4_LAT_GRID, lat_true1);
1040 atm_fields_compact.set_grid_name(GFIELD4_LON_GRID, "Longitude");
1041 atm_fields_compact.set_grid(GFIELD4_LON_GRID, lon_true1);
1042 }
1043}
1044
1045/* Workspace method: Doxygen documentation will be auto-generated */
1047 Matrix& iy,
1048 ArrayOfMatrix& iy_aux,
1049 Ppath& ppath,
1050 ArrayOfTensor3& diy_dx,
1051 const ArrayOfString& iy_aux_vars,
1052 const Index& iy_agenda_call1,
1053 const Tensor3& iy_transmittance,
1054 const Vector& rte_pos,
1055 const Vector& rte_los,
1056 const Vector& rte_pos2,
1057 const Index& stokes_dim,
1058 const Vector& f_grid,
1059 const Agenda& iy_loop_freqs_agenda,
1060 const Verbosity&) {
1061 // Throw error if unsupported features are requested
1062 ARTS_USER_ERROR_IF (!iy_agenda_call1,
1063 "Recursive usage not possible (iy_agenda_call1 must be 1).");
1064 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
1065 "*iy_transmittance* must be empty.");
1066
1067 const Index nf = f_grid.nelem();
1068
1069 for (Index i = 0; i < nf; i++) {
1070 // Variables for 1 frequency
1071 Matrix iy1;
1072 ArrayOfMatrix iy_aux1;
1073 ArrayOfTensor3 diy_dx1;
1074
1076 iy1,
1077 iy_aux1,
1078 ppath,
1079 diy_dx1,
1080 iy_agenda_call1,
1081 iy_transmittance,
1082 iy_aux_vars,
1083 0,
1084 Vector(1, f_grid[i]),
1085 rte_pos,
1086 rte_los,
1087 rte_pos2,
1088 iy_loop_freqs_agenda);
1089
1090 // After first frequency, give output its size
1091 if (i == 0) {
1092 iy.resize(nf, stokes_dim);
1093 //
1094 iy_aux.resize(iy_aux1.nelem());
1095 for (Index q = 0; q < iy_aux1.nelem(); q++) {
1096 iy_aux[q].resize(nf, stokes_dim);
1097 }
1098 //
1099 diy_dx.resize(diy_dx1.nelem());
1100 for (Index q = 0; q < diy_dx1.nelem(); q++) {
1101 diy_dx[q].resize(diy_dx1[q].npages(), nf, stokes_dim);
1102 }
1103 }
1104
1105 // Copy to output variables
1106 iy(i, joker) = iy1(0, joker);
1107 for (Index q = 0; q < iy_aux1.nelem(); q++) {
1108 iy_aux[q](i, joker) = iy_aux1[q](0, joker);
1109 }
1110 for (Index q = 0; q < diy_dx1.nelem(); q++) {
1111 diy_dx[q](joker, i, joker) = diy_dx1[q](joker, 0, joker);
1112 }
1113 }
1114}
1115
1116/* Workspace method: Doxygen documentation will be auto-generated */
1118 Matrix& iy,
1119 ArrayOfMatrix& iy_aux,
1120 ArrayOfTensor3& diy_dx,
1121 const Index& iy_agenda_call1,
1122 const Tensor3& iy_transmittance,
1123 const Vector& rte_pos,
1124 const Vector& rte_los,
1125 const ArrayOfString& iy_aux_vars,
1126 const Index& jacobian_do,
1127 const Index& atmosphere_dim,
1128 const Vector& p_grid,
1129 const Vector& lat_grid,
1130 const Vector& lon_grid,
1131 const Tensor3& z_field,
1132 const Tensor3& t_field,
1133 const Tensor4& vmr_field,
1134 const Vector& refellipsoid,
1135 const Matrix& z_surface,
1136 const Index& cloudbox_on,
1137 const ArrayOfIndex& cloudbox_limits,
1138 const Index& stokes_dim,
1139 const Vector& f_grid,
1140 const ArrayOfArrayOfSingleScatteringData& scat_data,
1141 const Agenda& iy_space_agenda,
1142 const Agenda& surface_rtprop_agenda,
1143 const Agenda& propmat_clearsky_agenda,
1144 const Agenda& ppath_step_agenda,
1145 const Numeric& ppath_lmax,
1146 const Numeric& ppath_lraytrace,
1147 const Tensor4& pnd_field,
1148 const String& iy_unit,
1149 const Numeric& mc_std_err,
1150 const Index& mc_max_time,
1151 const Index& mc_max_iter,
1152 const Index& mc_min_iter,
1153 const Numeric& mc_taustep_limit,
1154 const Index& t_interp_order,
1155 const Verbosity& verbosity) {
1156 // Throw error if unsupported features are requested
1157 ARTS_USER_ERROR_IF (atmosphere_dim != 3,
1158 "Only 3D atmospheres are allowed (atmosphere_dim must be 3)");
1159 ARTS_USER_ERROR_IF (!cloudbox_on,
1160 "The cloudbox must be activated (cloudbox_on must be 1)");
1161 ARTS_USER_ERROR_IF (jacobian_do,
1162 "This method does not provide any jacobians (jacobian_do must be 0)");
1163 ARTS_USER_ERROR_IF (!iy_agenda_call1,
1164 "Recursive usage not possible (iy_agenda_call1 must be 1)");
1165 ARTS_USER_ERROR_IF (iy_transmittance.ncols(),
1166 "*iy_transmittance* must be empty");
1167
1168 // Size output variables
1169 //
1170 const Index nf = f_grid.nelem();
1171 //
1172 iy.resize(nf, stokes_dim);
1173 diy_dx.resize(0);
1174 //
1175 //=== iy_aux part ===========================================================
1176 Index auxError = -1;
1177 {
1178 const Index naux = iy_aux_vars.nelem();
1179 iy_aux.resize(naux);
1180 //
1181 for (Index i = 0; i < naux; i++) {
1182 if (iy_aux_vars[i] == "Error (uncorrelated)") {
1183 auxError = i;
1184 iy_aux[i].resize(nf, stokes_dim);
1185 } else {
1187 "In *iy_aux_vars* you have included: \"", iy_aux_vars[i],
1188 "\"\nThis choice is not recognised.")
1189 }
1190 }
1191 }
1192 //===========================================================================
1193
1194 // Some MC variables are only local here
1195 //
1196 MCAntenna mc_antenna;
1197 mc_antenna.set_pencil_beam();
1198
1199 // Pos and los must be matrices
1200 Matrix pos(1, 3), los(1, 2);
1201 //
1202 pos(0, joker) = rte_pos;
1203 los(0, joker) = rte_los;
1204
1205 Workspace l_ws(ws);
1206 Agenda l_ppath_step_agenda(ppath_step_agenda);
1207 Agenda l_iy_space_agenda(iy_space_agenda);
1208 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
1209 Agenda l_surface_rtprop_agenda(surface_rtprop_agenda);
1210
1211 String fail_msg;
1212 bool failed = false;
1213
1214 if (nf)
1215#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) \
1216 firstprivate(l_ws, \
1217 l_ppath_step_agenda, \
1218 l_iy_space_agenda, \
1219 l_propmat_clearsky_agenda, \
1220 l_surface_rtprop_agenda)
1221 for (Index f_index = 0; f_index < nf; f_index++) {
1222 if (failed) continue;
1223
1224 try {
1225 // Seed reset for each loop. If not done, the errors
1226 // appear to be highly correlated.
1227 Index mc_seed;
1228 MCSetSeedFromTime(mc_seed, verbosity);
1229
1230 Vector y, mc_error;
1231 Index mc_iteration_count;
1232 Tensor3 mc_points;
1233 ArrayOfIndex mc_scat_order, mc_source_domain;
1234
1235 MCGeneral(l_ws,
1236 y,
1237 mc_iteration_count,
1238 mc_error,
1239 mc_points,
1240 mc_scat_order,
1241 mc_source_domain,
1242 mc_antenna,
1243 f_grid,
1244 f_index,
1245 pos,
1246 los,
1247 stokes_dim,
1248 atmosphere_dim,
1249 l_ppath_step_agenda,
1250 ppath_lmax,
1251 ppath_lraytrace,
1252 l_iy_space_agenda,
1253 l_surface_rtprop_agenda,
1254 l_propmat_clearsky_agenda,
1255 p_grid,
1256 lat_grid,
1257 lon_grid,
1258 z_field,
1259 refellipsoid,
1260 z_surface,
1261 t_field,
1262 vmr_field,
1263 cloudbox_on,
1264 cloudbox_limits,
1265 pnd_field,
1266 scat_data,
1267 1,
1268 1,
1269 1,
1270 1,
1271 iy_unit,
1272 mc_seed,
1273 mc_std_err,
1274 mc_max_time,
1275 mc_max_iter,
1276 mc_min_iter,
1277 mc_taustep_limit,
1278 1,
1279 t_interp_order,
1280 verbosity);
1281
1282 ARTS_ASSERT(y.nelem() == stokes_dim);
1283
1284 iy(f_index, joker) = y;
1285
1286 if (auxError >= 0) {
1287 iy_aux[auxError](f_index, joker) = mc_error;
1288 }
1289 } catch (const std::exception& e) {
1290 ostringstream os;
1291 os << "Error for f_index = " << f_index << " (" << f_grid[f_index]
1292 << ")" << endl
1293 << e.what();
1294#pragma omp critical(iyMC_fail)
1295 {
1296 failed = true;
1297 fail_msg = os.str();
1298 }
1299 continue;
1300 }
1301 }
1302
1303 ARTS_USER_ERROR_IF (failed, fail_msg);
1304}
1305
1306/* Workspace method: Doxygen documentation will be auto-generated */
1308 const ArrayOfMatrix& iy_aux,
1309 const ArrayOfString& iy_aux_vars,
1310 const Index& jacobian_do,
1311 const String& aux_var,
1312 const Verbosity&) {
1313 ARTS_USER_ERROR_IF (iy_aux.nelem() != iy_aux_vars.nelem(),
1314 "*iy_aux* and *iy_aux_vars* must have the same "
1315 "number of elements.");
1316
1317 ARTS_USER_ERROR_IF (jacobian_do,
1318 "This method can not provide any jacobians and "
1319 "*jacobian_do* must be 0.");
1320
1321 bool ready = false;
1322
1323 for (Index i = 0; i < iy_aux.nelem() && !ready; i++) {
1324 if (iy_aux_vars[i] == aux_var) {
1325 iy = iy_aux[i];
1326 ready = true;
1327 }
1328 }
1329
1330 ARTS_USER_ERROR_IF (!ready,
1331 "The selected auxiliary variable to insert in *iy* "
1332 "is either not defined at all or is not set.");
1333}
1334
1335/* Workspace method: Doxygen documentation will be auto-generated */
1337 Matrix& ppvar_optical_depth,
1338 const Tensor4& ppvar_trans_cumulat,
1339 const Verbosity&) {
1340 ppvar_optical_depth = ppvar_trans_cumulat(joker, joker, 0, 0);
1341 transform(ppvar_optical_depth, log, ppvar_optical_depth);
1342 ppvar_optical_depth *= -1;
1343}
1344
1345/* Workspace method: Doxygen documentation will be auto-generated */
1347 Vector& y,
1348 Vector& y_f,
1349 ArrayOfIndex& y_pol,
1350 Matrix& y_pos,
1351 Matrix& y_los,
1352 ArrayOfVector& y_aux,
1353 Matrix& y_geo,
1354 Matrix& jacobian,
1355 const Index& atmgeom_checked,
1356 const Index& atmfields_checked,
1357 const Index& atmosphere_dim,
1358 const EnergyLevelMap& nlte_field,
1359 const Index& cloudbox_on,
1360 const Index& cloudbox_checked,
1361 const Index& scat_data_checked,
1362 const Index& sensor_checked,
1363 const Index& stokes_dim,
1364 const Vector& f_grid,
1365 const Matrix& sensor_pos,
1366 const Matrix& sensor_los,
1367 const Matrix& transmitter_pos,
1368 const Matrix& mblock_dlos_grid,
1369 const Sparse& sensor_response,
1370 const Vector& sensor_response_f,
1371 const ArrayOfIndex& sensor_response_pol,
1372 const Matrix& sensor_response_dlos,
1373 const String& iy_unit,
1374 const Agenda& iy_main_agenda,
1375 const Agenda& geo_pos_agenda,
1376 const Agenda& jacobian_agenda,
1377 const Index& jacobian_do,
1378 const ArrayOfRetrievalQuantity& jacobian_quantities,
1379 const ArrayOfString& iy_aux_vars,
1380 const Verbosity& verbosity) {
1382
1383 // Basics
1384 //
1385 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1386 //
1387 ARTS_USER_ERROR_IF (f_grid.empty(),
1388 "The frequency grid is empty.");
1389 chk_if_increasing("f_grid", f_grid);
1390 ARTS_USER_ERROR_IF (f_grid[0] <= 0,
1391 "All frequencies in *f_grid* must be > 0.");
1392 //
1393 ARTS_USER_ERROR_IF (atmfields_checked != 1,
1394 "The atmospheric fields must be flagged to have\n"
1395 "passed a consistency check (atmfields_checked=1).");
1396 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
1397 "The atmospheric geometry must be flagged to have\n"
1398 "passed a consistency check (atmgeom_checked=1).");
1399 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
1400 "The cloudbox must be flagged to have\n"
1401 "passed a consistency check (cloudbox_checked=1).");
1402 if (cloudbox_on)
1403 ARTS_USER_ERROR_IF (scat_data_checked != 1,
1404 "The scattering data must be flagged to have\n"
1405 "passed a consistency check (scat_data_checked=1).");
1406 ARTS_USER_ERROR_IF (sensor_checked != 1,
1407 "The sensor variables must be flagged to have\n"
1408 "passed a consistency check (sensor_checked=1).");
1409
1410 // Some sizes
1411 const Index nf = f_grid.nelem();
1412 const Index nlos = mblock_dlos_grid.nrows();
1413 const Index n1y = sensor_response.nrows();
1414 const Index nmblock = sensor_pos.nrows();
1415 const Index niyb = nf * nlos * stokes_dim;
1416
1417 //---------------------------------------------------------------------------
1418 // Allocations and resizing
1419 //---------------------------------------------------------------------------
1420
1421 // Resize *y* and *y_XXX*
1422 //
1423 y.resize(nmblock * n1y);
1424 y_f.resize(nmblock * n1y);
1425 y_pol.resize(nmblock * n1y);
1426 y_pos.resize(nmblock * n1y, sensor_pos.ncols());
1427 y_los.resize(nmblock * n1y, sensor_los.ncols());
1428 y_geo.resize(nmblock * n1y, 5);
1429 y_geo = NAN; // Will be replaced if relavant data are provided (*geo_pos*)
1430
1431 // For y_aux we don't know the number of quantities, and we need to
1432 // store all output
1433 ArrayOfArrayOfVector iyb_aux_array(nmblock);
1434
1435 // Jacobian variables
1436 //
1437 Index j_analytical_do = 0;
1438 ArrayOfArrayOfIndex jacobian_indices;
1439 //
1440 if (jacobian_do) {
1441 bool any_affine;
1442 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
1443 jacobian.resize(nmblock * n1y,
1444 jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
1445 jacobian = 0;
1446 //
1447 FOR_ANALYTICAL_JACOBIANS_DO2(j_analytical_do = 1;)
1448 } else {
1449 jacobian.resize(0, 0);
1450 }
1451
1452 //---------------------------------------------------------------------------
1453 // The calculations
1454 //---------------------------------------------------------------------------
1455
1456 String fail_msg;
1457 bool failed = false;
1458
1459 if (nmblock >= arts_omp_get_max_threads() ||
1460 (nf <= nmblock && nmblock >= nlos)) {
1461 out3 << " Parallelizing mblock loop (" << nmblock << " iterations)\n";
1462
1463 // We have to make a local copy of the Workspace and the agendas because
1464 // only non-reference types can be declared firstprivate in OpenMP
1465 Workspace l_ws(ws);
1466 Agenda l_jacobian_agenda(jacobian_agenda);
1467 Agenda l_iy_main_agenda(iy_main_agenda);
1468 Agenda l_geo_pos_agenda(geo_pos_agenda);
1469
1470#pragma omp parallel for firstprivate( \
1471 l_ws, l_jacobian_agenda, l_iy_main_agenda, l_geo_pos_agenda)
1472 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
1473 // Skip remaining iterations if an error occurred
1474 if (failed) continue;
1475
1477 fail_msg,
1478 iyb_aux_array,
1479 l_ws,
1480 y,
1481 y_f,
1482 y_pol,
1483 y_pos,
1484 y_los,
1485 y_geo,
1486 jacobian,
1487 atmosphere_dim,
1488 nlte_field,
1489 cloudbox_on,
1490 stokes_dim,
1491 f_grid,
1492 sensor_pos,
1493 sensor_los,
1494 transmitter_pos,
1495 mblock_dlos_grid,
1496 sensor_response,
1497 sensor_response_f,
1498 sensor_response_pol,
1499 sensor_response_dlos,
1500 iy_unit,
1501 l_iy_main_agenda,
1502 l_geo_pos_agenda,
1503 l_jacobian_agenda,
1504 jacobian_do,
1505 jacobian_quantities,
1506 jacobian_indices,
1507 iy_aux_vars,
1508 verbosity,
1509 mblock_index,
1510 n1y,
1511 j_analytical_do);
1512 } // End mblock loop
1513 } else {
1514 out3 << " Not parallelizing mblock loop (" << nmblock << " iterations)\n";
1515
1516 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
1517 // Skip remaining iterations if an error occurred
1518 if (failed) continue;
1519
1521 fail_msg,
1522 iyb_aux_array,
1523 ws,
1524 y,
1525 y_f,
1526 y_pol,
1527 y_pos,
1528 y_los,
1529 y_geo,
1530 jacobian,
1531 atmosphere_dim,
1532 nlte_field,
1533 cloudbox_on,
1534 stokes_dim,
1535 f_grid,
1536 sensor_pos,
1537 sensor_los,
1538 transmitter_pos,
1539 mblock_dlos_grid,
1540 sensor_response,
1541 sensor_response_f,
1542 sensor_response_pol,
1543 sensor_response_dlos,
1544 iy_unit,
1545 iy_main_agenda,
1546 geo_pos_agenda,
1547 jacobian_agenda,
1548 jacobian_do,
1549 jacobian_quantities,
1550 jacobian_indices,
1551 iy_aux_vars,
1552 verbosity,
1553 mblock_index,
1554 n1y,
1555 j_analytical_do);
1556 } // End mblock loop
1557 }
1558
1559 // Rethrow exception if a runtime error occurred in the mblock loop
1560 ARTS_USER_ERROR_IF (failed, fail_msg);
1561
1562 // Compile y_aux
1563 //
1564 const Index nq = iyb_aux_array[0].nelem();
1565 y_aux.resize(nq);
1566 //
1567 for (Index q = 0; q < nq; q++) {
1568 y_aux[q].resize(nmblock * n1y);
1569 //
1570 for (Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
1571 const Range rowind =
1572 get_rowindex_for_mblock(sensor_response, mblock_index);
1573 const Index row0 = rowind.get_start();
1574
1575 // The sensor response must be applied in a special way for
1576 // uncorrelated errors. Schematically: sqrt( H.^2 * y.^2 )
1577 if (iy_aux_vars[q] == "Error (uncorrelated)") {
1578 for (Index i = 0; i < n1y; i++) {
1579 const Index row = row0 + i;
1580 y_aux[q][row] = 0;
1581 for (Index j = 0; j < niyb; j++) {
1582 y_aux[q][row] +=
1583 pow(sensor_response(i, j) * iyb_aux_array[mblock_index][q][j],
1584 (Numeric)2.0);
1585 }
1586 y_aux[q][row] = sqrt(y_aux[q][row]);
1587 }
1588 } else {
1589 mult(y_aux[q][rowind], sensor_response, iyb_aux_array[mblock_index][q]);
1590 }
1591 }
1592 }
1593}
1594
1595/* Workspace method: Doxygen documentation will be auto-generated */
1597 Vector& y,
1598 Vector& y_f,
1599 ArrayOfIndex& y_pol,
1600 Matrix& y_pos,
1601 Matrix& y_los,
1602 ArrayOfVector& y_aux,
1603 Matrix& y_geo,
1604 Matrix& jacobian,
1605 ArrayOfRetrievalQuantity& jacobian_quantities,
1606 const Index& atmfields_checked,
1607 const Index& atmgeom_checked,
1608 const Index& atmosphere_dim,
1609 const EnergyLevelMap& nlte_field,
1610 const Index& cloudbox_on,
1611 const Index& cloudbox_checked,
1612 const Index& scat_data_checked,
1613 const Index& sensor_checked,
1614 const Index& stokes_dim,
1615 const Vector& f_grid,
1616 const Matrix& sensor_pos,
1617 const Matrix& sensor_los,
1618 const Matrix& transmitter_pos,
1619 const Matrix& mblock_dlos_grid,
1620 const Sparse& sensor_response,
1621 const Vector& sensor_response_f,
1622 const ArrayOfIndex& sensor_response_pol,
1623 const Matrix& sensor_response_dlos,
1624 const String& iy_unit,
1625 const Agenda& iy_main_agenda,
1626 const Agenda& geo_pos_agenda,
1627 const Agenda& jacobian_agenda,
1628 const Index& jacobian_do,
1629 const ArrayOfString& iy_aux_vars,
1630 const ArrayOfRetrievalQuantity& jacobian_quantities_copy,
1631 const Index& append_instrument_wfs,
1632 const Verbosity& verbosity) {
1633 // The jacobian indices of old and new part (without transformations)
1634 ArrayOfArrayOfIndex jacobian_indices, jacobian_indices_copy;
1635 {
1636 bool any_affine;
1638 jacobian_indices_copy, any_affine, jacobian_quantities_copy, true);
1639 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
1640 }
1641
1642 // Check consistency of data representing first measurement
1643 const Index n1 = y.nelem();
1644 Index nrq1 = 0;
1646 "Input *y* is empty. Use *yCalc*");
1647 ARTS_USER_ERROR_IF (y_f.nelem() != n1,
1648 "Lengths of input *y* and *y_f* are inconsistent.");
1649 ARTS_USER_ERROR_IF (y_pol.nelem() != n1,
1650 "Lengths of input *y* and *y_pol* are inconsistent.");
1651 ARTS_USER_ERROR_IF (y_pos.nrows() != n1,
1652 "Sizes of input *y* and *y_pos* are inconsistent.");
1653 ARTS_USER_ERROR_IF (y_los.nrows() != n1,
1654 "Sizes of input *y* and *y_los* are inconsistent.");
1655 ARTS_USER_ERROR_IF (y_geo.nrows() != n1,
1656 "Sizes of input *y* and *y_geo* are inconsistent.");
1657 if (jacobian_do) {
1658 nrq1 = jacobian_quantities_copy.nelem();
1659 ARTS_USER_ERROR_IF (jacobian.nrows() != n1,
1660 "Sizes of *y* and *jacobian* are inconsistent.");
1661 ARTS_USER_ERROR_IF (jacobian.ncols() != jacobian_indices_copy[nrq1 - 1][1] + 1,
1662 "Size of input *jacobian* and size implied "
1663 "*jacobian_quantities_copy* are inconsistent.");
1664 }
1665
1666 // Calculate new measurement
1667 //
1668 Vector y2, y_f2;
1669 Matrix y_pos2, y_los2, y_geo2, jacobian2;
1670 ArrayOfIndex y_pol2;
1671 ArrayOfVector y_aux2;
1672 //
1673 yCalc(ws,
1674 y2,
1675 y_f2,
1676 y_pol2,
1677 y_pos2,
1678 y_los2,
1679 y_aux2,
1680 y_geo2,
1681 jacobian2,
1682 atmfields_checked,
1683 atmgeom_checked,
1684 atmosphere_dim,
1685 nlte_field,
1686 cloudbox_on,
1687 cloudbox_checked,
1688 scat_data_checked,
1689 sensor_checked,
1690 stokes_dim,
1691 f_grid,
1692 sensor_pos,
1693 sensor_los,
1694 transmitter_pos,
1695 mblock_dlos_grid,
1696 sensor_response,
1697 sensor_response_f,
1698 sensor_response_pol,
1699 sensor_response_dlos,
1700 iy_unit,
1701 iy_main_agenda,
1702 geo_pos_agenda,
1703 jacobian_agenda,
1704 jacobian_do,
1705 jacobian_quantities,
1706 iy_aux_vars,
1707 verbosity);
1708
1709 // Consistency checks
1710 ARTS_USER_ERROR_IF (y_pos.ncols() != y_pos2.ncols(),
1711 "Different number of columns in *y_pos* between the measurements.");
1712 ARTS_USER_ERROR_IF (y_los.ncols() != y_los2.ncols(),
1713 "Different number of columns in *y_los* between the measurements.");
1714
1715 // y and y_XXX
1716 //
1717 const Index n2 = y2.nelem();
1718 //
1719 {
1720 // Make copy of old measurement
1721 const Vector y1 = y, y_f1 = y_f;
1722 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
1723 const ArrayOfIndex y_pol1 = y_pol;
1724 const ArrayOfVector y_aux1 = y_aux;
1725 //
1726 y.resize(n1 + n2);
1727 y[Range(0, n1)] = y1;
1728 y[Range(n1, n2)] = y2;
1729 //
1730 y_f.resize(n1 + n2);
1731 y_f[Range(0, n1)] = y_f1;
1732 y_f[Range(n1, n2)] = y_f2;
1733 //
1734 y_pos.resize(n1 + n2, y_pos1.ncols());
1735 y_pos(Range(0, n1), joker) = y_pos1;
1736 y_pos(Range(n1, n2), joker) = y_pos2;
1737 //
1738 y_los.resize(n1 + n2, y_los1.ncols());
1739 y_los(Range(0, n1), joker) = y_los1;
1740 y_los(Range(n1, n2), joker) = y_los2;
1741 //
1742 y_geo.resize(n1 + n2, y_geo1.ncols());
1743 y_geo(Range(0, n1), joker) = y_geo1;
1744 y_geo(Range(n1, n2), joker) = y_geo2;
1745 //
1746 y_pol.resize(n1 + n2);
1747 for (Index i = 0; i < n1; i++) {
1748 y_pol[i] = y_pol1[i];
1749 }
1750 for (Index i = 0; i < n2; i++) {
1751 y_pol[n1 + i] = y_pol2[i];
1752 }
1753
1754 // y_aux
1755 const Index na1 = y_aux1.nelem();
1756 const Index na2 = y_aux2.nelem();
1757 const Index na = max(na1, na2);
1758 //
1759 y_aux.resize(na);
1760 //
1761 for (Index a = 0; a < na; a++) {
1762 y_aux[a].resize(n1 + n2);
1763 if (a < na1) {
1764 y_aux[a][Range(0, n1)] = y_aux1[a];
1765 } else {
1766 y_aux[a][Range(0, n1)] = 0;
1767 }
1768 if (a < na2) {
1769 y_aux[a][Range(n1, n2)] = y_aux2[a];
1770 } else {
1771 y_aux[a][Range(n1, n2)] = 0;
1772 }
1773 }
1774 }
1775
1776 // Jacobian and friends
1777 if (jacobian_do) {
1778 // Put in old jacobian_quantities and jacobian_indices as first part in
1779 // new version of these variables
1780 ArrayOfRetrievalQuantity jacobian_quantities2 = jacobian_quantities;
1781 ArrayOfArrayOfIndex jacobian_indices2 = jacobian_indices;
1782 //
1783 jacobian_quantities = jacobian_quantities_copy;
1784 jacobian_indices = jacobian_indices_copy;
1785
1786 // Loop new jacobian_quantities to determine how new jacobian data shall
1787 // be inserted
1788 //
1789 const Index nrq2 = jacobian_quantities2.nelem();
1790 ArrayOfIndex map_table(nrq2);
1791 //
1792 for (Index q2 = 0; q2 < nrq2; q2++) {
1793 Index pos = -1;
1794
1795 // Compare to old quantities, to determine if append shall be
1796 // considered. Some special checks performed here, grids checked later
1797 if (jacobian_quantities2[q2].Target().isSpeciesVMR() ||
1798 jacobian_quantities2[q2] == Jacobian::Atm::Temperature ||
1799 jacobian_quantities2[q2] == Jacobian::Special::ScatteringString ||
1800 jacobian_quantities2[q2].Target().isWind() ||
1801 jacobian_quantities2[q2] == Jacobian::Special::SurfaceString ||
1802 append_instrument_wfs) {
1803 for (Index q1 = 0; q1 < nrq1; q1++ && pos < 0) { // FIXME: What is with this "&& pos < 0"
1804 if (jacobian_quantities2[q2].Target().sameTargetType(jacobian_quantities_copy[q1].Target())) {
1805 if (jacobian_quantities2[q2].Target().isSpeciesVMR()) {
1806 if (jacobian_quantities2[q2].Subtag() ==
1807 jacobian_quantities_copy[q1].Subtag()) {
1808 if (jacobian_quantities2[q2].Mode() ==
1809 jacobian_quantities_copy[q1].Mode()) {
1810 pos = q1;
1811 } else {
1813 "Jacobians for ", jacobian_quantities2[q2],
1814 " shall be appended.\nThis requires "
1815 "that the same retrieval unit is used "
1816 "but it seems that this requirement is "
1817 "not met.")
1818 }
1819 }
1820 } else if (jacobian_quantities2[q2] == Jacobian::Atm::Temperature) {
1821 if (jacobian_quantities2[q2].Subtag() ==
1822 jacobian_quantities_copy[q1].Subtag()) {
1823 pos = q1;
1824 } else {
1826 "Jacobians for ", jacobian_quantities2[q2],
1827 " shall be appended.\nThis requires "
1828 "that HSE is either ON or OFF for both "
1829 "parts but it seems that this requirement "
1830 "is not met.")
1831 }
1832 } else if (jacobian_quantities[q2] == Jacobian::Special::ScatteringString) {
1833 if ((jacobian_quantities2[q2].Subtag() ==
1834 jacobian_quantities_copy[q1].Subtag()) &&
1835 (jacobian_quantities2[q2].SubSubtag() ==
1836 jacobian_quantities_copy[q1].SubSubtag())) {
1837 pos = q1;
1838 }
1839 } else if (jacobian_quantities2[q2].Subtag() == jacobian_quantities_copy[q1].Subtag()) {
1840 pos = q1;
1841 }
1842 }
1843 }
1844 }
1845
1846 // New quantity
1847 if (pos < 0) {
1848 map_table[q2] = jacobian_quantities.nelem();
1849 jacobian_quantities.push_back(jacobian_quantities2[q2]);
1850 ArrayOfIndex indices(2);
1851 indices[0] = jacobian_indices[jacobian_indices.nelem() - 1][1] + 1;
1852 indices[1] =
1853 indices[0] + jacobian_indices2[q2][1] - jacobian_indices2[q2][0];
1854 jacobian_indices.push_back(indices);
1855 }
1856 // Existing quantity
1857 else {
1858 map_table[q2] = pos;
1859 // Check if grids are equal
1860 ArrayOfVector grids1 = jacobian_quantities_copy[pos].Grids();
1861 ArrayOfVector grids2 = jacobian_quantities2[q2].Grids();
1862 bool any_wrong = false;
1863 if (grids1.nelem() != grids2.nelem()) {
1864 any_wrong = true;
1865 } else {
1866 for (Index g = 0; g < grids1.nelem(); g++) {
1867 if (grids1[g].nelem() != grids2[g].nelem()) {
1868 any_wrong = true;
1869 } else {
1870 for (Index e = 0; e < grids1[g].nelem(); e++) {
1871 const Numeric v1 = grids1[g][e];
1872 const Numeric v2 = grids2[g][e];
1873 if ((v1 == 0 && abs(v2) > 1e-9) || abs(v1 - v2) / v1 > 1e-6) {
1874 any_wrong = true;
1875 }
1876 }
1877 }
1878 }
1879 }
1880 if (any_wrong) {
1882 "Jacobians for ", jacobian_quantities2[q2],
1883 " shall be appended.\nThis requires that the "
1884 "same grids are used for both measurements,\nbut "
1885 "it seems that this requirement is not met.")
1886 }
1887 }
1888 }
1889
1890 // Create and fill *jacobian*
1891 //
1892 const Index nrq = jacobian_quantities.nelem();
1893 const Matrix jacobian1 = jacobian;
1894 //
1895 jacobian.resize(n1 + n2, jacobian_indices[nrq - 1][1] + 1);
1896 jacobian = 0;
1897 //
1898 // Put in old part in top-left corner
1899 jacobian(Range(0, n1), Range(0, jacobian_indices_copy[nrq1 - 1][1] + 1)) =
1900 jacobian1;
1901 // New parts
1902 for (Index q2 = 0; q2 < nrq2; q2++) {
1903 jacobian(Range(n1, n2),
1904 Range(jacobian_indices[map_table[q2]][0],
1905 jacobian_indices[map_table[q2]][1] -
1906 jacobian_indices[map_table[q2]][0] + 1)) =
1907 jacobian2(
1908 joker,
1909 Range(jacobian_indices2[q2][0],
1910 jacobian_indices2[q2][1] - jacobian_indices2[q2][0] + 1));
1911 }
1912 }
1913}
1914
1915/* Workspace method: Doxygen documentation will be auto-generated */
1917 Matrix& jacobian,
1918 const Vector& y_f,
1919 const ArrayOfIndex& y_pol,
1920 const String& iy_unit,
1921 const Verbosity&) {
1922 ARTS_USER_ERROR_IF (iy_unit == "1",
1923 "No need to use this method with *iy_unit* = \"1\".");
1924
1925 ARTS_USER_ERROR_IF (max(y) > 1e-3,
1926 "The spectrum vector *y* is required to have original radiance\n"
1927 "unit, but this seems not to be the case. This as a value above\n"
1928 "1e-3 is found in *y*.")
1929
1930 // Is jacobian set?
1931 //
1932 const Index ny = y.nelem();
1933 //
1934 const bool do_j = jacobian.nrows() == ny;
1935
1936 // Some jacobian quantities can not be handled
1937 ARTS_USER_ERROR_IF (do_j && max(jacobian) > 1e-3,
1938 "The method can not be used with jacobian quantities that are not\n"
1939 "obtained through radiative transfer calculations. One example on\n"
1940 "quantity that can not be handled is *jacobianAddPolyfit*.\n"
1941 "The maximum value of *jacobian* indicates that one or several\n"
1942 "such jacobian quantities are included.")
1943
1944 // Planck-Tb
1945 //--------------------------------------------------------------------------
1946 if (iy_unit == "PlanckBT") {
1947 // Hard to use telescoping here as the data are sorted differently in y
1948 // and jacobian, than what is expected apply_iy_unit. Copy to temporary
1949 // variables instead.
1950
1951 // Handle the elements in "frequency chunks"
1952
1953 Index i0 = 0; // Index of first element for present chunk
1954 //
1955 while (i0 < ny) {
1956 // Find number of values for this chunk
1957 Index n = 1;
1958 //
1959 while (i0 + n < ny && y_f[i0] == y_f[i0 + n]) {
1960 n++;
1961 }
1962
1963 Matrix yv(1, n);
1964 ArrayOfIndex i_pol(n);
1965 bool any_quv = false;
1966 //
1967 for (Index i = 0; i < n; i++) {
1968 const Index ix = i0 + i;
1969 yv(0, i) = y[ix];
1970 i_pol[i] = y_pol[ix];
1971 if (i_pol[i] > 1 && i_pol[i] < 5) {
1972 any_quv = true;
1973 }
1974 }
1975
1976 // Index of elements to convert
1977 Range ii(i0, n);
1978
1979 if (do_j) {
1980 ARTS_USER_ERROR_IF (any_quv && i_pol[0] != 1,
1981 "The conversion to PlanckBT, of the Jacobian and "
1982 "errors for Q, U and V, requires that I (first Stokes "
1983 "element) is at hand and that the data are sorted in "
1984 "such way that I comes first for each frequency.")
1985
1986 // Jacobian
1987 Tensor3 J(jacobian.ncols(), 1, n);
1988 J(joker, 0, joker) = transpose(jacobian(ii, joker));
1989 apply_iy_unit2(J, yv, iy_unit, y_f[i0], 1, i_pol);
1990 jacobian(ii, joker) = transpose(J(joker, 0, joker));
1991 }
1992
1993 // y (must be done last)
1994 apply_iy_unit(yv, iy_unit, y_f[i0], 1, i_pol);
1995 y[ii] = yv(0, joker);
1996
1997 i0 += n;
1998 }
1999 }
2000
2001 // Other conversions
2002 //--------------------------------------------------------------------------
2003 else {
2004 // Here we take each element of y separately.
2005
2006 Matrix yv(1, 1);
2007 ArrayOfIndex i_pol(1);
2008
2009 for (Index i = 0; i < ny; i++) {
2010 yv(0, 0) = y[i];
2011 i_pol[0] = y_pol[i];
2012
2013 // Jacobian
2014 if (do_j) {
2016 MatrixView(jacobian(i, joker)), yv, iy_unit, y_f[i], 1, i_pol);
2017 }
2018
2019 // y (must be done last)
2020 apply_iy_unit(yv, iy_unit, y_f[i], 1, i_pol);
2021 y[i] = yv(0, 0);
2022 }
2023 }
2024}
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:39
Array< Tensor3 > ArrayOfTensor3
An array of Tensor3.
Definition: array.h:73
The global header file for ARTS.
Index npages
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:23454
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:23520
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:24134
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:23337
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 &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 &l_mc_scat_order, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
Definition: m_montecarlo.cc:87
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &verbosity)
WORKSPACE METHOD: MCSetSeedFromTime.
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:44
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
bool empty() const
Check if variable is empty.
Definition: matpackIII.cc:36
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:58
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
bool empty() const
Check if variable is empty.
Definition: matpackIV.cc:50
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:46
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
const Tensor4 & Data() const noexcept
Energy level type.
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.
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:51
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
The MatrixView class.
Definition: matpackI.h:1125
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
Radiation Vector for Stokes dimension 1-4.
The range class.
Definition: matpackI.h:165
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
The Sparse class.
Definition: matpackII.h:67
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
The Tensor4 class.
Definition: matpackIV.h:421
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
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:1147
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
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:571
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:563
#define q1
#define abs(x)
#define q
#define ns
#define min(a, b)
#define max(a, b)
Header file for logic.cc.
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:1046
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:160
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:1596
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:613
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:1117
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:1916
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:1307
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:1346
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:1336
const Index GFIELD4_LON_GRID
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:1472
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
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:57
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
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_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_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:2248
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:1973
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:2206
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:2068
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
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.
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath.h:86
Index np
Number of points describing the ppath.
Definition: ppath.h:52
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath.h:64
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath.h:84
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
Vector lstep
The length between ppath points.
Definition: ppath.h:70
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
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