ARTS  2.4.0(git:4fb77825)
m_ppath.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
29 /*===========================================================================
30  === External declarations
31  ===========================================================================*/
32 
33 #include <cmath>
34 #include "arts.h"
35 #include "auto_md.h"
36 #include "check_input.h"
37 #include "geodetic.h"
38 #include "lin_alg.h"
39 #include "m_general.h"
40 #include "m_xml.h"
41 #include "math_funcs.h"
42 #include "messages.h"
43 #include "ppath.h"
44 #include "refraction.h"
45 #include "special_interp.h"
46 #include "xml_io.h"
47 
48 extern const Numeric RAD2DEG;
49 extern const Numeric DEG2RAD;
50 
51 /*===========================================================================
52  === The functions (in alphabetical order)
53  ===========================================================================*/
54 
55 /* Workspace method: Doxygen documentation will be auto-generated */
56 void AddZaAa(Matrix& new_los,
57  const Vector& ref_los,
58  const Matrix& dlos,
59  const Verbosity&) {
60  if (ref_los.nelem() != 2)
61  throw runtime_error("*ref_los* must have two columns.");
62  if (dlos.ncols() != 2) throw runtime_error("*dlos* must have two columns.");
63 
64  const Index nlos = dlos.nrows();
65 
66  new_los.resize(nlos, 2);
67 
68  for (Index i = 0; i < nlos; i++) {
69  add_za_aa(new_los(i, 0),
70  new_los(i, 1),
71  ref_los[0],
72  ref_los[1],
73  dlos(i, 0),
74  dlos(i, 1));
75  }
76 }
77 
78 /* Workspace method: Doxygen documentation will be auto-generated */
79 void DiffZaAa(Matrix& dlos,
80  const Vector& ref_los,
81  const Matrix& other_los,
82  const Verbosity&) {
83  if (ref_los.nelem() != 2)
84  throw runtime_error("*ref_los* must have two columns.");
85  if (other_los.ncols() != 2)
86  throw runtime_error("*other_los* must have two columns.");
87 
88  const Index nlos = other_los.nrows();
89 
90  dlos.resize(nlos, 2);
91 
92  for (Index i = 0; i < nlos; i++) {
93  diff_za_aa(dlos(i, 0),
94  dlos(i, 1),
95  ref_los[0],
96  ref_los[1],
97  other_los(i, 0),
98  other_los(i, 1));
99  }
100 }
101 
102 /* Workspace method: Doxygen documentation will be auto-generated */
104  const Ppath& ppath,
105  const Verbosity& verbosity) {
106  geo_pos.resize(5);
107  geo_pos = NAN;
108 
109  geo_pos[Range(0, ppath.pos.ncols())] =
110  ppath.pos(ppath.np - 1, Range(0, ppath.pos.ncols()));
111  geo_pos[Range(3, ppath.los.ncols())] =
112  ppath.los(ppath.np - 1, Range(0, ppath.los.ncols()));
113 
114  CREATE_OUT2;
115  out2 << " Sets geo-position to:\n" << geo_pos;
116 }
117 
118 /* Workspace method: Doxygen documentation will be auto-generated */
120  const Ppath& ppath,
121  const Verbosity& verbosity) {
122  geo_pos.resize(5);
123  geo_pos = NAN;
124 
125  // Take first point of ppath as first guess
126  geo_pos[Range(0, ppath.pos.ncols())] =
127  ppath.pos(0, Range(0, ppath.pos.ncols()));
128  geo_pos[Range(3, ppath.los.ncols())] =
129  ppath.los(0, Range(0, ppath.los.ncols()));
130 
131  for (Index i = 1; i < ppath.np; i++) {
132  if (ppath.pos(i, 0) < geo_pos[0]) {
133  geo_pos[Range(0, ppath.pos.ncols())] =
134  ppath.pos(i, Range(0, ppath.pos.ncols()));
135  geo_pos[Range(3, ppath.los.ncols())] =
136  ppath.los(i, Range(0, ppath.los.ncols()));
137  }
138  }
139 
140  CREATE_OUT2;
141  out2 << " Sets geo-position to:\n" << geo_pos;
142 }
143 
144 /* Workspace method: Doxygen documentation will be auto-generated */
146  const Ppath& ppath,
147  const Numeric& z_ref,
148  const Verbosity& verbosity) {
149  geo_pos.resize(5);
150  geo_pos = NAN;
151 
152  bool found = false;
153  Index ihit = 0;
154  bool above = false;
155 
156  if (ppath.pos(0, 0) >= z_ref) {
157  above = true;
158  }
159 
160  while (!found && ihit < ppath.np - 1) {
161  ihit += 1;
162  if (above && ppath.pos(ihit, 0) < z_ref) {
163  found = true;
164  } else if (!above && ppath.pos(ihit, 0) >= z_ref) {
165  found = true;
166  }
167  }
168 
169  if (found) {
170  geo_pos[0] = z_ref;
171 
172  // Make a simple linear interpolation to determine lat and lon
173  const Numeric w = (z_ref - ppath.pos(ihit - 1, 0)) /
174  (ppath.pos(ihit, 0) - ppath.pos(ihit - 1, 0));
175 
176  geo_pos[3] = w * ppath.los(ihit, 0) + (1 - w) * ppath.los(ihit - 1, 0);
177 
178  if (ppath.pos.ncols() > 1) {
179  geo_pos[1] = w * ppath.pos(ihit, 1) + (1 - w) * ppath.pos(ihit - 1, 1);
180 
181  if (ppath.pos.ncols() > 2) {
182  geo_pos[2] = w * ppath.pos(ihit, 2) + (1 - w) * ppath.pos(ihit - 1, 2);
183  geo_pos[4] = w * ppath.los(ihit, 1) + (1 - w) * ppath.los(ihit - 1, 1);
184  }
185  }
186  }
187 
188  CREATE_OUT2;
189  out2 << " Sets geo-position to:\n" << geo_pos;
190 }
191 
192 /* Workspace method: Doxygen documentation will be auto-generated */
194  Ppath& ppath,
195  const Agenda& ppath_agenda,
196  const Numeric& ppath_lmax,
197  const Numeric& ppath_lraytrace,
198  const Index& atmgeom_checked,
199  const Vector& f_grid,
200  const Index& cloudbox_on,
201  const Index& cloudbox_checked,
203  const Vector& rte_pos,
204  const Vector& rte_los,
205  const Vector& rte_pos2,
206  const Verbosity&) {
207  // Basics
208  //
209  if (atmgeom_checked != 1)
210  throw runtime_error(
211  "The atmospheric geometry must be flagged to have "
212  "passed a consistency check (atmgeom_checked=1).");
213  if (cloudbox_checked != 1)
214  throw runtime_error(
215  "The cloudbox must be flagged to have "
216  "passed a consistency check (cloudbox_checked=1).");
217 
219  ppath,
220  ppath_lmax,
222  rte_pos,
223  rte_los,
224  rte_pos2,
225  cloudbox_on,
227  f_grid,
228  ppath_agenda);
229 }
230 
231 /* Workspace method: Doxygen documentation will be auto-generated */
233  Ppath& ppath,
234  const Agenda& ppath_agenda,
235  const Numeric& ppath_lmax,
236  const Numeric& ppath_lraytrace,
237  const Index& atmgeom_checked,
238  const Vector& f_grid,
239  const Index& cloudbox_on,
240  const Index& cloudbox_checked,
242  const Vector& rte_pos,
243  const Vector& rte_los,
244  const Vector& rte_pos2,
245  const Numeric& altitude,
246  const Numeric& accuracy,
247  const Verbosity& verbosity) {
248  ppathCalc(ws,
249  ppath,
250  ppath_agenda,
251  ppath_lmax,
254  f_grid,
255  cloudbox_on,
258  rte_pos,
259  rte_los,
260  rte_pos2,
261  verbosity);
262 
263  // Iterate until converging at altitude of interest
264  Index pos = first_pos_before_altitude(ppath, altitude);
265  Numeric lmax = ppath_lmax;
266  while (true) {
267  lmax *= 0.5;
268 
269  ppathCalc(ws,
270  ppath,
271  ppath_agenda,
272  lmax,
275  f_grid,
276  cloudbox_on,
279  ppath.dim == 1 ? ppath.pos(pos, 0) : ppath.pos(pos, joker),
280  ppath.dim == 1 ? ppath.los(pos, 0) : ppath.los(pos, joker),
281  rte_pos2,
282  verbosity);
283  pos = first_pos_before_altitude(ppath, altitude);
284 
285  if (std::abs(ppath.pos(pos, 0) - altitude) < accuracy) {
286  ppathCalc(ws,
287  ppath,
288  ppath_agenda,
289  ppath_lmax,
292  f_grid,
293  cloudbox_on,
296  ppath.dim == 1 ? ppath.pos(pos, 0) : ppath.pos(pos, joker),
297  ppath.dim == 1 ? ppath.los(pos, 0) : ppath.los(pos, joker),
298  rte_pos2,
299  verbosity);
300  break;
301  }
302  }
303 }
304 
305 /* Workspace method: Doxygen documentation will be auto-generated */
307  Ppath& ppath,
308  Vector& rte_los,
310  const Agenda& ppath_step_agenda,
311  const Index& atmosphere_dim,
312  const Vector& p_grid,
313  const Vector& lat_grid,
314  const Vector& lon_grid,
315  const Tensor3& z_field,
316  const Vector& f_grid,
317  const Vector& refellipsoid,
318  const Matrix& z_surface,
319  const Vector& rte_pos,
320  const Vector& rte_pos2,
321  const Numeric& ppath_lmax,
322  const Numeric& za_accuracy,
323  const Numeric& pplrt_factor,
324  const Numeric& pplrt_lowest,
325  const Verbosity& verbosity) {
326  //--- Check input -----------------------------------------------------------
327  if (atmosphere_dim == 2)
328  throw runtime_error(
329  "2D atmospheres not yet handled. Support for negative"
330  " zenith angles needed. Remind me (Patrick) to fix this.");
331  //---------------------------------------------------------------------------
332 
333  // Geometric LOS from rte_pos to rte_pos2
334  Vector rte_los_geom;
337  lat_grid,
338  lon_grid,
339  refellipsoid,
340  rte_pos,
341  rte_pos2,
342  verbosity);
343 
344  // Radius of rte_pos and rte_pos2
345  const Numeric r1 =
347  rte_pos[0];
348  const Numeric r2 =
350  rte_pos2[0];
351 
352  // Geometric distance between rte_pos and rte_pos2, effective 2D-lat for
353  // rte_pos and and Cartesian coordinates of rte_pos:
354  Numeric l12, lat1 = 0, x1, y1 = 0, z1;
355  if (atmosphere_dim <= 2) {
356  if (atmosphere_dim == 2) {
357  lat1 = rte_pos[1];
358  }
359  distance2D(l12, r1, lat1, r2, rte_pos2[1]);
360  pol2cart(x1, z1, r1, lat1);
361  } else {
362  distance3D(l12, r1, rte_pos[1], rte_pos[2], r2, rte_pos2[1], rte_pos2[2]);
363  sph2cart(x1, y1, z1, r1, rte_pos[1], rte_pos[2]);
364  }
365 
366  // Define remaining variables used in the while-loop below
367  //
368  // Basic bookkeeping variables
369  Numeric za_upp_limit = 180;
370  Numeric za_low_limit = 0;
371  //
372  // Various variables associated with the ppath, and the point of the path
373  // closest to the transmitter
374  Ppath ppt; // "Test ppath"
375  Index ip = -999; // Index of closest ppath point
376  Numeric xip, yip = 0, zip; // Cartesian coords. of the closest ppath point
377  Numeric dxip, dyip = 0, dzip; // Cartesian LOS of the closest ppath point
378  //
379  // Data for the intersection of the l12-sphere
380  Vector posc(max(Index(2), atmosphere_dim));
381  Numeric rc, xc, yc = 0, zc;
382 
383  CREATE_OUT2;
384  CREATE_OUT3;
385 
386  const Index maxiter = 99;
387  Vector t_za(maxiter, -999), t_dza(maxiter, -999);
388  Index it = -1;
389 
390  // Keep trying until ready or ground intersetion determined
391  //
392  bool ground = false;
393  bool failed = false;
394  Index ntries = 0;
395  //
396  while (true) {
397  // Path for present rte_los (no cloudbox!)
398  ppath_calc(ws,
399  ppt,
402  p_grid,
403  lat_grid,
404  lon_grid,
405  z_field,
406  f_grid,
407  refellipsoid,
408  z_surface,
409  0,
410  ArrayOfIndex(0),
411  rte_pos,
412  rte_los,
413  ppath_lmax,
415  0,
416  verbosity);
417 
418  // Find the point closest to rte_pos2, on the side towards rte_pos.
419  // We do this by looking at the distance to rte_pos, that should be
420  // as close to l12 as possible, but not exceed it.
421  Numeric lip = 99e99;
422  ip = ppt.np;
423  //
424  while (lip >= l12 && ip > 0) {
425  ip--;
426  if (atmosphere_dim <= 2) {
427  distance2D(lip, r1, lat1, ppt.r[ip], ppt.pos(ip, 1));
428  } else {
429  distance3D(lip,
430  r1,
431  rte_pos[1],
432  rte_pos[2],
433  ppt.r[ip],
434  ppt.pos(ip, 1),
435  ppt.pos(ip, 2));
436  }
437  }
438 
439  Numeric za_new, daa = 0;
440 
441  // Surface intersection:
442  // Not OK if the ground position is too far from rte_pos2.
443  // (30 km selected to allow misses of smaller size when rte_pos2 is at
444  // surface level, but surface interference never OK if rte_pos above TOA)
445  if (ppath_what_background(ppt) == 2 && ip == ppt.np - 1 &&
446  l12 - lip > 30e3) {
447  za_new = rte_los[0] - 1;
448  za_upp_limit = rte_los[0];
449  }
450 
451  // Ppath OK
452  else {
453  // Estimate ppath at the distance of l12, and calculate size
454  // of "miss" (measured in diffference in geometric angles)
455  Vector los;
456  Numeric dza;
457  if (atmosphere_dim <= 2) {
458  // Convert pos and los for point ip to cartesian coordinates
459  poslos2cart(
460  xip, zip, dxip, dzip, ppt.r[ip], ppt.pos(ip, 1), ppt.los(ip, 0));
461  // Find where the extension from point ip crosses the l12
462  // sphere: point c
463  Numeric latc;
464  line_circle_intersect(xc, zc, xip, zip, dxip, dzip, x1, z1, l12);
465  cart2pol(rc, latc, xc, zc, ppt.pos(ip, 1), ppt.los(ip, 0));
466  posc[1] = latc;
467  posc[0] =
468  rc - pos2refell_r(
470  } else {
471  // Convert pos and los for point ip to cartesian coordinates
472  poslos2cart(xip,
473  yip,
474  zip,
475  dxip,
476  dyip,
477  dzip,
478  ppt.r[ip],
479  ppt.pos(ip, 1),
480  ppt.pos(ip, 2),
481  ppt.los(ip, 0),
482  ppt.los(ip, 1));
483  // Find where the extension from point ip crosses the l12
484  // sphere: point c
485  Numeric latc, lonc;
487  xc, yc, zc, xip, yip, zip, dxip, dyip, dzip, x1, y1, z1, l12);
488  cart2sph(rc,
489  latc,
490  lonc,
491  xc,
492  yc,
493  zc,
494  ppt.pos(ip, 1),
495  ppt.pos(ip, 2),
496  ppt.los(ip, 0),
497  ppt.los(ip, 1));
498  posc[1] = latc;
499  posc[2] = lonc;
500  posc[0] =
501  rc - pos2refell_r(
503  }
504  //
507  lat_grid,
508  lon_grid,
509  refellipsoid,
510  rte_pos,
511  posc,
512  verbosity);
513  //
514  dza = los[0] - rte_los_geom[0];
515 
516  // Update bookkeeping variables
517  it++;
518  t_za[it] = rte_los[0];
519  t_dza[it] = dza;
520  //
521  if (dza > 0 && rte_los[0] < za_upp_limit) {
522  za_upp_limit = rte_los[0];
523  } else if (dza < 0 && rte_los[0] > za_low_limit) {
524  za_low_limit = rte_los[0];
525  }
526 
527  // Ready ?
528  if (abs(dza) <= za_accuracy) {
529  break;
530  } else if (za_upp_limit - za_low_limit <= za_accuracy / 10) {
531  if (max(t_dza) < -10 * za_accuracy) {
532  ground = true;
533  out3 << " Ground intersection determined !!!\n";
534  break;
535  } else {
536  failed = true;
537  out3 << " Zenith angle search range closed !!!\n";
538  break;
539  }
540  }
541  // Catch non-convergence (just for extra safety, za-range should be
542  // closed quicker than this)
543  ntries += 1;
544  if (ntries >= maxiter) {
545  failed = true;
546  out3 << " Too many iterations !!!\n";
547  break;
548  }
549 
550  // Estimate new angle
551  if (it < 1) {
552  za_new = rte_los[0] - dza;
553  } else {
554  // Estimate new angle by linear regression over some of the
555  // last calculations
556  const Index nfit = min(it + 1, (Index)3);
557  const Index i0 = it - nfit + 1;
558  Vector p;
559  linreg(p, t_za[Range(i0, nfit)], t_dza[Range(i0, nfit)]);
560  za_new = -p[0] / p[1];
561  }
562  //
563  if (atmosphere_dim == 3) {
564  daa = los[1] - rte_los_geom[1];
565  }
566  }
567 
568  // Update rte_los. Use bisection of za_new is basically
569  // identical to old angle, or is outside lower or upper
570  // limit. Otherwise use reult of linear reg.
571  if (std::isinf(za_new) || std::isnan(za_new) ||
572  abs(za_new - rte_los[0]) < 0.99 * za_accuracy ||
573  za_new <= za_low_limit || za_new >= za_upp_limit) {
574  rte_los[0] = (za_low_limit + za_upp_limit) / 2;
575  } else {
576  rte_los[0] = za_new;
577  if (atmosphere_dim == 3) {
578  rte_los[1] -= daa;
579  if (rte_los[1] < -180) {
580  rte_los[1] += 360;
581  } else if (rte_los[1] > 180) {
582  rte_los[1] -= 360;
583  }
584  }
585  }
586  } // while
587  //--------------------------------------------------------------------------
588 
589  // If failed re-try with a shorter ppath_lraytrace, if not ending up with
590  // a too small value.
591  if (failed) {
592  ppath_lraytrace /= pplrt_factor;
593 
594  if (ppath_lraytrace >= pplrt_lowest) {
595  out2 << " Re-start with ppath_lraytrace = " << ppath_lraytrace;
596  ppathFromRtePos2(ws,
597  ppath,
598  rte_los,
602  p_grid,
603  lat_grid,
604  lon_grid,
605  z_field,
606  f_grid,
607  refellipsoid,
608  z_surface,
609  rte_pos,
610  rte_pos2,
611  ppath_lmax,
612  za_accuracy,
613  pplrt_factor,
614  pplrt_lowest,
615  verbosity);
616  } else {
619  }
620  return; // --->
621  }
622 
623  // Create final ppath.
624  // If ground intersection: Set to length 1 and ground background,
625  // to flag non-OK path
626  // Otherwise: Fill path and set background to transmitter
627 
628  if (ground) {
631  }
632 
633  else {
634  // Distance between point ip of ppt and posc
635  Numeric ll;
636  if (atmosphere_dim <= 2) {
637  distance2D(ll, rc, posc[1], ppt.r[ip], ppt.pos(ip, 1));
638  } else {
639  distance3D(
640  ll, rc, posc[1], posc[2], ppt.r[ip], ppt.pos(ip, 1), ppt.pos(ip, 2));
641  }
642 
643  // Last point of ppt closest to rte_pos2. No point to add, maybe
644  // calculate start_lstep and start_los:
645  if (ip == ppt.np - 1) {
647  ppath_copy(ppath, ppt, -1);
648  if (ppath_what_background(ppath) == 1) {
649  ppath.start_lstep = ll;
650  Numeric d1, d2 = 0, d3;
651  if (atmosphere_dim <= 2) {
652  cart2poslos(d1,
653  d3,
654  ppath.start_los[0],
655  xc,
656  zc,
657  dxip,
658  dzip,
659  ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
660  ppt.pos(ip, 1),
661  ppt.los(ip, 0));
662  } else {
663  cart2poslos(d1,
664  d2,
665  d3,
666  ppath.start_los[0],
667  ppath.start_los[1],
668  xc,
669  yc,
670  zc,
671  dxip,
672  dyip,
673  dzip,
674  ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
675  xip,
676  yip,
677  zip, // Added 161027,PE
678  ppt.pos(ip, 1),
679  ppt.pos(ip, 2),
680  ppt.los(ip, 0),
681  ppt.los(ip, 1));
682  }
683  }
684  }
685  // rte_pos2 inside the atmosphere (posc entered as end point)
686  else {
688  ppath_copy(ppath, ppt, ip + 1);
689  //
690  const Index i = ip + 1;
691  if (atmosphere_dim <= 2) {
692  cart2poslos(ppath.r[i],
693  ppath.pos(i, 1),
694  ppath.los(i, 0),
695  xc,
696  zc,
697  dxip,
698  dzip,
699  ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
700  ppt.pos(ip, 1),
701  ppt.los(ip, 0));
702  } else {
703  cart2poslos(ppath.r[i],
704  ppath.pos(i, 1),
705  ppath.pos(i, 2),
706  ppath.los(i, 0),
707  ppath.los(i, 1),
708  xc,
709  yc,
710  zc,
711  dxip,
712  dyip,
713  dzip,
714  ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
715  xip,
716  yip,
717  zip, // Added 161027,PE
718  ppt.pos(ip, 1),
719  ppt.pos(ip, 2),
720  ppt.los(ip, 0),
721  ppt.los(ip, 1));
722  }
723  //
724  ppath.pos(i, joker) = posc;
725  ppath.lstep[i - 1] = ll;
726  ppath.start_los = ppath.los(i, joker);
727 
728  // n by linear interpolation
729  // Gets tripped when ll is very close to (slightly greater than) lstep (ISA)
730  assert(ll < ppt.lstep[i - 1]);
731  const Numeric w = ll / ppt.lstep[i - 1];
732  ppath.nreal[i] = (1 - w) * ppt.nreal[i - 1] + w * ppt.nreal[i];
733  ppath.ngroup[i] = (1 - w) * ppt.ngroup[i - 1] + w * ppt.ngroup[i];
734 
735  // Grid positions
736  GridPos gp_lat, gp_lon;
737  rte_pos2gridpos(ppath.gp_p[i],
738  gp_lat,
739  gp_lon,
741  p_grid,
742  lat_grid,
743  lon_grid,
744  z_field,
745  ppath.pos(i, Range(0, atmosphere_dim)));
746  if (atmosphere_dim >= 2) {
747  gridpos_copy(ppath.gp_lat[i], gp_lat);
748  if (atmosphere_dim == 3) {
749  gridpos_copy(ppath.gp_lon[i], gp_lon);
750  }
751  }
752  }
753 
754  // Common stuff
756  ppath.start_pos = rte_pos2;
757  }
758 }
759 
760 /* Workspace method: Doxygen documentation will be auto-generated */
762  const Index& atmosphere_dim,
763  const Tensor3& z_field,
764  const Matrix& z_surface,
765  const Index& cloudbox_on,
768  const Vector& rte_pos,
769  const Vector& rte_los,
770  const Numeric& ppath_lmax,
771  const Verbosity&) {
772  // This function is a WSM but it is normally only called from yCalc.
773  // For that reason, this function does not repeat input checks that are
774  // performed in yCalc, it only performs checks regarding the sensor
775  // position and LOS.
776 
777  const Numeric z_sensor = rte_pos[0];
778  const Numeric za_sensor = rte_los[0];
779  const Index nz = z_field.npages();
780  const Numeric z_toa = z_field(nz - 1, 0, 0);
781  const bool above_toa = z_sensor > z_toa ? true : false;
782  const Numeric z_end = above_toa ? z_toa : z_sensor;
783  const Numeric dz2dl = abs(1 / cos(DEG2RAD * za_sensor));
784  Index background = -99;
785 
786  // Basics checks of input
787  if (atmosphere_dim != 1)
788  throw runtime_error("The function can only be used for 1D atmospheres.");
792  throw runtime_error(
793  "The WSV *ppath_inside_cloudbox_do* can only be set "
794  "to 1 if also *cloudbox_on* is 1.");
795  if (z_sensor < z_surface(0, 0)) {
796  ostringstream os;
797  os << "The sensor is below the surface."
798  << " altitude of sensor : " << z_sensor << endl
799  << " altitude of surface : " << z_surface(0, 0);
800  throw runtime_error(os.str());
801  }
802  if (abs(za_sensor - 90) < 0.1) {
803  ostringstream os;
804  os << "The zenith angle is " << za_sensor << endl
805  << "The method does not allow this. The zenith angle must deviate\n"
806  << "from 90 deg with at least 1 deg. That is, to be outside [89.9,90.1].";
807  throw runtime_error(os.str());
808  }
809 
810  // Find end grid position
811  GridPos gp_end;
812  // To avoid compiler warnings, start to assuming above_toa
813  gp_end.idx = nz - 2;
814  gp_end.fd[0] = 1;
815  gp_end.fd[1] = 0;
816  if (!above_toa) {
817  for (Index i = 0; i < nz - 1; i++) {
818  if (z_sensor < z_field(i + 1, 0, 0)) {
819  gp_end.idx = i;
820  gp_end.fd[0] = (z_sensor - z_field(i, 0, 0)) /
821  (z_field(i + 1, 0, 0) - z_field(i, 0, 0));
822  gp_end.fd[1] = 1 - gp_end.fd[0];
823  break;
824  }
825  }
826  }
827 
828  // Catch cases resulting in a ppath with 1 point
829  bool path_to_follow = true;
830  if (above_toa && za_sensor < 90) {
831  // Path fully in space
833  background = 1;
834  path_to_follow = false;
835  } else if (z_sensor == z_surface(0, 0) && za_sensor > 90) {
836  // On ground, looking down
838  background = 2;
839  path_to_follow = false;
840  } else if (cloudbox_on) {
842  z_sensor > z_field(cloudbox_limits[0], 0, 0) &&
843  z_sensor < z_field(cloudbox_limits[1], 0, 0)) {
844  // Inside cloud box
846  background = 4;
847  path_to_follow = false;
848  } else if ((z_sensor == z_field(cloudbox_limits[0], 0, 0) &&
849  za_sensor > 90) ||
850  (z_sensor == z_field(cloudbox_limits[1], 0, 0) &&
851  za_sensor < 90)) {
852  // Cloud box boundary
854  background = 3;
855  path_to_follow = false;
856  } else if (above_toa && cloudbox_limits[1] == nz - 1) {
857  // Cloud box boundary is at TOA
859  background = 3;
860  path_to_follow = false;
861  }
862  }
863 
864  // Determine ppath
865  if (path_to_follow) {
866  const Numeric max_dz = ppath_lmax > 0 ? ppath_lmax / dz2dl : 9e99;
867 
868  // Variables to describe each "break-point" of ppath. Point 0 is the end
869  // point. Not all nz points are necessarily passed.
870  ArrayOfIndex l_idx(nz);
871  ArrayOfVector l_fd0(nz);
872  ArrayOfVector l_z(nz);
873  Index nptot = 0;
874 
875  // Determine number of ppath points in each layer
876  {
877  Numeric z = z_end;
878  Index iout = -1;
879 
880  // Code similar, but for simplicity, we handle down- and
881  // up-ward separately:
882  if (za_sensor > 90) // Downward-looking
883  {
884  // Here we go down to next pressure level (or the surface) in each
885  // step. That is, if above surface, last point of step has fd[0]=0.
886 
887  // Put in end point
888  iout++;
889  nptot++;
890  l_fd0[0].resize(1);
891  l_z[0].resize(1);
892  l_idx[0] = gp_end.idx;
893  l_fd0[0][0] = gp_end.fd[0];
894  l_z[0][0] = z_end;
895 
896  for (Index i = gp_end.idx; i >= 0 && background < 0; i--) {
897  // Surface inside layer?
898  Numeric dz_step;
899  if (z_field(i, 0, 0) > z_surface(0, 0)) {
900  dz_step = z - z_field(i, 0, 0);
901  } else {
902  dz_step = z - z_surface(0, 0);
903  background = 2;
904  }
905 
906  const Index np =
907  dz_step <= max_dz ? 1 : Index(ceil(dz_step / max_dz));
908  const Numeric dz = dz_step / Numeric(np);
909  const Numeric dz_layer = z_field(i + 1, 0, 0) - z_field(i, 0, 0);
910 
911  // Update counters and resize
912  iout++;
913  nptot += np;
914  l_fd0[iout].resize(np);
915  l_z[iout].resize(np);
916 
917  // Intermediate points
918  for (Index j = 0; j < np - 1; j++) {
919  l_z[iout][j] = z - (Numeric(j) + 1) * dz;
920  l_fd0[iout][j] = (l_z[iout][j] - z_field(i, 0, 0)) / dz_layer;
921  }
922 
923  // End points handled seperately to avoid numerical problems
924  l_idx[iout] = i;
925  if (background == 2) // Surface is reached
926  {
927  l_z[iout][np - 1] = z_surface(0, 0);
928  l_fd0[iout][np - 1] =
929  (l_z[iout][np - 1] - z_field(i, 0, 0)) / dz_layer;
930  } else {
931  l_z[iout][np - 1] = z_field(i, 0, 0);
932  l_fd0[iout][np - 1] = 0;
933  //
934  if (cloudbox_on &&
935  (i == cloudbox_limits[1] || i == cloudbox_limits[0])) {
936  background = 3;
937  }
938  }
939 
940  // Update z
941  z = z_field(i, 0, 0);
942  }
943  } else // Upward-looking
944  {
945  // Here we have that first point of step has fd[0]=0, if not at
946  // sensor
947  for (Index i = gp_end.idx; i < nz && background < 0; i++) {
948  Numeric dz_layer;
949  Numeric dz_step;
950  if (cloudbox_on && i != gp_end.idx &&
951  (i == cloudbox_limits[0] ||
952  i == cloudbox_limits[1])) { // At an active cloudbox boundary
953  dz_step = 0;
954  dz_layer = 1;
955  background = 3;
956  } else if (i == nz - 1) { // At TOA
957  dz_step = 0;
958  dz_layer = 1;
959  background = 1;
960  } else {
961  dz_step = z_field(i + 1, 0, 0) - z;
962  dz_layer = z_field(i + 1, 0, 0) - z_field(i, 0, 0);
963  }
964 
965  const Index np =
966  dz_step <= max_dz ? 1 : Index(ceil(dz_step / max_dz));
967  const Numeric dz = dz_step / Numeric(np);
968 
969  // Update counters and resize
970  iout++;
971  nptot += np;
972  l_fd0[iout].resize(np);
973  l_z[iout].resize(np);
974 
975  // Start points handled seperately to avoid numerical problems
976  if (i == gp_end.idx) { // At sensor
977  l_idx[iout] = i;
978  l_z[iout][0] = z_sensor;
979  l_fd0[iout][0] = gp_end.fd[0];
980  } else if (i == nz - 1) { // At TOA
981  l_idx[iout] = i - 1;
982  l_z[iout][0] = z_field(i, 0, 0);
983  l_fd0[iout][0] = 1;
984  } else {
985  l_idx[iout] = i;
986  l_z[iout][0] = z_field(i, 0, 0);
987  l_fd0[iout][0] = 0;
988  }
989 
990  // Intermediate points
991  for (Index j = 1; j < np; j++) {
992  l_z[iout][j] = z + Numeric(j) * dz;
993  l_fd0[iout][j] = (l_z[iout][j] - z_field(i, 0, 0)) / dz_layer;
994  }
995 
996  // Update z
997  if (background < 0) {
998  z = z_field(i + 1, 0, 0);
999  }
1000  }
1001  }
1002  }
1003 
1005 
1006  // Fill ppath.pos(joker,0), ppath.gp_p and ppath.lstep
1007  Index iout = -1;
1008  Numeric z_last = -999;
1009  for (Index i = 0; i < nz; i++) {
1010  for (Index j = 0; j < l_z[i].nelem(); j++) {
1011  iout++;
1012  ppath.pos(iout, 0) = l_z[i][j];
1013  ppath.gp_p[iout].idx = l_idx[i];
1014  ppath.gp_p[iout].fd[0] = l_fd0[i][j];
1015  ppath.gp_p[iout].fd[1] = 1 - l_fd0[i][j];
1016  if (iout == 0) {
1017  z_last = ppath.pos(iout, 0);
1018  } else {
1019  ppath.lstep[iout - 1] = dz2dl * abs(z_last - l_z[i][j]);
1020  z_last = l_z[i][j];
1021  }
1022  }
1023  }
1024  }
1025 
1026  // Remaining data
1027  ppath_set_background(ppath, background);
1028  if (ppath.np == 1) {
1029  ppath.pos(0, 0) = z_end;
1030  ppath.gp_p[0] = gp_end;
1031  }
1032  ppath.pos(joker, 1) = 0;
1033  ppath.los(joker, 0) = za_sensor;
1034  ppath.constant = INFINITY; // Not defined here as r = Inf
1035  ppath.r = INFINITY;
1036  ppath.start_pos[0] = ppath.pos(ppath.np - 1, 0);
1037  ppath.start_pos[1] = 0;
1038  ppath.start_los[0] = za_sensor;
1039  ppath.end_pos[0] = z_sensor;
1040  ppath.end_pos[1] = 0;
1041  ppath.end_los[0] = za_sensor;
1042  if (above_toa) {
1043  ppath.end_lstep = dz2dl * (z_sensor - z_toa);
1044  }
1045  ppath.nreal = 1;
1046  ppath.ngroup = 1;
1047 }
1048 
1049 /* Workspace method: Doxygen documentation will be auto-generated */
1051  Ppath& ppath,
1052  const Agenda& ppath_step_agenda,
1054  const Index& atmosphere_dim,
1055  const Vector& p_grid,
1056  const Vector& lat_grid,
1057  const Vector& lon_grid,
1058  const Tensor3& z_field,
1059  const Vector& f_grid,
1060  const Vector& refellipsoid,
1061  const Matrix& z_surface,
1062  const Index& cloudbox_on,
1064  const Vector& rte_pos,
1065  const Vector& rte_los,
1066  const Numeric& ppath_lmax,
1067  const Numeric& ppath_lraytrace,
1068  const Verbosity& verbosity) {
1069  ppath_calc(ws,
1070  ppath,
1073  p_grid,
1074  lat_grid,
1075  lon_grid,
1076  z_field,
1077  f_grid,
1078  refellipsoid,
1079  z_surface,
1080  cloudbox_on,
1082  rte_pos,
1083  rte_los,
1084  ppath_lmax,
1087  verbosity);
1088 }
1089 
1090 /* Workspace method: Doxygen documentation will be auto-generated */
1091 void ppathWriteXMLPartial( //WS Input:
1092  const String& file_format,
1093  const Ppath& ppath,
1094  // WS Generic Input:
1095  const String& f,
1096  const Index& file_index,
1097  const Verbosity& verbosity) {
1098  String filename = f;
1099  Ppath ppath_partial = ppath;
1100  ArrayOfGridPos empty_gp;
1101  //Vector empty_v;
1102 
1103  ppath_partial.gp_p = empty_gp;
1104  ppath_partial.gp_lat = empty_gp;
1105  ppath_partial.gp_lon = empty_gp;
1106  //ppath_partial.nreal = empty_v;
1107  //ppath_partial.ngroup = empty_v;
1108 
1109  if (file_index >= 0) {
1110  // Create default filename if empty
1111  filename_xml_with_index(filename, file_index, "ppath");
1112  }
1113 
1114  WriteXML(file_format, ppath_partial, filename, 0, "ppath", "", "", verbosity);
1115 }
1116 
1117 // FIXMEDOC@Richard TRy to describe the meaning of ppath_field
1118 
1119 /* Workspace method: Doxygen documentation will be auto-generated */
1122  const Agenda& ppath_agenda,
1123  const Numeric& ppath_lmax,
1124  const Numeric& ppath_lraytrace,
1125  const Index& atmgeom_checked,
1126  const Tensor3& z_field,
1127  const Vector& f_grid,
1128  const Index& cloudbox_on,
1129  const Index& cloudbox_checked,
1131  const Vector& rte_pos,
1132  const Vector& rte_los,
1133  const Vector& rte_pos2,
1134  const Vector& refellipsoid,
1135  const Index& atmosphere_dim,
1136  const Index& zenith_angles_per_position,
1137  const Verbosity& verbosity) {
1138  if (atmosphere_dim not_eq 1)
1139  throw std::runtime_error("Only for 1D atmospheres");
1140  if (refellipsoid[1] not_eq 0.0)
1141  throw std::runtime_error("Not allowed for non-spherical planets");
1142  if (ppath_lmax >= 0)
1143  throw std::runtime_error("Only allowed for long paths (ppath_lmax < 0)");
1144 
1145  // Positions and angles of interest
1146  const Numeric zmin = z_field(0, 0, 0);
1147  const Numeric zmax = z_field(z_field.npages() - 1, 0, 0);
1148  const Numeric r = refellipsoid[0];
1149  const Numeric above_surface_tangent =
1150  90 - RAD2DEG * std::acos((r) / (r + zmax)) + 1e-4;
1151  const Numeric below_surface_tangent =
1152  90 - RAD2DEG * std::acos((r) / (r + zmax)) - 1e-4;
1153  const Numeric top_tangent = 90 - 1e-4;
1154 
1155  ppath_field.resize(3 * zenith_angles_per_position);
1156  Index ppath_field_pos = 0;
1157 
1158  Vector zenith_angles(zenith_angles_per_position);
1159 
1160  // Upwards:
1161  nlinspace(zenith_angles, 0, 90, zenith_angles_per_position);
1162  Vector rte_pos_true = rte_pos;
1163  rte_pos_true[0] = zmin;
1164  Vector rte_los_true = rte_los;
1165  for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1166  rte_los_true[0] = zenith_angles[iz];
1167 
1168  ppathCalc(ws,
1169  ppath_field[ppath_field_pos],
1170  ppath_agenda,
1171  ppath_lmax,
1174  f_grid,
1175  cloudbox_on,
1178  rte_pos_true,
1179  rte_los_true,
1180  rte_pos2,
1181  verbosity);
1182 
1183  ppath_field_pos++;
1184  }
1185 
1186  // Limb:
1187  nlinspace(zenith_angles,
1188  above_surface_tangent,
1189  top_tangent,
1190  zenith_angles_per_position);
1191  rte_pos_true[0] = zmax;
1192  for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1193  rte_los_true[0] = 180 - zenith_angles[iz];
1194 
1195  ppathCalc(ws,
1196  ppath_field[ppath_field_pos],
1197  ppath_agenda,
1198  ppath_lmax,
1201  f_grid,
1202  cloudbox_on,
1205  rte_pos_true,
1206  rte_los_true,
1207  rte_pos2,
1208  verbosity);
1209 
1210  ppath_field_pos++;
1211  }
1212 
1213  // Downwards:
1214  nlinspace(
1215  zenith_angles, 0, below_surface_tangent, zenith_angles_per_position);
1216  for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1217  rte_los_true[0] = 180 - zenith_angles[iz];
1218 
1219  ppathCalc(ws,
1220  ppath_field[ppath_field_pos],
1221  ppath_agenda,
1222  ppath_lmax,
1225  f_grid,
1226  cloudbox_on,
1229  rte_pos_true,
1230  rte_los_true,
1231  rte_pos2,
1232  verbosity);
1233 
1234  ppath_field_pos++;
1235  }
1236 }
1237 
1238 /* Workspace method: Doxygen documentation will be auto-generated */
1241  const Agenda& ppath_agenda,
1242  const Numeric& ppath_lmax,
1243  const Numeric& ppath_lraytrace,
1244  const Index& atmgeom_checked,
1245  const Vector& f_grid,
1246  const Index& cloudbox_on,
1247  const Index& cloudbox_checked,
1249  const Matrix& sensor_pos,
1250  const Matrix& sensor_los,
1251  const Vector& rte_pos2,
1252  const Verbosity& verbosity) {
1253  auto n = sensor_pos.nrows();
1254  ppath_field.resize(n);
1255 
1256  if (sensor_los.nrows() not_eq n)
1257  throw std::runtime_error(
1258  "Your sensor position matrix and sensor line of sight matrix do not match in size.\n");
1259 
1260  for (auto i = 0; i < n; i++)
1261  ppathCalc(ws,
1262  ppath_field[i],
1263  ppath_agenda,
1264  ppath_lmax,
1267  f_grid,
1268  cloudbox_on,
1271  sensor_pos(i, joker),
1272  sensor_los(i, joker),
1273  rte_pos2,
1274  verbosity);
1275 }
1276 
1277 /* Workspace method: Doxygen documentation will be auto-generated */
1278 void ppath_stepGeometric( // WS Output:
1279  Ppath& ppath_step,
1280  // WS Input:
1281  const Index& atmosphere_dim,
1282  const Vector& lat_grid,
1283  const Vector& lon_grid,
1284  const Tensor3& z_field,
1285  const Vector& refellipsoid,
1286  const Matrix& z_surface,
1287  const Numeric& ppath_lmax,
1288  const Verbosity&) {
1289  // Input checks here would be rather costly as this function is called
1290  // many times. So we perform asserts in the sub-functions, but no checks
1291  // here.
1292 
1293  // A call with background set, just wants to obtain the refractive index for
1294  // complete ppaths consistent of a single point.
1296  if (atmosphere_dim == 1) {
1298  z_field(joker, 0, 0),
1299  refellipsoid,
1300  z_surface(0, 0),
1301  ppath_lmax);
1302  }
1303 
1304  else if (atmosphere_dim == 2) {
1306  lat_grid,
1307  z_field(joker, joker, 0),
1308  refellipsoid,
1309  z_surface(joker, 0),
1310  ppath_lmax);
1311  }
1312 
1313  else if (atmosphere_dim == 3) {
1315  lat_grid,
1316  lon_grid,
1317  z_field,
1318  refellipsoid,
1319  z_surface,
1320  ppath_lmax);
1321  }
1322 
1323  else {
1324  throw runtime_error("The atmospheric dimensionality must be 1-3.");
1325  }
1326  }
1327 
1328  else {
1329  assert(ppath_step.np == 1);
1330  ppath_step.nreal[0] = 1;
1331  ppath_step.ngroup[0] = 1;
1332  }
1333 }
1334 
1335 /* Workspace method: Doxygen documentation will be auto-generated */
1337  Ppath& ppath_step,
1339  const Index& atmosphere_dim,
1340  const Vector& p_grid,
1341  const Vector& lat_grid,
1342  const Vector& lon_grid,
1343  const Tensor3& z_field,
1344  const Tensor3& t_field,
1345  const Tensor4& vmr_field,
1346  const Vector& refellipsoid,
1347  const Matrix& z_surface,
1348  const Vector& f_grid,
1349  const Numeric& ppath_lmax,
1350  const Numeric& ppath_lraytrace,
1351  const Verbosity&) {
1352  // Input checks here would be rather costly as this function is called
1353  // many times.
1354  assert(ppath_lraytrace > 0);
1355 
1356  // A call with background set, just wants to obtain the refractive index for
1357  // complete ppaths consistent of a single point.
1359  if (atmosphere_dim == 1) {
1360  ppath_step_refr_1d(ws,
1361  ppath_step,
1362  p_grid,
1363  z_field,
1364  t_field,
1365  vmr_field,
1366  f_grid,
1367  refellipsoid,
1368  z_surface(0, 0),
1369  ppath_lmax,
1371  "linear_basic",
1372  ppath_lraytrace);
1373  } else if (atmosphere_dim == 2) {
1374  ppath_step_refr_2d(ws,
1375  ppath_step,
1376  p_grid,
1377  lat_grid,
1378  z_field,
1379  t_field,
1380  vmr_field,
1381  f_grid,
1382  refellipsoid,
1383  z_surface(joker, 0),
1384  ppath_lmax,
1386  "linear_basic",
1387  ppath_lraytrace);
1388  } else if (atmosphere_dim == 3) {
1389  ppath_step_refr_3d(ws,
1390  ppath_step,
1391  p_grid,
1392  lat_grid,
1393  lon_grid,
1394  z_field,
1395  t_field,
1396  vmr_field,
1397  f_grid,
1398  refellipsoid,
1399  z_surface,
1400  ppath_lmax,
1402  "linear_basic",
1403  ppath_lraytrace);
1404  } else {
1405  throw runtime_error("The atmospheric dimensionality must be 1-3.");
1406  }
1407  }
1408 
1409  else {
1410  assert(ppath_step.np == 1);
1411  if (atmosphere_dim == 1) {
1412  get_refr_index_1d(ws,
1413  ppath_step.nreal[0],
1414  ppath_step.ngroup[0],
1416  p_grid,
1417  refellipsoid,
1418  z_field,
1419  t_field,
1420  vmr_field,
1421  f_grid,
1422  ppath_step.r[0]);
1423  } else if (atmosphere_dim == 2) {
1424  get_refr_index_2d(ws,
1425  ppath_step.nreal[0],
1426  ppath_step.ngroup[0],
1428  p_grid,
1429  lat_grid,
1430  refellipsoid,
1431  z_field,
1432  t_field,
1433  vmr_field,
1434  f_grid,
1435  ppath_step.r[0],
1436  ppath_step.pos(0, 1));
1437  } else {
1438  get_refr_index_3d(ws,
1439  ppath_step.nreal[0],
1440  ppath_step.ngroup[0],
1442  p_grid,
1443  lat_grid,
1444  lon_grid,
1445  refellipsoid,
1446  z_field,
1447  t_field,
1448  vmr_field,
1449  f_grid,
1450  ppath_step.r[0],
1451  ppath_step.pos(0, 1),
1452  ppath_step.pos(0, 2));
1453  }
1454  }
1455 }
1456 
1457 /* Workspace method: Doxygen documentation will be auto-generated */
1459  const Index& atmosphere_dim,
1460  const Numeric& za,
1461  const Numeric& aa,
1462  const Verbosity&) {
1463  // Check input
1464  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1465 
1466  if (atmosphere_dim == 1) {
1467  rte_los.resize(1);
1468  } else {
1469  rte_los.resize(2);
1470  rte_los[1] = aa;
1471  }
1472  rte_los[0] = za;
1473 }
1474 
1475 /* Workspace method: Doxygen documentation will be auto-generated */
1477  const Index& atmosphere_dim,
1478  const Vector& lat_grid,
1479  const Vector& lon_grid,
1480  const Vector& refellipsoid,
1481  const Vector& rte_pos,
1482  const Vector& rte_pos2,
1483  const Verbosity&) {
1484  // Check input
1487 
1488  // Radius of rte_pos and rte_pos2
1489  const Numeric r1 =
1491  rte_pos[0];
1492  const Numeric r2 =
1494  rte_pos2[0];
1495 
1496  // Remaining polar and cartesian coordinates of rte_pos
1497  Numeric lat1, lon1 = 0, x1, y1 = 0, z1;
1498  // Cartesian coordinates of rte_pos2
1499  Numeric x2, y2 = 0, z2;
1500  //
1501  if (atmosphere_dim == 1) {
1502  // Latitude distance implicitly checked by chk_rte_pos
1503  lat1 = 0;
1504  pol2cart(x1, z1, r1, lat1);
1505  pol2cart(x2, z2, r2, rte_pos2[1]);
1506  } else if (atmosphere_dim == 2) {
1507  lat1 = rte_pos[1];
1508  pol2cart(x1, z1, r1, lat1);
1509  pol2cart(x2, z2, r2, rte_pos2[1]);
1510  } else {
1511  lat1 = rte_pos[1];
1512  lon1 = rte_pos[2];
1513  sph2cart(x1, y1, z1, r1, lat1, lon1);
1514  sph2cart(x2, y2, z2, r2, rte_pos2[1], rte_pos2[2]);
1515  }
1516 
1517  // Geometrical LOS to transmitter
1518  Numeric za, aa;
1519  //
1520  los2xyz(za, aa, r1, lat1, lon1, x1, y1, z1, x2, y2, z2);
1521  //
1522  if (atmosphere_dim == 3) {
1523  rte_los.resize(2);
1524  rte_los[0] = za;
1525  rte_los[1] = aa;
1526  } else {
1527  rte_los.resize(1);
1528  rte_los[0] = za;
1529  if (atmosphere_dim == 2 && aa < 0) // Should 2D-za be negative?
1530  {
1531  rte_los[0] = -za;
1532  }
1533  }
1534 }
1535 
1536 /* Workspace method: Doxygen documentation will be auto-generated */
1538  const Index& atmosphere_dim,
1539  const Numeric& z,
1540  const Numeric& lat,
1541  const Numeric& lon,
1542  const Verbosity&) {
1543  // Check input
1544  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1545 
1546  rte_pos.resize(atmosphere_dim);
1547  rte_pos[0] = z;
1548  if (atmosphere_dim >= 2) {
1549  rte_pos[1] = lat;
1550  }
1551  if (atmosphere_dim == 3) {
1552  rte_pos[2] = lon;
1553  }
1554 }
1555 
1556 /* Workspace method: Doxygen documentation will be auto-generated */
1558  Vector& rte_los,
1559  const Index& atmosphere_dim,
1560  const Ppath& ppath,
1561  const Verbosity&) {
1562  const Index np = ppath.np;
1563 
1564  // Check input
1565  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1566  if (np == 0) throw runtime_error("The input *ppath* is empty.");
1567  if (ppath.pos.nrows() != np)
1568  throw runtime_error(
1569  "Internal inconsistency in *ppath* (size of data "
1570  "does not match np).");
1571 
1572  rte_pos = ppath.pos(np - 1, Range(0, atmosphere_dim));
1573  if (atmosphere_dim < 3) {
1574  rte_los = ppath.los(np - 1, Range(0, 1));
1575  } else {
1576  rte_los = ppath.los(np - 1, Range(0, 2));
1577  }
1578 }
1579 
1580 /* Workspace method: Doxygen documentation will be auto-generated */
1582  Matrix& sensor_los,
1583  const Index& atmosphere_dim,
1584  const Vector& lat_grid,
1585  const Vector& lon_grid,
1586  const Vector& refellipsoid,
1587  const Matrix& sensor_pos,
1588  const Matrix& target_pos,
1589  const Verbosity& verbosity) {
1590  const Index n = sensor_pos.nrows();
1591 
1592  if (sensor_pos.ncols() != atmosphere_dim)
1593  throw runtime_error(
1594  "The number of columns of sensor_pos must be "
1595  "equal to the atmospheric dimensionality.");
1596  if ((atmosphere_dim == 1 && target_pos.ncols() != 2) ||
1597  (atmosphere_dim >= 2 && target_pos.ncols() != atmosphere_dim))
1598  throw runtime_error(
1599  "The number of columns of targe_pos must be equal to "
1600  "the atmospheric dimensionality, except for 1D where "
1601  "two columns are demended (as for *rte_pos2*).");
1602  if (target_pos.nrows() != n)
1603  throw runtime_error(
1604  "*sensor_pos* and *target_pos* must have the same "
1605  "number of rows.");
1606 
1607  atmosphere_dim < 3 ? sensor_los.resize(n, 1) : sensor_los.resize(n, 2);
1608  Vector rte_los;
1609  for (Index i = 0; i < n; i++) {
1612  lat_grid,
1613  lon_grid,
1614  refellipsoid,
1615  sensor_pos(i, joker),
1616  target_pos(i, joker),
1617  verbosity);
1618  sensor_los(i, joker) = rte_los;
1619  }
1620 }
1621 
1622 /* Workspace method: Doxygen documentation will be auto-generated */
1624  const Ppath& ppath,
1625  const Verbosity&) {
1626  Index it;
1627  find_tanpoint(it, ppath);
1628 
1629  tan_pos.resize(ppath.pos.ncols());
1630 
1631  if (it < 0) {
1632  tan_pos = std::numeric_limits<Numeric>::quiet_NaN();
1633  } else {
1634  tan_pos[0] = ppath.pos(it, 0);
1635  tan_pos[1] = ppath.pos(it, 1);
1636  if (ppath.pos.ncols() == 3) {
1637  tan_pos[2] = ppath.pos(it, 2);
1638  }
1639  }
1640 }
1641 
1642 /* Workspace method: Doxygen documentation will be auto-generated */
1644  const Index& level,
1645  const Verbosity& verbosity) {
1646  Index it;
1647  find_tanpoint(it, ppath);
1648 
1649  ostringstream os;
1650 
1651  if (it < 0) {
1652  os << "Lowest altitude found at the end of the propagation path.\n"
1653  << "This indicates that the tangent point is either above the\n"
1654  << "top-of-the-atmosphere or below the planet's surface.";
1655  } else {
1656  os << "Tangent point position:\n-----------------------\n"
1657  << " z = " << ppath.pos(it, 0) / 1e3 << " km\n"
1658  << " lat = " << ppath.pos(it, 1) << " deg";
1659  if (ppath.pos.ncols() == 3)
1660  os << "\n lon: " << ppath.pos(it, 2) << " deg";
1661  }
1662 
1663  CREATE_OUTS;
1664  SWITCH_OUTPUT(level, os.str());
1665 }
1666 
1667 /* Workspace method: Doxygen documentation will be auto-generated */
1669  Vector& za_vector,
1671  const Matrix& sensor_pos,
1672  const Vector& p_grid,
1673  const Tensor3& t_field,
1674  const Tensor3& z_field,
1675  const Tensor4& vmr_field,
1676  const Vector& refellipsoid,
1677  const Index& atmosphere_dim,
1678  const Vector& f_grid,
1679  const Vector& ztan_vector,
1680  const Verbosity&) {
1681  if (atmosphere_dim != 1) {
1682  throw runtime_error("The function can only be used for 1D atmospheres.");
1683  }
1684 
1685  if (ztan_vector.nelem() != sensor_pos.nrows()) {
1686  ostringstream os;
1687  os << "The number of altitudes in true tangent altitude vector must\n"
1688  << "match the number of positions in *sensor_pos*.";
1689  throw runtime_error(os.str());
1690  }
1691 
1692  // Set za_vector's size equal to ztan_vector
1693  za_vector.resize(ztan_vector.nelem());
1694 
1695  // Define refraction variables
1697 
1698  // Calculate refractive index for the tangential altitudes
1699  for (Index i = 0; i < ztan_vector.nelem(); i++) {
1700  if (ztan_vector[i] > sensor_pos(i, 0)) {
1701  ostringstream os;
1702  os << "Invalid observation geometry: sensor (at z=" << sensor_pos(i, 0)
1703  << "m) is located below the requested tangent altitude (tanh="
1704  << ztan_vector[i] << "m)";
1705  throw runtime_error(os.str());
1706  } else {
1707  get_refr_index_1d(ws,
1711  p_grid,
1712  refellipsoid[0],
1713  z_field,
1714  t_field,
1715  vmr_field,
1716  f_grid,
1717  ztan_vector[i] + refellipsoid[0]);
1718 
1719  // Calculate zenith angle
1720  za_vector[i] = 180 - RAD2DEG * asin(refr_index_air *
1721  (refellipsoid[0] + ztan_vector[i]) /
1722  (refellipsoid[0] + sensor_pos(i, 0)));
1723  }
1724  }
1725 }
1726 
1727 /* Workspace method: Doxygen documentation will be auto-generated */
1728 void VectorZtanToZa1D(Vector& za_vector,
1729  const Matrix& sensor_pos,
1730  const Vector& refellipsoid,
1731  const Index& atmosphere_dim,
1732  const Vector& ztan_vector,
1733  const Verbosity&) {
1734  if (atmosphere_dim != 1) {
1735  throw runtime_error("The function can only be used for 1D atmospheres.");
1736  }
1737 
1738  const Index npos = sensor_pos.nrows();
1739 
1740  if (ztan_vector.nelem() != npos) {
1741  ostringstream os;
1742  os << "The number of altitudes in the geometric tangent altitude vector\n"
1743  << "must match the number of positions in *sensor_pos*.";
1744  throw runtime_error(os.str());
1745  }
1746 
1747  za_vector.resize(npos);
1748 
1749  for (Index i = 0; i < npos; i++) {
1750  if (ztan_vector[i] > sensor_pos(i, 0)) {
1751  ostringstream os;
1752  os << "Invalid observation geometry: sensor (at z=" << sensor_pos(i, 0)
1753  << "m) is located below the requested tangent altitude (tanh="
1754  << ztan_vector[i] << "m)";
1755  throw runtime_error(os.str());
1756  } else {
1757  za_vector[i] = geompath_za_at_r(refellipsoid[0] + ztan_vector[i],
1758  100,
1759  refellipsoid[0] + sensor_pos(i, 0));
1760  }
1761  }
1762 }
1763 
Matrix
The Matrix class.
Definition: matpackI.h:1193
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
sph2cart
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition: geodetic.cc:1236
ARTS::Var::ppath_step
Ppath ppath_step(Workspace &ws) noexcept
Definition: autoarts.h:5213
first_pos_before_altitude
Index first_pos_before_altitude(const Ppath &p, const Numeric &alt)
Determines ppath position just below an altitude.
Definition: ppath.cc:537
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
DiffZaAa
void DiffZaAa(Matrix &dlos, const Vector &ref_los, const Matrix &other_los, const Verbosity &)
WORKSPACE METHOD: DiffZaAa.
Definition: m_ppath.cc:79
ppath_init_structure
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
Initiates a Ppath structure to hold the given number of points.
Definition: ppath.cc:1426
ppath_stepRefractionBasic
void ppath_stepRefractionBasic(Workspace &ws, Ppath &ppath_step, const Agenda &refr_index_air_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Vector &f_grid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Verbosity &)
WORKSPACE METHOD: ppath_stepRefractionBasic.
Definition: m_ppath.cc:1336
ARTS::Var::ppath_step_agenda
Agenda ppath_step_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5220
ARTS::Var::lat
Numeric lat(Workspace &ws) noexcept
Definition: autoarts.h:3933
auto_md.h
gridpos_copy
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Definition: interpolation.cc:473
ARTS::Var::atmgeom_checked
Index atmgeom_checked(Workspace &ws) noexcept
Definition: autoarts.h:2489
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
ppath_agendaExecute
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:24775
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
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
TangentPointExtract
void TangentPointExtract(Vector &tan_pos, const Ppath &ppath, const Verbosity &)
WORKSPACE METHOD: TangentPointExtract.
Definition: m_ppath.cc:1623
rte_posSet
void rte_posSet(Vector &rte_pos, const Index &atmosphere_dim, const Numeric &z, const Numeric &lat, const Numeric &lon, const Verbosity &)
WORKSPACE METHOD: rte_posSet.
Definition: m_ppath.cc:1537
linreg
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:2302
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
ARTS::Var::ppath_inside_cloudbox_do
Index ppath_inside_cloudbox_do(Workspace &ws) noexcept
Definition: autoarts.h:5172
ppath_step_refr_1d
void ppath_step_refr_1d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 1D propagation path steps including effects of refraction.
Definition: ppath.cc:3576
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
ARTS::Var::refr_index_air_agenda
Agenda refr_index_air_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5539
ppath_step_geom_2d
void ppath_step_geom_2d(Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax)
Calculates 2D geometrical propagation path steps.
Definition: ppath.cc:2736
ppath_step_geom_1d
void ppath_step_geom_1d(Ppath &ppath, ConstVectorView z_field, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax)
Calculates 1D geometrical propagation path steps.
Definition: ppath.cc:2372
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ppath_fieldFromDownUpLimbGeoms
void ppath_fieldFromDownUpLimbGeoms(Workspace &ws, ArrayOfPpath &ppath_field, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Tensor3 &z_field, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Vector &refellipsoid, const Index &atmosphere_dim, const Index &zenith_angles_per_position, const Verbosity &verbosity)
WORKSPACE METHOD: ppath_fieldFromDownUpLimbGeoms.
Definition: m_ppath.cc:1120
ARTS::Var::refr_index_air
Numeric refr_index_air(Workspace &ws) noexcept
Definition: autoarts.h:5532
WriteXML
void WriteXML(const String &file_format, const T &v, const String &f, const Index &no_clobber, const String &v_name, const String &f_name, const String &no_clobber_name, const Verbosity &verbosity)
WORKSPACE METHOD: WriteXML.
Definition: m_xml.h:118
SWITCH_OUTPUT
#define SWITCH_OUTPUT(x, y)
Definition: m_general.h:46
get_refr_index_1d
void get_refr_index_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
get_refr_index_1d
Definition: refraction.cc:184
geo_posWherePpathPassesZref
void geo_posWherePpathPassesZref(Vector &geo_pos, const Ppath &ppath, const Numeric &z_ref, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posWherePpathPassesZref.
Definition: m_ppath.cc:145
geo_posLowestAltitudeOfPpath
void geo_posLowestAltitudeOfPpath(Vector &geo_pos, const Ppath &ppath, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posLowestAltitudeOfPpath.
Definition: m_ppath.cc:119
VectorZtanToZa1D
void VectorZtanToZa1D(Vector &za_vector, const Matrix &sensor_pos, const Vector &refellipsoid, const Index &atmosphere_dim, const Vector &ztan_vector, const Verbosity &)
WORKSPACE METHOD: VectorZtanToZa1D.
Definition: m_ppath.cc:1728
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
line_sphere_intersect
void line_sphere_intersect(Numeric &x, Numeric &y, Numeric &z, const Numeric &xl, const Numeric &yl, const Numeric &zl, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &xc, const Numeric &yc, const Numeric &zc, const Numeric &r)
geomtanpoint
Definition: geodetic.cc:823
ppath_fieldCalc
void ppath_fieldCalc(Workspace &ws, ArrayOfPpath &ppath_field, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Matrix &sensor_pos, const Matrix &sensor_los, const Vector &rte_pos2, const Verbosity &verbosity)
WORKSPACE METHOD: ppath_fieldCalc.
Definition: m_ppath.cc:1239
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
chk_rte_los
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
Definition: check_input.cc:1839
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
DEG2RAD
const Numeric DEG2RAD
get_refr_index_2d
void get_refr_index_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
get_refr_index_2d
Definition: refraction.cc:264
ppathPlaneParallel
void ppathPlaneParallel(Ppath &ppath, const Index &atmosphere_dim, const Tensor3 &z_field, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Verbosity &)
WORKSPACE METHOD: ppathPlaneParallel.
Definition: m_ppath.cc:761
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
GridPos::fd
Numeric fd[2]
Definition: interpolation.h:75
TangentPointPrint
void TangentPointPrint(const Ppath &ppath, const Index &level, const Verbosity &verbosity)
WORKSPACE METHOD: TangentPointPrint.
Definition: m_ppath.cc:1643
sensor_losGeometricFromSensorPosToOtherPositions
void sensor_losGeometricFromSensorPosToOtherPositions(Matrix &sensor_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &sensor_pos, const Matrix &target_pos, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_losGeometricFromSensorPosToOtherPositions.
Definition: m_ppath.cc:1581
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Ppath::gp_p
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
Agenda
The Agenda class.
Definition: agenda_class.h:44
ARTS::Var::lon
Numeric lon(Workspace &ws) noexcept
Definition: autoarts.h:4064
ppath_step_refr_3d
void ppath_step_refr_3d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 3D propagation path steps, with refraction, using a simple and fast ray tracing scheme.
Definition: ppath.cc:4335
cart2poslos
void cart2poslos(Numeric &r, Numeric &lat, Numeric &za, const Numeric &x, const Numeric &z, const Numeric &dx, const Numeric &dz, const Numeric &ppc, const Numeric &lat0, const Numeric &za0)
cart2poslos
Definition: geodetic.cc:113
Array
This can be used to make arrays out of anything.
Definition: array.h:108
ppathCalcFromAltitude
void ppathCalcFromAltitude(Workspace &ws, Ppath &ppath, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Numeric &altitude, const Numeric &accuracy, const Verbosity &verbosity)
WORKSPACE METHOD: ppathCalcFromAltitude.
Definition: m_ppath.cc:232
rte_pos_losMoveToStartOfPpath
void rte_pos_losMoveToStartOfPpath(Vector &rte_pos, Vector &rte_los, const Index &atmosphere_dim, const Ppath &ppath, const Verbosity &)
WORKSPACE METHOD: rte_pos_losMoveToStartOfPpath.
Definition: m_ppath.cc:1557
cart2pol
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:74
min
#define min(a, b)
Definition: legacy_continua.cc:20628
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:207
Ppath::gp_lat
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath.h:84
ppathCalc
void ppathCalc(Workspace &ws, Ppath &ppath, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Verbosity &)
WORKSPACE METHOD: ppathCalc.
Definition: m_ppath.cc:193
messages.h
Declarations having to do with the four output streams.
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
Ppath::ngroup
Vector ngroup
The group index of refraction.
Definition: ppath.h:80
pos2refell_r
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1081
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
my_basic_string< char >
ppath_stepGeometric
void ppath_stepGeometric(Ppath &ppath_step, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Numeric &ppath_lmax, const Verbosity &)
WORKSPACE METHOD: ppath_stepGeometric.
Definition: m_ppath.cc:1278
rte_losSet
void rte_losSet(Vector &rte_los, const Index &atmosphere_dim, const Numeric &za, const Numeric &aa, const Verbosity &)
WORKSPACE METHOD: rte_losSet.
Definition: m_ppath.cc:1458
Ppath::los
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
Ppath::gp_lon
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath.h:86
line_circle_intersect
void line_circle_intersect(Numeric &x, Numeric &z, const Numeric &xl, const Numeric &zl, const Numeric &dx, const Numeric &dz, const Numeric &xc, const Numeric &zc, const Numeric &r)
geomtanpoint2d
Definition: geodetic.cc:278
geompath_za_at_r
Numeric geompath_za_at_r(const Numeric &ppc, const Numeric &a_za, const Numeric &r)
Calculates the zenith angle for a given radius along a geometrical propagation path.
Definition: ppath.cc:103
Ppath::np
Index np
Number of points describing the ppath.
Definition: ppath.h:52
VectorZtanToZaRefr1D
void VectorZtanToZaRefr1D(Workspace &ws, Vector &za_vector, const Agenda &refr_index_air_agenda, const Matrix &sensor_pos, const Vector &p_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Index &atmosphere_dim, const Vector &f_grid, const Vector &ztan_vector, const Verbosity &)
WORKSPACE METHOD: VectorZtanToZaRefr1D.
Definition: m_ppath.cc:1668
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
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
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
los2xyz
void los2xyz(Numeric &za, Numeric &aa, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &x1, const Numeric &y1, const Numeric &z1, const Numeric &x2, const Numeric &y2, const Numeric &z2)
los2xyz
Definition: geodetic.cc:929
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::ppath_field
ArrayOfPpath ppath_field(Workspace &ws) noexcept
Definition: autoarts.h:5158
ARTS::Var::refr_index_air_group
Numeric refr_index_air_group(Workspace &ws) noexcept
Definition: autoarts.h:5553
math_funcs.h
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::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
ARTS::Var::sensor_pos
Matrix sensor_pos(Workspace &ws) noexcept
Definition: autoarts.h:6279
max
#define max(a, b)
Definition: legacy_continua.cc:20629
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
chk_rte_pos
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
Definition: check_input.cc:1770
ppath_step_refr_2d
void ppath_step_refr_2d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 2D propagation path steps, with refraction, using a simple and fast ray tracing scheme.
Definition: ppath.cc:3920
m_xml.h
Workspace methods and template functions for supergeneric XML IO.
ppathWriteXMLPartial
void ppathWriteXMLPartial(const String &file_format, const Ppath &ppath, const String &f, const Index &file_index, const Verbosity &verbosity)
WORKSPACE METHOD: ppathWriteXMLPartial.
Definition: m_ppath.cc:1091
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
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ARTS::Var::ppath_lmax
Numeric ppath_lmax(Workspace &ws) noexcept
Definition: autoarts.h:5183
ARTS::Var::ppath_lraytrace
Numeric ppath_lraytrace(Workspace &ws) noexcept
Definition: autoarts.h:5195
ppathStepByStep
void ppathStepByStep(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &ppath_inside_cloudbox_do, 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 Verbosity &verbosity)
WORKSPACE METHOD: ppathStepByStep.
Definition: m_ppath.cc:1050
Range
The range class.
Definition: matpackI.h:160
ppath_copy
void ppath_copy(Ppath &ppath1, const Ppath &ppath2, const Index &ncopy)
Copy the content in ppath2 to ppath1.
Definition: ppath.cc:1515
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ARTS::Var::ppath_agenda
Agenda ppath_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5146
ppath.h
Propagation path structure and functions.
find_tanpoint
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
Definition: ppath.cc:525
ppathFromRtePos2
void ppathFromRtePos2(Workspace &ws, Ppath &ppath, Vector &rte_los, Numeric &ppath_lraytrace, 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 Vector &rte_pos, const Vector &rte_pos2, const Numeric &ppath_lmax, const Numeric &za_accuracy, const Numeric &pplrt_factor, const Numeric &pplrt_lowest, const Verbosity &verbosity)
WORKSPACE METHOD: ppathFromRtePos2.
Definition: m_ppath.cc:306
ARTS::Var::rte_pos2
Vector rte_pos2(Workspace &ws) noexcept
Definition: autoarts.h:5696
geo_posEndOfPpath
void geo_posEndOfPpath(Vector &geo_pos, const Ppath &ppath, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posEndOfPpath.
Definition: m_ppath.cc:103
refraction.h
Refraction functions.
pol2cart
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition: geodetic.cc:331
x1
#define x1
Definition: lineshapemodel.h:414
ppath_step_geom_3d
void ppath_step_geom_3d(Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax)
Calculates 3D geometrical propagation path steps.
Definition: ppath.cc:3270
ll
#define ll
Definition: legacy_continua.cc:21959
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
Ppath::r
Vector r
Radius of each ppath point.
Definition: ppath.h:68
special_interp.h
Header file for special_interp.cc.
get_refr_index_3d
void get_refr_index_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
Definition: refraction.cc:357
RAD2DEG
const Numeric RAD2DEG
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
rte_pos2gridpos
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
Definition: special_interp.cc:788
ARTS::Var::file_index
Index file_index(Workspace &ws) noexcept
Definition: autoarts.h:3475
GridPos::idx
Index idx
Definition: interpolation.h:74
diff_za_aa
void diff_za_aa(Numeric &dza, Numeric &daa, const Numeric &za0, const Numeric &aa0, const Numeric &za, const Numeric &aa)
Takes the difference of zenith and azimuth angles.
Definition: ppath.cc:444
AddZaAa
void AddZaAa(Matrix &new_los, const Vector &ref_los, const Matrix &dlos, const Verbosity &)
WORKSPACE METHOD: AddZaAa.
Definition: m_ppath.cc:56
ppath_set_background
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
Definition: ppath.cc:1467
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
check_input.h
rte_losGeometricFromRtePosToRtePos2
void rte_losGeometricFromRtePosToRtePos2(Vector &rte_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Vector &rte_pos, const Vector &rte_pos2, const Verbosity &)
WORKSPACE METHOD: rte_losGeometricFromRtePosToRtePos2.
Definition: m_ppath.cc:1476
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
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
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
x2
#define x2
Definition: lineshapemodel.h:415
ARTS::Var::sensor_los
Matrix sensor_los(Workspace &ws) noexcept
Definition: autoarts.h:6201
ARTS::Var::cloudbox_checked
Index cloudbox_checked(Workspace &ws) noexcept
Definition: autoarts.h:2650
ARTS::Var::z_surface
Matrix z_surface(Workspace &ws) noexcept
Definition: autoarts.h:7754
filename_xml_with_index
void filename_xml_with_index(String &filename, const Index &file_index, const String &varname, const Index &digits)
Gives the default filename, with file index, for the XML formats.
Definition: xml_io.cc:450
distance2D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:182
CREATE_OUTS
#define CREATE_OUTS
Definition: messages.h:209
arts.h
The global header file for ARTS.
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
xml_io.h
This file contains basic functions to handle XML data files.
m_general.h
Template functions for general supergeneric ws methods.