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