ARTS 2.5.11 (git: 725533f0)
ppath.cc
Go to the documentation of this file.
1
14/*===========================================================================
15 === External declarations
16 ===========================================================================*/
17
18#include "ppath.h"
19#include <cmath>
20#include <stdexcept>
21#include "agenda_class.h"
22#include "array.h"
23#include "arts_conversions.h"
24#include "arts_omp.h"
25#include "auto_md.h"
26#include "check_input.h"
27#include "geodetic.h"
28#include "logic.h"
29#include "math_funcs.h"
30#include "messages.h"
31#include "mystring.h"
32#include "poly_roots.h"
33#include "refraction.h"
34#include "rte.h"
35#include "special_interp.h"
36
37inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
38inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
39
40/*===========================================================================
41 === Precision variables
42 ===========================================================================*/
43
44// This variable defines the maximum allowed error tolerance for radius.
45// The variable is, for example, used to check that a given a radius is
46// consistent with the specified grid cell.
47//
48const Numeric RTOL = 1e-3;
49
50// As RTOL but for latitudes and longitudes.
51//
52DEBUG_ONLY(const Numeric LATLONTOL = 1e-8;)
53
54// Accuarcy for length comparisons.
55//
56const Numeric LACC = 1e-5;
57
58// Values to apply if some calculation does not provide a solution
59const Numeric R_NOT_FOUND = -1; // A value below zero
60const Numeric L_NOT_FOUND = 99e99; // Some very large value for l/lat/lon
61const Numeric LAT_NOT_FOUND = 99e99;
62const Numeric LON_NOT_FOUND = 99e99;
63
64/*===========================================================================
65 === Functions related to geometrical propagation paths
66 ===========================================================================*/
67
80Numeric geometrical_ppc(const Numeric& r, const Numeric& za) {
81 ARTS_ASSERT(r > 0);
82 ARTS_ASSERT(abs(za) <= 180);
83
84 return r * sin(DEG2RAD * abs(za));
85}
86
87Numeric geompath_za_at_r(const Numeric& ppc,
88 const Numeric& a_za,
89 const Numeric& r) {
90 ARTS_ASSERT(ppc >= 0);
91 ARTS_ASSERT(abs(a_za) <= 180);
92 ARTS_ASSERT(r >= ppc - RTOL);
93
94 if (r > ppc) {
95 Numeric za = RAD2DEG * asin(ppc / r);
96 if (abs(a_za) > 90) {
97 za = 180 - za;
98 }
99 if (a_za < 0) {
100 za = -za;
101 }
102 return za;
103 } else {
104 if (a_za > 0) {
105 return 90;
106 } else {
107 return -90;
108 }
109 }
110}
111
125Numeric geompath_r_at_za(const Numeric& ppc, const Numeric& za) {
126 ARTS_ASSERT(ppc >= 0);
127 ARTS_ASSERT(abs(za) <= 180);
128
129 return ppc / sin(DEG2RAD * abs(za));
130}
131
132Numeric geompath_lat_at_za(const Numeric& za0,
133 const Numeric& lat0,
134 const Numeric& za) {
135 ARTS_ASSERT(abs(za0) <= 180);
136 ARTS_ASSERT(abs(za) <= 180);
137 ARTS_ASSERT((za0 >= 0 && za >= 0) || (za0 < 0 && za < 0));
138
139 return lat0 + za0 - za;
140}
141
142Numeric geompath_l_at_r(const Numeric& ppc, const Numeric& r) {
143 ARTS_ASSERT(ppc >= 0);
144 ARTS_ASSERT(r >= ppc - RTOL);
145
146 if (r > ppc) {
147 return sqrt(r * r - ppc * ppc);
148 } else {
149 return 0;
150 }
151}
152
153Numeric geompath_r_at_l(const Numeric& ppc, const Numeric& l) {
154 ARTS_ASSERT(ppc >= 0);
155
156 return sqrt(l * l + ppc * ppc);
157}
158
171Numeric geompath_r_at_lat(const Numeric& ppc,
172 const Numeric& lat0,
173 const Numeric& za0,
174 const Numeric& lat) {
175 ARTS_ASSERT(ppc >= 0);
176 ARTS_ASSERT(abs(za0) <= 180);
177 ARTS_ASSERT((za0 >= 0 && lat >= lat0) || (za0 <= 0 && lat <= lat0));
178
179 // Zenith angle at the new latitude
180 const Numeric za = za0 + lat0 - lat;
181
182 return geompath_r_at_za(ppc, za);
183}
184
208 Vector& lat,
209 Vector& za,
210 Numeric& lstep,
211 const Numeric& ppc,
212 const Numeric& r1,
213 const Numeric& lat1,
214 const Numeric& za1,
215 const Numeric& r2,
216 const bool& tanpoint,
217 const Numeric& lmax) {
218 // Calculate length from tangent point, along the path for point 1 and 2.
219 // Length defined in the viewing direction, ie. negative at end of sensor
220 Numeric l1 = geompath_l_at_r(ppc, r1);
221 if (abs(za1) > 90) {
222 l1 *= -1;
223 }
224 Numeric l2 = geompath_l_at_r(ppc, r2);
225 if (l1 < 0) {
226 l2 *= -1;
227 }
228 if (tanpoint) {
229 l2 *= -1;
230 }
231
232 // Calculate needed number of steps, considering a possible length criterion
233 Index n;
234 if (lmax > 0) {
235 // The absolute value of the length distance is needed here
236 // We can't accept n=0, which is the case if l1 = l2
237 n = max(Index(1), Index(ceil(abs(l2 - l1) / lmax)));
238 } else {
239 n = 1;
240 }
241
242 // Length of path steps (note that lstep here can become negative)
243 lstep = (l2 - l1) / (Numeric)n;
244
245 // Allocate vectors and put in point 1
246 r.resize(n + 1);
247 lat.resize(n + 1);
248 za.resize(n + 1);
249 r[0] = r1;
250 lat[0] = lat1;
251 za[0] = za1;
252
253 // Loop steps (beside last) and calculate radius and zenith angle
254 for (Index i = 1; i < n; i++) {
255 const Numeric l = l1 + lstep * (Numeric)i;
256 r[i] = geompath_r_at_l(ppc, l); // Sign of l does not matter here
257 // Set a zenith angle to 80 or 100 depending on sign of l
258 za[i] = geompath_za_at_r(ppc, sign(za1) * (90 - sign(l) * 10), r[i]);
259 }
260
261 // For maximum accuracy, set last radius to be exactly r2.
262 r[n] = r2; // Don't use r[n] below
263 za[n] = geompath_za_at_r(ppc, sign(za1) * (90 - sign(l2) * 10), r2);
264
265 // Ensure that zenith and nadir observations keep their zenith angle
266 if (abs(za1) < ANGTOL || abs(za1) > 180 - ANGTOL) {
267 za = za1;
268 }
269
270 // Calculate latitudes
271 for (Index i = 1; i <= n; i++) {
272 lat[i] = geompath_lat_at_za(za1, lat1, za[i]);
273 }
274
275 // Take absolute value of lstep
276 lstep = abs(lstep);
277}
278
279/*===========================================================================
280 === Functions focusing on zenith and azimuth angles
281 ===========================================================================*/
282
283void cart2zaaa(Numeric& za,
284 Numeric& aa,
285 const Numeric& dx,
286 const Numeric& dy,
287 const Numeric& dz) {
288 const Numeric r = sqrt(dx * dx + dy * dy + dz * dz);
289
290 ARTS_ASSERT(r > 0);
291
292 za = RAD2DEG * acos(dz / r);
293 aa = RAD2DEG * atan2(dy, dx);
294}
295
296void zaaa2cart(Numeric& dx,
297 Numeric& dy,
298 Numeric& dz,
299 const Numeric& za,
300 const Numeric& aa) {
301 const Numeric zarad = DEG2RAD * za;
302 const Numeric aarad = DEG2RAD * aa;
303
304 dz = cos(zarad);
305 dx = sin(zarad);
306 dy = sin(aarad) * dx;
307 dx = cos(aarad) * dx;
308}
309
310void enu2zaaa(Numeric& za,
311 Numeric& aa,
312 const Numeric& de,
313 const Numeric& dn,
314 const Numeric& du) {
315 za = RAD2DEG * acos( du );
316 aa = RAD2DEG * atan2( de, dn );
317}
318
319void zaaa2enu(Numeric& de,
320 Numeric& dn,
321 Numeric& du,
322 const Numeric& za,
323 const Numeric& aa) {
324 const Numeric st = sin(DEG2RAD*za);
325 de = st * sin(DEG2RAD*aa);
326 dn = st * cos(DEG2RAD*aa);
327 du = cos(DEG2RAD*za);
328}
329
346void rotationmat3D(Matrix& R, ConstVectorView vrot, const Numeric& a) {
347 ARTS_ASSERT(R.ncols() == 3);
348 ARTS_ASSERT(R.nrows() == 3);
349 ARTS_ASSERT(vrot.nelem() == 3);
350
351 const Numeric l =
352 sqrt(vrot[0] * vrot[0] + vrot[1] * vrot[1] + vrot[2] * vrot[2]);
353 ARTS_ASSERT(l > 1 - 9);
354 const Numeric u = vrot[0] / l;
355 const Numeric v = vrot[1] / l;
356 const Numeric w = vrot[2] / l;
357 const Numeric u2 = u * u;
358 const Numeric v2 = v * v;
359 const Numeric w2 = w * w;
360 const Numeric c = cos(DEG2RAD * a);
361 const Numeric s = sin(DEG2RAD * a);
362
363 // Fill R
364 R(0, 0) = u2 + (v2 + w2) * c;
365 R(0, 1) = u * v * (1 - c) - w * s;
366 R(0, 2) = u * w * (1 - c) + v * s;
367 R(1, 0) = u * v * (1 - c) + w * s;
368 R(1, 1) = v2 + (u2 + w2) * c;
369 R(1, 2) = v * w * (1 - c) - u * s;
370 R(2, 0) = u * w * (1 - c) - v * s;
371 R(2, 1) = v * w * (1 - c) + u * s;
372 R(2, 2) = w2 + (u2 + v2) * c;
373}
374
375void add_za_aa(Numeric& za,
376 Numeric& aa,
377 const Numeric& za0,
378 const Numeric& aa0,
379 const Numeric& dza,
380 const Numeric& daa) {
381 Vector xyz(3);
382 Vector vrot(3);
383 Vector u(3);
384
385 // Unit vector towards aa0 at za=90
386 //
387 zaaa2cart(xyz[0], xyz[1], xyz[2], 90, aa0);
388
389 // Find vector around which rotation shall be performed
390 //
391 // We can write this as cross([0 0 1],xyz). It turns out that the result
392 // of this operation is just [-y,x,0].
393 //
394 vrot[0] = -xyz[1];
395 vrot[1] = xyz[0];
396 vrot[2] = 0;
397
398 // Unit vector towards aa0+aa at za=90
399 //
400 zaaa2cart(xyz[0], xyz[1], xyz[2], 90 + dza, aa0 + daa);
401
402 // Apply rotation
403 //
404 Matrix R(3, 3);
405 rotationmat3D(R, vrot, za0 - 90);
406 mult(u, R, xyz);
407
408 // Calculate za and aa for rotated u
409 //
410 cart2zaaa(za, aa, u[0], u[1], u[2]);
411}
412
413void diff_za_aa(Numeric& dza,
414 Numeric& daa,
415 const Numeric& za0,
416 const Numeric& aa0,
417 const Numeric& za,
418 const Numeric& aa) {
419 Vector xyz(3);
420 Vector vrot(3);
421 Vector u(3);
422
423 ARTS_ASSERT(za != 0 && za != 180);
424
425 // Unit vector towards aa0 at za=90
426 //
427 zaaa2cart(xyz[0], xyz[1], xyz[2], za0, aa0);
428
429 // Find vector around which rotation shall be performed
430 //
431 // We can write this as cross([0 0 1],xyz). It turns out that the result
432 // of this operation is just [-y,x,0].
433 //
434 vrot[0] = -xyz[1];
435 vrot[1] = xyz[0];
436 vrot[2] = 0;
437
438 // Unit vector towards aa0+aa at za=90
439 //
440 zaaa2cart(xyz[0], xyz[1], xyz[2], za, aa);
441
442 // Apply rotation
443 //
444 Matrix R(3, 3);
445 rotationmat3D(R, vrot, -(za0 - 90));
446 mult(u, R, xyz);
447
448 // Calculate za and aa for rotated u
449 //
450 Numeric za_tmp, aa_tmp;
451 cart2zaaa(za_tmp, aa_tmp, u[0], u[1], u[2]);
452
453 // Calculate dza and daa
454 dza = za_tmp - 90;
455 daa = aa_tmp - aa0;
456}
457
458/*===========================================================================
459 === Various functions
460 ===========================================================================*/
461
475Numeric refraction_ppc(const Numeric& r,
476 const Numeric& za,
477 const Numeric& refr_index_air) {
478 ARTS_ASSERT(r > 0);
479 ARTS_ASSERT(abs(za) <= 180);
480
481 return r * refr_index_air * sin(DEG2RAD * abs(za));
482}
483
484void resolve_lon(Numeric& lon, const Numeric& lon5, const Numeric& lon6) {
485 ARTS_ASSERT(lon6 >= lon5);
486
487 if (lon < lon5 && lon + 180 <= lon6) {
488 lon += 360;
489 } else if (lon > lon6 && lon - 180 >= lon5) {
490 lon -= 360;
491 }
492}
493
494void find_tanpoint(Index& it, const Ppath& ppath) {
495 Numeric zmin = 99e99;
496 it = -1;
497 while (it < ppath.np - 1 && ppath.pos(it + 1, 0) < zmin) {
498 it++;
499 zmin = ppath.pos(it, 0);
500 }
501 if (it == 0 || it == ppath.np - 1) {
502 it = -1;
503 }
504}
505
506Index first_pos_before_altitude(const Ppath& p, const Numeric& alt) {
507 // Checker flags
508 bool below = false, above = false;
509
510 // Check first pos before altitude
511 for (Index i = 0; i < p.np; i++) {
512 if (p.pos(i, 0) < alt) {
513 if (above) return i - 1;
514 below = true;
515 } else {
516 if (below) return i - 1;
517 above = true;
518 }
519 }
520
521 return -1;
522}
523
524void error_if_limb_ppath(const Ppath& ppath) {
525 if (ppath.np > 2) {
526 Numeric signfac = sign(ppath.pos(1, 0) - ppath.pos(0, 0));
527 for (Index i = 2; i < ppath.np; i++) {
528 ARTS_USER_ERROR_IF (signfac * (ppath.pos(i, 0) - ppath.pos(i - 1, 0)) < 0,
529 "A propagation path of limb character found. Such "
530 "viewing geometries are not supported by this method. "
531 "Propagation paths must result in strictly "
532 "increasing or decreasing altitudes.");
533 }
534 }
535}
536
537/*===========================================================================
538 = 2D functions for surface and pressure level slope and tilt
539 ===========================================================================*/
540
555Numeric rsurf_at_lat(const Numeric& lat1,
556 const Numeric& lat3,
557 const Numeric& r1,
558 const Numeric& r3,
559 const Numeric& lat) {
560 return r1 + (lat - lat1) * (r3 - r1) / (lat3 - lat1);
561}
562
563void plevel_slope_2d(Numeric& c1,
564 ConstVectorView lat_grid,
565 ConstVectorView refellipsoid,
566 ConstVectorView z_surf,
567 const GridPos& gp,
568 const Numeric& za) {
569 Index i1 = gridpos2gridrange(gp, za >= 0);
570 const Numeric r1 = refell2r(refellipsoid, lat_grid[i1]) + z_surf[i1];
571 const Numeric r2 = refell2r(refellipsoid, lat_grid[i1 + 1]) + z_surf[i1 + 1];
572 //
573 c1 = (r2 - r1) / (lat_grid[i1 + 1] - lat_grid[i1]);
574}
575
592void plevel_slope_2d(Numeric& c1,
593 const Numeric& lat1,
594 const Numeric& lat2,
595 const Numeric& r1,
596 const Numeric& r2) {
597 c1 = (r2 - r1) / (lat2 - lat1);
598}
599
600Numeric plevel_angletilt(const Numeric& r, const Numeric& c1) {
601 // The tilt (in radians) is c1/r if c1 is converted to m/radian. So we get
602 // conversion RAD2DEG twice
603 return RAD2DEG * RAD2DEG * c1 / r;
604}
605
606bool is_los_downwards(const Numeric& za, const Numeric& tilt) {
607 ARTS_ASSERT(abs(za) <= 180);
608
609 // Yes, it shall be -tilt in both cases, if you wonder.
610 if (za > (90 - tilt) || za < (-90 - tilt)) {
611 return true;
612 } else {
613 return false;
614 }
615}
616
638void r_crossing_2d(Numeric& lat,
639 Numeric& l,
640 const Numeric& r_hit,
641 const Numeric& r_start,
642 const Numeric& lat_start,
643 const Numeric& za_start,
644 const Numeric& ppc) {
645 ARTS_ASSERT(abs(za_start) <= 180);
646 ARTS_ASSERT(r_start >= ppc);
647
648 const Numeric absza = abs(za_start);
649
650 // If above and looking upwards or r_hit below tangent point,
651 // we have no crossing:
652 if ((r_start >= r_hit && absza <= 90) || ppc > r_hit) {
653 lat = LAT_NOT_FOUND;
654 l = L_NOT_FOUND;
655 }
656
657 else {
658 // Passages of tangent point
659 if (absza > 90 && r_start <= r_hit) {
660 Numeric za = geompath_za_at_r(ppc, sign(za_start) * 89, r_hit);
661 lat = geompath_lat_at_za(za_start, lat_start, za);
662 l = geompath_l_at_r(ppc, r_start) + geompath_l_at_r(ppc, r_hit);
663 }
664
665 else {
666 Numeric za = geompath_za_at_r(ppc, za_start, r_hit);
667 lat = geompath_lat_at_za(za_start, lat_start, za);
668 l = abs(geompath_l_at_r(ppc, r_start) - geompath_l_at_r(ppc, r_hit));
669 ARTS_ASSERT(l > 0);
670 }
671 }
672}
673
709Numeric rslope_crossing2d(const Numeric& rp,
710 const Numeric& za,
711 const Numeric& r0,
712 Numeric c1) {
713 const Numeric zaabs = abs(za);
714
715 ARTS_ASSERT(za != 0);
716 ARTS_ASSERT(zaabs != 180);
717 ARTS_ASSERT(abs(c1) > 0); // c1=0 should work, but unnecessary to use this func.
718
719 // Convert slope to m/radian and consider viewing direction
720 c1 *= RAD2DEG;
721 if (za < 0) {
722 c1 = -c1;
723 }
724
725 // The nadir angle in radians, and cosine and sine of that angle
726 const Numeric beta = DEG2RAD * (180 - zaabs);
727 const Numeric cv = cos(beta);
728 const Numeric sv = sin(beta);
729
730 // Some repeated terms
731 const Numeric r0s = r0 * sv;
732 const Numeric r0c = r0 * cv;
733 const Numeric cs = c1 * sv;
734 const Numeric cc = c1 * cv;
735
736 // The vector of polynomial coefficients
737 //
738 Index n = 6;
739 Vector p0(n + 1);
740 //
741 p0[0] = r0s - rp * sv;
742 p0[1] = r0c + cs;
743 p0[2] = -r0s / 2 + cc;
744 p0[3] = -r0c / 6 - cs / 2;
745 p0[4] = r0s / 24 - cc / 6;
746 p0[5] = r0c / 120 + cs / 24;
747 p0[6] = -r0s / 720 + cc / 120;
748 //
749 // The accuracy when solving the polynomial equation gets worse when
750 // approaching 0 and 180 deg. The solution is to let the start polynomial
751 // order decrease when approaching these angles. The values below based on
752 // practical experience, don't change without making extremely careful tests.
753 //
754 if (abs(90 - zaabs) > 89.9) {
755 n = 1;
756 } else if (abs(90 - zaabs) > 75) {
757 n = 4;
758 }
759
760 // Calculate roots of the polynomial
761 Matrix roots;
762 int solutionfailure = 1;
763 //
764 while (solutionfailure) {
765 roots.resize(n, 2);
766 Vector p;
767 p = p0[Range(0, n + 1)];
768 solutionfailure = poly_root_solve(roots, p);
769 if (solutionfailure) {
770 n -= 1;
771 ARTS_ASSERT(n > 0);
772 }
773 }
774
775 // If r0=rp, numerical inaccuracy can give a false solution, very close
776 // to 0, that we must throw away.
777 Numeric dmin = 0;
778 if (abs(r0 - rp) < 1e-9) // 1 nm set based on practical experience.
779 {
780 dmin = 5e-12;
781 }
782
783 // Find the smallest root with imaginary part = 0, and real part > 0.
784 Numeric dlat = 1.571; // Not interested in solutions above 90 deg!
785 //
786 for (Index i = 0; i < n; i++) {
787 if (roots(i, 1) == 0 && roots(i, 0) > dmin && roots(i, 0) < dlat) {
788 dlat = roots(i, 0);
789 }
790 }
791
792 if (dlat < 1.57) // A somewhat smaller value than start one
793 {
794 // Convert back to degrees
795 dlat = RAD2DEG * dlat;
796
797 // Change sign if zenith angle is negative
798 if (za < 0) {
799 dlat = -dlat;
800 }
801 } else {
802 dlat = LAT_NOT_FOUND;
803 }
804
805 return dlat;
806}
807
838void plevel_crossing_2d(Numeric& r,
839 Numeric& lat,
840 Numeric& l,
841 const Numeric& r_start0,
842 const Numeric& lat_start,
843 const Numeric& za_start,
844 const Numeric& ppc,
845 const Numeric& lat1,
846 const Numeric& lat3,
847 const Numeric& r1,
848 const Numeric& r3,
849 const bool& above) {
850 const Numeric absza = abs(za_start);
851
852 ARTS_ASSERT(absza <= 180);
853 ARTS_ASSERT(lat_start >= lat1 && lat_start <= lat3);
854
855 // Zenith case
856 if (absza < ANGTOL) {
857 if (above) {
858 r = R_NOT_FOUND;
859 lat = LAT_NOT_FOUND;
860 l = L_NOT_FOUND;
861 } else {
862 lat = lat_start;
863 r = rsurf_at_lat(lat1, lat3, r1, r3, lat);
864 l = max(1e-9, r - r_start0); // Max to ensure a small positive
865 } // step, to handle numerical issues
866 }
867
868 // Nadir case
869 else if (absza > 180 - ANGTOL) {
870 if (above) {
871 lat = lat_start;
872 r = rsurf_at_lat(lat1, lat3, r1, r3, lat);
873 l = max(1e-9, r_start0 - r); // As above
874 } else {
875 r = R_NOT_FOUND;
876 lat = LAT_NOT_FOUND;
877 l = L_NOT_FOUND;
878 }
879 }
880
881 // The general case
882 else {
883 const Numeric rmin = min(r1, r3);
884 const Numeric rmax = max(r1, r3);
885
886 // The case of negligible slope
887 if (rmax - rmin < 1e-6) {
888 // Set r_start and r, considering impact of numerical problems
889 Numeric r_start = r_start0;
890 r = r1;
891 if (above) {
892 if (r_start < rmax) {
893 r_start = r = rmax;
894 }
895 } else {
896 if (r_start > rmin) {
897 r_start = r = rmin;
898 }
899 }
900
901 r_crossing_2d(lat, l, r, r_start, lat_start, za_start, ppc);
902
903 // Check if inside [lat1,lat3]
904 if (lat > lat3 || lat < lat1) {
905 r = R_NOT_FOUND;
906 lat = LAT_NOT_FOUND;
907 }
908 }
909
910 // With slope
911 else {
912 // Set r_start, considering impact of numerical problems
913 Numeric r_start = r_start0;
914 if (above) {
915 if (r_start < rmin) {
916 r_start = rmin;
917 }
918 } else {
919 if (r_start > rmax) {
920 r_start = rmax;
921 }
922 }
923
924 Numeric za = 999;
925
926 // Calculate crossing with closest radius
927 if (r_start > rmax) {
928 r = rmax;
929 r_crossing_2d(lat, l, r, r_start, lat_start, za_start, ppc);
930 } else if (r_start < rmin) {
931 r = rmin;
932 r_crossing_2d(lat, l, r, r_start, lat_start, za_start, ppc);
933 } else {
934 r = r_start;
935 lat = lat_start;
936 l = 0;
937 za = za_start;
938 }
939
940 // lat must be inside [lat1,lat3] if relevant to continue
941 if (lat < lat1 || lat > lat3) {
942 r = R_NOT_FOUND;
943 } // lat already set by r_crossing_2d
944
945 // Otherwise continue from found point, considering the level slope
946 else {
947 // Level slope and radius at lat
948 Numeric cpl;
949 plevel_slope_2d(cpl, lat1, lat3, r1, r3);
950 const Numeric rpl = r1 + cpl * (lat - lat1);
951
952 // Make adjustment if numerical problems
953 if (above) {
954 if (r < rpl) {
955 r = rpl;
956 }
957 } else {
958 if (r > rpl) {
959 r = rpl;
960 }
961 }
962
963 // Calculate zenith angle at (r,lat) (if <180, already set above)
964 if (za > 180) // lat+za preserved (also with negative za)
965 {
966 za = lat_start + za_start - lat;
967 };
968
969 // Latitude distance from present point to actual crossing
970 const Numeric dlat = rslope_crossing2d(r, za, rpl, cpl);
971
972 // Update lat and check if still inside [lat1,lat3].
973 // If yes, determine r and l
974 lat += dlat;
975 if (lat < lat1 || lat > lat3) {
976 r = R_NOT_FOUND;
977 lat = LAT_NOT_FOUND;
978 l = L_NOT_FOUND;
979 } else {
980 // It was tested to calculate r from geompath functions, but
981 // appeared to give poorer accuracy around zenith/nadir
982 r = rpl + cpl * dlat;
983
984 // Zenith angle needed to check tangent point
985 za = lat_start + za_start - lat;
986
987 // Passage of tangent point requires special attention
988 if (absza > 90 && abs(za) < 90) {
989 l = geompath_l_at_r(ppc, r_start) + geompath_l_at_r(ppc, r);
990 } else {
991 l = abs(geompath_l_at_r(ppc, r_start) - geompath_l_at_r(ppc, r));
992 }
993 }
994 }
995 }
996 }
997}
998
999/*===========================================================================
1000 = 3D functions for level slope and tilt, and lat/lon crossings
1001 ===========================================================================*/
1002
1022Numeric rsurf_at_latlon(const Numeric& lat1,
1023 const Numeric& lat3,
1024 const Numeric& lon5,
1025 const Numeric& lon6,
1026 const Numeric& r15,
1027 const Numeric& r35,
1028 const Numeric& r36,
1029 const Numeric& r16,
1030 const Numeric& lat,
1031 const Numeric& lon) {
1032 // We can't have any ARTS_ASSERT of *lat* and *lon* here as we can go outside
1033 // the ranges when called from *plevel_slope_3d*.
1034
1035 if (lat == lat1) {
1036 return r15 + (lon - lon5) * (r16 - r15) / (lon6 - lon5);
1037 } else if (lat == lat3) {
1038 return r35 + (lon - lon5) * (r36 - r35) / (lon6 - lon5);
1039 } else if (lon == lon5) {
1040 return r15 + (lat - lat1) * (r35 - r15) / (lat3 - lat1);
1041 } else if (lon == lon6) {
1042 return r16 + (lat - lat1) * (r36 - r16) / (lat3 - lat1);
1043 } else {
1044 const Numeric fdlat = (lat - lat1) / (lat3 - lat1);
1045 const Numeric fdlon = (lon - lon5) / (lon6 - lon5);
1046 return (1 - fdlat) * (1 - fdlon) * r15 + fdlat * (1 - fdlon) * r35 +
1047 (1 - fdlat) * fdlon * r16 + fdlat * fdlon * r36;
1048 }
1049}
1050
1051void plevel_slope_3d(Numeric& c1,
1052 Numeric& c2,
1053 const Numeric& lat1,
1054 const Numeric& lat3,
1055 const Numeric& lon5,
1056 const Numeric& lon6,
1057 const Numeric& r15,
1058 const Numeric& r35,
1059 const Numeric& r36,
1060 const Numeric& r16,
1061 const Numeric& lat,
1062 const Numeric& lon,
1063 const Numeric& aa) {
1064 // Save time and avoid numerical problems if all r are equal
1065 if (r15 == r35 && r15 == r36 && r15 == r16 && r35 == r36 && r35 == r16 &&
1066 r36 == r16) {
1067 c1 = 0;
1068 c2 = 0;
1069 }
1070
1071 else {
1072 // Radius at point of interest
1073 const Numeric r0 =
1074 rsurf_at_latlon(lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat, lon);
1075
1076 // Size of test angular distance. Unit is degrees.
1077 // dang = 1e-4 corresponds to about 10 m shift horisontally.
1078 // Smaller values should probably be avoided. For example, 1e-5 gave
1079 // significant c2 when it should be zero.
1080 const Numeric dang = 1e-4;
1081
1082 Numeric lat2, lon2;
1083 latlon_at_aa(lat2, lon2, lat, lon, aa, dang);
1084 resolve_lon(lon2, lon5, lon6);
1085 const Numeric dr1 =
1087 lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat2, lon2) -
1088 r0;
1089
1090 latlon_at_aa(lat2, lon2, lat, lon, aa, 2 * dang);
1091 resolve_lon(lon2, lon5, lon6);
1092 const Numeric dr2 =
1094 lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat2, lon2) -
1095 r0;
1096
1097 // Derive linear and quadratic coefficient
1098 c1 = 0.5 * (4 * dr1 - dr2);
1099 c2 = (dr1 - c1) / (dang * dang);
1100 c1 /= dang;
1101 }
1102}
1103
1129void plevel_slope_3d(Numeric& c1,
1130 Numeric& c2,
1131 ConstVectorView lat_grid,
1132 ConstVectorView lon_grid,
1133 ConstVectorView refellipsoid,
1134 ConstMatrixView z_surf,
1135 const GridPos& gp_lat,
1136 const GridPos& gp_lon,
1137 const Numeric& aa) {
1138 Index ilat = gridpos2gridrange(gp_lat, abs(aa) >= 0);
1139 Index ilon = gridpos2gridrange(gp_lon, aa >= 0);
1140
1141 // Allow that we are at the upper edge of the grids only for special cases:
1142 const Index llat = lat_grid.nelem() - 1;
1143 // At North pole:
1144 if (ilat >= llat) {
1145 if (lat_grid[llat] > POLELAT) {
1146 ilat = llat - 1;
1147 } else {
1149 "The upper latitude end of the atmosphere "
1150 "reached, that is not allowed.");
1151 }
1152 }
1153 // Complete 360deg lon coverage:
1154 if (ilon >= lon_grid.nelem() - 1) {
1155 if (is_lon_cyclic(lon_grid)) {
1156 ilon = 0;
1157 } else {
1159 "The upper longitude end of the atmosphere "
1160 "reached, that is not allowed.");
1161 }
1162 }
1163
1164 // Restore latitude and longitude values
1165 Vector itw(2);
1166 Numeric lat, lon;
1167 interpweights(itw, gp_lat);
1168 lat = interp(itw, lat_grid, gp_lat);
1169 interpweights(itw, gp_lon);
1170 lon = interp(itw, lon_grid, gp_lon);
1171
1172 // Extract values that defines the grid cell
1173 const Numeric lat1 = lat_grid[ilat];
1174 const Numeric lat3 = lat_grid[ilat + 1];
1175 const Numeric lon5 = lon_grid[ilon];
1176 const Numeric lon6 = lon_grid[ilon + 1];
1177 const Numeric re1 = refell2r(refellipsoid, lat1);
1178 const Numeric re3 = refell2r(refellipsoid, lat3);
1179 const Numeric r15 = re1 + z_surf(ilat, ilon);
1180 const Numeric r35 = re3 + z_surf(ilat + 1, ilon);
1181 const Numeric r36 = re3 + z_surf(ilat + 1, ilon + 1);
1182 const Numeric r16 = re1 + z_surf(ilat, ilon + 1);
1183
1185 c1, c2, lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat, lon, aa);
1186}
1187
1203Numeric rslope_crossing3d(const Numeric& rp,
1204 const Numeric& za,
1205 const Numeric& r0,
1206 Numeric c1,
1207 Numeric c2) {
1208 // Even if the four corner radii of the grid cell differ, both c1 and c2 can
1209 // turn up to be 0. So in contrast to the 2D version, here we accept zero c.
1210
1211 // Convert c1 and c2 from degrees to radians
1212 c1 *= RAD2DEG;
1213 c2 *= RAD2DEG * RAD2DEG;
1214
1215 // The nadir angle in radians, and cosine and sine of that angle
1216 const Numeric beta = DEG2RAD * (180 - za);
1217 const Numeric cv = cos(beta);
1218 const Numeric sv = sin(beta);
1219
1220 // Some repeated terms
1221 const Numeric r0s = r0 * sv;
1222 const Numeric r0c = r0 * cv;
1223 const Numeric c1s = c1 * sv;
1224 const Numeric c1c = c1 * cv;
1225 const Numeric c2s = c2 * sv;
1226 const Numeric c2c = c2 * cv;
1227
1228 // The vector of polynomial coefficients
1229 //
1230 Index n = 6;
1231 Vector p0(n + 1);
1232 //
1233 p0[0] = r0s - rp * sv;
1234 p0[1] = r0c + c1s;
1235 p0[2] = -r0s / 2 + c1c + c2s;
1236 p0[3] = -r0c / 6 - c1s / 2 + c2c;
1237 p0[4] = r0s / 24 - c1c / 6 - c2s / 2;
1238 p0[5] = r0c / 120 + c1s / 24 - c2c / 6;
1239 p0[6] = -r0s / 720 + c1c / 120 + c2s / 24;
1240 //
1241 // The accuracy when solving the polynomial equation gets worse when
1242 // approaching 0 and 180 deg. The solution is to let the start polynomial
1243 // order decrease when approaching these angles. The values below based
1244 // on practical experience, don't change without making extremly careful
1245 // tests.
1246 //
1247 if (abs(90 - za) > 89.9) {
1248 n = 1;
1249 } else if (abs(90 - za) > 75) {
1250 n = 4;
1251 }
1252
1253 // Calculate roots of the polynomial
1254 Matrix roots;
1255 int solutionfailure = 1;
1256 //
1257 while (solutionfailure) {
1258 roots.resize(n, 2);
1259 Vector p;
1260 p = p0[Range(0, n + 1)];
1261 solutionfailure = poly_root_solve(roots, p);
1262 if (solutionfailure) {
1263 n -= 1;
1264 ARTS_ASSERT(n > 0);
1265 }
1266 }
1267
1268 // If r0=rp, numerical inaccuracy can give a false solution, very close
1269 // to 0, that we must throw away.
1270 Numeric dmin = 0;
1271 if (r0 == rp) {
1272 dmin = 1e-6;
1273 }
1274
1275 // Find the smallest root with imaginary part = 0, and real part > 0.
1276 //
1277 Numeric dlat = 1.571; // Not interested in solutions above 90 deg!
1278 //
1279 for (Index i = 0; i < n; i++) {
1280 if (roots(i, 1) == 0 && roots(i, 0) > dmin && roots(i, 0) < dlat) {
1281 dlat = roots(i, 0);
1282 }
1283 }
1284
1285 if (dlat < 1.57) // A somewhat smaller value than start one
1286 {
1287 // Convert back to degrees
1288 dlat = RAD2DEG * dlat;
1289 } else {
1290 dlat = LAT_NOT_FOUND;
1291 }
1292
1293 return dlat;
1294}
1295
1327void r_crossing_3d(Numeric& lat,
1328 Numeric& lon,
1329 Numeric& l,
1330 const Numeric& r_hit,
1331 const Numeric& r_start,
1332 const Numeric& lat_start,
1333 const Numeric& lon_start,
1334 const Numeric& za_start,
1335 const Numeric& ppc,
1336 const Numeric& x,
1337 const Numeric& y,
1338 const Numeric& z,
1339 const Numeric& dx,
1340 const Numeric& dy,
1341 const Numeric& dz) {
1342 ARTS_ASSERT(za_start >= 0);
1343 ARTS_ASSERT(za_start <= 180);
1344
1345 // If above and looking upwards or r_hit below tangent point,
1346 // we have no crossing:
1347 if ((r_start >= r_hit && za_start <= 90) || ppc > r_hit) {
1348 lat = LAT_NOT_FOUND;
1349 lon = LON_NOT_FOUND;
1350 l = L_NOT_FOUND;
1351 }
1352
1353 else {
1354 // Zenith/nadir
1355 if (za_start < ANGTOL || za_start > 180 - ANGTOL) {
1356 l = abs(r_hit - r_start);
1357 lat = lat_start;
1358 lon = lon_start;
1359 }
1360
1361 else {
1362 const Numeric p = x * dx + y * dy + z * dz;
1363 const Numeric pp = p * p;
1364 const Numeric q = x * x + y * y + z * z - r_hit * r_hit;
1365 const Numeric sq = sqrt(pp - q);
1366 const Numeric l1 = -p + sq;
1367 const Numeric l2 = -p - sq;
1368
1369 const Numeric lmin = min(l1, l2);
1370 const Numeric lmax = max(l1, l2);
1371
1372 // If r_start equals r_hit solutions just above zero can appear (that
1373 // shall be rejected). So we pick the biggest solution if lmin is
1374 // negative or just above zero.
1375 // (Tried to use "if( r_start != r_hit )", but failed occasionally)
1376 if (lmin < 1e-6) {
1377 l = lmax;
1378 } else {
1379 l = lmin;
1380 }
1381 ARTS_ASSERT(l > 0);
1382
1383 lat = RAD2DEG * asin((z + dz * l) / r_hit);
1384 lon = RAD2DEG * atan2(y + dy * l, x + dx * l);
1385 }
1386 }
1387}
1388
1389/*===========================================================================
1390 === Basic functions for the Ppath structure
1391 ===========================================================================*/
1392
1394 const Index& atmosphere_dim,
1395 const Index& np) {
1396 ARTS_ASSERT(np > 0);
1397 ARTS_ASSERT(atmosphere_dim >= 1);
1398 ARTS_ASSERT(atmosphere_dim <= 3);
1399
1400 ppath.dim = atmosphere_dim;
1401 ppath.np = np;
1402 ppath.constant = -1;
1403
1404 const Index npos = max(Index(2), atmosphere_dim);
1405 const Index nlos = max(Index(1), atmosphere_dim - 1);
1406
1407 ppath.start_pos.resize(npos);
1408 ppath.start_pos = -999;
1409 ppath.start_los.resize(nlos);
1410 ppath.start_los = -999;
1411 ppath.start_lstep = 0;
1412 ppath.end_pos.resize(npos);
1413 ppath.end_los.resize(nlos);
1414 ppath.end_lstep = 0;
1415
1416 ppath.pos.resize(np, npos);
1417 ppath.los.resize(np, nlos);
1418 ppath.r.resize(np);
1419 ppath.lstep.resize(np - 1);
1420
1421 ppath.gp_p.resize(np);
1422 if (atmosphere_dim >= 2) {
1423 ppath.gp_lat.resize(np);
1424 if (atmosphere_dim == 3) {
1425 ppath.gp_lon.resize(np);
1426 }
1427 }
1428
1429 ppath_set_background(ppath, 0);
1430 ppath.nreal.resize(np);
1431 ppath.ngroup.resize(np);
1432}
1433
1434void ppath_set_background(Ppath& ppath, const Index& case_nr) {
1435 switch (case_nr) {
1436 case 0:
1437 ppath.background = "unvalid";
1438 break;
1439 case 1:
1440 ppath.background = "space";
1441 break;
1442 case 2:
1443 ppath.background = "surface";
1444 break;
1445 case 3:
1446 ppath.background = "cloud box level";
1447 break;
1448 case 4:
1449 ppath.background = "cloud box interior";
1450 break;
1451 case 9:
1452 ppath.background = "transmitter";
1453 break;
1454 default:
1456 "Case number ", case_nr, " is not defined.")
1457 }
1458}
1459
1460Index ppath_what_background(const Ppath& ppath) {
1461 if (ppath.background == "unvalid") {
1462 return 0;
1463 } else if (ppath.background == "space") {
1464 return 1;
1465 } else if (ppath.background == "surface") {
1466 return 2;
1467 } else if (ppath.background == "cloud box level") {
1468 return 3;
1469 } else if (ppath.background == "cloud box interior") {
1470 return 4;
1471 } else if (ppath.background == "transmitter") {
1472 return 9;
1473 } else {
1475 "The string ", ppath.background,
1476 " is not a valid background case.")
1477 }
1478}
1479
1480void ppath_copy(Ppath& ppath1, const Ppath& ppath2, const Index& ncopy) {
1481 Index n;
1482 if (ncopy < 0) {
1483 n = ppath2.np;
1484 } else {
1485 n = ncopy;
1486 }
1487
1488 ARTS_ASSERT(ppath1.np >= n);
1489
1490 // The field np shall not be copied !
1491
1492 ppath1.dim = ppath2.dim;
1493 ppath1.constant = ppath2.constant;
1494 ppath1.background = ppath2.background;
1495
1496 // As index 0 always is included in the copying, the end point is covered
1497 ppath1.end_pos = ppath2.end_pos;
1498 ppath1.end_los = ppath2.end_los;
1499 ppath1.end_lstep = ppath2.end_lstep;
1500
1501 // Copy start point only if copying fills up ppath1
1502 if (n == ppath1.np) {
1503 ppath1.start_pos = ppath2.start_pos;
1504 ppath1.start_los = ppath2.start_los;
1505 ppath1.start_lstep = ppath2.start_lstep;
1506 }
1507
1508 ppath1.pos(Range(0, n), joker) = ppath2.pos(Range(0, n), joker);
1509 ppath1.los(Range(0, n), joker) = ppath2.los(Range(0, n), joker);
1510 ppath1.r[Range(0, n)] = ppath2.r[Range(0, n)];
1511 ppath1.nreal[Range(0, n)] = ppath2.nreal[Range(0, n)];
1512 ppath1.ngroup[Range(0, n)] = ppath2.ngroup[Range(0, n)];
1513 if (n > 1) {
1514 ppath1.lstep[Range(0, n - 1)] = ppath2.lstep[Range(0, n - 1)];
1515 }
1516
1517 for (Index i = 0; i < n; i++) {
1518 gridpos_copy(ppath1.gp_p[i], ppath2.gp_p[i]);
1519
1520 if (ppath1.dim >= 2) {
1521 gridpos_copy(ppath1.gp_lat[i], ppath2.gp_lat[i]);
1522 }
1523
1524 if (ppath1.dim == 3) {
1525 gridpos_copy(ppath1.gp_lon[i], ppath2.gp_lon[i]);
1526 }
1527 }
1528}
1529
1546void ppath_append(Ppath& ppath1, const Ppath& ppath2) {
1547 const Index n1 = ppath1.np;
1548 const Index n2 = ppath2.np;
1549
1550 Ppath ppath;
1551 ppath_init_structure(ppath, ppath1.dim, n1);
1552 ppath_copy(ppath, ppath1, -1);
1553
1554 ppath_init_structure(ppath1, ppath1.dim, n1 + n2 - 1);
1555 ppath_copy(ppath1, ppath, -1);
1556
1557 // Append data from ppath2
1558 Index i1;
1559 for (Index i = 1; i < n2; i++) {
1560 i1 = n1 + i - 1;
1561
1562 ppath1.pos(i1, 0) = ppath2.pos(i, 0);
1563 ppath1.pos(i1, 1) = ppath2.pos(i, 1);
1564 ppath1.los(i1, 0) = ppath2.los(i, 0);
1565 ppath1.r[i1] = ppath2.r[i];
1566 ppath1.nreal[i1] = ppath2.nreal[i];
1567 ppath1.ngroup[i1] = ppath2.ngroup[i];
1568 gridpos_copy(ppath1.gp_p[i1], ppath2.gp_p[i]);
1569
1570 if (ppath1.dim >= 2) {
1571 gridpos_copy(ppath1.gp_lat[i1], ppath2.gp_lat[i]);
1572 }
1573
1574 if (ppath1.dim == 3) {
1575 ppath1.pos(i1, 2) = ppath2.pos(i, 2);
1576 ppath1.los(i1, 1) = ppath2.los(i, 1);
1577 gridpos_copy(ppath1.gp_lon[i1], ppath2.gp_lon[i]);
1578 }
1579
1580 ppath1.lstep[i1 - 1] = ppath2.lstep[i - 1];
1581 }
1582
1583 if (ppath_what_background(ppath2)) {
1584 ppath1.background = ppath2.background;
1585 }
1586
1587 ppath.start_pos = ppath2.start_pos;
1588 ppath.start_los = ppath2.start_los;
1589 ppath.start_lstep = ppath2.start_lstep;
1590}
1591
1592/*===========================================================================
1593 === 1D/2D/3D start and end ppath functions
1594 ===========================================================================*/
1595
1606void ppath_start_1d(Numeric& r_start,
1607 Numeric& lat_start,
1608 Numeric& za_start,
1609 Index& ip,
1610 const Ppath& ppath) {
1611 // Number of points in the incoming ppath
1612 const Index imax = ppath.np - 1;
1613
1614 // Extract starting radius, zenith angle and latitude
1615 r_start = ppath.r[imax];
1616 lat_start = ppath.pos(imax, 1);
1617 za_start = ppath.los(imax, 0);
1618
1619 // Determine index of the pressure level being the lower limit for the
1620 // grid range of interest.
1621 //
1622 ip = gridpos2gridrange(ppath.gp_p[imax], za_start <= 90);
1623}
1624
1636 ConstVectorView r_v,
1637 ConstVectorView lat_v,
1638 ConstVectorView za_v,
1639 ConstVectorView lstep,
1640 ConstVectorView n_v,
1641 ConstVectorView ng_v,
1642 ConstVectorView z_field,
1643 ConstVectorView refellipsoid,
1644 const Index& ip,
1645 const Index& endface,
1646 const Numeric& ppc) {
1647 // Number of path points
1648 const Index np = r_v.nelem();
1649
1650 // Re-allocate ppath for return results and fill the structure
1651 ppath_init_structure(ppath, 1, np);
1652
1653 ppath.constant = ppc;
1654
1655 // Help variables that are common for all points.
1656 const Numeric r1 = refellipsoid[0] + z_field[ip];
1657 const Numeric dr = z_field[ip + 1] - z_field[ip];
1658
1659 for (Index i = 0; i < np; i++) {
1660 ppath.r[i] = r_v[i];
1661 ppath.pos(i, 0) = r_v[i] - refellipsoid[0];
1662 ppath.pos(i, 1) = lat_v[i];
1663 ppath.los(i, 0) = za_v[i];
1664 ppath.nreal[i] = n_v[i];
1665 ppath.ngroup[i] = ng_v[i];
1666
1667 ppath.gp_p[i].idx = ip;
1668 ppath.gp_p[i].fd[0] = (r_v[i] - r1) / dr;
1669 ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
1670 gridpos_check_fd(ppath.gp_p[i]);
1671
1672 if (i > 0) {
1673 ppath.lstep[i - 1] = lstep[i - 1];
1674 }
1675 }
1676 gridpos_check_fd(ppath.gp_p[np - 1]);
1677
1678 //--- End point is the surface
1679 if (endface == 7) {
1680 ppath_set_background(ppath, 2);
1681 }
1682
1683 //--- End point is on top of a pressure level
1684 else if (endface <= 4) {
1685 gridpos_force_end_fd(ppath.gp_p[np - 1], z_field.nelem());
1686 }
1687}
1688
1699void ppath_start_2d(Numeric& r_start,
1700 Numeric& lat_start,
1701 Numeric& za_start,
1702 Index& ip,
1703 Index& ilat,
1704 Numeric& lat1,
1705 Numeric& lat3,
1706 Numeric& r1a,
1707 Numeric& r3a,
1708 Numeric& r3b,
1709 Numeric& r1b,
1710 Numeric& rsurface1,
1711 Numeric& rsurface3,
1712 Ppath& ppath,
1713 ConstVectorView lat_grid,
1714 ConstMatrixView z_field,
1715 ConstVectorView refellipsoid,
1716 ConstVectorView z_surface) {
1717 // Number of points in the incoming ppath
1718 const Index imax = ppath.np - 1;
1719
1720 // Extract starting radius, zenith angle and latitude
1721 r_start = ppath.r[imax];
1722 lat_start = ppath.pos(imax, 1);
1723 za_start = ppath.los(imax, 0);
1724
1725 // Determine interesting latitude grid range and latitude end points of
1726 // the range.
1727 //
1728 ilat = gridpos2gridrange(ppath.gp_lat[imax], za_start >= 0);
1729 //
1730 lat1 = lat_grid[ilat];
1731 lat3 = lat_grid[ilat + 1];
1732
1733 // Determine interesting pressure grid range. Do this first assuming that
1734 // the pressure levels are not tilted (that is, abs(za_start<=90) always
1735 // mean upward observation).
1736 // Set radius for the corners of the grid cell and the radial slope of
1737 // pressure level limits of the grid cell to match the found ip.
1738 //
1739 ip = gridpos2gridrange(ppath.gp_p[imax], abs(za_start) <= 90);
1740 //
1741 const Numeric re1 = refell2r(refellipsoid, lat_grid[ilat]);
1742 const Numeric re3 = refell2r(refellipsoid, lat_grid[ilat + 1]);
1743 //
1744 r1a = re1 + z_field(ip, ilat); // lower-left
1745 r3a = re3 + z_field(ip, ilat + 1); // lower-right
1746 r3b = re3 + z_field(ip + 1, ilat + 1); // upper-right
1747 r1b = re1 + z_field(ip + 1, ilat); // upper-left
1748
1749 // This part is a fix to catch start postions on top of a pressure level
1750 // that does not have an end fractional distance for the first step.
1751 {
1752 // Radius of lower and upper pressure level at the start position
1753 const Numeric rlow = rsurf_at_lat(lat1, lat3, r1a, r3a, lat_start);
1754 const Numeric rupp = rsurf_at_lat(lat1, lat3, r1b, r3b, lat_start);
1755 if (abs(r_start - rlow) < RTOL || abs(r_start - rupp) < RTOL) {
1756 gridpos_force_end_fd(ppath.gp_p[imax], z_field.nrows());
1757 }
1758 }
1759
1760 // Slopes of pressure levels
1761 Numeric c2, c4;
1762 plevel_slope_2d(c2, lat1, lat3, r1a, r3a);
1763 plevel_slope_2d(c4, lat1, lat3, r1b, r3b);
1764
1765 // Check if the LOS zenith angle happen to be between 90 and the zenith angle
1766 // of the pressure level (that is, 90 + tilt of pressure level), and in
1767 // that case if ip must be changed. This check is only needed when the
1768 // start point is on a pressure level.
1769 //
1770 if (is_gridpos_at_index_i(ppath.gp_p[imax], ip)) {
1771 Numeric tilt = plevel_angletilt(r_start, c2);
1772
1773 if (is_los_downwards(za_start, tilt)) {
1774 ip--;
1775 r1b = r1a;
1776 r3b = r3a;
1777 r1a = re1 + z_field(ip, ilat);
1778 r3a = re3 + z_field(ip, ilat + 1);
1779 plevel_slope_2d(c2, lat1, lat3, r1a, r3a);
1780 }
1781 } else if (is_gridpos_at_index_i(ppath.gp_p[imax], ip + 1)) {
1782 Numeric tilt = plevel_angletilt(r_start, c4);
1783
1784 if (!is_los_downwards(za_start, tilt)) {
1785 ip++;
1786 r1a = r1b;
1787 r3a = r3b;
1788 r3b = re3 + z_field(ip + 1, ilat + 1);
1789 r1b = re1 + z_field(ip + 1, ilat);
1790 plevel_slope_2d(c4, lat1, lat3, r1b, r3b);
1791 }
1792 }
1793
1794 // Surface radius at latitude end points
1795 rsurface1 = re1 + z_surface[ilat];
1796 rsurface3 = re3 + z_surface[ilat + 1];
1797}
1798
1810 ConstVectorView r_v,
1811 ConstVectorView lat_v,
1812 ConstVectorView za_v,
1813 ConstVectorView lstep,
1814 ConstVectorView n_v,
1815 ConstVectorView ng_v,
1816 ConstVectorView lat_grid,
1817 ConstMatrixView z_field,
1818 ConstVectorView refellipsoid,
1819 const Index& ip,
1820 const Index& ilat,
1821 const Index& endface,
1822 const Numeric& ppc) {
1823 // Number of path points
1824 const Index np = r_v.nelem();
1825 const Index imax = np - 1;
1826
1827 // Re-allocate ppath for return results and fill the structure
1828 //
1829 ppath_init_structure(ppath, 2, np);
1830
1831 ppath.constant = ppc;
1832
1833 // Help variables that are common for all points.
1834 const Numeric dlat = lat_grid[ilat + 1] - lat_grid[ilat];
1835 const Numeric z1low = z_field(ip, ilat);
1836 const Numeric z1upp = z_field(ip + 1, ilat);
1837 const Numeric dzlow = z_field(ip, ilat + 1) - z1low;
1838 const Numeric dzupp = z_field(ip + 1, ilat + 1) - z1upp;
1839 Numeric re = refell2r(refellipsoid, lat_grid[ilat]);
1840 const Numeric r1low = re + z1low;
1841 const Numeric r1upp = re + z1upp;
1842 re = refell2r(refellipsoid, lat_grid[ilat + 1]);
1843 const Numeric drlow = re + z_field(ip, ilat + 1) - r1low;
1844 const Numeric drupp = re + z_field(ip + 1, ilat + 1) - r1upp;
1845
1846 for (Index i = 0; i < np; i++) {
1847 ppath.r[i] = r_v[i];
1848 ppath.pos(i, 1) = lat_v[i];
1849 ppath.los(i, 0) = za_v[i];
1850 ppath.nreal[i] = n_v[i];
1851 ppath.ngroup[i] = ng_v[i];
1852
1853 // Weight in the latitude direction
1854 Numeric w = (lat_v[i] - lat_grid[ilat]) / dlat;
1855
1856 // Radius of lower and upper face at present latitude
1857 const Numeric rlow = r1low + w * drlow;
1858 const Numeric rupp = r1upp + w * drupp;
1859
1860 // Geometrical altitude of lower and upper face at present latitude
1861 const Numeric zlow = z1low + w * dzlow;
1862 const Numeric zupp = z1upp + w * dzupp;
1863
1864 ppath.gp_p[i].idx = ip;
1865 ppath.gp_p[i].fd[0] = (r_v[i] - rlow) / (rupp - rlow);
1866 ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
1867 gridpos_check_fd(ppath.gp_p[i]);
1868
1869 ppath.pos(i, 0) = zlow + ppath.gp_p[i].fd[0] * (zupp - zlow);
1870
1871 ppath.gp_lat[i].idx = ilat;
1872 ppath.gp_lat[i].fd[0] = (lat_v[i] - lat_grid[ilat]) / dlat;
1873 ppath.gp_lat[i].fd[1] = 1 - ppath.gp_lat[i].fd[0];
1874 gridpos_check_fd(ppath.gp_lat[i]);
1875
1876 if (i > 0) {
1877 ppath.lstep[i - 1] = lstep[i - 1];
1878 }
1879 }
1880 gridpos_check_fd(ppath.gp_p[imax]);
1881 gridpos_check_fd(ppath.gp_lat[imax]);
1882
1883 // Do end-face specific tasks
1884 if (endface == 7) {
1885 ppath_set_background(ppath, 2);
1886 }
1887
1888 // Set fractional distance for end point
1889 //
1890 if (endface == 1 || endface == 3) {
1891 gridpos_force_end_fd(ppath.gp_lat[np - 1], lat_grid.nelem());
1892 } else if (endface == 2 || endface == 4) {
1893 gridpos_force_end_fd(ppath.gp_p[np - 1], z_field.nrows());
1894 }
1895
1896 // Handle cases where exactly a corner is hit, or when slipping outside of
1897 // the grid box due to numerical inaccuarcy
1898 if (ppath.gp_p[imax].fd[0] < 0 || ppath.gp_p[imax].fd[1] < 0) {
1899 gridpos_force_end_fd(ppath.gp_p[imax], z_field.nrows());
1900 }
1901 if (ppath.gp_lat[imax].fd[0] < 0 || ppath.gp_lat[imax].fd[1] < 0) {
1902 gridpos_force_end_fd(ppath.gp_lat[imax], lat_grid.nelem());
1903 }
1904}
1905
1916void ppath_start_3d(Numeric& r_start,
1917 Numeric& lat_start,
1918 Numeric& lon_start,
1919 Numeric& za_start,
1920 Numeric& aa_start,
1921 Index& ip,
1922 Index& ilat,
1923 Index& ilon,
1924 Numeric& lat1,
1925 Numeric& lat3,
1926 Numeric& lon5,
1927 Numeric& lon6,
1928 Numeric& r15a,
1929 Numeric& r35a,
1930 Numeric& r36a,
1931 Numeric& r16a,
1932 Numeric& r15b,
1933 Numeric& r35b,
1934 Numeric& r36b,
1935 Numeric& r16b,
1936 Numeric& rsurface15,
1937 Numeric& rsurface35,
1938 Numeric& rsurface36,
1939 Numeric& rsurface16,
1940 Ppath& ppath,
1941 ConstVectorView lat_grid,
1942 ConstVectorView lon_grid,
1943 ConstTensor3View z_field,
1944 ConstVectorView refellipsoid,
1945 ConstMatrixView z_surface) {
1946 // Index of last point in the incoming ppath
1947 const Index imax = ppath.np - 1;
1948
1949 // Extract starting radius, zenith angle and latitude
1950 r_start = ppath.r[imax];
1951 lat_start = ppath.pos(imax, 1);
1952 lon_start = ppath.pos(imax, 2);
1953 za_start = ppath.los(imax, 0);
1954 aa_start = ppath.los(imax, 1);
1955
1956 // Number of lat/lon
1957 const Index nlat = lat_grid.nelem();
1958 const Index nlon = lon_grid.nelem();
1959
1960 // Lower index of lat and lon ranges of interest
1961 //
1962 // The longitude is undefined at the poles and as the azimuth angle
1963 // is defined in other way at the poles.
1964 //
1965 if (lat_start == 90) {
1966 ilat = nlat - 2;
1967 GridPos gp_tmp;
1968 gridpos(gp_tmp, lon_grid, aa_start);
1969 if (aa_start < 180) {
1970 ilon = gridpos2gridrange(gp_tmp, 1);
1971 } else {
1972 ilon = gridpos2gridrange(gp_tmp, 0);
1973 }
1974 } else if (lat_start == -90) {
1975 ilat = 0;
1976 GridPos gp_tmp;
1977 gridpos(gp_tmp, lon_grid, aa_start);
1978 if (aa_start < 180) {
1979 ilon = gridpos2gridrange(gp_tmp, 1);
1980 } else {
1981 ilon = gridpos2gridrange(gp_tmp, 0);
1982 }
1983 } else {
1984 if (lat_start > 0) {
1985 ilat = gridpos2gridrange(ppath.gp_lat[imax], abs(aa_start) < 90);
1986 } else {
1987 ilat = gridpos2gridrange(ppath.gp_lat[imax], abs(aa_start) <= 90);
1988 }
1989 if (lon_start < lon_grid[nlon - 1]) {
1990 ilon = gridpos2gridrange(ppath.gp_lon[imax], aa_start >= 0);
1991 } else {
1992 ilon = nlon - 2;
1993 }
1994 }
1995 //
1996 lat1 = lat_grid[ilat];
1997 lat3 = lat_grid[ilat + 1];
1998 lon5 = lon_grid[ilon];
1999 lon6 = lon_grid[ilon + 1];
2000
2001 // Determine interesting pressure grid range. Do this first assuming that
2002 // the pressure levels are not tilted (that is, abs(za_start<=90) always
2003 // mean upward observation).
2004 // Set radius for the corners of the grid cell and the radial slope of
2005 // pressure level limits of the grid cell to match the found ip.
2006 //
2007 ip = gridpos2gridrange(ppath.gp_p[imax], za_start <= 90);
2008 //
2009 const Numeric re1 = refell2r(refellipsoid, lat_grid[ilat]);
2010 const Numeric re3 = refell2r(refellipsoid, lat_grid[ilat + 1]);
2011 //
2012 r15a = re1 + z_field(ip, ilat, ilon);
2013 r35a = re3 + z_field(ip, ilat + 1, ilon);
2014 r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2015 r16a = re1 + z_field(ip, ilat, ilon + 1);
2016 r15b = re1 + z_field(ip + 1, ilat, ilon);
2017 r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2018 r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2019 r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2020
2021 // Check if the LOS zenith angle happen to be between 90 and the zenith angle
2022 // of the pressure level (that is, 90 + tilt of pressure level), and in
2023 // that case if ip must be changed. This check is only needed when the
2024 // start point is on a pressure level.
2025 //
2026 if (fabs(za_start - 90) <= 10) // To save time. Ie. max tilt assumed =10 deg
2027 {
2028 if (is_gridpos_at_index_i(ppath.gp_p[imax], ip)) {
2029 // Slope and angular tilt of lower pressure level
2030 Numeric c2a, c2b;
2031 plevel_slope_3d(c2a,
2032 c2b,
2033 lat1,
2034 lat3,
2035 lon5,
2036 lon6,
2037 r15a,
2038 r35a,
2039 r36a,
2040 r16a,
2041 lat_start,
2042 lon_start,
2043 aa_start);
2044 Numeric tilt = plevel_angletilt(r_start, c2a);
2045 // Negelect very small tilts, likely caused by numerical problems
2046 if (abs(tilt) > 1e-4 && is_los_downwards(za_start, tilt)) {
2047 ip--;
2048 r15b = r15a;
2049 r35b = r35a;
2050 r36b = r36a;
2051 r16b = r16a;
2052 r15a = re1 + z_field(ip, ilat, ilon);
2053 r35a = re3 + z_field(ip, ilat + 1, ilon);
2054 r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2055 r16a = re1 + z_field(ip, ilat, ilon + 1);
2056 }
2057 } else if (is_gridpos_at_index_i(ppath.gp_p[imax], ip + 1)) {
2058 // Slope and angular tilt of upper pressure level
2059 Numeric c4a, c4b;
2060 plevel_slope_3d(c4a,
2061 c4b,
2062 lat1,
2063 lat3,
2064 lon5,
2065 lon6,
2066 r15b,
2067 r35b,
2068 r36b,
2069 r16b,
2070 lat_start,
2071 lon_start,
2072 aa_start);
2073 Numeric tilt = plevel_angletilt(r_start, c4a);
2074 //
2075 if (!is_los_downwards(za_start, tilt)) {
2076 ip++;
2077 r15a = r15b;
2078 r35a = r35b;
2079 r36a = r36b;
2080 r16a = r16b;
2081 r15b = re1 + z_field(ip + 1, ilat, ilon);
2082 r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2083 r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2084 r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2085 }
2086 }
2087 }
2088
2089 // Surface radius at latitude/longitude corner points
2090 rsurface15 = re1 + z_surface(ilat, ilon);
2091 rsurface35 = re3 + z_surface(ilat + 1, ilon);
2092 rsurface36 = re3 + z_surface(ilat + 1, ilon + 1);
2093 rsurface16 = re1 + z_surface(ilat, ilon + 1);
2094}
2095
2107 ConstVectorView r_v,
2108 ConstVectorView lat_v,
2109 ConstVectorView lon_v,
2110 ConstVectorView za_v,
2111 ConstVectorView aa_v,
2112 ConstVectorView lstep,
2113 ConstVectorView n_v,
2114 ConstVectorView ng_v,
2115 ConstVectorView lat_grid,
2116 ConstVectorView lon_grid,
2117 ConstTensor3View z_field,
2118 ConstVectorView refellipsoid,
2119 const Index& ip,
2120 const Index& ilat,
2121 const Index& ilon,
2122 const Index& endface,
2123 const Numeric& ppc) {
2124 // Number of path points
2125 const Index np = r_v.nelem();
2126 const Index imax = np - 1;
2127
2128 // Re-allocate ppath for return results and fill the structure
2129 //
2130 ppath_init_structure(ppath, 3, np);
2131 //
2132 ppath.constant = ppc;
2133
2134 // Help variables that are common for all points.
2135 const Numeric lat1 = lat_grid[ilat];
2136 const Numeric lat3 = lat_grid[ilat + 1];
2137 const Numeric lon5 = lon_grid[ilon];
2138 const Numeric lon6 = lon_grid[ilon + 1];
2139 const Numeric re1 = refell2r(refellipsoid, lat_grid[ilat]);
2140 const Numeric re3 = refell2r(refellipsoid, lat_grid[ilat + 1]);
2141 const Numeric r15a = re1 + z_field(ip, ilat, ilon);
2142 const Numeric r35a = re3 + z_field(ip, ilat + 1, ilon);
2143 const Numeric r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2144 const Numeric r16a = re1 + z_field(ip, ilat, ilon + 1);
2145 const Numeric r15b = re1 + z_field(ip + 1, ilat, ilon);
2146 const Numeric r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2147 const Numeric r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2148 const Numeric r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2149 const Numeric dlat = lat3 - lat1;
2150 const Numeric dlon = lon6 - lon5;
2151
2152 for (Index i = 0; i < np; i++) {
2153 // Radius of pressure levels at present lat and lon
2154 const Numeric rlow = rsurf_at_latlon(
2155 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_v[i], lon_v[i]);
2156 const Numeric rupp = rsurf_at_latlon(
2157 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_v[i], lon_v[i]);
2158
2159 // Position and LOS
2160 ppath.r[i] = r_v[i];
2161 ppath.pos(i, 1) = lat_v[i];
2162 ppath.pos(i, 2) = lon_v[i];
2163 ppath.los(i, 0) = za_v[i];
2164 ppath.los(i, 1) = aa_v[i];
2165 ppath.nreal[i] = n_v[i];
2166 ppath.ngroup[i] = ng_v[i];
2167
2168 // Pressure grid index
2169 ppath.gp_p[i].idx = ip;
2170 ppath.gp_p[i].fd[0] = (r_v[i] - rlow) / (rupp - rlow);
2171 ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
2172 gridpos_check_fd(ppath.gp_p[i]);
2173
2174 // Geometrical altitude
2175 const Numeric re = rsurf_at_latlon(
2176 lat1, lat3, lon5, lon6, re1, re3, re3, re1, lat_v[i], lon_v[i]);
2177 const Numeric zlow = rlow - re;
2178 const Numeric zupp = rupp - re;
2179 //
2180 ppath.pos(i, 0) = zlow + ppath.gp_p[i].fd[0] * (zupp - zlow);
2181
2182 // Latitude grid index
2183 ppath.gp_lat[i].idx = ilat;
2184 ppath.gp_lat[i].fd[0] = (lat_v[i] - lat1) / dlat;
2185 ppath.gp_lat[i].fd[1] = 1 - ppath.gp_lat[i].fd[0];
2186 gridpos_check_fd(ppath.gp_lat[i]);
2187
2188 // Longitude grid index
2189 //
2190 // The longitude is undefined at the poles. The grid index is set to
2191 // the start point.
2192 //
2193 if (abs(lat_v[i]) < POLELAT) {
2194 ppath.gp_lon[i].idx = ilon;
2195 ppath.gp_lon[i].fd[0] = (lon_v[i] - lon5) / dlon;
2196 ppath.gp_lon[i].fd[1] = 1 - ppath.gp_lon[i].fd[0];
2197 gridpos_check_fd(ppath.gp_lon[i]);
2198 } else {
2199 ppath.gp_lon[i].idx = 0;
2200 ppath.gp_lon[i].fd[0] = 0;
2201 ppath.gp_lon[i].fd[1] = 1;
2202 }
2203
2204 if (i > 0) {
2205 ppath.lstep[i - 1] = lstep[i - 1];
2206 }
2207 }
2208
2209 // Do end-face specific tasks
2210 if (endface == 7) {
2211 ppath_set_background(ppath, 2);
2212 }
2213
2214 // Set fractional distance for end point
2215 //
2216 if (endface == 1 || endface == 3) {
2217 gridpos_force_end_fd(ppath.gp_lat[imax], lat_grid.nelem());
2218 } else if (endface == 2 || endface == 4) {
2219 gridpos_force_end_fd(ppath.gp_p[imax], z_field.npages());
2220 } else if (endface == 5 || endface == 6) {
2221 gridpos_force_end_fd(ppath.gp_lon[imax], lon_grid.nelem());
2222 }
2223
2224 // Handle cases where exactly a corner is hit, or when slipping outside of
2225 // the grid box due to numerical inaccuarcy
2226 if (ppath.gp_p[imax].fd[0] < 0 || ppath.gp_p[imax].fd[1] < 0) {
2227 gridpos_force_end_fd(ppath.gp_p[imax], z_field.npages());
2228 }
2229 if (ppath.gp_lat[imax].fd[0] < 0 || ppath.gp_lat[imax].fd[1] < 0) {
2230 gridpos_force_end_fd(ppath.gp_lat[imax], lat_grid.nelem());
2231 }
2232 if (ppath.gp_lon[imax].fd[0] < 0 || ppath.gp_lon[imax].fd[1] < 0) {
2233 gridpos_force_end_fd(ppath.gp_lon[imax], lon_grid.nelem());
2234 }
2235}
2236
2237/*===========================================================================
2238 === Core functions for geometrical ppath_step calculations
2239 ===========================================================================*/
2240
2265void do_gridrange_1d(Vector& r_v,
2266 Vector& lat_v,
2267 Vector& za_v,
2268 Numeric& lstep,
2269 Index& endface,
2270 const Numeric& r_start0,
2271 const Numeric& lat_start,
2272 const Numeric& za_start,
2273 const Numeric& ppc,
2274 const Numeric& lmax,
2275 const Numeric& ra,
2276 const Numeric& rb,
2277 const Numeric& rsurface) {
2278 Numeric r_start = r_start0;
2279
2280 ARTS_ASSERT(rb > ra);
2281 ARTS_ASSERT(r_start >= ra - RTOL);
2282 ARTS_ASSERT(r_start <= rb + RTOL);
2283
2284 // Shift radius if outside
2285 if (r_start < ra) {
2286 r_start = ra;
2287 } else if (r_start > rb) {
2288 r_start = rb;
2289 }
2290
2291 Numeric r_end;
2292 bool tanpoint = false;
2293 endface = -1;
2294
2295 // If upward, end point radius is always rb
2296 if (za_start <= 90) {
2297 endface = 4;
2298 r_end = rb;
2299 }
2300
2301 else {
2302 // Path reaches ra:
2303 if (ra > rsurface && ra > ppc) {
2304 endface = 2;
2305 r_end = ra;
2306 }
2307
2308 // Path reaches the surface:
2309 else if (rsurface > ppc) {
2310 endface = 7;
2311 r_end = rsurface;
2312 }
2313
2314 // The remaining option is a tangent point and back to rb
2315 else {
2316 endface = 4;
2317 r_end = rb;
2318 tanpoint = true;
2319 }
2320 }
2321
2322 ARTS_ASSERT(endface > 0);
2323
2325 lat_v,
2326 za_v,
2327 lstep,
2328 ppc,
2329 r_start,
2330 lat_start,
2331 za_start,
2332 r_end,
2333 tanpoint,
2334 lmax);
2335}
2336
2338 ConstVectorView z_field,
2339 ConstVectorView refellipsoid,
2340 const Numeric& z_surface,
2341 const Numeric& lmax) {
2342 // Starting radius, zenith angle and latitude
2343 Numeric r_start, lat_start, za_start;
2344
2345 // Index of the pressure level being the lower limit for the
2346 // grid range of interest.
2347 Index ip;
2348
2349 // Determine the variables defined above, and make asserts of input
2350 ppath_start_1d(r_start, lat_start, za_start, ip, ppath);
2351
2352 // If the field "constant" is negative, this is the first call of the
2353 // function and the path constant shall be calculated.
2354 Numeric ppc;
2355 if (ppath.constant < 0) {
2356 ppc = geometrical_ppc(r_start, za_start);
2357 } else {
2358 ppc = ppath.constant;
2359 }
2360
2361 // The path is determined by another function. Determine some variables
2362 // needed by that function and call the function.
2363 //
2364 // Vars to hold found path points, path step length and coding for end face
2365 Vector r_v, lat_v, za_v;
2366 Numeric lstep;
2367 Index endface;
2368 //
2369 do_gridrange_1d(r_v,
2370 lat_v,
2371 za_v,
2372 lstep,
2373 endface,
2374 r_start,
2375 lat_start,
2376 za_start,
2377 ppc,
2378 lmax,
2379 refellipsoid[0] + z_field[ip],
2380 refellipsoid[0] + z_field[ip + 1],
2381 refellipsoid[0] + z_surface);
2382
2383 // Fill *ppath*
2384 const Index np = r_v.nelem();
2385 ppath_end_1d(ppath,
2386 r_v,
2387 lat_v,
2388 za_v,
2389 Vector(np - 1, lstep),
2390 Vector(np, 1),
2391 Vector(np, 1),
2392 z_field,
2393 refellipsoid,
2394 ip,
2395 endface,
2396 ppc);
2397}
2398
2404void do_gridcell_2d_byltest(Vector& r_v,
2405 Vector& lat_v,
2406 Vector& za_v,
2407 Numeric& lstep,
2408 Index& endface,
2409 const Numeric& r_start0,
2410 const Numeric& lat_start0,
2411 const Numeric& za_start,
2412 const Numeric& l_start,
2413 const Index& icall,
2414 const Numeric& ppc,
2415 const Numeric& lmax,
2416 const Numeric& lat1,
2417 const Numeric& lat3,
2418 const Numeric& r1a,
2419 const Numeric& r3a,
2420 const Numeric& r3b,
2421 const Numeric& r1b,
2422 const Numeric& rsurface1,
2423 const Numeric& rsurface3) {
2424 Numeric r_start = r_start0;
2425 Numeric lat_start = lat_start0;
2426
2427 ARTS_ASSERT(icall < 10);
2428
2429 // Assert latitude and longitude
2430 ARTS_ASSERT(lat_start >= lat1 - LATLONTOL);
2431 ARTS_ASSERT(lat_start <= lat3 + LATLONTOL);
2432
2433 // Shift latitude and longitude if outside
2434 if (lat_start < lat1) {
2435 lat_start = lat1;
2436 } else if (lat_start > lat3) {
2437 lat_start = lat3;
2438 }
2439
2440 // Radius of lower and upper pressure level at the start position
2441 Numeric rlow = rsurf_at_lat(lat1, lat3, r1a, r3a, lat_start);
2442 Numeric rupp = rsurf_at_lat(lat1, lat3, r1b, r3b, lat_start);
2443
2444 // Assert radius (some extra tolerance is needed for radius)
2445 ARTS_ASSERT(r_start >= rlow - RTOL);
2446 ARTS_ASSERT(r_start <= rupp + RTOL);
2447
2448 // Shift radius if outside
2449 if (r_start < rlow) {
2450 r_start = rlow;
2451 } else if (r_start > rupp) {
2452 r_start = rupp;
2453 }
2454
2455 // Position and LOS in cartesian coordinates
2456 Numeric x, z, dx, dz;
2457 poslos2cart(x, z, dx, dz, r_start, lat_start, za_start);
2458
2459 // Some booleans to check for recursive call
2460 bool unsafe = false;
2461 bool do_surface = false;
2462
2463 // Determine the position of the end point
2464 //
2465 endface = 0;
2466 //
2467 Numeric r_end, lat_end, l_end;
2468
2469 // Zenith and nadir looking are handled as special cases
2470 const Numeric absza = abs(za_start);
2471
2472 // Zenith looking
2473 if (absza < ANGTOL) {
2474 l_end = rupp - r_start;
2475 endface = 4;
2476 }
2477
2478 // Nadir looking
2479 else if (absza > 180 - ANGTOL) {
2480 const Numeric rsurface =
2481 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_start);
2482
2483 if (rlow > rsurface) {
2484 l_end = r_start - rlow;
2485 endface = 2;
2486 } else {
2487 l_end = r_start - rsurface;
2488 endface = 7;
2489 }
2490 }
2491
2492 else {
2493 unsafe = true;
2494
2495 // Calculate correction terms for the position to compensate for
2496 // numerical inaccuracy.
2497 //
2498 Numeric r_corr, lat_corr;
2499 //
2500 cart2pol(r_corr, lat_corr, x, z, lat_start, za_start);
2501 //
2502 r_corr -= r_start;
2503 lat_corr -= lat_start;
2504
2505 // The end point is found by testing different step lengths until the
2506 // step length has been determined by a precision of *l_acc*.
2507 //
2508 // The first step is to found a length that goes outside the grid cell,
2509 // to find an upper limit. The lower limit is set to 0. The upper
2510 // limit is either set to l_start or it is estimated from the size of
2511 // the grid cell.
2512 // The search algorith is bisection, the next length to test is the
2513 // mean value of the minimum and maximum length limits.
2514 //
2515 if (l_start > 0) {
2516 l_end = l_start;
2517 } else {
2518 l_end = 2 * (rupp - rlow);
2519 }
2520 //
2521 Numeric l_in = 0, l_out = l_end;
2522 bool ready = false, startup = true;
2523
2524 if (rsurface1 + RTOL >= r1a || rsurface3 + RTOL >= r3a) {
2525 do_surface = true;
2526 }
2527
2528 while (!ready) {
2529 cart2pol(
2530 r_end, lat_end, x + dx * l_end, z + dz * l_end, lat_start, za_start);
2531 r_end -= r_corr;
2532 lat_end -= lat_corr;
2533
2534 bool inside = true;
2535
2536 rlow = rsurf_at_lat(lat1, lat3, r1a, r3a, lat_end);
2537
2538 if (do_surface) {
2539 const Numeric r_surface =
2540 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_end);
2541 if (r_surface >= rlow && r_end <= r_surface) {
2542 inside = false;
2543 endface = 7;
2544 }
2545 }
2546
2547 if (inside) {
2548 if (lat_end < lat1) {
2549 inside = false;
2550 endface = 1;
2551 } else if (lat_end > lat3) {
2552 inside = false;
2553 endface = 3;
2554 } else if (r_end < rlow) {
2555 inside = false;
2556 endface = 2;
2557 } else {
2558 rupp = rsurf_at_lat(lat1, lat3, r1b, r3b, lat_end);
2559 if (r_end > rupp) {
2560 inside = false;
2561 endface = 4;
2562 }
2563 }
2564 }
2565
2566 if (startup) {
2567 if (inside) {
2568 l_in = l_end;
2569 l_end *= 5;
2570 } else {
2571 l_out = l_end;
2572 l_end = (l_out + l_in) / 2;
2573 startup = false;
2574 }
2575 } else {
2576 if (inside) {
2577 l_in = l_end;
2578 } else {
2579 l_out = l_end;
2580 }
2581
2582 if ((l_out - l_in) < LACC) {
2583 ready = true;
2584 } else {
2585 l_end = (l_out + l_in) / 2;
2586 }
2587 }
2588 }
2589 }
2590
2591 //--- Create return vectors
2592 //
2593 Index n = 1;
2594 //
2595 if (lmax > 0) {
2596 n = Index(ceil(abs(l_end / lmax)));
2597 if (n < 1) {
2598 n = 1;
2599 }
2600 }
2601 //
2602 r_v.resize(n + 1);
2603 lat_v.resize(n + 1);
2604 za_v.resize(n + 1);
2605 //
2606 r_v[0] = r_start;
2607 lat_v[0] = lat_start;
2608 za_v[0] = za_start;
2609 //
2610 lstep = l_end / (Numeric)n;
2611 Numeric l;
2612 bool ready = true;
2613 //
2614 for (Index j = 1; j <= n; j++) {
2615 l = lstep * (Numeric)j;
2616 cart2poslos(r_v[j],
2617 lat_v[j],
2618 za_v[j],
2619 x + dx * l,
2620 z + dz * l,
2621 dx,
2622 dz,
2623 ppc,
2624 lat_start,
2625 za_start);
2626
2627 if (j < n) {
2628 if (unsafe) {
2629 // Check that r_v[j] is above lower pressure level and the
2630 // surface. This can fail around tangent points. For p-levels
2631 // with constant r this is easy to handle analytically, but the
2632 // problem is tricky in the general case with a non-spherical
2633 // geometry, and this crude solution is used instead. Not the
2634 // most elegant solution, but it works! Added later the same
2635 // check for upper level, after getting assert in that direction.
2636 // The z_field was crazy, but still formerly correct.
2637 rlow = rsurf_at_lat(lat1, lat3, r1a, r3a, lat_v[j]);
2638 if (do_surface) {
2639 const Numeric r_surface =
2640 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_v[j]);
2641 const Numeric r_test = max(r_surface, rlow);
2642 if (r_v[j] < r_test) {
2643 ready = false;
2644 break;
2645 }
2646 } else if (r_v[j] < rlow) {
2647 ready = false;
2648 break;
2649 }
2650
2651 rupp = rsurf_at_lat(lat1, lat3, r1b, r3b, lat_v[j]);
2652 if (r_v[j] > rupp) {
2653 ready = false;
2654 break;
2655 }
2656 }
2657 } else // j==n
2658 {
2659 if (unsafe) {
2660 // Set end point to be consistent with found endface.
2661 //
2662 if (endface == 1) {
2663 lat_v[n] = lat1;
2664 } else if (endface == 2) {
2665 r_v[n] = rsurf_at_lat(lat1, lat3, r1a, r3a, lat_v[n]);
2666 } else if (endface == 3) {
2667 lat_v[n] = lat3;
2668 } else if (endface == 4) {
2669 r_v[n] = rsurf_at_lat(lat1, lat3, r1b, r3b, lat_v[n]);
2670 } else if (endface == 7) {
2671 r_v[n] = rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_v[n]);
2672 }
2673 }
2674 }
2675 }
2676
2677 if (!ready) { // If an "outside" point found, restart with l as start search length
2679 lat_v,
2680 za_v,
2681 lstep,
2682 endface,
2683 r_start,
2684 lat_start,
2685 za_start,
2686 l,
2687 icall + 1,
2688 ppc,
2689 lmax,
2690 lat1,
2691 lat3,
2692 r1a,
2693 r3a,
2694 r3b,
2695 r1b,
2696 rsurface1,
2697 rsurface3);
2698 }
2699}
2700
2702 ConstVectorView lat_grid,
2703 ConstMatrixView z_field,
2704 ConstVectorView refellipsoid,
2705 ConstVectorView z_surface,
2706 const Numeric& lmax) {
2707 // Radius, zenith angle and latitude of start point.
2708 Numeric r_start, lat_start, za_start;
2709
2710 // Lower grid index for the grid cell of interest.
2711 Index ip, ilat;
2712
2713 // Radii and latitudes set by *ppath_start_2d*.
2714 Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
2715
2716 // Determine the variables defined above and make all possible asserts
2717 ppath_start_2d(r_start,
2718 lat_start,
2719 za_start,
2720 ip,
2721 ilat,
2722 lat1,
2723 lat3,
2724 r1a,
2725 r3a,
2726 r3b,
2727 r1b,
2728 rsurface1,
2729 rsurface3,
2730 ppath,
2731 lat_grid,
2732 z_field,
2733 refellipsoid,
2734 z_surface);
2735
2736 // If the field "constant" is negative, this is the first call of the
2737 // function and the path constant shall be calculated.
2738 Numeric ppc;
2739 if (ppath.constant < 0) {
2740 ppc = geometrical_ppc(r_start, za_start);
2741 } else {
2742 ppc = ppath.constant;
2743 }
2744
2745 // Vars to hold found path points, path step length and coding for end face
2746 Vector r_v, lat_v, za_v;
2747 Numeric lstep;
2748 Index endface;
2749
2751 lat_v,
2752 za_v,
2753 lstep,
2754 endface,
2755 r_start,
2756 lat_start,
2757 za_start,
2758 -1,
2759 0,
2760 ppc,
2761 lmax,
2762 lat1,
2763 lat3,
2764 r1a,
2765 r3a,
2766 r3b,
2767 r1b,
2768 rsurface1,
2769 rsurface3);
2770 // Fill *ppath*
2771 const Index np = r_v.nelem();
2772 ppath_end_2d(ppath,
2773 r_v,
2774 lat_v,
2775 za_v,
2776 Vector(np - 1, lstep),
2777 Vector(np, 1),
2778 Vector(np, 1),
2779 lat_grid,
2780 z_field,
2781 refellipsoid,
2782 ip,
2783 ilat,
2784 endface,
2785 ppc);
2786}
2787
2793void do_gridcell_3d_byltest(Vector& r_v,
2794 Vector& lat_v,
2795 Vector& lon_v,
2796 Vector& za_v,
2797 Vector& aa_v,
2798 Numeric& lstep,
2799 Index& endface,
2800 const Numeric& r_start0,
2801 const Numeric& lat_start0,
2802 const Numeric& lon_start0,
2803 const Numeric& za_start,
2804 const Numeric& aa_start,
2805 const Numeric& l_start,
2806 const Index& icall,
2807 const Numeric& ppc,
2808 const Numeric& lmax,
2809 const Numeric& lat1,
2810 const Numeric& lat3,
2811 const Numeric& lon5,
2812 const Numeric& lon6,
2813 const Numeric& r15a,
2814 const Numeric& r35a,
2815 const Numeric& r36a,
2816 const Numeric& r16a,
2817 const Numeric& r15b,
2818 const Numeric& r35b,
2819 const Numeric& r36b,
2820 const Numeric& r16b,
2821 const Numeric& rsurface15,
2822 const Numeric& rsurface35,
2823 const Numeric& rsurface36,
2824 const Numeric& rsurface16) {
2825 Numeric r_start = r_start0;
2826 Numeric lat_start = lat_start0;
2827 Numeric lon_start = lon_start0;
2828
2829 ARTS_ASSERT(icall < 10);
2830
2831 // Assert latitude and longitude
2832 ARTS_ASSERT(lat_start >= lat1 - LATLONTOL);
2833 ARTS_ASSERT(lat_start <= lat3 + LATLONTOL);
2834 ARTS_ASSERT(!(abs(lat_start) < POLELAT && lon_start < lon5 - LATLONTOL));
2835 ARTS_ASSERT(!(abs(lat_start) < POLELAT && lon_start > lon6 + LATLONTOL));
2836
2837 // Assertthat ppc and path data are consustent
2838 ARTS_ASSERT(abs(ppc - r_start0 * sin(DEG2RAD * za_start)) < 0.1);
2839
2840 // Shift latitude and longitude if outside
2841 if (lat_start < lat1) {
2842 lat_start = lat1;
2843 } else if (lat_start > lat3) {
2844 lat_start = lat3;
2845 }
2846 if (lon_start < lon5) {
2847 lon_start = lon5;
2848 } else if (lon_start > lon6) {
2849 lon_start = lon6;
2850 }
2851
2852 // Radius of lower and upper pressure level at the start position
2853 Numeric rlow = rsurf_at_latlon(
2854 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_start, lon_start);
2855 Numeric rupp = rsurf_at_latlon(
2856 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_start, lon_start);
2857
2858 // Assert radius (some extra tolerance is needed for radius)
2859 ARTS_ASSERT(r_start >= rlow - RTOL);
2860 ARTS_ASSERT(r_start <= rupp + RTOL);
2861
2862 // Shift radius if outside
2863 if (r_start < rlow) {
2864 r_start = rlow;
2865 } else if (r_start > rupp) {
2866 r_start = rupp;
2867 }
2868
2869 // Position and LOS in cartesian coordinates
2870 Numeric x, y, z, dx, dy, dz;
2872 x, y, z, dx, dy, dz, r_start, lat_start, lon_start, za_start, aa_start);
2873
2874 // Some booleans to check for recursive call
2875 bool unsafe = false;
2876 bool do_surface = false;
2877
2878 // Determine the position of the end point
2879 //
2880 endface = 0;
2881 //
2882 Numeric r_end, lat_end, lon_end, l_end;
2883
2884 // Zenith and nadir looking are handled as special cases
2885
2886 // Zenith looking
2887 if (za_start < ANGTOL) {
2888 l_end = rupp - r_start;
2889 endface = 4;
2890 }
2891
2892 // Nadir looking
2893 else if (za_start > 180 - ANGTOL) {
2894 const Numeric rsurface = rsurf_at_latlon(lat1,
2895 lat3,
2896 lon5,
2897 lon6,
2898 rsurface15,
2899 rsurface35,
2900 rsurface36,
2901 rsurface16,
2902 lat_start,
2903 lon_start);
2904
2905 if (rlow > rsurface) {
2906 l_end = r_start - rlow;
2907 endface = 2;
2908 } else {
2909 l_end = r_start - rsurface;
2910 endface = 7;
2911 }
2912 }
2913
2914 else {
2915 unsafe = true;
2916
2917 // Calculate correction terms for the position to compensate for
2918 // numerical inaccuracy.
2919 //
2920 Numeric r_corr, lat_corr, lon_corr;
2921 //
2922 cart2sph(r_corr,
2923 lat_corr,
2924 lon_corr,
2925 x,
2926 y,
2927 z,
2928 lat_start,
2929 lon_start,
2930 za_start,
2931 aa_start);
2932 //
2933 r_corr -= r_start;
2934 lat_corr -= lat_start;
2935 lon_corr -= lon_start;
2936
2937 // The end point is found by testing different step lengths until the
2938 // step length has been determined by a precision of *l_acc*.
2939 //
2940 // The first step is to found a length that goes outside the grid cell,
2941 // to find an upper limit. The lower limit is set to 0. The upper
2942 // limit is either set to l_start or it is estimated from the size of
2943 // the grid cell.
2944 // The search algorith is bisection, the next length to test is the
2945 // mean value of the minimum and maximum length limits.
2946 //
2947 if (l_start > 0) {
2948 l_end = l_start;
2949 } else {
2950 l_end = 2 * (rupp - rlow);
2951 }
2952 //
2953 Numeric l_in = 0, l_out = l_end;
2954 bool ready = false, startup = true;
2955
2956 if (rsurface15 + RTOL >= r15a || rsurface35 + RTOL >= r35a ||
2957 rsurface36 + RTOL >= r36a || rsurface16 + RTOL >= r16a) {
2958 do_surface = true;
2959 }
2960
2961 while (!ready) {
2962 cart2sph(r_end,
2963 lat_end,
2964 lon_end,
2965 x + dx * l_end,
2966 y + dy * l_end,
2967 z + dz * l_end,
2968 lat_start,
2969 lon_start,
2970 za_start,
2971 aa_start);
2972 r_end -= r_corr;
2973 lat_end -= lat_corr;
2974 lon_end -= lon_corr;
2975 resolve_lon(lon_end, lon5, lon6);
2976
2977 // Special fix for north-south observations
2978 if (abs(lat_start) < POLELAT && abs(lat_end) < POLELAT &&
2979 (abs(aa_start) < ANGTOL || abs(aa_start) > 180 - ANGTOL)) {
2980 lon_end = lon_start;
2981 }
2982
2983 bool inside = true;
2984
2985 rlow = rsurf_at_latlon(
2986 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_end, lon_end);
2987
2988 if (do_surface) {
2989 const Numeric r_surface = rsurf_at_latlon(lat1,
2990 lat3,
2991 lon5,
2992 lon6,
2993 rsurface15,
2994 rsurface35,
2995 rsurface36,
2996 rsurface16,
2997 lat_end,
2998 lon_end);
2999 if (r_surface >= rlow && r_end <= r_surface) {
3000 inside = false;
3001 endface = 7;
3002 }
3003 }
3004
3005 if (inside) {
3006 if (lat_end < lat1) {
3007 inside = false;
3008 endface = 1;
3009 } else if (lat_end > lat3) {
3010 inside = false;
3011 endface = 3;
3012 } else if (lon_end < lon5) {
3013 inside = false;
3014 endface = 5;
3015 } else if (lon_end > lon6) {
3016 inside = false;
3017 endface = 6;
3018 } else if (r_end < rlow) {
3019 inside = false;
3020 endface = 2;
3021 } else {
3022 rupp = rsurf_at_latlon(
3023 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_end, lon_end);
3024 if (r_end > rupp) {
3025 inside = false;
3026 endface = 4;
3027 }
3028 }
3029 }
3030
3031 if (startup) {
3032 if (inside) {
3033 l_in = l_end;
3034 l_end *= 5;
3035 } else {
3036 l_out = l_end;
3037 l_end = (l_out + l_in) / 2;
3038 startup = false;
3039 }
3040 } else {
3041 if (inside) {
3042 l_in = l_end;
3043 } else {
3044 l_out = l_end;
3045 }
3046
3047 if ((l_out - l_in) < LACC) {
3048 ready = true;
3049 } else {
3050 l_end = (l_out + l_in) / 2;
3051 }
3052 }
3053 }
3054 }
3055
3056 //--- Create return vectors
3057 //
3058 Index n = 1;
3059 //
3060 if (lmax > 0) {
3061 n = Index(ceil(abs(l_end / lmax)));
3062 if (n < 1) {
3063 n = 1;
3064 }
3065 }
3066 //
3067 r_v.resize(n + 1);
3068 lat_v.resize(n + 1);
3069 lon_v.resize(n + 1);
3070 za_v.resize(n + 1);
3071 aa_v.resize(n + 1);
3072 //
3073 r_v[0] = r_start;
3074 lat_v[0] = lat_start;
3075 lon_v[0] = lon_start;
3076 za_v[0] = za_start;
3077 aa_v[0] = aa_start;
3078 //
3079 lstep = l_end / (Numeric)n;
3080 Numeric l;
3081 bool ready = true;
3082 //
3083 for (Index j = 1; j <= n; j++) {
3084 l = lstep * (Numeric)j;
3085 cart2poslos(r_v[j],
3086 lat_v[j],
3087 lon_v[j],
3088 za_v[j],
3089 aa_v[j],
3090 x + dx * l,
3091 y + dy * l,
3092 z + dz * l,
3093 dx,
3094 dy,
3095 dz,
3096 ppc,
3097 x,
3098 y,
3099 z,
3100 lat_start,
3101 lon_start,
3102 za_start,
3103 aa_start);
3104
3105 // Shall lon values be shifted (value 0 is already OK)?
3106 resolve_lon(lon_v[j], lon5, lon6);
3107
3108 if (j < n) {
3109 if (unsafe) {
3110 // Check that r_v[j] is above lower pressure level and the
3111 // surface. This can fail around tangent points. For p-levels
3112 // with constant r this is easy to handle analytically, but the
3113 // problem is tricky in the general case with a non-spherical
3114 // geometry, and this crude solution is used instead. Not the
3115 // most elegant solution, but it works! Added later the same
3116 // check for upper level, after getting assert in that direction.
3117 // The z_field was crazy, but still formerly correct.
3118 rlow = rsurf_at_latlon(
3119 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_v[j], lon_v[j]);
3120 if (do_surface) {
3121 const Numeric r_surface = rsurf_at_latlon(lat1,
3122 lat3,
3123 lon5,
3124 lon6,
3125 rsurface15,
3126 rsurface35,
3127 rsurface36,
3128 rsurface16,
3129 lat_v[j],
3130 lon_v[j]);
3131 const Numeric r_test = max(r_surface, rlow);
3132 if (r_v[j] < r_test) {
3133 ready = false;
3134 break;
3135 }
3136 } else if (r_v[j] < rlow) {
3137 ready = false;
3138 break;
3139 }
3140
3141 rupp = rsurf_at_latlon(
3142 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_v[j], lon_v[j]);
3143 if (r_v[j] > rupp) {
3144 ready = false;
3145 break;
3146 }
3147 }
3148 } else // j==n
3149 {
3150 if (unsafe) {
3151 // Set end point to be consistent with found endface.
3152 //
3153 if (endface == 1) {
3154 lat_v[n] = lat1;
3155 } else if (endface == 2) {
3156 r_v[n] = rsurf_at_latlon(lat1,
3157 lat3,
3158 lon5,
3159 lon6,
3160 r15a,
3161 r35a,
3162 r36a,
3163 r16a,
3164 lat_v[n],
3165 lon_v[n]);
3166 } else if (endface == 3) {
3167 lat_v[n] = lat3;
3168 } else if (endface == 4) {
3169 r_v[n] = rsurf_at_latlon(lat1,
3170 lat3,
3171 lon5,
3172 lon6,
3173 r15b,
3174 r35b,
3175 r36b,
3176 r16b,
3177 lat_v[n],
3178 lon_v[n]);
3179 } else if (endface == 5) {
3180 lon_v[n] = lon5;
3181 } else if (endface == 6) {
3182 lon_v[n] = lon6;
3183 } else if (endface == 7) {
3184 r_v[n] = rsurf_at_latlon(lat1,
3185 lat3,
3186 lon5,
3187 lon6,
3188 rsurface15,
3189 rsurface35,
3190 rsurface36,
3191 rsurface16,
3192 lat_v[n],
3193 lon_v[n]);
3194 }
3195 }
3196 }
3197 }
3198
3199 if (!ready) { // If an "outside" point found, restart with l as start search length
3201 lat_v,
3202 lon_v,
3203 za_v,
3204 aa_v,
3205 lstep,
3206 endface,
3207 r_start,
3208 lat_start,
3209 lon_start,
3210 za_start,
3211 aa_start,
3212 l,
3213 icall + 1,
3214 ppc,
3215 lmax,
3216 lat1,
3217 lat3,
3218 lon5,
3219 lon6,
3220 r15a,
3221 r35a,
3222 r36a,
3223 r16a,
3224 r15b,
3225 r35b,
3226 r36b,
3227 r16b,
3228 rsurface15,
3229 rsurface35,
3230 rsurface36,
3231 rsurface16);
3232 }
3233}
3234
3236 ConstVectorView lat_grid,
3237 ConstVectorView lon_grid,
3238 ConstTensor3View z_field,
3239 ConstVectorView refellipsoid,
3240 ConstMatrixView z_surface,
3241 const Numeric& lmax) {
3242 // Radius, zenith angle and latitude of start point.
3243 Numeric r_start, lat_start, lon_start, za_start, aa_start;
3244
3245 // Lower grid index for the grid cell of interest.
3246 Index ip, ilat, ilon;
3247
3248 // Radius for corner points, latitude and longitude of the grid cell
3249 //
3250 Numeric lat1, lat3, lon5, lon6;
3251 Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
3252 Numeric rsurface15, rsurface35, rsurface36, rsurface16;
3253
3254 // Determine the variables defined above and make all possible asserts
3255 ppath_start_3d(r_start,
3256 lat_start,
3257 lon_start,
3258 za_start,
3259 aa_start,
3260 ip,
3261 ilat,
3262 ilon,
3263 lat1,
3264 lat3,
3265 lon5,
3266 lon6,
3267 r15a,
3268 r35a,
3269 r36a,
3270 r16a,
3271 r15b,
3272 r35b,
3273 r36b,
3274 r16b,
3275 rsurface15,
3276 rsurface35,
3277 rsurface36,
3278 rsurface16,
3279 ppath,
3280 lat_grid,
3281 lon_grid,
3282 z_field,
3283 refellipsoid,
3284 z_surface);
3285
3286 // If the field "constant" is negative, this is the first call of the
3287 // function and the path constant shall be calculated.
3288 Numeric ppc;
3289 if (ppath.constant < 0) {
3290 ppc = geometrical_ppc(r_start, za_start);
3291 } else {
3292 ppc = ppath.constant;
3293 }
3294
3295 // Vars to hold found path points, path step length and coding for end face
3296 Vector r_v, lat_v, lon_v, za_v, aa_v;
3297 Numeric lstep;
3298 Index endface;
3299
3301 lat_v,
3302 lon_v,
3303 za_v,
3304 aa_v,
3305 lstep,
3306 endface,
3307 r_start,
3308 lat_start,
3309 lon_start,
3310 za_start,
3311 aa_start,
3312 -1,
3313 0,
3314 ppc,
3315 lmax,
3316 lat1,
3317 lat3,
3318 lon5,
3319 lon6,
3320 r15a,
3321 r35a,
3322 r36a,
3323 r16a,
3324 r15b,
3325 r35b,
3326 r36b,
3327 r16b,
3328 rsurface15,
3329 rsurface35,
3330 rsurface36,
3331 rsurface16);
3332
3333 // Fill *ppath*
3334 const Index np = r_v.nelem();
3335 ppath_end_3d(ppath,
3336 r_v,
3337 lat_v,
3338 lon_v,
3339 za_v,
3340 aa_v,
3341 Vector(np - 1, lstep),
3342 Vector(np, 1),
3343 Vector(np, 1),
3344 lat_grid,
3345 lon_grid,
3346 z_field,
3347 refellipsoid,
3348 ip,
3349 ilat,
3350 ilon,
3351 endface,
3352 ppc);
3353}
3354
3355/*===========================================================================
3356 === Core functions for refraction *ppath_step* functions
3357 ===========================================================================*/
3358
3398 Array<Numeric>& r_array,
3399 Array<Numeric>& lat_array,
3400 Array<Numeric>& za_array,
3401 Array<Numeric>& l_array,
3402 Array<Numeric>& n_array,
3403 Array<Numeric>& ng_array,
3404 Index& endface,
3405 ConstVectorView p_grid,
3406 ConstVectorView refellipsoid,
3407 ConstTensor3View z_field,
3408 ConstTensor3View t_field,
3409 ConstTensor4View vmr_field,
3410 ConstVectorView f_grid,
3411 const Numeric& lmax,
3412 const Agenda& refr_index_air_agenda,
3413 const Numeric& lraytrace,
3414 const Numeric& rsurface,
3415 const Numeric& r1,
3416 const Numeric& r3,
3417 Numeric r,
3418 Numeric lat,
3419 Numeric za) {
3420 // Loop boolean
3421 bool ready = false;
3422
3423 // Store first point
3424 Numeric refr_index_air, refr_index_air_group;
3426 refr_index_air,
3427 refr_index_air_group,
3428 refr_index_air_agenda,
3429 p_grid,
3430 refellipsoid,
3431 z_field,
3432 t_field,
3433 vmr_field,
3434 f_grid,
3435 r);
3436 r_array.push_back(r);
3437 lat_array.push_back(lat);
3438 za_array.push_back(za);
3439 n_array.push_back(refr_index_air);
3440 ng_array.push_back(refr_index_air_group);
3441
3442 // Variables for output from do_gridcell_2d
3443 Vector r_v, lat_v, za_v;
3444 Numeric lstep, lcum = 0, dlat;
3445
3446 while (!ready) {
3447 // Constant for the geometrical step to make
3448 const Numeric ppc_step = geometrical_ppc(r, za);
3449
3450 // Where will a geometric path exit the grid cell?
3451 do_gridrange_1d(r_v,
3452 lat_v,
3453 za_v,
3454 lstep,
3455 endface,
3456 r,
3457 lat,
3458 za,
3459 ppc_step,
3460 -1,
3461 r1,
3462 r3,
3463 rsurface);
3464 ARTS_ASSERT(r_v.nelem() == 2);
3465
3466 // If *lstep* is <= *lraytrace*, extract the found end point (if not
3467 // a tangent point, we are ready).
3468 // Otherwise, we make a geometrical step with length *lraytrace*.
3469
3470 Numeric za_flagside = za;
3471
3472 if (lstep <= lraytrace) {
3473 r = r_v[1];
3474 dlat = lat_v[1] - lat;
3475 lat = lat_v[1];
3476 lcum += lstep;
3477 ready = true;
3478 } else {
3479 Numeric l;
3480 if (za <= 90) {
3481 l = geompath_l_at_r(ppc_step, r) + lraytrace;
3482 } else {
3483 l = geompath_l_at_r(ppc_step, r) - lraytrace;
3484 if (l < 0) {
3485 za_flagside = 180 - za_flagside;
3486 } // Tangent point passed!
3487 }
3488
3489 r = geompath_r_at_l(ppc_step, l);
3490
3491 const Numeric lat_new = geompath_lat_at_za(
3492 za, lat, geompath_za_at_r(ppc_step, za_flagside, r));
3493 dlat = lat_new - lat;
3494 lat = lat_new;
3495 lstep = lraytrace;
3496 lcum += lraytrace;
3497 }
3498
3499 // Refractive index at new point
3500 Numeric dndr;
3502 refr_index_air,
3503 refr_index_air_group,
3504 dndr,
3505 refr_index_air_agenda,
3506 p_grid,
3507 refellipsoid,
3508 z_field,
3509 t_field,
3510 vmr_field,
3511 f_grid,
3512 r);
3513
3514 // Calculate LOS zenith angle at found point.
3515 const Numeric za_rad = DEG2RAD * za;
3516 //
3517 za += -dlat; // Pure geometrical effect
3518 //
3519 za += (RAD2DEG * lstep / refr_index_air) * (-sin(za_rad) * dndr);
3520
3521 // Make sure that obtained *za* is inside valid range
3522 if (za < 0) {
3523 za = -za;
3524 } else if (za > 180) {
3525 za -= 360;
3526 }
3527
3528 // Store found point?
3529 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
3530 r_array.push_back(r);
3531 lat_array.push_back(lat);
3532 za_array.push_back(za);
3533 n_array.push_back(refr_index_air);
3534 ng_array.push_back(refr_index_air_group);
3535 l_array.push_back(lcum);
3536 lcum = 0;
3537 }
3538 }
3539}
3540
3542 Ppath& ppath,
3543 ConstVectorView p_grid,
3544 ConstTensor3View z_field,
3545 ConstTensor3View t_field,
3546 ConstTensor4View vmr_field,
3547 ConstVectorView f_grid,
3548 ConstVectorView refellipsoid,
3549 const Numeric& z_surface,
3550 const Numeric& lmax,
3551 const Agenda& refr_index_air_agenda,
3552 const String& rtrace_method,
3553 const Numeric& lraytrace) {
3554 // Starting radius, zenith angle and latitude
3555 Numeric r_start, lat_start, za_start;
3556
3557 // Index of the pressure level being the lower limit for the
3558 // grid range of interest.
3559 Index ip;
3560
3561 // Determine the variables defined above, and make asserts of input
3562 ppath_start_1d(r_start, lat_start, za_start, ip, ppath);
3563
3564 // If the field "constant" is negative, this is the first call of the
3565 // function and the path constant shall be calculated.
3566 // If the sensor is placed outside the atmosphere, the constant is
3567 // already set.
3568 Numeric ppc;
3569 if (ppath.constant < 0) {
3570 Numeric refr_index_air, refr_index_air_group;
3572 refr_index_air,
3573 refr_index_air_group,
3574 refr_index_air_agenda,
3575 p_grid,
3576 refellipsoid,
3577 z_field,
3578 t_field,
3579 vmr_field,
3580 f_grid,
3581 r_start);
3582 ppc = refraction_ppc(r_start, za_start, refr_index_air);
3583 } else {
3584 ppc = ppath.constant;
3585 }
3586
3587 // Perform the ray tracing
3588 //
3589 // Arrays to store found ray tracing points
3590 // (Vectors don't work here as we don't know how many points there will be)
3591 Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
3592 Index endface;
3593 //
3594 if (rtrace_method == "linear_basic") {
3595 /*
3596 raytrace_1d_linear_basic( ws, r_array, lat_array, za_array, l_array,
3597 n_array, ng_array, endface, refellipsoid, p_grid, z_field, t_field,
3598 vmr_field, f_grid, lmax, refr_index_air_agenda,
3599 lraytrace, ppc, refellipsoid[0] + z_surface,
3600 refellipsoid[0]+z_field[ip], refellipsoid[0]+z_field[ip+1],
3601 r_start, lat_start, za_start );
3602 */
3604 r_array,
3605 lat_array,
3606 za_array,
3607 l_array,
3608 n_array,
3609 ng_array,
3610 endface,
3611 p_grid,
3612 refellipsoid,
3613 z_field,
3614 t_field,
3615 vmr_field,
3616 f_grid,
3617 lmax,
3618 refr_index_air_agenda,
3619 lraytrace,
3620 refellipsoid[0] + z_surface,
3621 refellipsoid[0] + z_field(ip, 0, 0),
3622 refellipsoid[0] + z_field(ip + 1, 0, 0),
3623 r_start,
3624 lat_start,
3625 za_start);
3626 } else {
3627 // Make sure we fail if called with an invalid rtrace_method.
3628 ARTS_ASSERT(false);
3629 }
3630
3631 // Fill *ppath*
3632 //
3633 const Index np = r_array.nelem();
3634 Vector r_v(np), lat_v(np), za_v(np), l_v(np - 1), n_v(np), ng_v(np);
3635 for (Index i = 0; i < np; i++) {
3636 r_v[i] = r_array[i];
3637 lat_v[i] = lat_array[i];
3638 za_v[i] = za_array[i];
3639 n_v[i] = n_array[i];
3640 ng_v[i] = ng_array[i];
3641 if (i < np - 1) {
3642 l_v[i] = l_array[i];
3643 }
3644 }
3645 //
3646 ppath_end_1d(ppath,
3647 r_v,
3648 lat_v,
3649 za_v,
3650 l_v,
3651 n_v,
3652 ng_v,
3653 z_field(joker, 0, 0),
3654 refellipsoid,
3655 ip,
3656 endface,
3657 ppc);
3658}
3659
3704 Array<Numeric>& r_array,
3705 Array<Numeric>& lat_array,
3706 Array<Numeric>& za_array,
3707 Array<Numeric>& l_array,
3708 Array<Numeric>& n_array,
3709 Array<Numeric>& ng_array,
3710 Index& endface,
3711 ConstVectorView p_grid,
3712 ConstVectorView lat_grid,
3713 ConstVectorView refellipsoid,
3714 ConstTensor3View z_field,
3715 ConstTensor3View t_field,
3716 ConstTensor4View vmr_field,
3717 ConstVectorView f_grid,
3718 const Numeric& lmax,
3719 const Agenda& refr_index_air_agenda,
3720 const Numeric& lraytrace,
3721 const Numeric& lat1,
3722 const Numeric& lat3,
3723 const Numeric& rsurface1,
3724 const Numeric& rsurface3,
3725 const Numeric& r1a,
3726 const Numeric& r3a,
3727 const Numeric& r3b,
3728 const Numeric& r1b,
3729 Numeric r,
3730 Numeric lat,
3731 Numeric za) {
3732 // Loop boolean
3733 bool ready = false;
3734
3735 // Store first point
3736 Numeric refr_index_air, refr_index_air_group;
3738 refr_index_air,
3739 refr_index_air_group,
3740 refr_index_air_agenda,
3741 p_grid,
3742 lat_grid,
3743 refellipsoid,
3744 z_field,
3745 t_field,
3746 vmr_field,
3747 f_grid,
3748 r,
3749 lat);
3750 r_array.push_back(r);
3751 lat_array.push_back(lat);
3752 za_array.push_back(za);
3753 n_array.push_back(refr_index_air);
3754 ng_array.push_back(refr_index_air_group);
3755
3756 // Variables for output from do_gridcell_2d
3757 Vector r_v, lat_v, za_v;
3758 Numeric lstep, lcum = 0, dlat;
3759
3760 while (!ready) {
3761 // Constant for the geometrical step to make
3762 const Numeric ppc_step = geometrical_ppc(r, za);
3763
3764 // Where will a geometric path exit the grid cell?
3766 lat_v,
3767 za_v,
3768 lstep,
3769 endface,
3770 r,
3771 lat,
3772 za,
3773 lraytrace,
3774 0,
3775 ppc_step,
3776 -1,
3777 lat1,
3778 lat3,
3779 r1a,
3780 r3a,
3781 r3b,
3782 r1b,
3783 rsurface1,
3784 rsurface3);
3785 ARTS_ASSERT(r_v.nelem() == 2);
3786
3787 // If *lstep* is <= *lraytrace*, extract the found end point (if not
3788 // a tangent point, we are ready).
3789 // Otherwise, we make a geometrical step with length *lraytrace*.
3790
3791 Numeric za_flagside = za;
3792
3793 if (lstep <= lraytrace) {
3794 r = r_v[1];
3795 dlat = lat_v[1] - lat;
3796 lat = lat_v[1];
3797 lcum += lstep;
3798 ready = true;
3799 } else {
3800 Numeric l;
3801 if (abs(za) <= 90) {
3802 l = geompath_l_at_r(ppc_step, r) + lraytrace;
3803 } else {
3804 l = geompath_l_at_r(ppc_step, r) - lraytrace;
3805 if (l < 0) // Tangent point passed!
3806 {
3807 za_flagside = sign(za) * 180 - za_flagside;
3808 }
3809 }
3810
3811 r = geompath_r_at_l(ppc_step, l);
3812
3813 const Numeric lat_new = geompath_lat_at_za(
3814 za, lat, geompath_za_at_r(ppc_step, za_flagside, r));
3815 dlat = lat_new - lat;
3816 lat = lat_new;
3817 lstep = lraytrace;
3818 lcum += lraytrace;
3819
3820 // For paths along the latitude end faces we can end up outside the
3821 // grid cell. We simply look for points outisde the grid cell.
3822 if (lat < lat1) {
3823 lat = lat1;
3824 } else if (lat > lat3) {
3825 lat = lat3;
3826 }
3827 }
3828
3829 // Refractive index at new point
3830 Numeric dndr, dndlat;
3832 refr_index_air,
3833 refr_index_air_group,
3834 dndr,
3835 dndlat,
3836 refr_index_air_agenda,
3837 p_grid,
3838 lat_grid,
3839 refellipsoid,
3840 z_field,
3841 t_field,
3842 vmr_field,
3843 f_grid,
3844 r,
3845 lat);
3846
3847 // Calculate LOS zenith angle at found point.
3848 const Numeric za_rad = DEG2RAD * za;
3849 //
3850 za += -dlat; // Pure geometrical effect
3851 //
3852 za += (RAD2DEG * lstep / refr_index_air) *
3853 (-sin(za_rad) * dndr + cos(za_rad) * dndlat);
3854
3855 // Make sure that obtained *za* is inside valid range
3856 if (za < -180) {
3857 za += 360;
3858 } else if (za > 180) {
3859 za -= 360;
3860 }
3861
3862 // If the path is zenith/nadir along a latitude end face, we must check
3863 // that the path does not exit with new *za*.
3864 if (lat == lat1 && za < 0) {
3865 endface = 1;
3866 ready = 1;
3867 } else if (lat == lat3 && za > 0) {
3868 endface = 3;
3869 ready = 1;
3870 }
3871
3872 // Store found point?
3873 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
3874 r_array.push_back(r);
3875 lat_array.push_back(lat);
3876 za_array.push_back(za);
3877 n_array.push_back(refr_index_air);
3878 ng_array.push_back(refr_index_air_group);
3879 l_array.push_back(lcum);
3880 lcum = 0;
3881 }
3882 }
3883}
3884
3886 Ppath& ppath,
3887 ConstVectorView p_grid,
3888 ConstVectorView lat_grid,
3889 ConstTensor3View z_field,
3890 ConstTensor3View t_field,
3891 ConstTensor4View vmr_field,
3892 ConstVectorView f_grid,
3893 ConstVectorView refellipsoid,
3894 ConstVectorView z_surface,
3895 const Numeric& lmax,
3896 const Agenda& refr_index_air_agenda,
3897 const String& rtrace_method,
3898 const Numeric& lraytrace) {
3899 // Radius, zenith angle and latitude of start point.
3900 Numeric r_start, lat_start, za_start;
3901
3902 // Lower grid index for the grid cell of interest.
3903 Index ip, ilat;
3904
3905 // Radii and latitudes set by *ppath_start_2d*.
3906 Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
3907
3908 // Determine the variables defined above and make all possible asserts
3909 ppath_start_2d(r_start,
3910 lat_start,
3911 za_start,
3912 ip,
3913 ilat,
3914 lat1,
3915 lat3,
3916 r1a,
3917 r3a,
3918 r3b,
3919 r1b,
3920 rsurface1,
3921 rsurface3,
3922 ppath,
3923 lat_grid,
3924 z_field(joker, joker, 0),
3925 refellipsoid,
3926 z_surface);
3927
3928 // Perform the ray tracing
3929 //
3930 // No constant for the path is valid here.
3931 //
3932 // Arrays to store found ray tracing points
3933 // (Vectors don't work here as we don't know how many points there will be)
3934 Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
3935 Index endface;
3936 //
3937 if (rtrace_method == "linear_basic") {
3939 r_array,
3940 lat_array,
3941 za_array,
3942 l_array,
3943 n_array,
3944 ng_array,
3945 endface,
3946 p_grid,
3947 lat_grid,
3948 refellipsoid,
3949 z_field,
3950 t_field,
3951 vmr_field,
3952 f_grid,
3953 lmax,
3954 refr_index_air_agenda,
3955 lraytrace,
3956 lat1,
3957 lat3,
3958 rsurface1,
3959 rsurface3,
3960 r1a,
3961 r3a,
3962 r3b,
3963 r1b,
3964 r_start,
3965 lat_start,
3966 za_start);
3967 } else {
3968 // Make sure we fail if called with an invalid rtrace_method.
3969 ARTS_ASSERT(false);
3970 }
3971
3972 // Fill *ppath*
3973 //
3974 const Index np = r_array.nelem();
3975 Vector r_v(np), lat_v(np), za_v(np), l_v(np - 1), n_v(np), ng_v(np);
3976 for (Index i = 0; i < np; i++) {
3977 r_v[i] = r_array[i];
3978 lat_v[i] = lat_array[i];
3979 za_v[i] = za_array[i];
3980 n_v[i] = n_array[i];
3981 ng_v[i] = ng_array[i];
3982 if (i < np - 1) {
3983 l_v[i] = l_array[i];
3984 }
3985 }
3986 //
3987 ppath_end_2d(ppath,
3988 r_v,
3989 lat_v,
3990 za_v,
3991 l_v,
3992 n_v,
3993 ng_v,
3994 lat_grid,
3995 z_field(joker, joker, 0),
3996 refellipsoid,
3997 ip,
3998 ilat,
3999 endface,
4000 -1);
4001}
4002
4061 Array<Numeric>& r_array,
4062 Array<Numeric>& lat_array,
4063 Array<Numeric>& lon_array,
4064 Array<Numeric>& za_array,
4065 Array<Numeric>& aa_array,
4066 Array<Numeric>& l_array,
4067 Array<Numeric>& n_array,
4068 Array<Numeric>& ng_array,
4069 Index& endface,
4070 ConstVectorView refellipsoid,
4071 ConstVectorView p_grid,
4072 ConstVectorView lat_grid,
4073 ConstVectorView lon_grid,
4074 ConstTensor3View z_field,
4075 ConstTensor3View t_field,
4076 ConstTensor4View vmr_field,
4077 ConstVectorView f_grid,
4078 const Numeric& lmax,
4079 const Agenda& refr_index_air_agenda,
4080 const Numeric& lraytrace,
4081 const Numeric& lat1,
4082 const Numeric& lat3,
4083 const Numeric& lon5,
4084 const Numeric& lon6,
4085 const Numeric& rsurface15,
4086 const Numeric& rsurface35,
4087 const Numeric& rsurface36,
4088 const Numeric& rsurface16,
4089 const Numeric& r15a,
4090 const Numeric& r35a,
4091 const Numeric& r36a,
4092 const Numeric& r16a,
4093 const Numeric& r15b,
4094 const Numeric& r35b,
4095 const Numeric& r36b,
4096 const Numeric& r16b,
4097 Numeric r,
4098 Numeric lat,
4099 Numeric lon,
4100 Numeric za,
4101 Numeric aa) {
4102 // Loop boolean
4103 bool ready = false;
4104
4105 // Store first point
4106 Numeric refr_index_air, refr_index_air_group;
4108 refr_index_air,
4109 refr_index_air_group,
4110 refr_index_air_agenda,
4111 p_grid,
4112 lat_grid,
4113 lon_grid,
4114 refellipsoid,
4115 z_field,
4116 t_field,
4117 vmr_field,
4118 f_grid,
4119 r,
4120 lat,
4121 lon);
4122 r_array.push_back(r);
4123 lat_array.push_back(lat);
4124 lon_array.push_back(lon);
4125 za_array.push_back(za);
4126 aa_array.push_back(aa);
4127 n_array.push_back(refr_index_air);
4128 ng_array.push_back(refr_index_air_group);
4129
4130 // Variables for output from do_gridcell_3d
4131 Vector r_v, lat_v, lon_v, za_v, aa_v;
4132 Numeric lstep, lcum = 0;
4133 Numeric za_new, aa_new;
4134
4135 while (!ready) {
4136 // Constant for the geometrical step to make
4137 const Numeric ppc_step = geometrical_ppc(r, za);
4138
4139 // Where will a geometric path exit the grid cell?
4141 lat_v,
4142 lon_v,
4143 za_v,
4144 aa_v,
4145 lstep,
4146 endface,
4147 r,
4148 lat,
4149 lon,
4150 za,
4151 aa,
4152 lraytrace,
4153 0,
4154 ppc_step,
4155 -1,
4156 lat1,
4157 lat3,
4158 lon5,
4159 lon6,
4160 r15a,
4161 r35a,
4162 r36a,
4163 r16a,
4164 r15b,
4165 r35b,
4166 r36b,
4167 r16b,
4168 rsurface15,
4169 rsurface35,
4170 rsurface36,
4171 rsurface16);
4172 ARTS_ASSERT(r_v.nelem() == 2);
4173
4174 // If *lstep* is <= *lraytrace*, extract the found end point.
4175 // Otherwise, we make a geometrical step with length *lraytrace*.
4176
4177 if (lstep <= lraytrace) {
4178 r = r_v[1];
4179 lat = lat_v[1];
4180 lon = lon_v[1];
4181 za_new = za_v[1];
4182 aa_new = aa_v[1];
4183 lcum += lstep;
4184 ready = true;
4185 } else {
4186 // Sensor pos and LOS in cartesian coordinates
4187 Numeric x, y, z, dx, dy, dz, lat_new, lon_new;
4188 //
4189 poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
4190 lstep = lraytrace;
4191 cart2poslos(r,
4192 lat_new,
4193 lon_new,
4194 za_new,
4195 aa_new,
4196 x + dx * lstep,
4197 y + dy * lstep,
4198 z + dz * lstep,
4199 dx,
4200 dy,
4201 dz,
4202 ppc_step,
4203 x,
4204 y,
4205 z,
4206 lat,
4207 lon,
4208 za,
4209 aa);
4210 lcum += lstep;
4211
4212 // Shall lon values be shifted?
4213 resolve_lon(lon_new, lon5, lon6);
4214
4215 lat = lat_new;
4216 lon = lon_new;
4217 }
4218
4219 // Refractive index at new point
4220 Numeric dndr, dndlat, dndlon;
4222 refr_index_air,
4223 refr_index_air_group,
4224 dndr,
4225 dndlat,
4226 dndlon,
4227 refr_index_air_agenda,
4228 p_grid,
4229 lat_grid,
4230 lon_grid,
4231 refellipsoid,
4232 z_field,
4233 t_field,
4234 vmr_field,
4235 f_grid,
4236 r,
4237 lat,
4238 lon);
4239
4240 // Calculate LOS zenith angle at found point.
4241 const Numeric aterm = RAD2DEG * lstep / refr_index_air;
4242 const Numeric za_rad = DEG2RAD * za;
4243 const Numeric aa_rad = DEG2RAD * aa;
4244 const Numeric sinza = sin(za_rad);
4245 const Numeric sinaa = sin(aa_rad);
4246 const Numeric cosaa = cos(aa_rad);
4247 //
4248 Vector los(2);
4249 los[0] = za_new;
4250 los[1] = aa_new;
4251 //
4252 if (za < ANGTOL || za > 180 - ANGTOL) {
4253 los[0] += aterm * (cos(za_rad) * (cosaa * dndlat + sinaa * dndlon));
4254 los[1] = RAD2DEG * atan2(dndlon, dndlat);
4255 } else {
4256 los[0] += aterm * (-sinza * dndr +
4257 cos(za_rad) * (cosaa * dndlat + sinaa * dndlon));
4258 los[1] += aterm * sinza * (cosaa * dndlon - sinaa * dndlat);
4259 }
4260 //
4261 adjust_los(los, 3);
4262 //
4263 za = los[0];
4264 aa = los[1];
4265
4266 // For some cases where the path goes along an end face,
4267 // it could be the case that the refraction bends the path out
4268 // of the grid cell.
4269 if (za > 0 && za < 180) {
4270 if (lon == lon5 && aa < 0) {
4271 endface = 5;
4272 ready = 1;
4273 } else if (lon == lon6 && aa > 0) {
4274 endface = 6;
4275 ready = 1;
4276 } else if (lat == lat1 && lat != -90 && abs(aa) > 90) {
4277 endface = 1;
4278 ready = 1;
4279 } else if (lat == lat3 && lat != 90 && abs(aa) < 90) {
4280 endface = 3;
4281 ready = 1;
4282 }
4283 }
4284
4285 // Store found point?
4286 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
4287 r_array.push_back(r);
4288 lat_array.push_back(lat);
4289 lon_array.push_back(lon);
4290 za_array.push_back(za);
4291 aa_array.push_back(aa);
4292 n_array.push_back(refr_index_air);
4293 ng_array.push_back(refr_index_air_group);
4294 l_array.push_back(lcum);
4295 lcum = 0;
4296 }
4297 }
4298}
4299
4301 Ppath& ppath,
4302 ConstVectorView p_grid,
4303 ConstVectorView lat_grid,
4304 ConstVectorView lon_grid,
4305 ConstTensor3View z_field,
4306 ConstTensor3View t_field,
4307 ConstTensor4View vmr_field,
4308 ConstVectorView f_grid,
4309 ConstVectorView refellipsoid,
4310 ConstMatrixView z_surface,
4311 const Numeric& lmax,
4312 const Agenda& refr_index_air_agenda,
4313 const String& rtrace_method,
4314 const Numeric& lraytrace) {
4315 // Radius, zenith angle and latitude of start point.
4316 Numeric r_start, lat_start, lon_start, za_start, aa_start;
4317
4318 // Lower grid index for the grid cell of interest.
4319 Index ip, ilat, ilon;
4320
4321 // Radius for corner points, latitude and longitude of the grid cell
4322 //
4323 Numeric lat1, lat3, lon5, lon6;
4324 Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
4325 Numeric rsurface15, rsurface35, rsurface36, rsurface16;
4326
4327 // Determine the variables defined above and make all possible asserts
4328 ppath_start_3d(r_start,
4329 lat_start,
4330 lon_start,
4331 za_start,
4332 aa_start,
4333 ip,
4334 ilat,
4335 ilon,
4336 lat1,
4337 lat3,
4338 lon5,
4339 lon6,
4340 r15a,
4341 r35a,
4342 r36a,
4343 r16a,
4344 r15b,
4345 r35b,
4346 r36b,
4347 r16b,
4348 rsurface15,
4349 rsurface35,
4350 rsurface36,
4351 rsurface16,
4352 ppath,
4353 lat_grid,
4354 lon_grid,
4355 z_field,
4356 refellipsoid,
4357 z_surface);
4358
4359 // Perform the ray tracing
4360 //
4361 // No constant for the path is valid here.
4362 //
4363 // Arrays to store found ray tracing points
4364 // (Vectors don't work here as we don't know how many points there will be)
4365 Array<Numeric> r_array, lat_array, lon_array, za_array, aa_array;
4366 Array<Numeric> l_array, n_array, ng_array;
4367 Index endface;
4368 //
4369 if (rtrace_method == "linear_basic") {
4371 r_array,
4372 lat_array,
4373 lon_array,
4374 za_array,
4375 aa_array,
4376 l_array,
4377 n_array,
4378 ng_array,
4379 endface,
4380 refellipsoid,
4381 p_grid,
4382 lat_grid,
4383 lon_grid,
4384 z_field,
4385 t_field,
4386 vmr_field,
4387 f_grid,
4388 lmax,
4389 refr_index_air_agenda,
4390 lraytrace,
4391 lat1,
4392 lat3,
4393 lon5,
4394 lon6,
4395 rsurface15,
4396 rsurface35,
4397 rsurface36,
4398 rsurface16,
4399 r15a,
4400 r35a,
4401 r36a,
4402 r16a,
4403 r15b,
4404 r35b,
4405 r36b,
4406 r16b,
4407 r_start,
4408 lat_start,
4409 lon_start,
4410 za_start,
4411 aa_start);
4412 } else {
4413 // Make sure we fail if called with an invalid rtrace_method.
4414 ARTS_ASSERT(false);
4415 }
4416
4417 // Fill *ppath*
4418 //
4419 const Index np = r_array.nelem();
4420 Vector r_v(np), lat_v(np), lon_v(np), za_v(np), aa_v(np), l_v(np - 1);
4421 Vector n_v(np), ng_v(np);
4422 for (Index i = 0; i < np; i++) {
4423 r_v[i] = r_array[i];
4424 lat_v[i] = lat_array[i];
4425 lon_v[i] = lon_array[i];
4426 za_v[i] = za_array[i];
4427 aa_v[i] = aa_array[i];
4428 n_v[i] = n_array[i];
4429 ng_v[i] = ng_array[i];
4430 if (i < np - 1) {
4431 l_v[i] = l_array[i];
4432 }
4433 }
4434 //
4435 // Fill *ppath*
4436 ppath_end_3d(ppath,
4437 r_v,
4438 lat_v,
4439 lon_v,
4440 za_v,
4441 aa_v,
4442 l_v,
4443 n_v,
4444 ng_v,
4445 lat_grid,
4446 lon_grid,
4447 z_field,
4448 refellipsoid,
4449 ip,
4450 ilat,
4451 ilon,
4452 endface,
4453 -1);
4454}
4455
4456/*===========================================================================
4457 === Main functions
4458 ===========================================================================*/
4459
4461 const Index& atmosphere_dim,
4462 ConstVectorView p_grid,
4463 ConstVectorView lat_grid,
4464 ConstVectorView lon_grid,
4465 ConstTensor3View z_field,
4466 ConstVectorView refellipsoid,
4467 ConstMatrixView z_surface,
4468 const Index& cloudbox_on,
4469 const ArrayOfIndex& cloudbox_limits,
4470 const bool& ppath_inside_cloudbox_do,
4471 ConstVectorView rte_pos,
4472 ConstVectorView rte_los,
4473 const Verbosity& verbosity) {
4475
4476 // This function contains no checks or asserts as it is only a sub-function.
4477
4478 // Allocate the ppath structure
4479 ppath_init_structure(ppath, atmosphere_dim, 1);
4480
4481 // Index of last pressure level
4482 const Index lp = p_grid.nelem() - 1;
4483
4484 // The different atmospheric dimensionalities are handled seperately
4485
4486 //-- 1D ---------------------------------------------------------------------
4487 if (atmosphere_dim == 1) {
4488 // End position and LOS
4489 ppath.end_pos[0] = rte_pos[0];
4490 ppath.end_pos[1] = 0;
4491 ppath.end_los[0] = rte_los[0];
4492
4493 // Sensor is inside the model atmosphere:
4494 if (rte_pos[0] < z_field(lp, 0, 0)) {
4495 // Check that the sensor is above the surface
4496 ARTS_USER_ERROR_IF ((rte_pos[0] + RTOL) < z_surface(0, 0),
4497 "The ppath starting point is placed ",
4498 (z_surface(0, 0) - rte_pos[0]) / 1e3, " km below the surface.")
4499
4500 // Set ppath
4501 ppath.pos(0, joker) = ppath.end_pos;
4502 ppath.r[0] = refellipsoid[0] + rte_pos[0];
4503 ppath.los(0, joker) = ppath.end_los;
4504 //
4505 gridpos(ppath.gp_p, z_field(joker, 0, 0), ExhaustiveConstVectorView{ppath.pos(0, 0)});
4506 gridpos_check_fd(ppath.gp_p[0]);
4507
4508 // Is the sensor on the surface looking down?
4509 // If yes and the sensor is inside the cloudbox, the background will
4510 // be changed below.
4511 if (ppath.pos(0, 0) <= z_surface(0, 0) && ppath.los(0, 0) > 90) {
4512 ppath_set_background(ppath, 2);
4513 }
4514
4515 // Outside cloudbox:
4516 // Check sensor position with respect to cloud box.
4517 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4518 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4519 if (fgp >= (Numeric)cloudbox_limits[0] &&
4520 fgp <= (Numeric)cloudbox_limits[1]) {
4521 ppath_set_background(ppath, 4);
4522 }
4523 }
4524
4525 // Inside cloudbox:
4526 DEBUG_ONLY(if (ppath_inside_cloudbox_do) {
4527 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4528 ARTS_ASSERT(fgp >= (Numeric)cloudbox_limits[0] &&
4529 fgp <= (Numeric)cloudbox_limits[1]);
4530 })
4531 }
4532
4533 // Sensor is outside the model atmosphere:
4534 else {
4535 // We can here set ppc and n as we are outside the atmosphere
4536 ppath.nreal = 1.0;
4537 ppath.ngroup = 1.0;
4538 ppath.constant =
4539 geometrical_ppc(refellipsoid[0] + rte_pos[0], rte_los[0]);
4540
4541 // Totally outside
4542 if (rte_los[0] <= 90 ||
4543 ppath.constant >= refellipsoid[0] + z_field(lp, 0, 0)) {
4544 ppath.pos(0, 0) = rte_pos[0];
4545 ppath.pos(0, 1) = 0;
4546 ppath.r[0] = refellipsoid[0] + rte_pos[0];
4547 ppath.los(0, 0) = rte_los[0];
4548 //
4549 ppath_set_background(ppath, 1);
4550 out1 << " --- WARNING ---, path is totally outside of the "
4551 << "model atmosphere\n";
4552 }
4553
4554 // Path enters the atmosphere.
4555 else {
4556 ppath.r[0] = refellipsoid[0] + z_field(lp, 0, 0);
4557 ppath.pos(0, 0) = z_field(lp, 0, 0);
4558 ppath.los(0, 0) =
4559 geompath_za_at_r(ppath.constant, rte_los[0], ppath.r[0]);
4560 ppath.pos(0, 1) = geompath_lat_at_za(rte_los[0], 0, ppath.los(0, 0));
4561 ppath.end_lstep =
4562 geompath_l_at_r(ppath.constant, refellipsoid[0] + rte_pos[0]) -
4563 geompath_l_at_r(ppath.constant, ppath.r[0]);
4564
4565 // Here we know the grid position exactly
4566 ppath.gp_p[0].idx = lp - 1;
4567 ppath.gp_p[0].fd[0] = 1;
4568 ppath.gp_p[0].fd[1] = 0;
4569
4570 // If cloud box reaching TOA, we have also found the background
4571 if (cloudbox_on && cloudbox_limits[1] == lp) {
4572 ppath_set_background(ppath, 3);
4573 }
4574 }
4575 }
4576 } // End 1D
4577
4578 //-- 2D ---------------------------------------------------------------------
4579 else if (atmosphere_dim == 2) {
4580 // End position and LOS
4581 ppath.end_pos[0] = rte_pos[0];
4582 ppath.end_pos[1] = rte_pos[1];
4583 ppath.end_los[0] = rte_los[0];
4584
4585 // Index of last latitude
4586 const Index llat = lat_grid.nelem() - 1;
4587
4588 // Is sensor inside range of lat_grid?
4589 // If yes, determine TOA altitude at sensor position
4590 GridPos gp_lat;
4591 Vector itw(2);
4592 bool islatin = false;
4593 Numeric r_e; // Ellipsoid radius at sensor position
4594 Numeric z_toa = -99e99;
4595 if (rte_pos[1] > lat_grid[0] && rte_pos[1] < lat_grid[llat]) {
4596 islatin = true;
4597 gridpos(gp_lat, lat_grid, rte_pos[1]);
4598 interpweights(itw, gp_lat);
4599 z_toa = interp(itw, z_field(lp, joker, 0), gp_lat);
4600 r_e = refell2d(refellipsoid, lat_grid, gp_lat);
4601 } else {
4602 r_e = refell2r(refellipsoid, rte_pos[1]);
4603 }
4604
4605 // Sensor is inside the model atmosphere:
4606 if (islatin && rte_pos[0] < z_toa) {
4607 const Numeric z_s = interp(itw, z_surface(joker, 0), gp_lat);
4608
4609 // Check that the sensor is above the surface
4610 ARTS_USER_ERROR_IF ((rte_pos[0] + RTOL) < z_s,
4611 "The ppath starting point is placed ", (z_s - rte_pos[0]) / 1e3,
4612 " km below the surface.")
4613
4614 // Set ppath
4615 ppath.pos(0, joker) = ppath.end_pos;
4616 ppath.r[0] = r_e + rte_pos[0];
4617 ppath.los(0, joker) = ppath.end_los;
4618
4619 // Determine gp_p (gp_lat is recalculated, but ...)
4620 GridPos gp_lon_dummy;
4621 rte_pos2gridpos(ppath.gp_p[0],
4622 ppath.gp_lat[0],
4623 gp_lon_dummy,
4624 atmosphere_dim,
4625 p_grid,
4626 lat_grid,
4627 lon_grid,
4628 z_field,
4629 rte_pos);
4630 gridpos_check_fd(ppath.gp_p[0]);
4631 gridpos_check_fd(ppath.gp_lat[0]);
4632
4633 // Is the sensor on the surface looking down?
4634 // If yes and the sensor is inside the cloudbox, the background will
4635 // be changed below.
4636 if (ppath.pos(0, 0) <= z_s) {
4637 // Calculate radial slope of the surface
4638 Numeric rslope;
4639 plevel_slope_2d(rslope,
4640 lat_grid,
4641 refellipsoid,
4642 z_surface(joker, 0),
4643 gp_lat,
4644 ppath.los(0, 0));
4645
4646 // Calculate angular tilt of the surface
4647 const Numeric atilt = plevel_angletilt(r_e + z_s, rslope);
4648
4649 // Are we looking down into the surface?
4650 // If yes and the sensor is inside the cloudbox, the background
4651 // will be changed below.
4652 if (is_los_downwards(ppath.los(0, 0), atilt)) {
4653 ppath_set_background(ppath, 2);
4654 }
4655 }
4656
4657 // Outside cloudbox:
4658 // Check sensor position with respect to cloud box.
4659 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4660 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4661 const Numeric fgl = fractional_gp(ppath.gp_lat[0]);
4662 if (fgp >= (Numeric)cloudbox_limits[0] &&
4663 fgp <= (Numeric)cloudbox_limits[1] &&
4664 fgl >= (Numeric)cloudbox_limits[2] &&
4665 fgl <= (Numeric)cloudbox_limits[3]) {
4666 ppath_set_background(ppath, 4);
4667 }
4668 }
4669
4670 // Inside cloudbox:
4671 DEBUG_ONLY(if (ppath_inside_cloudbox_do) {
4672 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4673 const Numeric fgl = fractional_gp(ppath.gp_lat[0]);
4674 ARTS_ASSERT(fgp >= (Numeric)cloudbox_limits[0] &&
4675 fgp <= (Numeric)cloudbox_limits[1] &&
4676 fgl >= (Numeric)cloudbox_limits[2] &&
4677 fgl <= (Numeric)cloudbox_limits[3]);
4678 })
4679 }
4680
4681 // Sensor is outside the model atmosphere:
4682 else {
4683 // Handle cases when the sensor looks in the wrong way
4684 ARTS_USER_ERROR_IF ((rte_pos[1] <= lat_grid[0] && rte_los[0] <= 0) ||
4685 (rte_pos[1] >= lat_grid[llat] && rte_los[0] >= 0),
4686 "The sensor is outside (or at the limit) of the model "
4687 "atmosphere but\nlooks in the wrong direction (wrong sign "
4688 "for the zenith angle?).\nThis case includes nadir "
4689 "looking exactly at the latitude end points.")
4690
4691 // We can here set ppc and n as we are outside the atmosphere
4692 ppath.nreal = 1.0;
4693 ppath.ngroup = 1.0;
4694 const Numeric r_p = r_e + rte_pos[0];
4695 ppath.constant = geometrical_ppc(r_p, rte_los[0]);
4696
4697 // Determine TOA radii, and min and max value
4698 Vector r_toa(llat + 1);
4699 Numeric r_toa_min = 99e99, r_toa_max = -1;
4700 for (Index ilat = 0; ilat <= llat; ilat++) {
4701 r_toa[ilat] =
4702 refell2r(refellipsoid, lat_grid[ilat]) + z_field(lp, ilat, 0);
4703 if (r_toa[ilat] < r_toa_min) {
4704 r_toa_min = r_toa[ilat];
4705 }
4706 if (r_toa[ilat] > r_toa_max) {
4707 r_toa_max = r_toa[ilat];
4708 }
4709 }
4710 ARTS_USER_ERROR_IF (!islatin && r_p <= r_toa_max,
4711 "The sensor is horizontally outside (or at the limit) of "
4712 "the model\natmosphere, but is at a radius smaller than "
4713 "the maximum value of\nthe top-of-the-atmosphere radii. "
4714 "This is not allowed. Make the\nmodel atmosphere larger "
4715 "to also cover the sensor position?")
4716
4717 // Upward:
4718 if (abs(rte_los[0]) <= 90) {
4719 ppath.pos(0, 0) = rte_pos[0];
4720 ppath.pos(0, 1) = rte_pos[1];
4721 ppath.r[0] = r_e + rte_pos[0];
4722 ppath.los(0, 0) = rte_los[0];
4723 //
4724 ppath_set_background(ppath, 1);
4725 out1 << " ------- WARNING -------: path is totally outside of "
4726 << "the model atmosphere\n";
4727 }
4728
4729 // Downward:
4730 else {
4731 bool above = false, ready = false, failed = false;
4732 Numeric rt = -1, latt, lt, lt_old = L_NOT_FOUND;
4733 GridPos gp_latt;
4734 Vector itwt(2);
4735
4736 // Check if clearly above the model atmosphere
4737 if (ppath.constant >= r_toa_max) {
4738 above = true;
4739 ready = true;
4740 } else { // Otherwise pick out a suitable first test radius
4741 if (islatin || ppath.constant > r_toa_min) {
4742 rt = r_toa_max;
4743 } else {
4744 rt = r_toa_min;
4745 }
4746 }
4747
4748 // Iterate until solution found or moving out from model atm.
4749 //
4750 Index ntries = 0;
4751 //
4752 while (!ready && !failed) {
4753 // If rt < ppath.constant, ppath above atmosphere
4754 if (rt < ppath.constant) {
4755 above = true;
4756 ready = true;
4757 }
4758
4759 else {
4760 // Calculate length to entrance point at rt
4762 latt, lt, rt, r_p, rte_pos[1], rte_los[0], ppath.constant);
4763 ARTS_ASSERT(lt < 9e9);
4764
4765 // Entrance outside range of lat_grid = fail
4766 if (latt < lat_grid[0] || latt > lat_grid[llat]) {
4767 failed = true;
4768 }
4769
4770 // OK iteration
4771 else {
4772 // Converged?
4773 if (abs(lt - lt_old) < LACC || ntries == 10000) {
4774 ready = true;
4775 }
4776
4777 // Update rt
4778 lt_old = lt;
4779 gridpos(gp_latt, lat_grid, latt);
4780 interpweights(itwt, gp_latt);
4781 rt = interp(itwt, r_toa, gp_latt);
4782 }
4783 }
4784 ntries++;
4785 // For a case with dense lat_grid the calculations here ended up in
4786 // an endless loop, jumping between two lt. Tried different things,
4787 // but nothing worked besides this hard stop of the iteration.
4788 if (ntries == 10000) {
4789 ready = true;
4790 }
4791 } // while
4792
4793 ARTS_USER_ERROR_IF (failed,
4794 "The path does not enter the model atmosphere. It "
4795 "reaches the\ntop of the atmosphere "
4796 "altitude around latitude ", latt, " deg.")
4797 if (above) {
4798 ppath.pos(0, 0) = rte_pos[0];
4799 ppath.pos(0, 1) = rte_pos[1];
4800 ppath.r[0] = r_e + rte_pos[0];
4801 ppath.los(0, 0) = rte_los[0];
4802 //
4803 ppath_set_background(ppath, 1);
4804 out1 << " ------- WARNING -------: path is totally outside "
4805 << "of the model atmosphere\n";
4806 } else {
4807 ppath.r[0] = rt;
4808 ppath.pos(0, 0) = interp(itwt, z_field(lp, joker, 0), gp_latt);
4809 // Calculate za first and use to determine lat
4810 ppath.los(0, 0) = geompath_za_at_r(ppath.constant, rte_los[0], rt);
4811 ppath.pos(0, 1) = interp(itwt, lat_grid, gp_latt);
4812 // This was used before, but replaced when inroduced ntries
4813 //geompath_lat_at_za(rte_los[0], rte_pos[1], ppath.los(0, 0));
4814 ppath.end_lstep = lt;
4815
4816 // Here we know the pressure grid position exactly
4817 ppath.gp_p[0].idx = lp - 1;
4818 ppath.gp_p[0].fd[0] = 1;
4819 ppath.gp_p[0].fd[1] = 0;
4820
4821 // Latitude grid position already calculated
4822 gridpos_copy(ppath.gp_lat[0], gp_latt);
4823
4824 // Hit with cloudbox reaching TOA?
4825 if (cloudbox_on && cloudbox_limits[1] == lp) {
4826 Numeric fgp = fractional_gp(gp_latt);
4827 if (fgp >= (Numeric)cloudbox_limits[2] &&
4828 fgp <= (Numeric)cloudbox_limits[3]) {
4829 ppath_set_background(ppath, 3);
4830 }
4831 }
4832 }
4833 } // Downward
4834 } // Outside
4835 } // End 2D
4836
4837 //-- 3D ---------------------------------------------------------------------
4838 else {
4839 // Index of last latitude and longitude
4840 const Index llat = lat_grid.nelem() - 1;
4841 const Index llon = lon_grid.nelem() - 1;
4842
4843 // Adjust longitude of rte_pos to range used in lon_grid
4844 Numeric lon2use = rte_pos[2];
4845 resolve_lon(lon2use, lon_grid[0], lon_grid[llon]);
4846
4847 // End position and LOS
4848 ppath.end_pos[0] = rte_pos[0];
4849 ppath.end_pos[1] = rte_pos[1];
4850 ppath.end_pos[2] = lon2use;
4851 ppath.end_los[0] = rte_los[0];
4852 ppath.end_los[1] = rte_los[1];
4853
4854 // Is sensor inside range of lat_grid and lon_grid?
4855 // If yes, determine TOA altitude at sensor position
4856 GridPos gp_lat, gp_lon;
4857 Vector itw(4);
4858 bool islatlonin = false;
4859 Numeric r_e; // Ellipsoid radius at sensor position
4860 Numeric z_toa = -99e99;
4861
4862 if (rte_pos[1] >= lat_grid[0] && rte_pos[1] <= lat_grid[llat] &&
4863 lon2use >= lon_grid[0] && lon2use <= lon_grid[llon]) {
4864 islatlonin = true;
4865 gridpos(gp_lat, lat_grid, rte_pos[1]);
4866 gridpos(gp_lon, lon_grid, lon2use);
4867 interpweights(itw, gp_lat, gp_lon);
4868 z_toa = interp(itw, z_field(lp, joker, joker), gp_lat, gp_lon);
4869 r_e = refell2d(refellipsoid, lat_grid, gp_lat);
4870 } else {
4871 r_e = refell2r(refellipsoid, rte_pos[1]);
4872 }
4873
4874 // Sensor is inside the model atmosphere:
4875 if (islatlonin && rte_pos[0] < z_toa) {
4876 const Numeric z_s = interp(itw, z_surface, gp_lat, gp_lon);
4877
4878 // Check that the sensor is above the surface
4879 ARTS_USER_ERROR_IF ((rte_pos[0] + RTOL) < z_s,
4880 "The ppath starting point is placed ", (z_s - rte_pos[0]) / 1e3,
4881 " km below the surface.")
4882
4883 // Set ppath
4884 ppath.pos(0, joker) = ppath.end_pos;
4885 ppath.r[0] = r_e + rte_pos[0];
4886 ppath.los(0, joker) = ppath.end_los;
4887
4888 // Determine gp_p (gp_lat and gp_lon are recalculated, but ...)
4889 rte_pos2gridpos(ppath.gp_p[0],
4890 ppath.gp_lat[0],
4891 ppath.gp_lon[0],
4892 atmosphere_dim,
4893 p_grid,
4894 lat_grid,
4895 lon_grid,
4896 z_field,
4897 rte_pos);
4898 gridpos_check_fd(ppath.gp_p[0]);
4899 gridpos_check_fd(ppath.gp_lat[0]);
4900 gridpos_check_fd(ppath.gp_lon[0]);
4901
4902 // Is the sensor on the surface looking down?
4903 // If yes and the sensor is inside the cloudbox, the background will
4904 // be changed below.
4905 if (ppath.pos(0, 0) <= z_s) {
4906 // Calculate radial slope of the surface
4907 Numeric c1, c2;
4908 plevel_slope_3d(c1,
4909 c2,
4910 lat_grid,
4911 lon_grid,
4912 refellipsoid,
4913 z_surface,
4914 gp_lat,
4915 gp_lon,
4916 ppath.los(0, 1));
4917
4918 // Calculate angular tilt of the surface
4919 const Numeric atilt = plevel_angletilt(r_e + z_s, c1);
4920
4921 // Are we looking down into the surface?
4922 // If yes and the sensor is inside the cloudbox, the background
4923 // will be changed below.
4924 if (is_los_downwards(ppath.los(0, 0), atilt)) {
4925 ppath_set_background(ppath, 2);
4926 }
4927 }
4928
4929 // Outside cloudbox:
4930 // Check sensor position with respect to cloud box.
4931 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4932 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4933 const Numeric fgl = fractional_gp(ppath.gp_lat[0]);
4934 const Numeric fgo = fractional_gp(ppath.gp_lon[0]);
4935 if (fgp >= (Numeric)cloudbox_limits[0] &&
4936 fgp <= (Numeric)cloudbox_limits[1] &&
4937 fgl >= (Numeric)cloudbox_limits[2] &&
4938 fgl <= (Numeric)cloudbox_limits[3] &&
4939 fgo >= (Numeric)cloudbox_limits[4] &&
4940 fgo <= (Numeric)cloudbox_limits[5]) {
4941 ppath_set_background(ppath, 4);
4942 }
4943 }
4944
4945 // Inside cloudbox:
4946 DEBUG_ONLY(if (ppath_inside_cloudbox_do) {
4947 const Numeric fgp = fractional_gp(ppath.gp_p[0]);
4948 const Numeric fgl = fractional_gp(ppath.gp_lat[0]);
4949 const Numeric fgo = fractional_gp(ppath.gp_lon[0]);
4950 ARTS_ASSERT(fgp >= (Numeric)cloudbox_limits[0] &&
4951 fgp <= (Numeric)cloudbox_limits[1] &&
4952 fgl >= (Numeric)cloudbox_limits[2] &&
4953 fgl <= (Numeric)cloudbox_limits[3] &&
4954 fgo >= (Numeric)cloudbox_limits[4] &&
4955 fgo <= (Numeric)cloudbox_limits[5]);
4956 })
4957 }
4958
4959 // Sensor is outside the model atmosphere:
4960 else {
4961 // We can here set ppc and n as we are outside the atmosphere
4962 ppath.nreal = 1.0;
4963 ppath.ngroup = 1.0;
4964 const Numeric r_p = r_e + rte_pos[0];
4965 ppath.constant = geometrical_ppc(r_p, rte_los[0]);
4966
4967 // Determine TOA radii, and min and max value
4968 Matrix r_toa(llat + 1, llon + 1);
4969 Numeric r_toa_min = 99e99, r_toa_max = -1;
4970 for (Index ilat = 0; ilat <= llat; ilat++) {
4971 const Numeric r_lat = refell2r(refellipsoid, lat_grid[ilat]);
4972 for (Index ilon = 0; ilon <= llon; ilon++) {
4973 r_toa(ilat, ilon) = r_lat + z_field(lp, ilat, ilon);
4974 if (r_toa(ilat, ilon) < r_toa_min) {
4975 r_toa_min = r_toa(ilat, ilon);
4976 }
4977 if (r_toa(ilat, ilon) > r_toa_max) {
4978 r_toa_max = r_toa(ilat, ilon);
4979 }
4980 }
4981 }
4982
4983 ARTS_USER_ERROR_IF (r_p < r_toa_max - RTOL,
4984 "r_p=",setprecision(18),r_p,"\n"
4985 "r_toa_max=",r_toa_max,"\n"
4986 "The sensor is horizontally outside of "
4987 "the model\natmosphere, but is at a radius smaller than "
4988 "the maximum value of\nthe top-of-the-atmosphere radii. "
4989 "This is not allowed. Make the\nmodel atmosphere larger "
4990 "to also cover the sensor position?")
4991
4992 // Upward:
4993 if (rte_los[0] <= 90) {
4994 ppath.pos(0, 0) = rte_pos[0];
4995 ppath.pos(0, 1) = rte_pos[1];
4996 ppath.pos(0, 2) = rte_pos[2];
4997 ppath.r[0] = r_e + rte_pos[0];
4998 ppath.los(0, 0) = rte_los[0];
4999 ppath.los(0, 1) = rte_los[1];
5000 //
5001 ppath_set_background(ppath, 1);
5002 out1 << " ------- WARNING -------: path is totally outside of "
5003 << "the model atmosphere\n";
5004 }
5005
5006 // Downward:
5007 else {
5008 bool above = false, ready = false, failed = false;
5009 Numeric rt = -1, latt, lont, lt, lt_old = L_NOT_FOUND;
5010 GridPos gp_latt, gp_lont;
5011 Vector itwt(4);
5012
5013 // Check if clearly above the model atmosphere
5014 if (ppath.constant >= r_toa_max) {
5015 above = true;
5016 ready = true;
5017 } else { // Otherwise pick out a suitable first test radius
5018 if (islatlonin || ppath.constant > r_toa_min) {
5019 rt = r_toa_max;
5020 } else {
5021 rt = r_toa_min;
5022 }
5023 }
5024
5025 // Sensor pos and LOS in cartesian coordinates
5026 Numeric x, y, z, dx, dy, dz;
5027 poslos2cart(x,
5028 y,
5029 z,
5030 dx,
5031 dy,
5032 dz,
5033 r_p,
5034 rte_pos[1],
5035 lon2use,
5036 rte_los[0],
5037 rte_los[1]);
5038
5039 // Iterate until solution found or moving out from model atm.
5040 //
5041 while (!ready && !failed) {
5042 // If rt < ppath.constant, ppath above atmosphere
5043 if (rt < ppath.constant) {
5044 above = true;
5045 ready = true;
5046 }
5047
5048 else {
5049 // Calculate lat and lon for entrance point at rt
5050 r_crossing_3d(latt,
5051 lont,
5052 lt,
5053 rt,
5054 r_p,
5055 rte_pos[1],
5056 lon2use,
5057 rte_los[0],
5058 ppath.constant,
5059 x,
5060 y,
5061 z,
5062 dx,
5063 dy,
5064 dz);
5065 resolve_lon(lont, lon_grid[0], lon_grid[llon]);
5066
5067 // Entrance outside range of lat/lon_grids = fail
5068 if (latt < lat_grid[0] || latt > lat_grid[llat] ||
5069 lont < lon_grid[0] || lont > lon_grid[llon]) {
5070 failed = true;
5071 }
5072
5073 // OK iteration
5074 else {
5075 // Converged?
5076 if (abs(lt - lt_old) < LACC) {
5077 ready = true;
5078 }
5079
5080 // Update rt
5081 lt_old = lt;
5082 gridpos(gp_latt, lat_grid, latt);
5083 gridpos(gp_lont, lon_grid, lont);
5084 interpweights(itwt, gp_latt, gp_lont);
5085 rt = interp(itwt, r_toa, gp_latt, gp_lont);
5086 }
5087 }
5088 } // while
5089
5090 ARTS_USER_ERROR_IF (failed,
5091 "The path does not enter the model atmosphere. It\n"
5092 "reaches the top of the atmosphere altitude around:\n"
5093 " lat: ", latt, " deg.\n lon: ", lont, " deg.")
5094 if (above) {
5095 ppath.pos(0, 0) = rte_pos[0];
5096 ppath.pos(0, 1) = rte_pos[1];
5097 ppath.pos(0, 2) = rte_pos[2];
5098// ppath.pos(0, 1) = lon2use;
5099 ppath.r[0] = r_e + rte_pos[0];
5100 ppath.los(0, 0) = rte_los[0];
5101 ppath.los(0, 1) = rte_los[1];
5102 //
5103 ppath_set_background(ppath, 1);
5104 out1 << " ------- WARNING -------: path is totally outside "
5105 << "of the model atmosphere\n";
5106 } else {
5107 // Calculate lt for last rt, and use it to determine pos/los
5108 lt = geompath_l_at_r(ppath.constant, r_p) -
5109 geompath_l_at_r(ppath.constant, rt);
5110 cart2poslos(ppath.r[0],
5111 ppath.pos(0, 1),
5112 ppath.pos(0, 2),
5113 ppath.los(0, 0),
5114 ppath.los(0, 1),
5115 x + dx * lt,
5116 y + dy * lt,
5117 z + dz * lt,
5118 dx,
5119 dy,
5120 dz,
5121 ppath.constant,
5122 x,
5123 y,
5124 z,
5125 rte_pos[1],
5126 lon2use,
5127 rte_los[0],
5128 rte_los[1]);
5129 ARTS_ASSERT(abs(ppath.r[0] - rt) < RTOL);
5130 resolve_lon(ppath.pos(0, 2), lon_grid[0], lon_grid[llon]);
5131 //
5132 ppath.pos(0, 0) =
5133 interp(itwt, z_field(lp, joker, joker), gp_latt, gp_lont);
5134 ppath.end_lstep = lt;
5135
5136 // Here we know the pressure grid position exactly
5137 ppath.gp_p[0].idx = lp - 1;
5138 ppath.gp_p[0].fd[0] = 1;
5139 ppath.gp_p[0].fd[1] = 0;
5140
5141 // Lat and lon grid position already calculated
5142 gridpos_copy(ppath.gp_lat[0], gp_latt);
5143 gridpos_copy(ppath.gp_lon[0], gp_lont);
5144
5145 // Hit with cloudbox reaching TOA?
5146 if (cloudbox_on && cloudbox_limits[1] == lp) {
5147 Numeric fgp1 = fractional_gp(gp_latt);
5148 Numeric fgp2 = fractional_gp(gp_lont);
5149 if (fgp1 >= (Numeric)cloudbox_limits[2] &&
5150 fgp1 <= (Numeric)cloudbox_limits[3] &&
5151 fgp2 >= (Numeric)cloudbox_limits[4] &&
5152 fgp2 <= (Numeric)cloudbox_limits[5]) {
5153 ppath_set_background(ppath, 3);
5154 }
5155 }
5156 }
5157 } // Downward
5158 } // Outside
5159 } // End 3D
5160}
5161
5163 Ppath& ppath,
5164 const Agenda& ppath_step_agenda,
5165 const Index& atmosphere_dim,
5166 const Vector& p_grid,
5167 const Vector& lat_grid,
5168 const Vector& lon_grid,
5169 const Tensor3& z_field,
5170 const Vector& f_grid,
5171 const Vector& refellipsoid,
5172 const Matrix& z_surface,
5173 const Index& cloudbox_on,
5174 const ArrayOfIndex& cloudbox_limits,
5175 const Vector& rte_pos,
5176 const Vector& rte_los,
5177 const Numeric& ppath_lmax,
5178 const Numeric& ppath_lraytrace,
5179 const bool& ppath_inside_cloudbox_do,
5180 const Verbosity& verbosity) {
5181 // This function is a WSM but it is normally only called from yCalc.
5182 // For that reason, this function does not repeat input checks that are
5183 // performed in yCalc, it only performs checks regarding the sensor
5184 // position and LOS.
5185
5186 //--- Check input -----------------------------------------------------------
5187 chk_rte_pos(atmosphere_dim, rte_pos);
5188 chk_rte_los(atmosphere_dim, rte_los);
5189 ARTS_USER_ERROR_IF (ppath_inside_cloudbox_do && !cloudbox_on,
5190 "The WSV *ppath_inside_cloudbox_do* can only be set "
5191 "to 1 if also *cloudbox_on* is 1.");
5192 //--- End: Check input ------------------------------------------------------
5193
5194 // Initiate the partial Ppath structure.
5195 // The function doing the work sets ppath_step to the point of the path
5196 // inside the atmosphere closest to the sensor, if the path is at all inside
5197 // the atmosphere.
5198 // If the background field is set by the function this flags that there is no
5199 // path to follow (for example when the sensor is inside the cloud box).
5200 // The function checks also that the sensor position and LOS give an
5201 // allowed path.
5202 //
5203 Ppath ppath_step;
5204 //
5205 ppath_start_stepping(ppath_step,
5206 atmosphere_dim,
5207 p_grid,
5208 lat_grid,
5209 lon_grid,
5210 z_field,
5211 refellipsoid,
5212 z_surface,
5213 cloudbox_on,
5214 cloudbox_limits,
5215 ppath_inside_cloudbox_do,
5216 rte_pos,
5217 rte_los,
5218 verbosity);
5219 // For debugging:
5220 //Print( ppath_step, 0, verbosity );
5221
5222 // The only data we need to store from this initial ppath_step is
5223 // end_pos/los/lstep
5224 const Numeric end_lstep = ppath_step.end_lstep;
5225 const Vector end_pos = ppath_step.end_pos;
5226 const Vector end_los = ppath_step.end_los;
5227
5228 // Perform propagation path steps until the starting point is found, which
5229 // is flagged by ppath_step by setting the background field.
5230 //
5231 // The results of each step, returned by ppath_step_agenda as a new
5232 // ppath_step, are stored as an array of Ppath structures.
5233 //
5234 Array<Ppath> ppath_array(0);
5235 Index np = 1; // Counter for number of points of the path
5236 Index istep = 0; // Counter for number of steps
5237 //
5238 const Index imax_p = p_grid.nelem() - 1;
5239 const Index imax_lat = lat_grid.nelem() - 1;
5240 const Index imax_lon = lon_grid.nelem() - 1;
5241 //
5242 bool ready = ppath_what_background(ppath_step);
5243 //
5244 while (!ready) {
5245 // Call ppath_step agenda.
5246 // The new path step is added to *ppath_array* last in the while block
5247 //
5248 istep++;
5249 //
5251 ppath_step,
5252 ppath_lmax,
5253 ppath_lraytrace,
5254 f_grid,
5255 ppath_step_agenda);
5256 // For debugging:
5257 //Print( ppath_step, 0, verbosity );
5258
5259 // Number of points in returned path step
5260 const Index n = ppath_step.np;
5261
5262 // Increase the total number
5263 np += n - 1;
5264
5265 ARTS_USER_ERROR_IF (istep > 10000,
5266 "10 000 path points have been reached. Is this an infinite loop?");
5267
5268 //----------------------------------------------------------------------
5269 //--- Check if some boundary is reached
5270 //----------------------------------------------------------------------
5271
5272 //--- Outside cloud box ------------------------------------------------
5273 if (!ppath_inside_cloudbox_do) {
5274 // Check if the top of the atmosphere is reached
5275 // (Perform first a strict check, and if upward propagation also a
5276 // non-strict check. The later to avoid fatal "fixes" at TOA.)
5277 if (is_gridpos_at_index_i(ppath_step.gp_p[n - 1], imax_p) ||
5278 (abs(ppath_step.los(n - 1, 0)) < 90 &&
5279 is_gridpos_at_index_i(ppath_step.gp_p[n - 1], imax_p, false))) {
5280 ppath_set_background(ppath_step, 1);
5281 }
5282
5283 // Check that path does not exit at a latitude or longitude end face
5284 if (atmosphere_dim == 2) {
5285 // Latitude
5286 ARTS_USER_ERROR_IF (is_gridpos_at_index_i(ppath_step.gp_lat[n - 1], 0),
5287 "The path exits the atmosphere through the lower "
5288 "latitude end face.\nThe exit point is at an altitude"
5289 " of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5290 ARTS_USER_ERROR_IF (is_gridpos_at_index_i(ppath_step.gp_lat[n - 1], imax_lat),
5291 "The path exits the atmosphere through the upper "
5292 "latitude end face.\nThe exit point is at an altitude"
5293 " of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5294 }
5295 if (atmosphere_dim == 3) {
5296 // Latitude
5297 ARTS_USER_ERROR_IF (lat_grid[0] > -90 &&
5298 is_gridpos_at_index_i(ppath_step.gp_lat[n - 1], 0),
5299 "The path exits the atmosphere through the lower "
5300 "latitude end face.\nThe exit point is at an altitude"
5301 " of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5302 ARTS_USER_ERROR_IF (lat_grid[imax_lat] < 90 &&
5303 is_gridpos_at_index_i(ppath_step.gp_lat[n - 1], imax_lat),
5304 "The path exits the atmosphere through the upper "
5305 "latitude end face.\nThe exit point is at an altitude"
5306 " of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5307
5308 // Longitude
5309 // Note that it must be if and else if here. Otherwise e.g. -180
5310 // will be shifted to 180 and then later back to -180.
5311 if (is_gridpos_at_index_i(ppath_step.gp_lon[n - 1], 0) &&
5312 ppath_step.los(n - 1, 1) < 0 &&
5313 abs(ppath_step.pos(n - 1, 1)) < 90) {
5314 // Check if the longitude point can be shifted +360 degrees
5315 if (lon_grid[imax_lon] - lon_grid[0] >= 360) {
5316 ppath_step.pos(n - 1, 2) = ppath_step.pos(n - 1, 2) + 360;
5317 gridpos(
5318 ppath_step.gp_lon[n - 1], lon_grid, ppath_step.pos(n - 1, 2));
5319 } else {
5321 "The path exits the atmosphere through the lower "
5322 "longitude end face.\nThe exit point is at an "
5323 "altitude of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5324 }
5325 } else if (is_gridpos_at_index_i(ppath_step.gp_lon[n - 1], imax_lon) &&
5326 ppath_step.los(n - 1, 1) > 0 &&
5327 abs(ppath_step.pos(n - 1, 1)) < 90) {
5328 // Check if the longitude point can be shifted -360 degrees
5329 if (lon_grid[imax_lon] - lon_grid[0] >= 360) {
5330 ppath_step.pos(n - 1, 2) = ppath_step.pos(n - 1, 2) - 360;
5331 gridpos(
5332 ppath_step.gp_lon[n - 1], lon_grid, ppath_step.pos(n - 1, 2));
5333 } else {
5335 "The path exits the atmosphere through the upper "
5336 "longitude end face.\nThe exit point is at an "
5337 "altitude of ", ppath_step.pos(n - 1, 0) / 1e3, " km.")
5338 }
5339 }
5340 }
5341
5342 // Check if there is an intersection with an active cloud box
5343 if (cloudbox_on) {
5344 Numeric ipos = fractional_gp(ppath_step.gp_p[n - 1]);
5345
5346 if (ipos >= Numeric(cloudbox_limits[0]) &&
5347 ipos <= Numeric(cloudbox_limits[1])) {
5348 if (atmosphere_dim == 1) {
5349 ppath_set_background(ppath_step, 3);
5350 } else {
5351 ipos = fractional_gp(ppath_step.gp_lat[n - 1]);
5352
5353 if (ipos >= Numeric(cloudbox_limits[2]) &&
5354 ipos <= Numeric(cloudbox_limits[3])) {
5355 if (atmosphere_dim == 2) {
5356 ppath_set_background(ppath_step, 3);
5357 } else {
5358 ipos = fractional_gp(ppath_step.gp_lon[n - 1]);
5359
5360 if (ipos >= Numeric(cloudbox_limits[4]) &&
5361 ipos <= Numeric(cloudbox_limits[5])) {
5362 ppath_set_background(ppath_step, 3);
5363 }
5364 }
5365 }
5366 }
5367 }
5368 }
5369 }
5370
5371 //--- Inside cloud box -------------------------------------------------
5372 else {
5373 // A first version just checked if point was at or outside any
5374 // boundary but numerical problems could cause that the start point
5375 // was taken as the exit point. So check of ppath direction had to be
5376 // added. Fractional distances used for this.
5377
5378 // Pressure dimension
5379 Numeric ipos1 = fractional_gp(ppath_step.gp_p[n - 1]);
5380 Numeric ipos2 = fractional_gp(ppath_step.gp_p[n - 2]);
5381 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[0]);
5382 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[1]);
5383 if (ipos1 <= (Numeric)cloudbox_limits[0] && ipos1 < ipos2) {
5384 ppath_set_background(ppath_step, 3);
5385 }
5386
5387 else if (ipos1 >= Numeric(cloudbox_limits[1]) && ipos1 > ipos2) {
5388 ppath_set_background(ppath_step, 3);
5389 }
5390
5391 else if (atmosphere_dim > 1) {
5392 // Latitude dimension
5393 ipos1 = fractional_gp(ppath_step.gp_lat[n - 1]);
5394 ipos2 = fractional_gp(ppath_step.gp_lat[n - 2]);
5395 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[2]);
5396 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[3]);
5397 if (ipos1 <= Numeric(cloudbox_limits[2]) && ipos1 < ipos2) {
5398 ppath_set_background(ppath_step, 3);
5399 }
5400
5401 else if (ipos1 >= Numeric(cloudbox_limits[3]) && ipos1 > ipos2) {
5402 ppath_set_background(ppath_step, 3);
5403 }
5404
5405 else if (atmosphere_dim > 2) {
5406 // Longitude dimension
5407 ipos1 = fractional_gp(ppath_step.gp_lon[n - 1]);
5408 ipos2 = fractional_gp(ppath_step.gp_lon[n - 2]);
5409 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[4]);
5410 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[5]);
5411 if (ipos1 <= Numeric(cloudbox_limits[4]) && ipos1 < ipos2) {
5412 ppath_set_background(ppath_step, 3);
5413 }
5414
5415 else if (ipos1 >= Numeric(cloudbox_limits[5]) && ipos1 > ipos2) {
5416 ppath_set_background(ppath_step, 3);
5417 }
5418 }
5419 }
5420 }
5421 //----------------------------------------------------------------------
5422 //--- End of boundary check
5423 //----------------------------------------------------------------------
5424
5425 // Ready?
5426 if (ppath_what_background(ppath_step)) {
5427 ppath_step.start_pos = ppath_step.pos(n - 1, joker);
5428 ppath_step.start_los = ppath_step.los(n - 1, joker);
5429 ready = true;
5430 }
5431
5432 // Put new ppath_step in ppath_array
5433 ppath_array.push_back(ppath_step);
5434
5435 } // End path steps
5436
5437 // Combine all structures in ppath_array to form the return Ppath structure.
5438 //
5439 ppath_init_structure(ppath, atmosphere_dim, np);
5440 //
5441 const Index na = ppath_array.nelem();
5442 //
5443 if (na == 0) // No path, just the starting point
5444 {
5445 ppath_copy(ppath, ppath_step, 1);
5446 // To set n for positions inside the atmosphere, ppath_step_agenda
5447 // must be called once. The later is not always the case. A fix to handle
5448 // those cases:
5449 if (ppath_what_background(ppath_step) > 1) {
5450 //Print( ppath_step, 0, verbosity );
5452 ppath_step,
5453 ppath_lmax,
5454 ppath_lraytrace,
5455 f_grid,
5456 ppath_step_agenda);
5457 ppath.nreal[0] = ppath_step.nreal[0];
5458 ppath.ngroup[0] = ppath_step.ngroup[0];
5459 }
5460 }
5461
5462 else // Otherwise, merge the array elelments
5463 {
5464 np = 0; // Now used as counter for points moved to ppath
5465 //
5466 for (Index i = 0; i < na; i++) {
5467 // For the first structure, the first point shall be included.
5468 // For later structures, the first point shall not be included, but
5469 // there will always be at least two points.
5470
5471 Index n = ppath_array[i].np;
5472
5473 // First index to include
5474 Index i1 = 1;
5475 if (i == 0) {
5476 i1 = 0;
5477 }
5478
5479 // Vectors and matrices that can be handled by ranges.
5480 ppath.r[Range(np, n - i1)] = ppath_array[i].r[Range(i1, n - i1)];
5481 ppath.pos(Range(np, n - i1), joker) =
5482 ppath_array[i].pos(Range(i1, n - i1), joker);
5483 ppath.los(Range(np, n - i1), joker) =
5484 ppath_array[i].los(Range(i1, n - i1), joker);
5485 ppath.nreal[Range(np, n - i1)] = ppath_array[i].nreal[Range(i1, n - i1)];
5486 ppath.ngroup[Range(np, n - i1)] =
5487 ppath_array[i].ngroup[Range(i1, n - i1)];
5488 ppath.lstep[Range(np - i1, n - 1)] = ppath_array[i].lstep;
5489
5490 // Grid positions must be handled by a loop
5491 for (Index j = i1; j < n; j++) {
5492 ppath.gp_p[np + j - i1] = ppath_array[i].gp_p[j];
5493 }
5494 if (atmosphere_dim >= 2) {
5495 for (Index j = i1; j < n; j++) {
5496 ppath.gp_lat[np + j - i1] = ppath_array[i].gp_lat[j];
5497 }
5498 }
5499 if (atmosphere_dim == 3) {
5500 for (Index j = i1; j < n; j++) {
5501 ppath.gp_lon[np + j - i1] = ppath_array[i].gp_lon[j];
5502 }
5503 }
5504
5505 // Increase number of points done
5506 np += n - i1;
5507 }
5508
5509 // Field just included once:
5510 // end_pos/los/lspace from first path_step (extracted above):
5511 ppath.end_lstep = end_lstep;
5512 ppath.end_pos = end_pos;
5513 ppath.end_los = end_los;
5514 // Constant, background and start_pos/los from last path_step:
5515 ppath.constant = ppath_step.constant;
5516 ppath.background = ppath_step.background;
5517 ppath.start_pos = ppath_step.start_pos;
5518 ppath.start_los = ppath_step.start_los;
5519 ppath.start_lstep = ppath_step.start_lstep;
5520 }
5521}
5522
5523
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
Definition array.h:128
base min(const Array< base > &x)
Min function.
Definition array.h:144
Common ARTS conversions.
Header file for helper functions for OpenMP.
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
Definition auto_md.cc:27543
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
The Agenda class.
This can be used to make arrays out of anything.
Definition array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Workspace class.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define DEBUG_ONLY(...)
Definition debug.h:55
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
#define ARTS_USER_ERROR(...)
Definition debug.h:153
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
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:585
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition geodetic.cc:1248
void latlon_at_aa(Numeric &lat2, Numeric &lon2, const Numeric &lat1, const Numeric &lon1, const Numeric &aa, const Numeric &ddeg)
latlon_at_aa
Definition geodetic.cc:996
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition geodetic.cc:338
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition geodetic.cc:1287
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:96
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition geodetic.cc:57
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Numeric sign(const Numeric &x)
sign
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition messages.h:187
This file contains the definition of String, the ARTS string class.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition poly_roots.cc:73
Contains the code to determine roots of polynomials.
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:3885
void ppath_end_1d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &endface, const Numeric &ppc)
Internal help function for 1D path calculations.
Definition ppath.cc:1635
void ppath_end_3d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView lon_v, ConstVectorView za_v, ConstVectorView aa_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &ilon, const Index &endface, const Numeric &ppc)
Internal help function for 3D path calculations.
Definition ppath.cc:2106
Numeric rsurf_at_lat(const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const Numeric &lat)
Determines the radius of a pressure level or the surface given the radius at the corners of a 2D grid...
Definition ppath.cc:555
bool is_los_downwards(const Numeric &za, const Numeric &tilt)
Determines if a line-of-sight is downwards compared to the angular tilt of the surface or a pressure ...
Definition ppath.cc:606
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:4300
void ppath_start_3d(Numeric &r_start, Numeric &lat_start, Numeric &lon_start, Numeric &za_start, Numeric &aa_start, Index &ip, Index &ilat, Index &ilon, Numeric &lat1, Numeric &lat3, Numeric &lon5, Numeric &lon6, Numeric &r15a, Numeric &r35a, Numeric &r36a, Numeric &r16a, Numeric &r15b, Numeric &r35b, Numeric &r36b, Numeric &r16b, Numeric &rsurface15, Numeric &rsurface35, Numeric &rsurface36, Numeric &rsurface16, Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface)
Internal help function for 3D path calculations.
Definition ppath.cc:1916
const Numeric LON_NOT_FOUND
Definition ppath.cc:62
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:296
void do_gridrange_1d(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lmax, const Numeric &ra, const Numeric &rb, const Numeric &rsurface)
Calculates the geometrical path through a 1D grid range.
Definition ppath.cc:2265
Numeric geompath_r_at_l(const Numeric &ppc, const Numeric &l)
Calculates the radius for a distance from the tangent point.
Definition ppath.cc:153
void raytrace_3d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &lon_array, Array< Numeric > &za_array, Array< Numeric > &aa_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView refellipsoid, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, Numeric r, Numeric lat, Numeric lon, Numeric za, Numeric aa)
Performs ray tracing for 3D with linear steps.
Definition ppath.cc:4060
void ppath_start_1d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, const Ppath &ppath)
Internal help function for 1D path calculations.
Definition ppath.cc:1606
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition ppath.cc:524
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:413
const Numeric LAT_NOT_FOUND
Definition ppath.cc:61
void raytrace_1d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &rsurface, const Numeric &r1, const Numeric &r3, Numeric r, Numeric lat, Numeric za)
Performs ray tracing for 1D with linear steps.
Definition ppath.cc:3397
void plevel_crossing_2d(Numeric &r, Numeric &lat, Numeric &l, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const bool &above)
Handles the crossing with a geometric ppaths step and a atmospheric grid box level for 2D.
Definition ppath.cc:838
void do_gridcell_3d_byltest(Vector &r_v, Vector &lat_v, Vector &lon_v, Vector &za_v, Vector &aa_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &lon_start0, const Numeric &za_start, const Numeric &aa_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16)
See ATD for a description of the algorithm.
Definition ppath.cc:2793
constexpr Numeric DEG2RAD
Definition ppath.cc:37
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition ppath.cc:1460
Index first_pos_before_altitude(const Ppath &p, const Numeric &alt)
Determines ppath position just below an altitude.
Definition ppath.cc:506
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:3235
Numeric geompath_r_at_lat(const Numeric &ppc, const Numeric &lat0, const Numeric &za0, const Numeric &lat)
Calculates the radius for a given latitude.
Definition ppath.cc:171
void ppath_append(Ppath &ppath1, const Ppath &ppath2)
Combines two Ppath structures.
Definition ppath.cc:1546
void rotationmat3D(Matrix &R, ConstVectorView vrot, const Numeric &a)
Creates a 3D rotation matrix.
Definition ppath.cc:346
Numeric geometrical_ppc(const Numeric &r, const Numeric &za)
Calculates the propagation path constant for pure geometrical calculations.
Definition ppath.cc:80
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:375
void zaaa2enu(Numeric &de, Numeric &dn, Numeric &du, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to ENU unit vector.
Definition ppath.cc:319
void ppath_copy(Ppath &ppath1, const Ppath &ppath2, const Index &ncopy)
Copy the content in ppath2 to ppath1.
Definition ppath.cc:1480
void enu2zaaa(Numeric &za, Numeric &aa, const Numeric &de, const Numeric &dn, const Numeric &du)
Converts ENU unit vector vector to zenith and azimuth.
Definition ppath.cc:310
void r_crossing_3d(Numeric &lat, Numeric &lon, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &lon_start, const Numeric &za_start, const Numeric &ppc, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Calculates where a 3D LOS crosses the specified radius.
Definition ppath.cc:1327
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:2701
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:87
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:283
constexpr Numeric RAD2DEG
Definition ppath.cc:38
void ppath_end_2d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &endface, const Numeric &ppc)
Internal help function for 2D path calculations.
Definition ppath.cc:1809
Numeric rslope_crossing3d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1, Numeric c2)
3D version of rslope_crossing2d.
Definition ppath.cc:1203
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
Definition ppath.cc:494
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
Definition ppath.cc:1434
const Numeric LACC
Definition ppath.cc:56
void plevel_slope_3d(Numeric &c1, Numeric &c2, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon, const Numeric &aa)
Definition ppath.cc:1051
void r_crossing_2d(Numeric &lat, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc)
Calculates where a 2D LOS crosses the specified radius.
Definition ppath.cc:638
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
Resolves which longitude angle that shall be used.
Definition ppath.cc:484
Numeric refraction_ppc(const Numeric &r, const Numeric &za, const Numeric &refr_index_air)
Calculates the propagation path constant for cases with refraction.
Definition ppath.cc:475
Numeric geompath_lat_at_za(const Numeric &za0, const Numeric &lat0, const Numeric &za)
Calculates the latitude for a given zenith angle along a geometrical propagation path.
Definition ppath.cc:132
void do_gridcell_2d_byltest(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &za_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, const Numeric &rsurface1, const Numeric &rsurface3)
Works as do_gridcell_3d_byltest, but downscaled to 2D.
Definition ppath.cc:2404
Numeric geompath_r_at_za(const Numeric &ppc, const Numeric &za)
Calculates the zenith angle for a given radius along a geometrical propagation path.
Definition ppath.cc:125
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:5162
const Numeric R_NOT_FOUND
Definition ppath.cc:59
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:2337
Numeric rsurf_at_latlon(const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon)
Determines the radius of a pressure level or the surface given the radius at the corners of a 3D grid...
Definition ppath.cc:1022
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:1393
Numeric plevel_angletilt(const Numeric &r, const Numeric &c1)
Calculates the angular tilt of the surface or a pressure level.
Definition ppath.cc:600
const Numeric L_NOT_FOUND
Definition ppath.cc:60
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:3541
void ppath_start_stepping(Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const bool &ppath_inside_cloudbox_do, ConstVectorView rte_pos, ConstVectorView rte_los, const Verbosity &verbosity)
Initiates a Ppath structure for calculation of a path with ppath_step.
Definition ppath.cc:4460
void ppath_start_2d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, Index &ilat, Numeric &lat1, Numeric &lat3, Numeric &r1a, Numeric &r3a, Numeric &r3b, Numeric &r1b, Numeric &rsurface1, Numeric &rsurface3, Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface)
Internal help function for 2D path calculations.
Definition ppath.cc:1699
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
Calculates the length from the tangent point for the given radius.
Definition ppath.cc:142
void geompath_from_r1_to_r2(Vector &r, Vector &lat, Vector &za, Numeric &lstep, const Numeric &ppc, const Numeric &r1, const Numeric &lat1, const Numeric &za1, const Numeric &r2, const bool &tanpoint, const Numeric &lmax)
Determines radii, latitudes and zenith angles between two points of a propagation path.
Definition ppath.cc:207
void plevel_slope_2d(Numeric &c1, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstVectorView z_surf, const GridPos &gp, const Numeric &za)
Calculates the radial slope of the surface or a pressure level for 2D.
Definition ppath.cc:563
void raytrace_2d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &rsurface1, const Numeric &rsurface3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, Numeric r, Numeric lat, Numeric za)
Performs ray tracing for 2D with linear steps.
Definition ppath.cc:3703
const Numeric RTOL
Definition ppath.cc:48
Numeric rslope_crossing2d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1)
Calculates the angular distance to a crossing with a level having a radial slope.
Definition ppath.cc:709
Propagation path structure and functions.
constexpr Numeric POLELAT
Size of north and south poles.
Definition ppath.h:28
constexpr Numeric ANGTOL
Width of zenith and nadir directions.
Definition ppath.h:39
void refr_gradients_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, Numeric &dndlon, 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)
refr_gradients_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)
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
void refr_gradients_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, 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)
refr_gradients_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
void refr_gradients_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, 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)
refr_gradients_1d
Refraction functions.
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:66
Declaration of functions in rte.cc.
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.
Header file for special_interp.cc.
Structure to store a grid position.
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
String background
Radiative background.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Numeric end_lstep
The distance between end pos and the first position in pos.
Vector start_pos
Start position.
Vector lstep
The length between ppath points.
Numeric start_lstep
Length between sensor and atmospheric boundary.
Numeric constant
The propagation path constant (only used for 1D)
Vector r
Radius of each ppath point.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Vector ngroup
The group index of refraction.
Index dim
Atmospheric dimensionality.
Vector end_pos
End position.
Vector start_los
Start line-of-sight.
Vector nreal
The real part of the refractive index at each path position.
Vector end_los
End line-of-sight.
#define u
#define v
#define w
#define a
#define c