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