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