ARTS  2.4.0(git:4fb77825)
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 "rte.h"
33 #include <cmath>
34 #include <stdexcept>
35 #include "auto_md.h"
36 #include "check_input.h"
37 #include "legacy_continua.h"
38 #include "geodetic.h"
39 #include "lin_alg.h"
40 #include "logic.h"
41 #include "math_funcs.h"
42 #include "montecarlo.h"
43 #include "physics_funcs.h"
44 #include "ppath.h"
45 #include "refraction.h"
46 #include "special_interp.h"
47 
48 extern const String SURFACE_MAINTAG;
49 extern const String PROPMAT_SUBSUBTAG;
50 extern const String TEMPERATURE_MAINTAG;
51 extern const String SCATSPECIES_MAINTAG;
52 extern const Numeric SPEED_OF_LIGHT;
53 extern const Numeric TEMP_0_C;
54 
55 /*===========================================================================
56  === The functions in alphabetical order
57  ===========================================================================*/
58 
61  ArrayOfStokesVector& dS_dx,
63  ConstVectorView ppath_f_grid,
64  ConstVectorView ppath_line_of_sight,
65  ConstVectorView ppath_vmrs,
66  const Numeric& ppath_temperature,
67  const Numeric& ppath_pressure,
68  const ArrayOfIndex& jacobian_species,
69  const ArrayOfIndex& jacobian_wind,
70  const Index& lte,
71  const Index& atmosphere_dim,
72  const bool& jacobian_do) {
73  if (not jacobian_do) return;
74 
75  // All relevant quantities are extracted first
76  const Index nq = jacobian_quantities.nelem();
77 
78  // Computational temporary vector
79  Vector a;
80 
81  for (Index i = 0; i < nq; i++) {
82  if (not jacobian_quantities[i].Analytical()) continue;
83 
84  Index component;
85 
86  switch (jacobian_wind[i]) {
87  case Index(JacobianType::AbsWind):
88  component = 0;
89  break;
90  case Index(JacobianType::WindFieldU):
91  component = 1;
92  break;
93  case Index(JacobianType::WindFieldV):
94  component = 2;
95  break;
96  case Index(JacobianType::WindFieldW):
97  component = 3;
98  break;
99  default:
100  component = -1; // This could be frequency derivative...
101  }
102  if (component not_eq -1) {
104  a, component, ppath_line_of_sight, ppath_f_grid, atmosphere_dim);
105 
106  // Apply conversion to K-matrix partial derivative
107  dK_dx[i] *= a;
108 
109  // Apply conversion to source vector partial derivative
110  if (not lte) dS_dx[i] *= a;
111  } else if (jacobian_species[i] > -1) {
112  const bool from_propmat =
113  jacobian_quantities[i].SubSubtag() == PROPMAT_SUBSUBTAG;
114  const Index& isp = jacobian_species[i];
115 
116  // Computational factor
117  Numeric factor;
118 
119  // Scaling factors to handle retrieval unit
120  if (not from_propmat) {
121  //vmrunitscf(factor, jacobian_quantities[i].Mode(),
122  vmrunitscf(
123  factor, "vmr", ppath_vmrs[isp], ppath_pressure, ppath_temperature);
124  } else {
125  //dxdvmrscf(factor, jacobian_quantities[i].Mode(),
126  dxdvmrscf(
127  factor, "vmr", ppath_vmrs[isp], ppath_pressure, ppath_temperature);
128  }
129  // Apply conversion to K-matrix partial derivative
130  dK_dx[i] *= factor;
131 
132  // Apply conversion to source vector partial derivative
133  if (not lte) dS_dx[i] *= factor;
134  } else {
135  // All other partial derivatives should already be adapted...
136  }
137  }
138 }
139 
141  if (atmosphere_dim == 1) {
142  if (los[0] < 0) {
143  los[0] = -los[0];
144  } else if (los[0] > 180) {
145  los[0] = 360 - los[0];
146  }
147  } else if (atmosphere_dim == 2) {
148  if (los[0] < -180) {
149  los[0] = los[0] + 360;
150  } else if (los[0] > 180) {
151  los[0] = los[0] - 360;
152  }
153  } else {
154  // If any of the angles out-of-bounds, use cart2zaaa to resolve
155  if (abs(los[0] - 90) > 90 || abs(los[1]) > 180) {
156  Numeric dx, dy, dz;
157  zaaa2cart(dx, dy, dz, los[0], los[1]);
158  cart2zaaa(los[0], los[1], dx, dy, dz);
159  }
160  }
161 }
162 
164  const String& iy_unit,
166  const Numeric& n,
167  const ArrayOfIndex& i_pol) {
168  // The code is largely identical between the two apply_iy_unit functions.
169  // If any change here, remember to update the other function.
170 
171  const Index nf = iy.nrows();
172  const Index ns = iy.ncols();
173 
174  assert(f_grid.nelem() == nf);
175  assert(i_pol.nelem() == ns);
176 
177  if (iy_unit == "1") {
178  if (n != 1) {
179  iy *= (n * n);
180  }
181  }
182 
183  else if (iy_unit == "RJBT") {
184  for (Index iv = 0; iv < nf; iv++) {
185  const Numeric scfac = invrayjean(1, f_grid[iv]);
186  for (Index is = 0; is < ns; is++) {
187  if (i_pol[is] < 5) // Stokes components
188  {
189  iy(iv, is) *= scfac;
190  } else // Measuement single pols
191  {
192  iy(iv, is) *= 2 * scfac;
193  }
194  }
195  }
196  }
197 
198  else if (iy_unit == "PlanckBT") {
199  for (Index iv = 0; iv < nf; iv++) {
200  for (Index is = ns - 1; is >= 0; is--) // Order must here be reversed
201  {
202  if (i_pol[is] == 1) {
203  iy(iv, is) = invplanck(iy(iv, is), f_grid[iv]);
204  } else if (i_pol[is] < 5) {
205  assert(i_pol[0] == 1);
206  iy(iv, is) = invplanck(0.5 * (iy(iv, 0) + iy(iv, is)), f_grid[iv]) -
207  invplanck(0.5 * (iy(iv, 0) - iy(iv, is)), f_grid[iv]);
208  } else {
209  iy(iv, is) = invplanck(2 * iy(iv, is), f_grid[iv]);
210  }
211  }
212  }
213  }
214 
215  else if (iy_unit == "W/(m^2 m sr)") {
216  for (Index iv = 0; iv < nf; iv++) {
217  const Numeric scfac = n * n * f_grid[iv] * (f_grid[iv] / SPEED_OF_LIGHT);
218  for (Index is = 0; is < ns; is++) {
219  iy(iv, is) *= scfac;
220  }
221  }
222  }
223 
224  else if (iy_unit == "W/(m^2 m-1 sr)") {
225  iy *= (n * n * SPEED_OF_LIGHT);
226  }
227 
228  else {
229  ostringstream os;
230  os << "Unknown option: iy_unit = \"" << iy_unit << "\"\n"
231  << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
232  << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"";
233 
234  throw runtime_error(os.str());
235  }
236 }
237 
240  const String& iy_unit,
242  const Numeric& n,
243  const ArrayOfIndex& i_pol) {
244  // The code is largely identical between the two apply_iy_unit functions.
245  // If any change here, remember to update the other function.
246 
247  const Index nf = iy.nrows();
248  const Index ns = iy.ncols();
249  const Index np = J.npages();
250 
251  assert(J.nrows() == nf);
252  assert(J.ncols() == ns);
253  assert(f_grid.nelem() == nf);
254  assert(i_pol.nelem() == ns);
255 
256  if (iy_unit == "1") {
257  if (n != 1) {
258  J *= (n * n);
259  }
260  }
261 
262  else if (iy_unit == "RJBT") {
263  for (Index iv = 0; iv < nf; iv++) {
264  const Numeric scfac = invrayjean(1, f_grid[iv]);
265  for (Index is = 0; is < ns; is++) {
266  if (i_pol[is] < 5) // Stokes componenets
267  {
268  for (Index ip = 0; ip < np; ip++) {
269  J(ip, iv, is) *= scfac;
270  }
271  } else // Measuement single pols
272  {
273  for (Index ip = 0; ip < np; ip++) {
274  J(ip, iv, is) *= 2 * scfac;
275  }
276  }
277  }
278  }
279  }
280 
281  else if (iy_unit == "PlanckBT") {
282  for (Index iv = 0; iv < f_grid.nelem(); iv++) {
283  for (Index is = ns - 1; is >= 0; is--) {
284  Numeric scfac = 1;
285  if (i_pol[is] == 1) {
286  scfac = dinvplanckdI(iy(iv, is), f_grid[iv]);
287  } else if (i_pol[is] < 5) {
288  assert(i_pol[0] == 1);
289  scfac = dinvplanckdI(0.5 * (iy(iv, 0) + iy(iv, is)), f_grid[iv]) +
290  dinvplanckdI(0.5 * (iy(iv, 0) - iy(iv, is)), f_grid[iv]);
291  } else {
292  scfac = dinvplanckdI(2 * iy(iv, is), f_grid[iv]);
293  }
294  //
295  for (Index ip = 0; ip < np; ip++) {
296  J(ip, iv, is) *= scfac;
297  }
298  }
299  }
300  }
301 
302  else if (iy_unit == "W/(m^2 m sr)") {
303  for (Index iv = 0; iv < nf; iv++) {
304  const Numeric scfac = n * n * f_grid[iv] * (f_grid[iv] / SPEED_OF_LIGHT);
305  for (Index ip = 0; ip < np; ip++) {
306  for (Index is = 0; is < ns; is++) {
307  J(ip, iv, is) *= scfac;
308  }
309  }
310  }
311  }
312 
313  else if (iy_unit == "W/(m^2 m-1 sr)") {
314  J *= (n * n * SPEED_OF_LIGHT);
315  }
316 
317  else {
318  ostringstream os;
319  os << "Unknown option: iy_unit = \"" << iy_unit << "\"\n"
320  << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
321  << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\"";
322 
323  throw runtime_error(os.str());
324  }
325 }
326 
327 void bending_angle1d(Numeric& alpha, const Ppath& ppath) {
328  Numeric theta;
329  if (ppath.dim < 3) {
330  theta = abs(ppath.start_pos[1] - ppath.end_pos[1]);
331  } else {
332  theta = sphdist(ppath.start_pos[1],
333  ppath.start_pos[2],
334  ppath.end_pos[1],
335  ppath.end_pos[2]);
336  }
337 
338  // Eq 17 in Kursinski et al., TAO, 2000:
339  alpha = ppath.start_los[0] - ppath.end_los[0] + theta;
340 
341  // This as
342  // phi_r = 180 - ppath.end_los[0]
343  // phi_t = ppath.start_los[0]
344 }
345 
373  Vector& pos,
374  Vector& rte_los,
375  Index& background,
377  const Numeric& lo0,
378  const Agenda& ppath_step_agenda,
379  const Numeric& ppath_lmax,
380  const Numeric& ppath_lraytrace,
381  const Index& atmosphere_dim,
389  const Verbosity& verbosity) {
390  // Special treatment of 1D around zenith/nadir
391  // (zenith angles outside [0,180] are changed by *adjust_los*)
392  bool invert_lat = false;
393  if (atmosphere_dim == 1 && (rte_los[0] < 0 || rte_los[0] > 180)) {
394  invert_lat = true;
395  }
396 
397  // Handle cases where angles have moved out-of-bounds due to disturbance
399 
400  // Calculate the ppath for disturbed rte_los
401  Ppath ppx;
402  //
403  ppath_calc(ws,
404  ppx,
407  p_grid,
408  lat_grid,
409  lon_grid,
410  z_field,
411  f_grid,
412  refellipsoid,
413  z_surface,
414  0,
415  ArrayOfIndex(0),
416  rte_pos,
417  rte_los,
418  ppath_lmax,
420  0,
421  verbosity);
422  //
423  background = ppath_what_background(ppx);
424 
425  // Calcualte cumulative optical path for ppx
426  Vector lox(ppx.np);
427  Index ilast = ppx.np - 1;
428  lox[0] = ppx.end_lstep;
429  for (Index i = 1; i <= ilast; i++) {
430  lox[i] =
431  lox[i - 1] + ppx.lstep[i - 1] * (ppx.nreal[i - 1] + ppx.nreal[i]) / 2.0;
432  }
433 
434  pos.resize(max(Index(2), atmosphere_dim));
435 
436  // Reciever at a longer distance (most likely out in space):
437  if (lox[ilast] < lo0) {
438  const Numeric dl = lo0 - lox[ilast];
439  if (atmosphere_dim < 3) {
440  Numeric x, z, dx, dz;
441  poslos2cart(
442  x, z, dx, dz, ppx.r[ilast], ppx.pos(ilast, 1), ppx.los(ilast, 0));
443  cart2pol(pos[0],
444  pos[1],
445  x + dl * dx,
446  z + dl * dz,
447  ppx.pos(ilast, 1),
448  ppx.los(ilast, 0));
449  } else {
450  Numeric x, y, z, dx, dy, dz;
451  poslos2cart(x,
452  y,
453  z,
454  dx,
455  dy,
456  dz,
457  ppx.r[ilast],
458  ppx.pos(ilast, 1),
459  ppx.pos(ilast, 2),
460  ppx.los(ilast, 0),
461  ppx.los(ilast, 1));
462  cart2sph(pos[0],
463  pos[1],
464  pos[2],
465  x + dl * dx,
466  y + dl * dy,
467  z + dl * dz,
468  ppx.pos(ilast, 1),
469  ppx.pos(ilast, 2),
470  ppx.los(ilast, 0),
471  ppx.los(ilast, 1));
472  }
473  }
474 
475  // Interpolate to lo0
476  else {
477  GridPos gp;
478  Vector itw(2);
479  gridpos(gp, lox, lo0);
480  interpweights(itw, gp);
481  //
482  pos[0] = interp(itw, ppx.r, gp);
483  pos[1] = interp(itw, ppx.pos(joker, 1), gp);
484  if (atmosphere_dim == 3) {
485  pos[2] = interp(itw, ppx.pos(joker, 2), gp);
486  }
487  }
488 
489  if (invert_lat) {
490  pos[1] = -pos[1];
491  }
492 }
493 
495  Numeric& dlf,
496  const Agenda& ppath_step_agenda,
497  const Index& atmosphere_dim,
505  const Ppath& ppath,
506  const Numeric& ppath_lmax,
507  const Numeric& ppath_lraytrace,
508  const Numeric& dza,
509  const Verbosity& verbosity) {
510  // Optical and physical path between transmitter and reciver
511  Numeric lo = ppath.start_lstep + ppath.end_lstep;
512  Numeric lp = lo;
513  for (Index i = 0; i <= ppath.np - 2; i++) {
514  lp += ppath.lstep[i];
515  lo += ppath.lstep[i] * (ppath.nreal[i] + ppath.nreal[i + 1]) / 2.0;
516  }
517  // Extract rte_pos and rte_los
518  const Vector rte_pos = ppath.start_pos[Range(0, atmosphere_dim)];
519  //
520  Vector rte_los0(max(Index(1), atmosphere_dim - 1)), rte_los;
521  mirror_los(rte_los, ppath.start_los, atmosphere_dim);
522  rte_los0 = rte_los[Range(0, max(Index(1), atmosphere_dim - 1))];
523 
524  // A new ppath with positive zenith angle off-set
525  //
526  Vector pos1;
527  Index backg1;
528  //
529  rte_los = rte_los0;
530  rte_los[0] += dza;
531  //
533  pos1,
534  rte_los,
535  backg1,
536  rte_pos,
537  lo,
539  ppath_lmax,
542  p_grid,
543  lat_grid,
544  lon_grid,
545  z_field,
546  f_grid,
547  refellipsoid,
548  z_surface,
549  verbosity);
550 
551  // Same thing with negative zenit angle off-set
552  Vector pos2;
553  Index backg2;
554  //
555  rte_los = rte_los0; // Use rte_los0 as rte_los can have been "adjusted"
556  rte_los[0] -= dza;
557  //
559  pos2,
560  rte_los,
561  backg2,
562  rte_pos,
563  lo,
565  ppath_lmax,
568  p_grid,
569  lat_grid,
570  lon_grid,
571  z_field,
572  f_grid,
573  refellipsoid,
574  z_surface,
575  verbosity);
576 
577  // Calculate distance between pos1 and 2, and derive the loss factor
578  // All appears OK:
579  if (backg1 == backg2) {
580  Numeric l12;
581  if (atmosphere_dim < 3) {
582  distance2D(l12, pos1[0], pos1[1], pos2[0], pos2[1]);
583  } else {
584  distance3D(l12, pos1[0], pos1[1], pos1[2], pos2[0], pos2[1], pos2[2]);
585  }
586  //
587  dlf = lp * 2 * DEG2RAD * dza / l12;
588  }
589  // If different backgrounds, then only use the second calculation
590  else {
591  Numeric l12;
592  if (atmosphere_dim == 1) {
593  const Numeric r = refellipsoid[0];
594  distance2D(l12, r + ppath.end_pos[0], 0, pos2[0], pos2[1]);
595  } else if (atmosphere_dim == 2) {
596  const Numeric r = refell2r(refellipsoid, ppath.end_pos[1]);
597  distance2D(l12, r + ppath.end_pos[0], ppath.end_pos[1], pos2[0], pos2[1]);
598  } else {
599  const Numeric r = refell2r(refellipsoid, ppath.end_pos[1]);
600  distance3D(l12,
601  r + ppath.end_pos[0],
602  ppath.end_pos[1],
603  ppath.end_pos[2],
604  pos2[0],
605  pos2[1],
606  pos2[2]);
607  }
608  //
609  dlf = lp * DEG2RAD * dza / l12;
610  }
611 }
612 
614  Numeric& dlf,
615  const Agenda& ppath_step_agenda,
616  const Index& atmosphere_dim,
624  const Ppath& ppath,
625  const Numeric& ppath_lmax,
626  const Numeric& ppath_lraytrace,
627  const Numeric& dza,
628  const Verbosity& verbosity) {
629  if (ppath.end_los[0] < 90 || ppath.start_los[0] > 90)
630  throw runtime_error(
631  "The function *defocusing_sat2sat* can only be used "
632  "for limb sounding geometry.");
633 
634  // Index of tangent point
635  Index it;
636  find_tanpoint(it, ppath);
637  assert(it >= 0);
638 
639  // Length between tangent point and transmitter/reciver
640  Numeric lt = ppath.start_lstep, lr = ppath.end_lstep;
641  for (Index i = it; i <= ppath.np - 2; i++) {
642  lt += ppath.lstep[i];
643  }
644  for (Index i = 0; i < it; i++) {
645  lr += ppath.lstep[i];
646  }
647 
648  // Bending angle and impact parameter for centre ray
649  Numeric alpha0, a0;
650  bending_angle1d(alpha0, ppath);
651  alpha0 *= DEG2RAD;
652  a0 = ppath.constant;
653 
654  // Azimuth loss term (Eq 18.5 in Kursinski et al.)
655  const Numeric lf = lr * lt / (lr + lt);
656  const Numeric alt = 1 / (1 - alpha0 * lf / refellipsoid[0]);
657 
658  // Calculate two new ppaths to get dalpha/da
659  Numeric alpha1, a1, alpha2, a2, dada;
660  Ppath ppt;
661  Vector rte_pos = ppath.end_pos[Range(0, atmosphere_dim)];
662  Vector rte_los = ppath.end_los;
663  //
664  rte_los[0] -= dza;
666  ppath_calc(ws,
667  ppt,
670  p_grid,
671  lat_grid,
672  lon_grid,
673  z_field,
674  f_grid,
675  refellipsoid,
676  z_surface,
677  0,
678  ArrayOfIndex(0),
679  rte_pos,
680  rte_los,
681  ppath_lmax,
683  0,
684  verbosity);
685  bending_angle1d(alpha2, ppt);
686  alpha2 *= DEG2RAD;
687  a2 = ppt.constant;
688  //
689  rte_los[0] += 2 * dza;
691  ppath_calc(ws,
692  ppt,
695  p_grid,
696  lat_grid,
697  lon_grid,
698  z_field,
699  f_grid,
700  refellipsoid,
701  z_surface,
702  0,
703  ArrayOfIndex(0),
704  rte_pos,
705  rte_los,
706  ppath_lmax,
708  0,
709  verbosity);
710  // This path can hit the surface. And we need to check if ppt is OK.
711  // (remember this function only deals with sat-to-sat links and OK
712  // background here is be space)
713  // Otherwise use the centre ray as the second one.
714  if (ppath_what_background(ppt) == 1) {
715  bending_angle1d(alpha1, ppt);
716  alpha1 *= DEG2RAD;
717  a1 = ppt.constant;
718  dada = (alpha2 - alpha1) / (a2 - a1);
719  } else {
720  dada = (alpha2 - alpha0) / (a2 - a0);
721  }
722 
723  // Zenith loss term (Eq 18 in Kursinski et al.)
724  const Numeric zlt = 1 / (1 - dada * lf);
725 
726  // Total defocusing loss
727  dlf = zlt * alt;
728 }
729 
731  const Numeric& u,
732  const Numeric& v,
733  const Numeric& w,
734  const Index& atmosphere_dim) {
735  // Strength of field
736  const Numeric f = sqrt(u * u + v * v + w * w);
737 
738  // Zenith and azimuth angle for field (in radians)
739  const Numeric za_f = acos(w / f);
740  const Numeric aa_f = atan2(u, v);
741 
742  // Zenith and azimuth angle for photon direction (in radians)
743  Vector los_p;
744  mirror_los(los_p, los, atmosphere_dim);
745  const Numeric za_p = DEG2RAD * los_p[0];
746  const Numeric aa_p = DEG2RAD * los_p[1];
747 
748  return f * (cos(za_f) * cos(za_p) + sin(za_f) * sin(za_p) * cos(aa_f - aa_p));
749 }
750 
751 void ext_mat_case(Index& icase,
753  const Index stokes_dim) {
754  if (icase == 0) {
755  icase = 1; // Start guess is diagonal
756 
757  //--- Scalar case ----------------------------------------------------------
758  if (stokes_dim == 1) {
759  }
760 
761  //--- Vector RT ------------------------------------------------------------
762  else {
763  // Check symmetries and analyse structure of exp_mat:
764  assert(ext_mat(1, 1) == ext_mat(0, 0));
765  assert(ext_mat(1, 0) == ext_mat(0, 1));
766 
767  if (ext_mat(1, 0) != 0) {
768  icase = 2;
769  }
770 
771  if (stokes_dim >= 3) {
772  assert(ext_mat(2, 2) == ext_mat(0, 0));
773  assert(ext_mat(2, 1) == -ext_mat(1, 2));
774  assert(ext_mat(2, 0) == ext_mat(0, 2));
775 
776  if (ext_mat(2, 0) != 0 || ext_mat(2, 1) != 0) {
777  icase = 3;
778  }
779 
780  if (stokes_dim > 3) {
781  assert(ext_mat(3, 3) == ext_mat(0, 0));
782  assert(ext_mat(3, 2) == -ext_mat(2, 3));
783  assert(ext_mat(3, 1) == -ext_mat(1, 3));
784  assert(ext_mat(3, 0) == ext_mat(0, 3));
785 
786  if (icase < 3) // if icase==3, already at most complex case
787  {
788  if (ext_mat(3, 0) != 0 || ext_mat(3, 1) != 0) {
789  icase = 3;
790  } else if (ext_mat(3, 2) != 0) {
791  icase = 2;
792  }
793  }
794  }
795  }
796  }
797  }
798 }
799 
800 void ext2trans(MatrixView trans_mat,
801  Index& icase,
803  const Numeric& lstep) {
804  const Index stokes_dim = ext_mat.ncols();
805 
806  assert(ext_mat.nrows() == stokes_dim);
807  assert(trans_mat.nrows() == stokes_dim && trans_mat.ncols() == stokes_dim);
808 
809  // Theoretically ext_mat(0,0) >= 0, but to demand this can cause problems for
810  // iterative retrievals, and the assert is skipped. Negative should be a
811  // result of negative vmr, and this issue is checked in basics_checkedCalc.
812  //assert( ext_mat(0,0) >= 0 );
813 
814  assert(icase >= 0 && icase <= 3);
815  assert(!is_singular(ext_mat));
816  assert(lstep >= 0);
817 
818  // Analyse ext_mat?
820 
821  // Calculation options:
822  if (icase == 1) {
823  trans_mat = 0;
824  trans_mat(0, 0) = exp(-ext_mat(0, 0) * lstep);
825  for (Index i = 1; i < stokes_dim; i++) {
826  trans_mat(i, i) = trans_mat(0, 0);
827  }
828  }
829 
830  else if (icase == 2 && stokes_dim < 3) {
831  // Expressions below are found in "Polarization in Spectral Lines" by
832  // Landi Degl'Innocenti and Landolfi (2004).
833  const Numeric tI = exp(-ext_mat(0, 0) * lstep);
834  const Numeric HQ = ext_mat(0, 1) * lstep;
835  trans_mat(0, 0) = tI * cosh(HQ);
836  trans_mat(1, 1) = trans_mat(0, 0);
837  trans_mat(1, 0) = -tI * sinh(HQ);
838  trans_mat(0, 1) = trans_mat(1, 0);
839  /* Does not work for stokes_dim==3, and commnted out 180502 by PE:
840  if( stokes_dim >= 3 )
841  {
842  trans_mat(2,0) = 0;
843  trans_mat(2,1) = 0;
844  trans_mat(0,2) = 0;
845  trans_mat(1,2) = 0;
846  const Numeric RQ = ext_mat(2,3) * lstep;
847  trans_mat(2,2) = tI * cos( RQ );
848  if( stokes_dim > 3 )
849  {
850  trans_mat(3,0) = 0;
851  trans_mat(3,1) = 0;
852  trans_mat(0,3) = 0;
853  trans_mat(1,3) = 0;
854  trans_mat(3,3) = trans_mat(2,2);
855  trans_mat(3,2) = tI * sin( RQ );
856  trans_mat(2,3) = -trans_mat(3,2);
857  }
858  }
859  */
860  } else {
861  Matrix ext_mat_ds = ext_mat;
862  ext_mat_ds *= -lstep;
863  //
864  Index q = 10; // index for the precision of the matrix exp function
865  //
866  switch (stokes_dim) {
867  case 4:
869  trans_mat, ext_mat_ds);
870  break;
871  default:
872  matrix_exp(trans_mat, ext_mat_ds, q);
873  }
874  }
875 }
876 
877 void get_iy(Workspace& ws,
878  Matrix& iy,
879  const Index& cloudbox_on,
881  const EnergyLevelMap& nlte_field,
885  const String& iy_unit,
886  const Agenda& iy_main_agenda) {
889  Ppath ppath;
890  Tensor3 iy_transmission(0, 0, 0);
891  const Index iy_agenda_call1 = 1;
892  const ArrayOfString iy_aux_vars(0);
893  const Index iy_id = 0;
894  const Index jacobian_do = 0;
895 
897  iy,
898  iy_aux,
899  ppath,
900  diy_dx,
903  iy_aux_vars,
904  iy_id,
905  iy_unit,
906  cloudbox_on,
907  jacobian_do,
908  f_grid,
909  nlte_field,
910  rte_pos,
911  rte_los,
912  rte_pos2,
914 }
915 
917  Matrix& iy,
920  const Index& iy_id,
921  const Index& jacobian_do,
923  const Ppath& ppath,
925  const Index& atmosphere_dim,
926  const EnergyLevelMap& nlte_field,
927  const Index& cloudbox_on,
928  const Index& stokes_dim,
930  const String& iy_unit,
932  const Agenda& iy_main_agenda,
933  const Agenda& iy_space_agenda,
934  const Agenda& iy_surface_agenda,
935  const Agenda& iy_cloudbox_agenda,
936  const Index& iy_agenda_call1,
937  const Verbosity& verbosity) {
938  CREATE_OUT3;
939 
940  // Some sizes
941  const Index nf = f_grid.nelem();
942  const Index np = ppath.np;
943 
944  // Set rtp_pos and rtp_los to match the last point in ppath.
945  //
946  // Note that the Ppath positions (ppath.pos) for 1D have one column more
947  // than expected by most functions. Only the first atmosphere_dim values
948  // shall be copied.
949  //
951  rtp_pos.resize(atmosphere_dim);
952  rtp_pos = ppath.pos(np - 1, Range(0, atmosphere_dim));
953  rtp_los.resize(ppath.los.ncols());
954  rtp_los = ppath.los(np - 1, joker);
955 
956  out3 << "Radiative background: " << ppath.background << "\n";
957 
958  // Handle the different background cases
959  //
960  String agenda_name;
961  //
962  switch (ppath_what_background(ppath)) {
963  case 1: //--- Space ----------------------------------------------------
964  {
965  agenda_name = "iy_space_agenda";
966  chk_not_empty(agenda_name, iy_space_agenda);
968  } break;
969 
970  case 2: //--- The surface -----------------------------------------------
971  {
972  agenda_name = "iy_surface_agenda";
973  chk_not_empty(agenda_name, iy_surface_agenda);
974  //
975  const Index los_id = iy_id % (Index)1000;
976  Index iy_id_new = iy_id + (Index)9 * los_id;
977  //
978  // Surface jacobian stuff:
980  if (jacobian_do && iy_agenda_call1) {
981  for (Index i = 0; i < jacobian_quantities.nelem(); i++) {
982  if (jacobian_quantities[i].MainTag() == SURFACE_MAINTAG) {
983  dsurface_names.push_back(jacobian_quantities[i].Subtag());
984  }
985  }
986  }
989  //
991  iy,
992  diy_dx,
995  iy_unit,
997  iy_id_new,
998  cloudbox_on,
999  jacobian_do,
1001  f_grid,
1002  nlte_field,
1003  rtp_pos,
1004  rtp_los,
1005  rte_pos2,
1009  } break;
1010 
1011  case 3: //--- Cloudbox boundary or interior ------------------------------
1012  case 4: {
1013  agenda_name = "iy_cloudbox_agenda";
1014  chk_not_empty(agenda_name, iy_cloudbox_agenda);
1017  } break;
1018 
1019  default: //--- ????? ----------------------------------------------------
1020  // Are we here, the coding is wrong somewhere
1021  assert(false);
1022  }
1023 
1024  if (iy.ncols() != stokes_dim || iy.nrows() != nf) {
1025  ostringstream os;
1026  os << "The size of *iy* returned from *" << agenda_name << "* is\n"
1027  << "not correct:\n"
1028  << " expected size = [" << nf << "," << stokes_dim << "]\n"
1029  << " size of iy = [" << iy.nrows() << "," << iy.ncols() << "]\n";
1030  throw runtime_error(os.str());
1031  }
1032 }
1033 
1035  Vector& ppath_t,
1036  EnergyLevelMap& ppath_nlte,
1037  Matrix& ppath_vmr,
1038  Matrix& ppath_wind,
1039  Matrix& ppath_mag,
1040  const Ppath& ppath,
1041  const Index& atmosphere_dim,
1044  const EnergyLevelMap& nlte_field,
1052  const Index np = ppath.np;
1053  // Pressure:
1054  ppath_p.resize(np);
1055  Matrix itw_p(np, 2);
1056  interpweights(itw_p, ppath.gp_p);
1057  itw2p(ppath_p, p_grid, ppath.gp_p, itw_p);
1058 
1059  // Temperature:
1060  ppath_t.resize(np);
1061  Matrix itw_field;
1063  itw_field, atmosphere_dim, ppath.gp_p, ppath.gp_lat, ppath.gp_lon);
1064  interp_atmfield_by_itw(ppath_t,
1066  t_field,
1067  ppath.gp_p,
1068  ppath.gp_lat,
1069  ppath.gp_lon,
1070  itw_field);
1071 
1072  // VMR fields:
1073  const Index ns = vmr_field.nbooks();
1074  ppath_vmr.resize(ns, np);
1075  for (Index is = 0; is < ns; is++) {
1076  interp_atmfield_by_itw(ppath_vmr(is, joker),
1078  vmr_field(is, joker, joker, joker),
1079  ppath.gp_p,
1080  ppath.gp_lat,
1081  ppath.gp_lon,
1082  itw_field);
1083  }
1084 
1085  // NLTE temperatures
1086  ppath_nlte = nlte_field.InterpToGridPos(atmosphere_dim, ppath.gp_p, ppath.gp_lat, ppath.gp_lon);
1087 
1088  // Winds:
1089  ppath_wind.resize(3, np);
1090  ppath_wind = 0;
1091  //
1092  if (wind_u_field.npages() > 0) {
1093  interp_atmfield_by_itw(ppath_wind(0, joker),
1095  wind_u_field,
1096  ppath.gp_p,
1097  ppath.gp_lat,
1098  ppath.gp_lon,
1099  itw_field);
1100  }
1101  if (wind_v_field.npages() > 0) {
1102  interp_atmfield_by_itw(ppath_wind(1, joker),
1104  wind_v_field,
1105  ppath.gp_p,
1106  ppath.gp_lat,
1107  ppath.gp_lon,
1108  itw_field);
1109  }
1110  if (wind_w_field.npages() > 0) {
1111  interp_atmfield_by_itw(ppath_wind(2, joker),
1113  wind_w_field,
1114  ppath.gp_p,
1115  ppath.gp_lat,
1116  ppath.gp_lon,
1117  itw_field);
1118  }
1119 
1120  // Magnetic field:
1121  ppath_mag.resize(3, np);
1122  ppath_mag = 0;
1123  //
1124  if (mag_u_field.npages() > 0) {
1125  interp_atmfield_by_itw(ppath_mag(0, joker),
1127  mag_u_field,
1128  ppath.gp_p,
1129  ppath.gp_lat,
1130  ppath.gp_lon,
1131  itw_field);
1132  }
1133  if (mag_v_field.npages() > 0) {
1134  interp_atmfield_by_itw(ppath_mag(1, joker),
1136  mag_v_field,
1137  ppath.gp_p,
1138  ppath.gp_lat,
1139  ppath.gp_lon,
1140  itw_field);
1141  }
1142  if (mag_w_field.npages() > 0) {
1143  interp_atmfield_by_itw(ppath_mag(2, joker),
1145  mag_w_field,
1146  ppath.gp_p,
1147  ppath.gp_lat,
1148  ppath.gp_lon,
1149  itw_field);
1150  }
1151 }
1152 
1154  Matrix& ppath_pnd,
1155  ArrayOfMatrix& ppath_dpnd_dx,
1156  const Ppath& ppath,
1157  const Index& atmosphere_dim,
1159  const Tensor4& pnd_field,
1160  const ArrayOfTensor4& dpnd_field_dx) {
1161  const Index np = ppath.np;
1162 
1163  // Pnd along the ppath
1164  ppath_pnd.resize(pnd_field.nbooks(), np);
1165  ppath_pnd = 0;
1166  ppath_dpnd_dx.resize(dpnd_field_dx.nelem());
1167 
1168  bool any_dpnd = false;
1169  for (Index iq = 0; iq < dpnd_field_dx.nelem(); iq++) {
1170  if (dpnd_field_dx[iq].empty()) {
1171  ppath_dpnd_dx[iq].resize(0, 0);
1172  } else {
1173  any_dpnd = true;
1174  ppath_dpnd_dx[iq].resize(pnd_field.nbooks(), np);
1175  }
1176  }
1177 
1178  // A variable that can map from ppath to particle containers.
1179  // If outside cloudbox or all (d)pnd=0, this variable holds -1.
1180  clear2cloudy.resize(np);
1181 
1182  // Determine ppath_pnd and ppath_dpnd_dx
1183  Index nin = 0;
1184  for (Index ip = 0; ip < np; ip++) // PPath point
1185  {
1186  Matrix itw(1, Index(pow(2.0, Numeric(atmosphere_dim))));
1187 
1188  ArrayOfGridPos gpc_p(1), gpc_lat(1), gpc_lon(1);
1189  GridPos gp_lat, gp_lon;
1190  if (atmosphere_dim >= 2) {
1191  gridpos_copy(gp_lat, ppath.gp_lat[ip]);
1192  }
1193  if (atmosphere_dim == 3) {
1194  gridpos_copy(gp_lon, ppath.gp_lon[ip]);
1195  }
1196 
1197  if (is_gp_inside_cloudbox(ppath.gp_p[ip],
1198  gp_lat,
1199  gp_lon,
1201  true,
1202  atmosphere_dim)) {
1204  gpc_p[0],
1205  gpc_lat[0],
1206  gpc_lon[0],
1207  ppath.gp_p[ip],
1208  gp_lat,
1209  gp_lon,
1211  cloudbox_limits);
1212  for (Index i = 0; i < pnd_field.nbooks(); i++) {
1213  interp_atmfield_by_itw(ppath_pnd(i, ip),
1215  pnd_field(i, joker, joker, joker),
1216  gpc_p,
1217  gpc_lat,
1218  gpc_lon,
1219  itw);
1220  }
1221  bool any_ppath_dpnd = false;
1222  if (any_dpnd) {
1223  for (Index iq = 0; iq < dpnd_field_dx.nelem();
1224  iq++) // Jacobian parameter
1225  {
1226  if (!dpnd_field_dx[iq].empty()) {
1227  for (Index i = 0; i < pnd_field.nbooks();
1228  i++) // Scattering element
1229  {
1230  interp_atmfield_by_itw(ppath_dpnd_dx[iq](i, ip),
1232  dpnd_field_dx[iq](i, joker, joker, joker),
1233  gpc_p,
1234  gpc_lat,
1235  gpc_lon,
1236  itw);
1237  }
1238  if (max(ppath_dpnd_dx[iq](joker, ip)) > 0. ||
1239  min(ppath_dpnd_dx[iq](joker, ip)) < 0.)
1240  any_ppath_dpnd = true;
1241  }
1242  }
1243  }
1244  if (max(ppath_pnd(joker, ip)) > 0. || min(ppath_pnd(joker, ip)) < 0. ||
1245  any_ppath_dpnd) {
1246  clear2cloudy[ip] = nin;
1247  nin++;
1248  } else {
1249  clear2cloudy[ip] = -1;
1250  }
1251  } else {
1252  clear2cloudy[ip] = -1;
1253  }
1254  }
1255 }
1256 
1257 void get_ppath_f(Matrix& ppath_f,
1258  const Ppath& ppath,
1260  const Index& atmosphere_dim,
1261  const Numeric& rte_alonglos_v,
1262  ConstMatrixView ppath_wind) {
1263  // Sizes
1264  const Index nf = f_grid.nelem();
1265  const Index np = ppath.np;
1266 
1267  ppath_f.resize(nf, np);
1268 
1269  // Doppler relevant velocity
1270  //
1271  for (Index ip = 0; ip < np; ip++) {
1272  // Start by adding rte_alonglos_v (most likely sensor effects)
1273  Numeric v_doppler = rte_alonglos_v;
1274 
1275  // Include wind
1276  if (ppath_wind(1, ip) != 0 || ppath_wind(0, ip) != 0 ||
1277  ppath_wind(2, ip) != 0) {
1278  // The dot product below is valid for the photon direction. Winds
1279  // along this direction gives a positive contribution.
1280  v_doppler += dotprod_with_los(ppath.los(ip, joker),
1281  ppath_wind(0, ip),
1282  ppath_wind(1, ip),
1283  ppath_wind(2, ip),
1284  atmosphere_dim);
1285  }
1286 
1287  // Determine frequency grid
1288  if (v_doppler == 0) {
1289  ppath_f(joker, ip) = f_grid;
1290  } else {
1291  // Positive v_doppler means that sensor measures lower rest
1292  // frequencies
1293  const Numeric a = 1 - v_doppler / SPEED_OF_LIGHT;
1294  for (Index iv = 0; iv < nf; iv++) {
1295  ppath_f(iv, ip) = a * f_grid[iv];
1296  }
1297  }
1298  }
1299 }
1300 
1302  const Index& mblock_index) {
1303  const Index n1y = sensor_response.nrows();
1304  return Range(n1y * mblock_index, n1y);
1305 }
1306 
1308  VectorView dB_dT,
1309  ConstVectorView ppath_f_grid,
1310  const Numeric& ppath_temperature,
1311  const bool& do_temperature_derivative) {
1312  const Index nf = ppath_f_grid.nelem();
1313 
1314  for (Index i = 0; i < nf; i++)
1315  B[i] = planck(ppath_f_grid[i], ppath_temperature);
1316 
1317  if (do_temperature_derivative)
1318  for (Index i = 0; i < nf; i++)
1319  dB_dT[i] = dplanck_dt(ppath_f_grid[i], ppath_temperature);
1320 }
1321 
1323  Workspace& ws,
1324  PropagationMatrix& K,
1325  StokesVector& S,
1326  Index& lte,
1327  ArrayOfPropagationMatrix& dK_dx,
1328  ArrayOfStokesVector& dS_dx,
1331  ConstVectorView ppath_f_grid,
1332  ConstVectorView ppath_magnetic_field,
1333  ConstVectorView ppath_line_of_sight,
1334  const EnergyLevelMap& ppath_nlte,
1335  ConstVectorView ppath_vmrs,
1336  const Numeric& ppath_temperature,
1337  const Numeric& ppath_pressure,
1338  const ArrayOfIndex& jacobian_species,
1339  const bool& jacobian_do) {
1340  // All relevant quantities are extracted first
1341  const Index nq = jacobian_quantities.nelem();
1342 
1343  // Local variables inside Agenda
1345  ArrayOfStokesVector nlte_source, dnlte_dx_source, nlte_dx_dsource_dx;
1346 
1347  // Perform the propagation matrix computations
1350  nlte_source,
1353  nlte_dx_dsource_dx,
1355  ppath_f_grid,
1356  ppath_magnetic_field,
1357  ppath_line_of_sight,
1358  ppath_pressure,
1359  ppath_temperature,
1360  ppath_nlte,
1361  ppath_vmrs,
1363 
1364  // We only now know how large the propagation matrix will be!
1365  const Index npmat = propmat_clearsky.nelem();
1366 
1367  // If therea re no NLTE elements, then set the LTE flag
1368  lte = nlte_source.nelem() ? 0 : 1;
1369 
1370  // Sum the propagation matrix
1371  K = propmat_clearsky[0];
1372  for (Index i = 1; i < npmat; i++) K += propmat_clearsky[i];
1373 
1374  // Set NLTE source term if applicable
1375  if (not lte) {
1376  // Add all source terms up
1377  S = nlte_source[0];
1378  for (Index i = 1; i < npmat; i++) S += nlte_source[i];
1379  } else {
1380  S.SetZero();
1381  }
1382 
1383  // Set the partial derivatives
1384  if (jacobian_do) {
1385  for (Index i = 0; i < nq; i++) {
1386  if (jacobian_quantities[i].MainTag() == SCATSPECIES_MAINTAG) {
1387  dK_dx[i].SetZero();
1388  dS_dx[i].SetZero();
1389  } else if (jacobian_quantities[i].SubSubtag() == PROPMAT_SUBSUBTAG) {
1390  // Find position of index in ppd
1392 
1393  dK_dx[i] = dpropmat_clearsky_dx[j];
1394  if (lte) {
1395  dS_dx[i].SetZero();
1396  } else {
1397  dS_dx[i] = dnlte_dx_source[j];
1398  dS_dx[i] += nlte_dx_dsource_dx[j];
1399 
1400  // TEST: Old routine applied unit conversion on only the last
1401  // part of the equation. Was this correct or wrong? If correct,
1402  // this is an issue. Otherwise, the old version was incorrect.
1403  // Have to setup perturbation test-case to study which is most
1404  // reasonable...
1405  }
1406  } else if (jacobian_species[i] > -1) // Did not compute values in Agenda
1407  {
1408  dK_dx[i] = propmat_clearsky[jacobian_species[i]];
1409 
1410  // We cannot know the NLTE jacobian if this method was used
1411  // because that information is thrown away. It is still faster
1412  // to retain this method since it requires less computations
1413  // when we do not need NLTE, which is most of the time...
1414  if (not lte) {
1415  ostringstream os;
1416 
1417  os << "We do not yet support species"
1418  << " tag and NLTE Jacobians.\n";
1419  throw std::runtime_error(os.str());
1420  }
1421  dS_dx[i].SetZero();
1422  }
1423  }
1424  }
1425 }
1426 
1428  MatrixView J,
1429  Tensor3View dJ_dx,
1430  const PropagationMatrix& K,
1431  const StokesVector& a,
1432  const StokesVector& S,
1433  const ArrayOfPropagationMatrix& dK_dx,
1434  const ArrayOfStokesVector& da_dx,
1435  const ArrayOfStokesVector& dS_dx,
1436  ConstVectorView B,
1437  ConstVectorView dB_dT,
1439  const bool& jacobian_do) {
1440  const Index nq = jacobian_quantities.nelem();
1441  const Index nf = K.NumberOfFrequencies();
1442  const Index ns = K.StokesDimensions();
1443 
1444  dJ_dx = 0;
1445  for (Index i1 = 0; i1 < nf; i1++) {
1446  Matrix invK(ns, ns);
1447  Vector j(ns);
1448 
1449  if (K.IsRotational(i1)) {
1450  dJ_dx(joker, i1, joker) = 0;
1451  J(i1, joker) = 0;
1452  continue;
1453  }
1454 
1455  K.MatrixInverseAtPosition(invK, i1);
1456 
1457  // Set a B to j
1458  j = a.VectorAtPosition(i1);
1459  ;
1460  j *= B[i1];
1461 
1462  // Add S to j
1463  if (not S.IsEmpty()) j += S.VectorAtPosition(i1);
1464  ;
1465 
1466  // Compute J = K^-1 (a B + S)
1467  mult(J(i1, joker), invK, j);
1468 
1469  // Compute dJ = K^-1((da B + a dB + dS) - dK K^-1(a B + S))
1470  if (jacobian_do)
1471  //FOR_ANALYTICAL_JACOBIANS_DO
1472  //(
1473  for (Index iq = 0; iq < nq; iq++) {
1474  if (jacobian_quantities[iq].Analytical()) {
1475  Matrix dk(ns, ns), tmp_matrix(ns, ns);
1476  Vector dj(ns, 0), tmp(ns);
1477 
1478  // Control parameters for special jacobians that are computed elsewhere
1479  //const bool has_dk = (not dK_dx[iq].IsEmpty()); // currently always
1480  //const bool has_ds = (not dS_dx[iq].IsEmpty()); // evaluate as true
1481  const bool has_dt =
1482  (jacobian_quantities[iq].MainTag() == TEMPERATURE_MAINTAG);
1483 
1484  // Sets the -K^-1 dK/dx K^-1 (a B + S) term
1485  //if(has_dk)
1486  //{
1487  dK_dx[iq].MatrixAtPosition(dk, i1);
1488  mult(tmp, dk, J(i1, joker));
1489 
1490  dj = da_dx[iq].VectorAtPosition(i1);
1491  dj *= B[i1];
1492 
1493  dj -= tmp;
1494 
1495  // Adds a dB to dj
1496  if (has_dt) {
1497  tmp = a.VectorAtPosition(i1);
1498  tmp *= dB_dT[i1];
1499  dj += tmp;
1500  }
1501 
1502  // Adds dS to dj
1503  //if(has_ds)
1504  dj += dS_dx[iq].VectorAtPosition(i1);
1505 
1506  mult(dJ_dx(iq, i1, joker), invK, dj);
1507  //}
1508  }
1509  // don't need that anymore now that we zero dJ_dx at the very beginning
1510  //else
1511  //{
1512  // dJ_dx(iq, i1, joker) = 0;
1513  //}
1514  }
1515  //)
1516  }
1517 
1518  if (nq) dJ_dx *= 0.5;
1519 }
1520 
1523  ConstVectorView ppath_wind,
1524  ConstVectorView ppath_line_of_sight,
1525  const Numeric& rte_alonglos_v,
1526  const Index& atmosphere_dim) {
1527  Numeric v_doppler = rte_alonglos_v;
1528 
1529  if (ppath_wind[0] not_eq 0 or ppath_wind[1] not_eq 0 or
1530  ppath_wind[2] not_eq 0)
1531  v_doppler += dotprod_with_los(ppath_line_of_sight,
1532  ppath_wind[0],
1533  ppath_wind[1],
1534  ppath_wind[2],
1535  atmosphere_dim);
1536  ppath_f_grid = f_grid;
1537 
1538  if (v_doppler not_eq 0) ppath_f_grid *= 1 - v_doppler / SPEED_OF_LIGHT;
1539 }
1540 
1542  const Index& component,
1543  ConstVectorView& line_of_sight,
1545  const Index& atmosphere_dim) {
1546  // component 0 means total speed
1547  // component 1 means u speed
1548  // component 2 means v speed
1549  // component 3 means w speed
1550 
1551  // Sizes
1552  const Index nf = f_grid.nelem();
1553 
1554  // Doppler relevant velocity
1555  //
1556  // initialize
1557  Numeric dv_doppler_dx = 0.0;
1558 
1559  switch (component) {
1560  case 0: // this is total and is already initialized to avoid compiler warnings
1561  dv_doppler_dx = 1.0;
1562  break;
1563  case 1: // this is the u-component
1564  dv_doppler_dx =
1565  (dotprod_with_los(line_of_sight, 1, 0, 0, atmosphere_dim));
1566  break;
1567  case 2: // this is v-component
1568  dv_doppler_dx =
1569  (dotprod_with_los(line_of_sight, 0, 1, 0, atmosphere_dim));
1570  break;
1571  case 3: // this is w-component
1572  dv_doppler_dx =
1573  (dotprod_with_los(line_of_sight, 0, 0, 1, atmosphere_dim));
1574  break;
1575  default:
1576  throw std::runtime_error(
1577  "This being seen means that there is a development bug in interactions with get_ppath_df_dW.\n");
1578  break;
1579  }
1580 
1581  // Determine frequency grid
1582  if (dv_doppler_dx == 0.0) {
1583  f_partials.resize(nf);
1584  f_partials = 0.0;
1585  } else {
1586  f_partials = f_grid;
1587  f_partials *= -dv_doppler_dx / SPEED_OF_LIGHT;
1588  }
1589 }
1590 
1592  StokesVector& ap,
1593  PropagationMatrix& Kp,
1594  ArrayOfStokesVector& dap_dx,
1595  ArrayOfPropagationMatrix& dKp_dx,
1597  ConstMatrixView ppath_1p_pnd, // the ppath_pnd at this ppath point
1598  const ArrayOfMatrix&
1599  ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1600  const Index ppath_1p_id,
1602  ConstVectorView ppath_line_of_sight,
1603  ConstVectorView ppath_temperature,
1604  const Index& atmosphere_dim,
1605  const bool& jacobian_do) {
1606  const Index nf = Kp.NumberOfFrequencies(), stokes_dim = Kp.StokesDimensions();
1607 
1608  //StokesVector da_aux(nf, stokes_dim);
1609  //PropagationMatrix dK_aux(nf, stokes_dim);
1610 
1612 
1613  // Direction of outgoing scattered radiation (which is reversed to
1614  // LOS). Only used for extracting scattering properties.
1615  Vector dir;
1616  mirror_los(dir, ppath_line_of_sight, atmosphere_dim);
1617  Matrix dir_array(1, 2, 0.);
1618  dir_array(0, joker) = dir;
1619 
1620  ArrayOfArrayOfTensor5 ext_mat_Nse;
1621  ArrayOfArrayOfTensor4 abs_vec_Nse;
1622  ArrayOfArrayOfIndex ptypes_Nse;
1623  Matrix t_ok;
1624  ArrayOfTensor5 ext_mat_ssbulk;
1625  ArrayOfTensor4 abs_vec_ssbulk;
1626  ArrayOfIndex ptype_ssbulk;
1627  Tensor5 ext_mat_bulk;
1628  Tensor4 abs_vec_bulk;
1629  Index ptype_bulk;
1630 
1631  // get per-scat-elem data here. and fold with pnd.
1632  // keep per-scat-elem data to fold with dpnd_dx further down in
1633  // analyt-jac-loop.
1634  opt_prop_NScatElems(ext_mat_Nse,
1635  abs_vec_Nse,
1636  ptypes_Nse,
1637  t_ok,
1638  scat_data,
1639  stokes_dim,
1640  ppath_temperature,
1641  dir_array,
1642  -1);
1643 
1644  opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1645  abs_vec_ssbulk,
1646  ptype_ssbulk,
1647  ext_mat_Nse,
1648  abs_vec_Nse,
1649  ptypes_Nse,
1650  ppath_1p_pnd,
1651  t_ok);
1652  opt_prop_Bulk(ext_mat_bulk,
1653  abs_vec_bulk,
1654  ptype_bulk,
1655  ext_mat_ssbulk,
1656  abs_vec_ssbulk,
1657  ptype_ssbulk);
1658 
1659  const Index nf_ssd = abs_vec_bulk.nbooks(); // number of freqs in extracted
1660  // optprops. if 1, we need to
1661  // duplicate the ext/abs output.
1662 
1663  for (Index iv = 0; iv < nf; iv++) {
1664  if (nf_ssd > 1) {
1665  ap.SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1666  Kp.SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker), iv);
1667  } else {
1668  ap.SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1669  Kp.SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1670  }
1671  }
1672 
1673  if (jacobian_do)
1675  if (ppath_dpnd_dx[iq].empty()) {
1676  dap_dx[iq].SetZero();
1677  dKp_dx[iq].SetZero();
1678  } else {
1679  // check, whether we have any non-zero ppath_dpnd_dx in this
1680  // pnd-affecting x? might speed up things a little bit.
1681  opt_prop_ScatSpecBulk(ext_mat_ssbulk,
1682  abs_vec_ssbulk,
1683  ptype_ssbulk,
1684  ext_mat_Nse,
1685  abs_vec_Nse,
1686  ptypes_Nse,
1687  ppath_dpnd_dx[iq](joker, Range(ppath_1p_id, 1)),
1688  t_ok);
1689  opt_prop_Bulk(ext_mat_bulk,
1690  abs_vec_bulk,
1691  ptype_bulk,
1692  ext_mat_ssbulk,
1693  abs_vec_ssbulk,
1694  ptype_ssbulk);
1695  for (Index iv = 0; iv < nf; iv++) {
1696  if (nf_ssd > 1) {
1697  dap_dx[iq].SetAtPosition(abs_vec_bulk(iv, 0, 0, joker), iv);
1698  dKp_dx[iq].SetAtPosition(ext_mat_bulk(iv, 0, 0, joker, joker),
1699  iv);
1700  } else {
1701  dap_dx[iq].SetAtPosition(abs_vec_bulk(0, 0, 0, joker), iv);
1702  dKp_dx[iq].SetAtPosition(ext_mat_bulk(0, 0, 0, joker, joker), iv);
1703  }
1704  }
1705  })
1706 }
1707 
1709  StokesVector& Sp,
1710  ArrayOfStokesVector& dSp_dx,
1712  ConstVectorView ppath_1p_pnd, // the ppath_pnd at this ppath point
1713  const ArrayOfMatrix&
1714  ppath_dpnd_dx, // the full ppath_dpnd_dx, ie all ppath points
1715  const Index ppath_1p_id,
1720  ConstMatrixView ppath_line_of_sight,
1721  const GridPos& ppath_pressure,
1722  const Vector& temperature,
1723  const Index& atmosphere_dim,
1724  const bool& jacobian_do,
1725  const Index& t_interp_order) {
1726  if (atmosphere_dim != 1)
1727  throw runtime_error("This function handles so far only 1D atmospheres.");
1728 
1729  const Index nf = Sp.NumberOfFrequencies();
1730  const Index stokes_dim = Sp.StokesDimensions();
1731  const Index ne = ppath_1p_pnd.nelem();
1732  assert(TotalNumberOfElements(scat_data) == ne);
1733  const Index nza = za_grid.nelem();
1734  const Index naa = aa_grid.nelem();
1735  const Index nq = jacobian_do ? jacobian_quantities.nelem() : 0;
1736 
1737  // interpolate incident field to this ppath point (no need to do this
1738  // separately per scatelem)
1739  GridPos gp_p;
1740  gridpos_copy(gp_p, ppath_pressure);
1741  Vector itw_p(2);
1742  interpweights(itw_p, gp_p);
1743  Tensor3 inc_field(nf, nza, stokes_dim, 0.);
1744  for (Index iv = 0; iv < nf; iv++) {
1745  for (Index iza = 0; iza < nza; iza++) {
1746  for (Index i = 0; i < stokes_dim; i++) {
1747  inc_field(iv, iza, i) =
1748  interp(itw_p, cloudbox_field(iv, joker, 0, 0, iza, 0, i), gp_p);
1749  }
1750  }
1751  }
1752 
1753  // create matrix of incident directions (flat representation of the
1754  // za_grid * aa_grid matrix)
1755  Matrix idir(nza * naa, 2);
1756  Index ia = 0;
1757  for (Index iza = 0; iza < nza; iza++) {
1758  for (Index iaa = 0; iaa < naa; iaa++) {
1759  idir(ia, 0) = za_grid[iza];
1760  idir(ia, 1) = aa_grid[iaa];
1761  ia++;
1762  }
1763  }
1764 
1765  // setting prop (aka scattered) direction
1766  Matrix pdir(1, 2);
1767  //if( ppath_line_of_sight.ncols()==2 )
1768  // pdir(0,joker) = ppath_line_of_sight;
1769  //else // 1D case only (currently the only handled case). azimuth not defined.
1770  {
1771  pdir(0, 0) = ppath_line_of_sight(0, 0);
1772  pdir(0, 1) = 0.;
1773  }
1774 
1775  // some further variables needed for pha_mat extraction
1776  Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
1777  Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
1778  Tensor6 pha_mat_1se(nf_ssd, 1, 1, nza * naa, stokes_dim, stokes_dim);
1779  Vector t_ok(1);
1780  Index ptype;
1781  Tensor3 scat_source_1se(ne, nf, stokes_dim, 0.);
1782 
1783  Index ise_flat = 0;
1784  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++) {
1785  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
1786  // determine whether we have some valid pnd for this
1787  // scatelem (in pnd or dpnd)
1788  Index val_pnd = 0;
1789  if (ppath_1p_pnd[ise_flat] != 0) {
1790  val_pnd = 1;
1791  } else if (jacobian_do) {
1792  for (Index iq = 0; (!val_pnd) && (iq < nq); iq++) {
1793  if (jacobian_quantities[iq].Analytical() &&
1794  !ppath_dpnd_dx[iq].empty() &&
1795  ppath_dpnd_dx[iq](ise_flat, ppath_1p_id) != 0) {
1796  val_pnd = 1;
1797  }
1798  }
1799  }
1800 
1801  if (val_pnd) {
1802  pha_mat_1ScatElem(pha_mat_1se,
1803  ptype,
1804  t_ok,
1805  scat_data[i_ss][i_se],
1806  temperature,
1807  pdir,
1808  idir,
1809  0,
1810  t_interp_order);
1811  if (t_ok[0] == 0) {
1812  ostringstream os;
1813  os << "Interpolation error for (flat-array) scattering "
1814  << "element #" << ise_flat << "\n"
1815  << "at location/temperature point #" << ppath_1p_id << "\n";
1816  throw runtime_error(os.str());
1817  }
1818 
1819  Index this_iv = 0;
1820  for (Index iv = 0; iv < nf; iv++) {
1821  if (!duplicate_freqs) {
1822  this_iv = iv;
1823  }
1824  Tensor3 product_fields(nza, naa, stokes_dim, 0.);
1825 
1826  ia = 0;
1827  for (Index iza = 0; iza < nza; iza++) {
1828  for (Index iaa = 0; iaa < naa; iaa++) {
1829  for (Index i = 0; i < stokes_dim; i++) {
1830  for (Index j = 0; j < stokes_dim; j++) {
1831  product_fields(iza, iaa, i) +=
1832  pha_mat_1se(this_iv, 0, 0, ia, i, j) *
1833  inc_field(iv, iza, j);
1834  }
1835  }
1836  ia++;
1837  }
1838  }
1839 
1840  for (Index i = 0; i < stokes_dim; i++) {
1841  scat_source_1se(ise_flat, iv, i) = AngIntegrate_trapezoid(
1842  product_fields(joker, joker, i), za_grid, aa_grid);
1843  }
1844  } // for iv
1845  } // if val_pnd
1846 
1847  ise_flat++;
1848 
1849  } // for i_se
1850  } // for i_ss
1851 
1852  for (Index iv = 0; iv < nf; iv++) {
1853  Vector scat_source(stokes_dim, 0.);
1854  for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1855  for (Index i = 0; i < stokes_dim; i++) {
1856  scat_source[i] +=
1857  scat_source_1se(ise_flat, iv, i) * ppath_1p_pnd[ise_flat];
1858  }
1859  }
1860 
1861  Sp.SetAtPosition(scat_source, iv);
1862 
1863  if (jacobian_do) {
1865  if (ppath_dpnd_dx[iq].empty()) { dSp_dx[iq].SetZero(); } else {
1866  scat_source = 0.;
1867  for (ise_flat = 0; ise_flat < ne; ise_flat++) {
1868  for (Index i = 0; i < stokes_dim; i++) {
1869  scat_source[i] += scat_source_1se(ise_flat, iv, i) *
1870  ppath_dpnd_dx[iq](ise_flat, ppath_1p_id);
1871  dSp_dx[iq].SetAtPosition(scat_source, iv);
1872  }
1873  }
1874  })
1875  }
1876  } // for iv
1877 }
1878 
1881  Tensor3View T,
1882  Tensor4View dT_close_dx,
1883  Tensor4View dT_far_dx,
1884  ConstTensor3View cumulative_transmission_close,
1885  const PropagationMatrix& K_close,
1886  const PropagationMatrix& K_far,
1887  const ArrayOfPropagationMatrix& dK_close_dx,
1888  const ArrayOfPropagationMatrix& dK_far_dx,
1889  const Numeric& ppath_distance,
1890  const bool& first_level,
1891  const Numeric& dppath_distance_dT_HSE_close,
1892  const Numeric& dppath_distance_dT_HSE_far,
1893  const Index& temperature_derivative_position_if_hse_is_active) {
1894  // Frequency counter
1895  const Index nf = K_close.NumberOfFrequencies();
1896  const Index stokes_dim = T.ncols();
1897 
1898  if (first_level) {
1899  if (stokes_dim > 1)
1900  for (Index iv = 0; iv < nf; iv++)
1902  else
1904  } else {
1905  // Compute the transmission of the layer between close and far
1906  if (not dK_close_dx.nelem())
1907  compute_transmission_matrix(T, ppath_distance, K_close, K_far);
1908  else
1910  T,
1911  dT_close_dx,
1912  dT_far_dx,
1913  ppath_distance,
1914  K_close,
1915  K_far,
1916  dK_close_dx,
1917  dK_far_dx,
1918  dppath_distance_dT_HSE_close,
1919  dppath_distance_dT_HSE_far,
1920  temperature_derivative_position_if_hse_is_active);
1921 
1922  // Cumulate transmission
1923  if (stokes_dim > 1)
1924  for (Index iv = 0; iv < nf; iv++)
1926  cumulative_transmission_close(iv, joker, joker),
1927  T(iv, joker, joker));
1928  else {
1929  cumulative_transmission = cumulative_transmission_close;
1931  }
1932  }
1933 }
1934 
1935 void iyb_calc_body(bool& failed,
1936  String& fail_msg,
1937  ArrayOfArrayOfMatrix& iy_aux_array,
1938  Workspace& ws,
1939  Ppath& ppath,
1940  Vector& iyb,
1941  ArrayOfMatrix& diyb_dx,
1942  const Index& mblock_index,
1943  const Index& atmosphere_dim,
1944  const EnergyLevelMap& nlte_field,
1945  const Index& cloudbox_on,
1946  const Index& stokes_dim,
1951  const String& iy_unit,
1952  const Agenda& iy_main_agenda,
1953  const Index& j_analytical_do,
1955  const ArrayOfArrayOfIndex& jacobian_indices,
1957  const ArrayOfString& iy_aux_vars,
1958  const Index& ilos,
1959  const Index& nf) {
1960  // The try block here is necessary to correctly handle
1961  // exceptions inside the parallel region.
1962  try {
1963  //--- LOS of interest
1964  //
1965  Vector los(sensor_los.ncols());
1966  //
1967  los = sensor_los(mblock_index, joker);
1968  if (mblock_dlos_grid.ncols() == 1) {
1969  los[0] += mblock_dlos_grid(ilos, 0);
1970  adjust_los(los, atmosphere_dim);
1971  } else {
1972  add_za_aa(los[0],
1973  los[1],
1974  los[0],
1975  los[1],
1976  mblock_dlos_grid(ilos, 0),
1977  mblock_dlos_grid(ilos, 1));
1978  }
1979 
1980  //--- rtp_pos 1 and 2
1981  //
1982  Vector rtp_pos, rtp_pos2(0);
1983  //
1985  if (transmitter_pos.nrows()) {
1986  rtp_pos2 = transmitter_pos(mblock_index, joker);
1987  }
1988 
1989  // Calculate iy and associated variables
1990  //
1991  Matrix iy;
1993  Tensor3 iy_transmission(0, 0, 0);
1994  const Index iy_agenda_call1 = 1;
1995  const Index iy_id =
1996  (Index)1e6 * (mblock_index + 1) + (Index)1e3 * (ilos + 1);
1997  //
1999  iy,
2000  iy_aux_array[ilos],
2001  ppath,
2002  diy_dx,
2005  iy_aux_vars,
2006  iy_id,
2007  iy_unit,
2008  cloudbox_on,
2009  j_analytical_do,
2010  f_grid,
2011  nlte_field,
2012  rtp_pos,
2013  los,
2014  rtp_pos2,
2015  iy_main_agenda);
2016 
2017  // Start row in iyb etc. for present LOS
2018  //
2019  const Index row0 = ilos * nf * stokes_dim;
2020 
2021  // Jacobian part
2022  //
2023  if (j_analytical_do) {
2025  for (Index ip = 0;
2026  ip < jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1;
2027  ip++) {
2028  for (Index is = 0; is < stokes_dim; is++) {
2029  diyb_dx[iq](Range(row0 + is, nf, stokes_dim), ip) =
2030  diy_dx[iq](ip, joker, is);
2031  }
2032  })
2033  }
2034 
2035  // iy : copy to iyb
2036  for (Index is = 0; is < stokes_dim; is++) {
2037  iyb[Range(row0 + is, nf, stokes_dim)] = iy(joker, is);
2038  }
2039 
2040  } // End try
2041 
2042  catch (const std::exception& e) {
2043 #pragma omp critical(iyb_calc_fail)
2044  {
2045  fail_msg = e.what();
2046  failed = true;
2047  }
2048  }
2049 }
2050 
2052  Vector& iyb,
2053  ArrayOfVector& iyb_aux,
2054  ArrayOfMatrix& diyb_dx,
2055  Matrix& geo_pos_matrix,
2056  const Index& mblock_index,
2057  const Index& atmosphere_dim,
2058  const EnergyLevelMap& nlte_field,
2059  const Index& cloudbox_on,
2060  const Index& stokes_dim,
2066  const String& iy_unit,
2067  const Agenda& iy_main_agenda,
2068  const Agenda& geo_pos_agenda,
2069  const Index& j_analytical_do,
2071  const ArrayOfArrayOfIndex& jacobian_indices,
2072  const ArrayOfString& iy_aux_vars,
2073  const Verbosity& verbosity) {
2074  CREATE_OUT3;
2075 
2076  // Sizes
2077  const Index nf = f_grid.nelem();
2078  const Index nlos = mblock_dlos_grid.nrows();
2079  const Index niyb = nf * nlos * stokes_dim;
2080  // Set up size of containers for data of 1 measurement block.
2081  // (can not be made below due to parallalisation)
2082  iyb.resize(niyb);
2083  //
2084  if (j_analytical_do) {
2085  diyb_dx.resize(jacobian_indices.nelem());
2086  FOR_ANALYTICAL_JACOBIANS_DO2(diyb_dx[iq].resize(
2087  niyb, jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1);)
2088  } else {
2089  diyb_dx.resize(0);
2090  }
2091  // Assume that geo_pos_agenda returns empty geo_pos.
2092  geo_pos_matrix.resize(nlos, 5);
2093  geo_pos_matrix = NAN;
2094 
2095  // For iy_aux we don't know the number of quantities, and we have to store
2096  // all outout
2097  ArrayOfArrayOfMatrix iy_aux_array(nlos);
2098 
2099  // We have to make a local copy of the Workspace and the agendas because
2100  // only non-reference types can be declared firstprivate in OpenMP
2101  Workspace l_ws(ws);
2102  Agenda l_iy_main_agenda(iy_main_agenda);
2103  Agenda l_geo_pos_agenda(geo_pos_agenda);
2104 
2105  String fail_msg;
2106  bool failed = false;
2107  if (nlos >= arts_omp_get_max_threads() || nlos * 10 >= nf) {
2108  out3 << " Parallelizing los loop (" << nlos << " iterations, " << nf
2109  << " frequencies)\n";
2110 
2111  // Start of actual calculations
2112 #pragma omp parallel for if (!arts_omp_in_parallel()) \
2113  firstprivate(l_ws, l_iy_main_agenda, l_geo_pos_agenda)
2114  for (Index ilos = 0; ilos < nlos; ilos++) {
2115  // Skip remaining iterations if an error occurred
2116  if (failed) continue;
2117 
2118  Ppath ppath;
2119  iyb_calc_body(failed,
2120  fail_msg,
2121  iy_aux_array,
2122  l_ws,
2123  ppath,
2124  iyb,
2125  diyb_dx,
2126  mblock_index,
2128  nlte_field,
2129  cloudbox_on,
2130  stokes_dim,
2131  sensor_pos,
2132  sensor_los,
2135  iy_unit,
2136  l_iy_main_agenda,
2137  j_analytical_do,
2139  jacobian_indices,
2140  f_grid,
2141  iy_aux_vars,
2142  ilos,
2143  nf);
2144 
2145  // Skip remaining iterations if an error occurred
2146  if (failed) continue;
2147 
2148  // Note that this code is found in two places inside the function
2149  Vector geo_pos;
2150  try {
2151  geo_pos_agendaExecute(l_ws, geo_pos, ppath, l_geo_pos_agenda);
2152  if (geo_pos.nelem()) {
2153  if (geo_pos.nelem() != 5)
2154  throw runtime_error(
2155  "Wrong size of *geo_pos* obtained from *geo_pos_agenda*.\n"
2156  "The length of *geo_pos* must be zero or five.");
2157 
2158  geo_pos_matrix(ilos, joker) = geo_pos;
2159  }
2160  } catch (const std::exception& e) {
2161 #pragma omp critical(iyb_calc_fail)
2162  {
2163  fail_msg = e.what();
2164  failed = true;
2165  }
2166  }
2167  }
2168  } else {
2169  out3 << " Not parallelizing los loop (" << nlos << " iterations, " << nf
2170  << " frequencies)\n";
2171 
2172  for (Index ilos = 0; ilos < nlos; ilos++) {
2173  // Skip remaining iterations if an error occurred
2174  if (failed) continue;
2175 
2176  Ppath ppath;
2177  iyb_calc_body(failed,
2178  fail_msg,
2179  iy_aux_array,
2180  l_ws,
2181  ppath,
2182  iyb,
2183  diyb_dx,
2184  mblock_index,
2186  nlte_field,
2187  cloudbox_on,
2188  stokes_dim,
2189  sensor_pos,
2190  sensor_los,
2193  iy_unit,
2194  l_iy_main_agenda,
2195  j_analytical_do,
2197  jacobian_indices,
2198  f_grid,
2199  iy_aux_vars,
2200  ilos,
2201  nf);
2202 
2203  // Skip remaining iterations if an error occurred
2204  if (failed) continue;
2205 
2206  // Note that this code is found in two places inside the function
2207  Vector geo_pos;
2208  try {
2209  geo_pos_agendaExecute(l_ws, geo_pos, ppath, l_geo_pos_agenda);
2210  if (geo_pos.nelem()) {
2211  if (geo_pos.nelem() != 5)
2212  throw runtime_error(
2213  "Wrong size of *geo_pos* obtained from *geo_pos_agenda*.\n"
2214  "The length of *geo_pos* must be zero or five.");
2215 
2216  geo_pos_matrix(ilos, joker) = geo_pos;
2217  }
2218  } catch (const std::exception& e) {
2219 #pragma omp critical(iyb_calc_fail)
2220  {
2221  fail_msg = e.what();
2222  failed = true;
2223  }
2224  }
2225  }
2226  }
2227 
2228  if (failed)
2229  throw runtime_error("Run-time error in function: iyb_calc\n" + fail_msg);
2230 
2231  // Compile iyb_aux
2232  //
2233  const Index nq = iy_aux_array[0].nelem();
2234  iyb_aux.resize(nq);
2235  //
2236  for (Index q = 0; q < nq; q++) {
2237  iyb_aux[q].resize(niyb);
2238  //
2239  for (Index ilos = 0; ilos < nlos; ilos++) {
2240  const Index row0 = ilos * nf * stokes_dim;
2241  for (Index iv = 0; iv < nf; iv++) {
2242  const Index row1 = row0 + iv * stokes_dim;
2243  const Index i1 = min(iv, iy_aux_array[ilos][q].nrows() - 1);
2244  for (Index is = 0; is < stokes_dim; is++) {
2245  Index i2 = min(is, iy_aux_array[ilos][q].ncols() - 1);
2246  iyb_aux[q][row1 + is] = iy_aux_array[ilos][q](i1, i2);
2247  }
2248  }
2249  }
2250  }
2251 }
2252 
2253 void iy_transmission_mult(Tensor3& iy_trans_total,
2254  ConstTensor3View iy_trans_old,
2255  ConstTensor3View iy_trans_new) {
2256  const Index nf = iy_trans_old.npages();
2257  const Index ns = iy_trans_old.ncols();
2258 
2259  assert(ns == iy_trans_old.nrows());
2260  assert(nf == iy_trans_new.npages());
2261  assert(ns == iy_trans_new.nrows());
2262  assert(ns == iy_trans_new.ncols());
2263 
2264  iy_trans_total.resize(nf, ns, ns);
2265 
2266  for (Index iv = 0; iv < nf; iv++) {
2267  mult(iy_trans_total(iv, joker, joker),
2268  iy_trans_old(iv, joker, joker),
2269  iy_trans_new(iv, joker, joker));
2270  }
2271 }
2272 
2274  ConstTensor3View iy_trans,
2275  ConstMatrixView iy_old) {
2276  const Index nf = iy_trans.npages();
2277  const Index ns = iy_trans.ncols();
2278 
2279  assert(ns == iy_trans.nrows());
2280  assert(nf == iy_old.nrows());
2281  assert(ns == iy_old.ncols());
2282 
2283  iy_new.resize(nf, ns);
2284 
2285  for (Index iv = 0; iv < nf; iv++) {
2286  mult(iy_new(iv, joker), iy_trans(iv, joker, joker), iy_old(iv, joker));
2287  }
2288 }
2289 
2290 void mirror_los(Vector& los_mirrored,
2291  ConstVectorView los,
2292  const Index& atmosphere_dim) {
2293  los_mirrored.resize(2);
2294  //
2295  if (atmosphere_dim == 1) {
2296  los_mirrored[0] = 180 - los[0];
2297  los_mirrored[1] = 180;
2298  } else if (atmosphere_dim == 2) {
2299  los_mirrored[0] = 180 - fabs(los[0]);
2300  if (los[0] >= 0) {
2301  los_mirrored[1] = 180;
2302  } else {
2303  los_mirrored[1] = 0;
2304  }
2305  } else if (atmosphere_dim == 3) {
2306  los_mirrored[0] = 180 - los[0];
2307  los_mirrored[1] = los[1] + 180;
2308  if (los_mirrored[1] > 180) {
2309  los_mirrored[1] -= 360;
2310  }
2311  }
2312 }
2313 
2315  Numeric& lon,
2316  const Index& atmosphere_dim,
2320  ConstVectorView pos) {
2321  assert(pos.nelem() == atmosphere_dim);
2322 
2323  if (atmosphere_dim == 1) {
2324  assert(lat_true.nelem() == 1);
2325  assert(lon_true.nelem() == 1);
2326  //
2327  lat = lat_true[0];
2328  lon = lon_true[0];
2329  }
2330 
2331  else if (atmosphere_dim == 2) {
2332  assert(lat_true.nelem() == lat_grid.nelem());
2333  assert(lon_true.nelem() == lat_grid.nelem());
2334  GridPos gp;
2335  Vector itw(2);
2336  gridpos(gp, lat_grid, pos[1]);
2337  interpweights(itw, gp);
2338  lat = interp(itw, lat_true, gp);
2339  lon = interp(itw, lon_true, gp);
2340  }
2341 
2342  else {
2343  lat = pos[1];
2344  lon = pos[2];
2345  }
2346 }
2347 
2349  ArrayOfIndex& jac_species_i,
2350  ArrayOfIndex& jac_scat_i,
2351  ArrayOfIndex& jac_is_t,
2352  ArrayOfIndex& jac_wind_i,
2353  ArrayOfIndex& jac_mag_i,
2354  ArrayOfIndex& jac_other,
2356  ArrayOfTensor3& diy_dpath,
2357  const Index& ns,
2358  const Index& nf,
2359  const Index& np,
2360  const Index& nq,
2362  const Index& cloudbox_on,
2363  const ArrayOfString& scat_species,
2366  const Index& iy_agenda_call1,
2367  const bool is_active) {
2368  const Index nn = is_active ? nf * np : nf;
2369 
2370  FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq].resize(np, nn, ns);
2371  diy_dpath[iq] = 0.0;)
2372 
2374  jac_scat_i,
2375  jac_is_t,
2376  jac_wind_i,
2377  jac_mag_i,
2379  abs_species,
2380  cloudbox_on,
2381  scat_species);
2382 
2384  jac_other[iq] = (jacobian_quantities[iq].PropMatType() ==
2387  : Index(JacobianType::None);
2388 
2389  if (jac_scat_i[iq] + 1) {
2390  if (dpnd_field_dx[iq].empty())
2391  throw runtime_error(
2392  "*dpnd_field_dx* not allowed to be empty for "
2393  "scattering Jacobian species.");
2394  }
2395  // FIXME: should we indeed check for that? remove if it causes issues.
2396  else {
2397  if (!dpnd_field_dx[iq].empty())
2398  throw runtime_error(
2399  "*dpnd_field_dx* must be empty for "
2400  "non-scattering Jacobian species.");
2401  })
2402 
2403  if (iy_agenda_call1) {
2404  diy_dx.resize(nq);
2405  //
2406  bool any_affine;
2407  ArrayOfArrayOfIndex jacobian_indices;
2408  jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
2409  //
2411  jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1, nn, ns);
2412  diy_dx[iq] = 0.0;)
2413  }
2414 }
2415 
2417  Workspace& ws,
2419  ArrayOfTensor3& diy_dpath,
2420  const Index& ns,
2421  const Index& nf,
2422  const Index& np,
2423  const Index& atmosphere_dim,
2424  const Ppath& ppath,
2425  const Vector& ppvar_p,
2426  const Vector& ppvar_t,
2427  const Matrix& ppvar_vmr,
2428  const Index& iy_agenda_call1,
2429  const Tensor3& iy_transmission,
2430  const Agenda& water_p_eq_agenda,
2432  const ArrayOfIndex jac_species_i,
2433  const ArrayOfIndex jac_is_t) {
2434  // Weight with iy_transmission
2435  if (!iy_agenda_call1) {
2436  Matrix X, Y;
2437  //
2439  Y.resize(ns, diy_dpath[iq].npages());
2440  for (Index iv = 0; iv < nf; iv++) {
2441  X = transpose(diy_dpath[iq](joker, iv, joker));
2442  mult(Y, iy_transmission(iv, joker, joker), X);
2443  diy_dpath[iq](joker, iv, joker) = transpose(Y);
2444  })
2445  }
2446 
2447  // Handle abs species retrieval units, both internally and impact on T-jacobian
2448  //
2449  Tensor3 water_p_eq(0, 0, 0);
2450  //
2451  // Conversion for abs species itself
2452  for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2453  // Let x be VMR, and z the selected retrieval unit.
2454  // We have then that diy/dz = diy/dx * dx/dz
2455  //
2456  if (jacobian_quantities[iq].Analytical() && jac_species_i[iq] >= 0) {
2457  if (jacobian_quantities[iq].Mode() == "vmr") {
2458  }
2459 
2460  else if (jacobian_quantities[iq].Mode() == "rel") {
2461  // Here x = vmr*z
2462  for (Index ip = 0; ip < np; ip++) {
2463  diy_dpath[iq](ip, joker, joker) *= ppvar_vmr(jac_species_i[iq], ip);
2464  }
2465  }
2466 
2467  else if (jacobian_quantities[iq].Mode() == "nd") {
2468  // Here x = z/nd_tot
2469  for (Index ip = 0; ip < np; ip++) {
2470  diy_dpath[iq](ip, joker, joker) /=
2471  number_density(ppvar_p[ip], ppvar_t[ip]);
2472  }
2473  }
2474 
2475  else if (jacobian_quantities[iq].Mode() == "rh") {
2476  // Here x = (p_sat/p) * z
2477  Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2478  t_data(joker, 0, 0) = ppvar_t;
2479  water_p_eq_agendaExecute(ws, water_p_eq, t_data, water_p_eq_agenda);
2480  for (Index ip = 0; ip < np; ip++) {
2481  diy_dpath[iq](ip, joker, joker) *= water_p_eq(ip, 0, 0) / ppvar_p[ip];
2482  }
2483  }
2484 
2485  else if (jacobian_quantities[iq].Mode() == "q") {
2486  // Here we use the approximation of x = z/0.622
2487  diy_dpath[iq](joker, joker, joker) /= 0.622;
2488  }
2489 
2490  else {
2491  assert(0);
2492  }
2493  }
2494  }
2495 
2496  // Correction of temperature Jacobian
2497  for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2498  // Let a be unit for abs species, and iy = f(T,a(T))
2499  // We have then that diy/dT = df/dT + df/da*da/dT
2500  // diy_dpath holds already df/dT. Remains is to add
2501  // df/da*da/dT for which abs species having da/dT != 0
2502  // This is only true for "nd" and "rh"
2503  //
2504  if (jac_is_t[iq] != Index(JacobianType::None)) {
2505  // Loop abs species, again
2506  for (Index ia = 0; ia < jacobian_quantities.nelem(); ia++) {
2507  if (jacobian_quantities[iq].Analytical() && jac_species_i[ia] >= 0) {
2508  if (jacobian_quantities[ia].Mode() == "nd") {
2509  for (Index ip = 0; ip < np; ip++) {
2510  Matrix ddterm = diy_dpath[ia](ip, joker, joker);
2511  ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2512  (number_density(ppvar_p[ip], ppvar_t[ip] + 1) -
2513  number_density(ppvar_p[ip], ppvar_t[ip]));
2514  diy_dpath[iq](ip, joker, joker) += ddterm;
2515  }
2516  } else if (jacobian_quantities[ia].Mode() == "rh") {
2517  Tensor3 t_data(ppvar_t.nelem(), 1, 1);
2518  t_data(joker, 0, 0) = ppvar_t;
2519  // Calculate water sat. pressure if not already done
2520  if (water_p_eq.npages() == 0) {
2522  ws, water_p_eq, t_data, water_p_eq_agenda);
2523  }
2524  // Sat.pressure for +1K
2525  Tensor3 water_p_eq1K;
2526  t_data(joker, 0, 0) += 1;
2528  ws, water_p_eq1K, t_data, water_p_eq_agenda);
2529 
2530  for (Index ip = 0; ip < np; ip++) {
2531  const Numeric p_eq = water_p_eq(ip, 0, 0);
2532  const Numeric p_eq1K = water_p_eq1K(ip, 0, 0);
2533  Matrix ddterm = diy_dpath[ia](ip, joker, joker);
2534  ddterm *= ppvar_vmr(jac_species_i[ia], ip) *
2535  (ppvar_p[ip] / pow(p_eq, 2.0)) * (p_eq1K - p_eq);
2536  diy_dpath[iq](ip, joker, joker) += ddterm;
2537  }
2538  }
2539  }
2540  }
2541  }
2542  }
2543 
2544  // Map to retrieval grids
2546  jacobian_quantities[iq],
2547  diy_dpath[iq],
2549  ppath,
2550  ppvar_p);)
2551 }
2552 
2554  Matrix& iy,
2556  Tensor3& ppvar_iy,
2557  const Index& ns,
2558  const Index& np,
2559  const Vector& f_grid,
2560  const Ppath& ppath,
2562  const Index& j_analytical_do,
2563  const String& iy_unit) {
2564  // Determine refractive index to use for the n2 radiance law
2565  Numeric n = 1.0; // First guess is that sensor is in space
2566  //
2567  if (ppath.end_lstep == 0) // If true, sensor inside the atmosphere
2568  {
2569  n = ppath.nreal[np - 1];
2570  }
2571 
2572  // Polarisation index variable
2573  ArrayOfIndex i_pol(ns);
2574  for (Index is = 0; is < ns; is++) {
2575  i_pol[is] = is + 1;
2576  }
2577 
2578  // Jacobian part (must be converted to Tb before iy for PlanckBT)
2579  //
2580  if (j_analytical_do) {
2582  apply_iy_unit2(diy_dx[iq], iy, iy_unit, f_grid, n, i_pol);)
2583  }
2584 
2585  // iy
2586  apply_iy_unit(iy, iy_unit, f_grid, n, i_pol);
2587 
2588  // ppvar_iy
2589  for (Index ip = 0; ip < ppath.np; ip++) {
2590  apply_iy_unit(
2591  ppvar_iy(joker, joker, ip), iy_unit, f_grid, ppath.nreal[ip], i_pol);
2592  }
2593 }
2594 
2595 void yCalc_mblock_loop_body(bool& failed,
2596  String& fail_msg,
2597  ArrayOfArrayOfVector& iyb_aux_array,
2598  Workspace& ws,
2599  Vector& y,
2600  Vector& y_f,
2602  Matrix& y_pos,
2603  Matrix& y_los,
2604  Matrix& y_geo,
2605  Matrix& jacobian,
2606  const Index& atmosphere_dim,
2607  const EnergyLevelMap& nlte_field,
2608  const Index& cloudbox_on,
2609  const Index& stokes_dim,
2610  const Vector& f_grid,
2611  const Matrix& sensor_pos,
2612  const Matrix& sensor_los,
2613  const Matrix& transmitter_pos,
2614  const Matrix& mblock_dlos_grid,
2615  const Sparse& sensor_response,
2616  const Vector& sensor_response_f,
2619  const String& iy_unit,
2620  const Agenda& iy_main_agenda,
2621  const Agenda& geo_pos_agenda,
2622  const Agenda& jacobian_agenda,
2623  const Index& jacobian_do,
2625  const ArrayOfArrayOfIndex& jacobian_indices,
2626  const ArrayOfString& iy_aux_vars,
2627  const Verbosity& verbosity,
2628  const Index& mblock_index,
2629  const Index& n1y,
2630  const Index& j_analytical_do) {
2631  try {
2632  // Calculate monochromatic pencil beam data for 1 measurement block
2633  //
2634  Vector iyb, iyb_error, yb(n1y);
2635  ArrayOfMatrix diyb_dx;
2636  Matrix geo_pos_matrix;
2637  //
2638  iyb_calc(ws,
2639  iyb,
2640  iyb_aux_array[mblock_index],
2641  diyb_dx,
2642  geo_pos_matrix,
2643  mblock_index,
2645  nlte_field,
2646  cloudbox_on,
2647  stokes_dim,
2648  f_grid,
2649  sensor_pos,
2650  sensor_los,
2653  iy_unit,
2656  j_analytical_do,
2658  jacobian_indices,
2659  iy_aux_vars,
2660  verbosity);
2661 
2662  // Apply sensor response matrix on iyb, and put into y
2663  //
2665  const Index row0 = rowind.get_start();
2666  //
2668  //
2669  y[rowind] = yb; // *yb* also used below, as input to jacobian_agenda
2670 
2671  // Fill information variables. And search for NaNs in *y*.
2672  //
2673  for (Index i = 0; i < n1y; i++) {
2674  const Index ii = row0 + i;
2675  if (std::isnan(y[ii]))
2676  throw runtime_error("One or several NaNs found in *y*.");
2677  y_f[ii] = sensor_response_f[i];
2678  y_pol[ii] = sensor_response_pol[i];
2681  y_los(ii, 0) += sensor_response_dlos(i, 0);
2682  if (sensor_response_dlos.ncols() > 1) {
2683  y_los(ii, 1) += sensor_response_dlos(i, 1);
2684  }
2685  }
2686 
2687  // Apply sensor response matrix on diyb_dx, and put into jacobian
2688  // (that is, analytical jacobian part)
2689  //
2690  if (j_analytical_do) {
2692  mult(jacobian(rowind,
2693  Range(jacobian_indices[iq][0],
2694  jacobian_indices[iq][1] -
2695  jacobian_indices[iq][0] + 1)),
2697  diyb_dx[iq]);)
2698  }
2699 
2700  // Calculate remaining parts of *jacobian*
2701  //
2702  if (jacobian_do) {
2705  }
2706 
2707  // Handle geo-positioning
2708  if (!std::isnan(geo_pos_matrix(0, 0))) // No data are flagged as NaN
2709  {
2710  // We set geo_pos based on the max value in sensor_response
2711  const Index nfs = f_grid.nelem() * stokes_dim;
2712  for (Index i = 0; i < n1y; i++) {
2713  Index jmax = -1;
2714  Numeric rmax = -99e99;
2715  for (Index j = 0; j < sensor_response.ncols(); j++) {
2716  if (sensor_response(i, j) > rmax) {
2717  rmax = sensor_response(i, j);
2718  jmax = j;
2719  }
2720  }
2721  const Index jhit = Index(floor(jmax / nfs));
2722  y_geo(row0 + i, joker) = geo_pos_matrix(jhit, joker);
2723  }
2724  }
2725  }
2726 
2727  catch (const std::exception& e) {
2728 #pragma omp critical(yCalc_fail)
2729  {
2730  fail_msg = e.what();
2731  failed = true;
2732  }
2733  }
2734 }
2735 
2737  const Vector& f_grid,
2738  const Numeric& ze_tref,
2739  const Numeric& k2) {
2740  const Index nf = f_grid.nelem();
2741 
2742  assert(fac.nelem() == nf);
2743 
2744  // Refractive index for water (if needed)
2745  Matrix complex_n(0, 0);
2746  if (k2 <= 0) {
2747  complex_n_water_liebe93(complex_n, f_grid, ze_tref);
2748  }
2749 
2750  // Common conversion factor
2751  const Numeric a = 4e18 / (PI * PI * PI * PI);
2752 
2753  for (Index iv = 0; iv < nf; iv++) {
2754  // Reference dielectric factor.
2755  Numeric K2;
2756  if (k2 >= 0) {
2757  K2 = k2;
2758  } else {
2759  Complex n(complex_n(iv, 0), complex_n(iv, 1));
2760  Complex n2 = n * n;
2761  Complex K = (n2 - Numeric(1.0)) / (n2 + Numeric(2.0));
2762  Numeric absK = abs(K);
2763  K2 = absK * absK;
2764  }
2765 
2766  // Wavelength
2767  Numeric la = SPEED_OF_LIGHT / f_grid[iv];
2768 
2769  fac[iv] = a * la * la * la * la / K2;
2770  }
2771 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::y_pos
Matrix y_pos(Workspace &ws) noexcept
Definition: autoarts.h:7531
legacy_continua.h
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::Ppath::pos
std::size_t pos() const noexcept
Definition: autoarts.h:1300
ARTS::Var::wind_v_field
Tensor3 wind_v_field(Workspace &ws) noexcept
Definition: autoarts.h:7246
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
pha_mat_1ScatElem
void pha_mat_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
Definition: optproperties.cc:881
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
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
ARTS::Var::geo_pos
Vector geo_pos(Workspace &ws) noexcept
Definition: autoarts.h:3544
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
SURFACE_MAINTAG
const String SURFACE_MAINTAG
geo_pos_agendaExecute
void geo_pos_agendaExecute(Workspace &ws, Vector &geo_pos, const Ppath &ppath, const Agenda &input_agenda)
Definition: auto_md.cc:23915
ARTS::Var::rtp_los
Vector rtp_los(Workspace &ws) noexcept
Definition: autoarts.h:5715
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
ARTS::Var::ppath_step_agenda
Agenda ppath_step_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5220
PROPMAT_SUBSUBTAG
const String PROPMAT_SUBSUBTAG
ARTS::Var::lat
Numeric lat(Workspace &ws) noexcept
Definition: autoarts.h:3933
bending_angle1d
void bending_angle1d(Numeric &alpha, const Ppath &ppath)
Calculates the bending angle for a 1D atmosphere.
Definition: rte.cc:327
auto_md.h
ARTS::Var::iy_cloudbox_agenda
Agenda iy_cloudbox_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3742
gridpos_copy
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Definition: interpolation.cc:473
PropagationMatrix::IsEmpty
bool IsEmpty() const
Asks if the class is empty.
Definition: propagationmatrix.h:221
iyb_calc_body
void iyb_calc_body(bool &failed, String &fail_msg, ArrayOfArrayOfMatrix &iy_aux_array, Workspace &ws, Ppath &ppath, Vector &iyb, ArrayOfMatrix &diyb_dx, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, ConstMatrixView sensor_pos, ConstMatrixView sensor_los, ConstMatrixView transmitter_pos, ConstMatrixView mblock_dlos_grid, const String &iy_unit, const Agenda &iy_main_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, ConstVectorView f_grid, const ArrayOfString &iy_aux_vars, const Index &ilos, const Index &nf)
Definition: rte.cc:1935
iy_space_agendaExecute
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:24299
ARTS::Var::lo
Numeric lo(Workspace &ws) noexcept
Definition: autoarts.h:4039
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
get_stepwise_frequency_grid
void get_stepwise_frequency_grid(VectorView ppath_f_grid, ConstVectorView f_grid, ConstVectorView ppath_wind, ConstVectorView ppath_line_of_sight, const Numeric &rte_alonglos_v, const Index &atmosphere_dim)
Inverse of get_stepwise_f_partials.
Definition: rte.cc:1521
ARTS::Var::transmitter_pos
Matrix transmitter_pos(Workspace &ws) noexcept
Definition: autoarts.h:7093
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
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
add_za_aa
void add_za_aa(Numeric &za, Numeric &aa, const Numeric &za0, const Numeric &aa0, const Numeric &dza, const Numeric &daa)
Adds up zenith and azimuth angles.
Definition: ppath.cc:406
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
ARTS::Var::ppvar_vmr
Matrix ppvar_vmr(Workspace &ws) noexcept
Definition: autoarts.h:5372
itw2p
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Definition: special_interp.cc:718
ARTS::Var::iy_transmission
Tensor3 iy_transmission(Workspace &ws) noexcept
Definition: autoarts.h:3836
ARTS::Var::sensor_response
Sparse sensor_response(Workspace &ws) noexcept
Definition: autoarts.h:6300
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
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
ConstTensor7View
A constant view of a Tensor7.
Definition: matpackVII.h:147
joker
const Joker joker
ze_cfac
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
Definition: rte.cc:2736
get_ppath_cloudvars
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Definition: rte.cc:1153
ARTS::Var::pnd_field
Tensor4 pnd_field(Workspace &ws) noexcept
Definition: autoarts.h:5081
ARTS::Var::iy_id
Index iy_id(Workspace &ws) noexcept
Definition: autoarts.h:3775
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
dplanck_dt
Numeric dplanck_dt(const Numeric &f, const Numeric &t)
dplanck_dt
Definition: physics_funcs.cc:337
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
sphdist
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1205
ARTS::Var::scat_data_mono
ArrayOfArrayOfSingleScatteringData scat_data_mono(Workspace &ws) noexcept
Definition: autoarts.h:5862
ARTS::Var::dnlte_dx_source
ArrayOfStokesVector dnlte_dx_source(Workspace &ws) noexcept
Definition: autoarts.h:3021
PropagationMatrix::MatrixInverseAtPosition
void MatrixInverseAtPosition(MatrixView ret, const Index iv=0, const Index iz=0, const Index ia=0) const
Return the matrix inverse at the position.
Definition: propagationmatrix.cc:1506
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
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
ARTS::Var::ppvar_iy
Tensor3 ppvar_iy(Workspace &ws) noexcept
Definition: autoarts.h:5246
Sparse
The Sparse class.
Definition: matpackII.h:60
chk_not_empty
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:694
vmrunitscf
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion between gas species units.
Definition: jacobian.cc:1001
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
ARTS::Var::jacobian_do
Index jacobian_do(Workspace &ws) noexcept
Definition: autoarts.h:3912
dinvplanckdI
Numeric dinvplanckdI(const Numeric &i, const Numeric &f)
dinvplanckdI
Definition: physics_funcs.cc:98
Ppath::end_lstep
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
ARTS::Var::wind_u_field
Tensor3 wind_u_field(Workspace &ws) noexcept
Definition: autoarts.h:7211
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
PropagationMatrix
Definition: propagationmatrix.h:87
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ARTS::Var::diy_dx
ArrayOfTensor3 diy_dx(Workspace &ws) noexcept
Definition: autoarts.h:3007
diy_from_path_to_rgrids
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Maps jacobian data for points along the propagation path, to jacobian retrieval grid data.
Definition: jacobian.cc:337
ARTS::Var::y_los
Matrix y_los(Workspace &ws) noexcept
Definition: autoarts.h:7498
is_singular
bool is_singular(ConstMatrixView A)
Checks if a square matrix is singular.
Definition: logic.cc:295
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
opt_prop_ScatSpecBulk
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
Definition: optproperties.cc:148
arts_omp_get_max_threads
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Definition: arts_omp.cc:46
ARTS::Var::mag_w_field
Tensor3 mag_w_field(Workspace &ws) noexcept
Definition: autoarts.h:4197
PropagationMatrix::IsRotational
bool IsRotational(const Index iv=0, const Index iz=0, const Index ia=0) const
False if diagonal element is non-zero in internal Matrix representation.
Definition: propagationmatrix.h:247
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:405
ARTS::Var::iy_unit
String iy_unit(Workspace &ws) noexcept
Definition: autoarts.h:3856
DEG2RAD
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
JacobianType::None
@ None
ARTS::Var::iy
Matrix iy(Workspace &ws) noexcept
Definition: autoarts.h:3690
cart2sph
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition: geodetic.cc:563
complex_n_water_liebe93
void complex_n_water_liebe93(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_water_liebe93
Definition: refraction.cc:71
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
montecarlo.h
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
compute_transmission_matrix
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Definition: propagationmatrix.cc:33
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
matrix_exp
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:391
ARTS::Var::iy_space_agenda
Agenda iy_space_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3801
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
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
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
Agenda
The Agenda class.
Definition: agenda_class.h:44
ARTS::Var::lon
Numeric lon(Workspace &ws) noexcept
Definition: autoarts.h:4064
ARTS::Var::dsurface_names
ArrayOfString dsurface_names(Workspace &ws) noexcept
Definition: autoarts.h:3344
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
dxdvmrscf
void dxdvmrscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion of derivatives with respect to VMR.
Definition: jacobian.cc:1029
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
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
Range::get_start
constexpr Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:327
ARTS::Var::lon_true
Vector lon_true(Workspace &ws) noexcept
Definition: autoarts.h:4113
get_stepwise_scattersky_propmat
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, ConstMatrixView ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstVectorView ppath_line_of_sight, ConstVectorView ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
Definition: rte.cc:1591
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
invrayjean
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
Definition: physics_funcs.cc:192
opt_prop_NScatElems
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
Definition: optproperties.cc:248
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
ARTS::Var::rte_alonglos_v
Numeric rte_alonglos_v(Workspace &ws) noexcept
Definition: autoarts.h:5628
ARTS::Var::nlte_field
EnergyLevelMap nlte_field(Workspace &ws) noexcept
Definition: autoarts.h:4604
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
number_density
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
Definition: physics_funcs.cc:219
ARTS::Var::sensor_response_f
Vector sensor_response_f(Workspace &ws) noexcept
Definition: autoarts.h:6376
cart2pol
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:74
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
min
#define min(a, b)
Definition: legacy_continua.cc:20628
get_stepwise_effective_source
void get_stepwise_effective_source(MatrixView J, Tensor3View dJ_dx, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK_dx, const ArrayOfStokesVector &da_dx, const ArrayOfStokesVector &dS_dx, ConstVectorView B, ConstVectorView dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Gets the effective source at propagation path point.
Definition: rte.cc:1427
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:207
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
ARTS::Var::y_f
Vector y_f(Workspace &ws) noexcept
Definition: autoarts.h:7466
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
my_basic_string< char >
ARTS::Var::y_pol
ArrayOfIndex y_pol(Workspace &ws) noexcept
Definition: autoarts.h:7515
Ppath::los
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
ARTS::Var::nlte_source
ArrayOfStokesVector nlte_source(Workspace &ws) noexcept
Definition: autoarts.h:4648
Ppath::np
Index np
Number of points describing the ppath.
Definition: ppath.h:52
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
ArrayOfIndex
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:40
VectorView
The VectorView class.
Definition: matpackI.h:610
iy_cloudbox_agendaExecute
void iy_cloudbox_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:23996
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::rte_los
Vector rte_los(Workspace &ws) noexcept
Definition: autoarts.h:5651
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
ARTS::Var::ppvar_t
Vector ppvar_t(Workspace &ws) noexcept
Definition: autoarts.h:5330
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
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::ppath
Ppath ppath(Workspace &ws) noexcept
Definition: autoarts.h:5139
PropagationMatrix::NumberOfFrequencies
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
Definition: propagationmatrix.h:204
ARTS::Var::dsurface_emission_dx
ArrayOfMatrix dsurface_emission_dx(Workspace &ws) noexcept
Definition: autoarts.h:3333
ARTS::Var::surface_props_data
Tensor3 surface_props_data(Workspace &ws) noexcept
Definition: autoarts.h:6745
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:49
ARTS::Var::iyb
Vector iyb(Workspace &ws) noexcept
Definition: autoarts.h:3873
iyb_calc
void iyb_calc(Workspace &ws, Vector &iyb, ArrayOfVector &iyb_aux, ArrayOfMatrix &diyb_dx, Matrix &geo_pos_matrix, const Index &mblock_index, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, ConstVectorView f_grid, ConstMatrixView sensor_pos, ConstMatrixView sensor_los, ConstMatrixView transmitter_pos, ConstMatrixView mblock_dlos_grid, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const Index &j_analytical_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
Performs calculations for one measurement block, on iy-level.
Definition: rte.cc:2051
propmat_clearsky_agendaExecute
void propmat_clearsky_agendaExecute(Workspace &ws, ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfStokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_dx_source, ArrayOfStokesVector &nlte_dsource_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:23513
ARTS::Var::iy_agenda_call1
Index iy_agenda_call1(Workspace &ws) noexcept
Definition: autoarts.h:3704
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
math_funcs.h
invplanck
Numeric invplanck(const Numeric &i, const Numeric &f)
invplanck
Definition: physics_funcs.cc:165
ppath_calc
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
This is the core for the WSM ppathStepByStep.
Definition: ppath.cc:5206
Ppath::lstep
Vector lstep
The length between ppath points.
Definition: ppath.h:70
geodetic.h
ARTS::Var::refellipsoid
Vector refellipsoid(Workspace &ws) noexcept
Definition: autoarts.h:5519
ARTS::Var::iy_aux_vars
ArrayOfString iy_aux_vars(Workspace &ws) noexcept
Definition: autoarts.h:3735
zaaa2cart
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to a cartesian unit vector.
Definition: ppath.cc:347
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
get_stepwise_scattersky_source
void get_stepwise_scattersky_source(StokesVector &Sp, ArrayOfStokesVector &dSp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, ConstVectorView ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstTensor7View cloudbox_field, ConstVectorView za_grid, ConstVectorView aa_grid, ConstMatrixView ppath_line_of_sight, const GridPos &ppath_pressure, const Vector &temperature, const Index &atmosphere_dim, const bool &jacobian_do, const Index &t_interp_order)
Calculates the stepwise scattering source terms.
Definition: rte.cc:1708
ARTS::Var::lat_true
Vector lat_true(Workspace &ws) noexcept
Definition: autoarts.h:3985
ARTS::Var::sensor_pos
Matrix sensor_pos(Workspace &ws) noexcept
Definition: autoarts.h:6279
ARTS::Var::ext_mat
PropagationMatrix ext_mat(Workspace &ws) noexcept
Definition: autoarts.h:3390
EnergyLevelMap
Definition: energylevelmap.h:60
max
#define max(a, b)
Definition: legacy_continua.cc:20629
defocusing_sat2sat
void defocusing_sat2sat(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Ppath &ppath, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
Calculates defocusing for limb measurements between two satellites.
Definition: rte.cc:613
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
AngIntegrate_trapezoid
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:296
SCATSPECIES_MAINTAG
const String SCATSPECIES_MAINTAG
a1
#define a1
Definition: complex.h:56
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Ppath::nreal
Vector nreal
The real part of the refractive index at each path position.
Definition: ppath.h:78
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
pos2true_latlon
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true, ConstVectorView pos)
Determines the true alt and lon for an "ARTS position".
Definition: rte.cc:2314
ARTS::Var::dpnd_field_dx
ArrayOfTensor4 dpnd_field_dx(Workspace &ws) noexcept
Definition: autoarts.h:3251
ARTS::Var::ppath_lmax
Numeric ppath_lmax(Workspace &ws) noexcept
Definition: autoarts.h:5183
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
ARTS::Var::ppath_lraytrace
Numeric ppath_lraytrace(Workspace &ws) noexcept
Definition: autoarts.h:5195
ARTS::Var::iy_surface_agenda
Agenda iy_surface_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3808
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
ARTS::Var::nrows
Index nrows(Workspace &ws) noexcept
Definition: autoarts.h:4682
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:1527
PropagationMatrix::SetZero
void SetZero()
Sets all data to zero.
Definition: propagationmatrix.h:934
water_p_eq_agendaExecute
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:25260
planck
Numeric planck(const Numeric &f, const Numeric &t)
planck
Definition: physics_funcs.cc:269
ARTS::Var::mblock_dlos_grid
Matrix mblock_dlos_grid(Workspace &ws) noexcept
Definition: autoarts.h:4247
Range
The range class.
Definition: matpackI.h:160
ARTS::Var::sensor_response_pol
ArrayOfIndex sensor_response_pol(Workspace &ws) noexcept
Definition: autoarts.h:6411
get_stepwise_f_partials
void get_stepwise_f_partials(Vector &f_partials, const Index &component, ConstVectorView &line_of_sight, ConstVectorView f_grid, const Index &atmosphere_dim)
Computes the ratio that a partial derivative with regards to frequency relates to the wind of come co...
Definition: rte.cc:1541
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
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
FOR_ANALYTICAL_JACOBIANS_DO2
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:412
ppath.h
Propagation path structure and functions.
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
q
#define q
Definition: legacy_continua.cc:21712
ARTS::Var::iy_aux
ArrayOfMatrix iy_aux(Workspace &ws) noexcept
Definition: autoarts.h:3719
ARTS::Var::geo_pos_agenda
Agenda geo_pos_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3551
find_tanpoint
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
Definition: ppath.cc:525
ARTS::Var::dpropmat_clearsky_dx
ArrayOfPropagationMatrix dpropmat_clearsky_dx(Workspace &ws) noexcept
Definition: autoarts.h:3269
Ppath::constant
Numeric constant
The propagation path constant (only used for 1D)
Definition: ppath.h:54
logic.h
Header file for logic.cc.
ARTS::Var::iy_main_agenda
Agenda iy_main_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3794
ARTS::Var::rte_pos2
Vector rte_pos2(Workspace &ws) noexcept
Definition: autoarts.h:5696
opt_prop_Bulk
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
Definition: optproperties.cc:102
refraction.h
Refraction functions.
mirror_los
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:2290
PropagationMatrix::StokesDimensions
Index StokesDimensions() const
The stokes dimension of the propagation matrix.
Definition: propagationmatrix.h:201
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
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
dotprod_with_los
Numeric dotprod_with_los(ConstVectorView los, const Numeric &u, const Numeric &v, const Numeric &w, const Index &atmosphere_dim)
Calculates the dot product between a field and a LOS.
Definition: rte.cc:730
TEMP_0_C
const Numeric TEMP_0_C
ARTS::Var::yb
Vector yb(Workspace &ws) noexcept
Definition: autoarts.h:7542
rte.h
Declaration of functions in rte.cc.
ext2trans
void ext2trans(MatrixView trans_mat, Index &icase, ConstMatrixView ext_mat, const Numeric &lstep)
Converts an extinction matrix to a transmission matrix.
Definition: rte.cc:800
ARTS::Var::mag_v_field
Tensor3 mag_v_field(Workspace &ws) noexcept
Definition: autoarts.h:4164
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
cart2zaaa
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Converts a cartesian directional vector to zenith and azimuth.
Definition: ppath.cc:312
get_iy
void get_iy(Workspace &ws, Matrix &iy, const Index &cloudbox_on, ConstVectorView f_grid, const EnergyLevelMap &nlte_field, ConstVectorView rte_pos, ConstVectorView rte_los, ConstVectorView rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda)
Basic call of iy_main_agenda.
Definition: rte.cc:877
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
Ppath::r
Vector r
Radius of each ppath point.
Definition: ppath.h:68
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:284
defocusing_general
void defocusing_general(Workspace &ws, Numeric &dlf, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Ppath &ppath, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Numeric &dza, const Verbosity &verbosity)
Defocusing for arbitrary geometry (zenith angle part only)
Definition: rte.cc:494
dx
#define dx
Definition: legacy_continua.cc:22180
special_interp.h
Header file for special_interp.cc.
equivalent_propmattype_index
Index equivalent_propmattype_index(const ArrayOfRetrievalQuantity &js, const Index i) noexcept
Returns a list of positions for the derivatives in Propagation Matrix calculations.
Definition: jacobian.cc:1107
ARTS::Var::wind_w_field
Tensor3 wind_w_field(Workspace &ws) noexcept
Definition: autoarts.h:7279
compute_transmission_matrix_and_derivative
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_dx_upper_level, Tensor4View dT_dx_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const ArrayOfPropagationMatrix &dupper_level_dx, const ArrayOfPropagationMatrix &dlower_level_dx, const Numeric &dr_dTu, const Numeric &dr_dTl, const Index it, const Index iz, const Index ia)
Definition: propagationmatrix.cc:478
poslos2cart
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:355
ARTS::Var::cloudbox_field
Tensor7 cloudbox_field(Workspace &ws) noexcept
Definition: autoarts.h:2676
get_stepwise_transmission_matrix
void get_stepwise_transmission_matrix(Tensor3View cumulative_transmission, Tensor3View T, Tensor4View dT_close_dx, Tensor4View dT_far_dx, ConstTensor3View cumulative_transmission_close, const PropagationMatrix &K_close, const PropagationMatrix &K_far, const ArrayOfPropagationMatrix &dK_close_dx, const ArrayOfPropagationMatrix &dK_far_dx, const Numeric &ppath_distance, const bool &first_level, const Numeric &dppath_distance_dT_HSE_close, const Numeric &dppath_distance_dT_HSE_far, const Index &temperature_derivative_position_if_hse_is_active)
Computes layer transmission matrix and cumulative transmission.
Definition: rte.cc:1879
jacobian_agendaExecute
void jacobian_agendaExecute(Workspace &ws, Matrix &jacobian, const Index mblock_index, const Vector &iyb, const Vector &yb, const Agenda &input_agenda)
Definition: auto_md.cc:24546
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
a2
#define a2
Definition: complex.h:58
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
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:24221
iy_surface_agendaExecute
void iy_surface_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, ArrayOfTensor4 &dsurface_rmatrix_dx, ArrayOfMatrix &dsurface_emission_dx, const String &iy_unit, const Tensor3 &iy_transmission, const Index iy_id, const Index cloudbox_on, const Index jacobian_do, const Agenda &iy_main_agenda, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rtp_pos, const Vector &rtp_los, const Vector &rte_pos2, const Tensor3 &surface_props_data, const ArrayOfString &dsurface_names, const Agenda &input_agenda)
Definition: auto_md.cc:24341
ARTS::Var::jacobian
Matrix jacobian(Workspace &ws) noexcept
Definition: autoarts.h:3892
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ARTS::Var::sensor_response_dlos
Matrix sensor_response_dlos(Workspace &ws) noexcept
Definition: autoarts.h:6340
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
ARTS::Var::scat_species
ArrayOfString scat_species(Workspace &ws) noexcept
Definition: autoarts.h:6086
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
check_input.h
ARTS::Var::mag_u_field
Tensor3 mag_u_field(Workspace &ws) noexcept
Definition: autoarts.h:4130
defocusing_general_sub
void defocusing_general_sub(Workspace &ws, Vector &pos, Vector &rte_los, Index &background, ConstVectorView rte_pos, const Numeric &lo0, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Verbosity &verbosity)
Just to avoid duplicatuion of code in defocusing_general.
Definition: rte.cc:372
ARTS::Var::ncols
Index ncols(Workspace &ws) noexcept
Definition: autoarts.h:4546
cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(MatrixView F, ConstMatrixView A)
Definition: lin_alg.cc:963
adjust_los
void adjust_los(VectorView los, const Index &atmosphere_dim)
Ensures that the zenith and azimuth angles of a line-of-sight vector are inside defined ranges.
Definition: rte.cc:140
ARTS::Var::mblock_index
Index mblock_index(Workspace &ws) noexcept
Definition: autoarts.h:4258
TEMPERATURE_MAINTAG
const String TEMPERATURE_MAINTAG
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1494
Vector
The Vector class.
Definition: matpackI.h:860
ARTS::Var::jacobian_agenda
Agenda jacobian_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3899
ARTS::Var::propmat_clearsky
ArrayOfPropagationMatrix propmat_clearsky(Workspace &ws) noexcept
Definition: autoarts.h:5398
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
distance3D
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition: geodetic.cc:618
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
ARTS::Var::ppvar_p
Vector ppvar_p(Workspace &ws) noexcept
Definition: autoarts.h:5304
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
StokesVector::SetAtPosition
void SetAtPosition(ConstVectorView x, const Index iv=0, const Index iz=0, const Index ia=0)
Definition: propagationmatrix.h:1285
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
ARTS::Var::water_p_eq_agenda
Agenda water_p_eq_agenda(Workspace &ws) noexcept
Definition: autoarts.h:7165
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Ppath::pos
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath.h:64
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Var::sensor_los
Matrix sensor_los(Workspace &ws) noexcept
Definition: autoarts.h:6201
ARTS::Var::y_geo
Matrix y_geo(Workspace &ws) noexcept
Definition: autoarts.h:7480
ARTS::Var::dsurface_rmatrix_dx
ArrayOfTensor4 dsurface_rmatrix_dx(Workspace &ws) noexcept
Definition: autoarts.h:3356
ext_mat_case
void ext_mat_case(Index &icase, ConstMatrixView ext_mat, const Index stokes_dim)
Definition: rte.cc:751
get_pointers_for_analytical_jacobians
void get_pointers_for_analytical_jacobians(ArrayOfIndex &abs_species_i, ArrayOfIndex &scat_species_i, ArrayOfIndex &is_t, ArrayOfIndex &wind_i, ArrayOfIndex &magfield_i, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &cloudbox_on, const ArrayOfString &scat_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:596
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
JacPropMatType::VMR
@ VMR
ARTS::Var::z_surface
Matrix z_surface(Workspace &ws) noexcept
Definition: autoarts.h:7754
StokesVector::VectorAtPosition
VectorView VectorAtPosition(const Index iv=0, const Index iz=0, const Index ia=0)
Get a vectorview to the position.
Definition: propagationmatrix.h:1266
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
refell2r
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1135
PI
const Numeric PI
Global constant, pi.
ARTS::Var::rtp_pos
Vector rtp_pos(Workspace &ws) noexcept
Definition: autoarts.h:5774
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
ns
#define ns
Definition: legacy_continua.cc:22183
distance2D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:182
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
PropagationMatrix::SetAtPosition
void SetAtPosition(const PropagationMatrix &x, const Index iv=0, const Index iz=0, const Index ia=0)
Set the At Position object.
Definition: propagationmatrix.h:392