ARTS 2.5.11 (git: 725533f0)
geodetic.cc
Go to the documentation of this file.
1/*===========================================================================
2 === File description
3 ===========================================================================*/
4
14/*===========================================================================
15 === External declarations
16 ===========================================================================*/
17
18#include "geodetic.h"
19#include <cmath>
20#include <stdexcept>
21#include "arts_conversions.h"
22#include "math_funcs.h"
23#include "ppath.h"
24
25inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
26inline constexpr Numeric RAD2DEG=Conversion::rad2deg(1);
27
28/*===========================================================================
29 === 2D functions
30 ===========================================================================*/
31
32// The 2D case is treated as being the 3D x/z-plane. That is, y-coordinate is
33// skipped. For simplicity, the angle coordinate is denoted as latitude.
34// However, the latitude is here not limited to [-90,90]. It is cyclic and can
35// have any value. The input *lat0* is used to shift the output from atan2 with
36// n*360 to return what should be the expected latitude. That is, it is assumed
37// that no operation moves the latitude more than 180 degrees from the initial
38// value *lat0*. Negative zeniath angles are handled, following ARTS definition
39// of 2D geometry.
40
42
57void cart2pol(Numeric& r,
58 Numeric& lat,
59 const Numeric& x,
60 const Numeric& z,
61 const Numeric& lat0,
62 const Numeric& za0) {
63 r = sqrt(x * x + z * z);
64
65 // Zenith and nadir cases
66 const Numeric absza = abs(za0);
67 if (absza < ANGTOL || absza > 180 - ANGTOL) {
68 lat = lat0;
69 }
70
71 else { // Latitude inside [0,360]
72 lat = RAD2DEG * atan2(z, x);
73 // Shift with n*360 to get as close as possible to lat0
74 lat = lat - 360.0 * Numeric(round((lat - lat0) / 360.0));
75 }
76}
77
79
96void cart2poslos(Numeric& r,
97 Numeric& lat,
98 Numeric& za,
99 const Numeric& x,
100 const Numeric& z,
101 const Numeric& dx,
102 const Numeric& dz,
103 const Numeric& ppc,
104 const Numeric& lat0,
105 const Numeric& za0) {
106 r = sqrt(x * x + z * z);
107
108 // Zenith and nadir cases
109 const Numeric absza = abs(za0);
110 if (absza < ANGTOL || absza > 180 - ANGTOL) {
111 lat = lat0;
112 za = za0;
113 }
114
115 else {
116 lat = RAD2DEG * atan2(z, x);
117
118 const Numeric latrad = DEG2RAD * lat;
119 const Numeric coslat = cos(latrad);
120 const Numeric sinlat = sin(latrad);
121 const Numeric dr = coslat * dx + sinlat * dz;
122
123 // Use ppc for max accuracy, but dr required to resolve if up-
124 // and downward cases.
125
126 // Another possibility to obtain (absolute value of) za is
127 // RAD2DEG*acos(dr).
128 // It is checked that the two ways give consistent results, but
129 // occasionally deviate with 1e-4 deg (due to numerical issues).
130
131 za = RAD2DEG * asin(ppc / r);
132 if (za0 > 0) {
133 if (std::isnan(za)) {
134 za = 90;
135 } else if (dr < 0) {
136 za = 180.0 - za;
137 }
138 } else {
139 if (std::isnan(za)) {
140 za = -90;
141 } else if (dr < 0) {
142 za = -180.0 + za;
143 } else {
144 za = -za;
145 }
146 }
147 }
148}
149
151
165void distance2D(Numeric& l,
166 const Numeric& r1,
167 const Numeric& lat1,
168 const Numeric& r2,
169 const Numeric& lat2) {
170 ARTS_ASSERT(abs(lat2 - lat1) <= 180);
171
172 Numeric x1, z1, x2, z2;
173 pol2cart(x1, z1, r1, lat1);
174 pol2cart(x2, z2, r2, lat2);
175
176 const Numeric dx = x2 - x1;
177 const Numeric dz = z2 - z1;
178 l = sqrt(dx * dx + dz * dz);
179}
180
182
212/*
213void geomtanpoint2d(
214 Numeric& r_tan,
215 Numeric& lat_tan,
216 ConstVectorView refellipsoid,
217 const Numeric& r,
218 const Numeric& lat,
219 const Numeric& za )
220{
221 ARTS_ASSERT( refellipsoid.nelem() == 2 );
222 ARTS_ASSERT( refellipsoid[0] > 0 );
223 ARTS_ASSERT( refellipsoid[1] >= 0 );
224 ARTS_ASSERT( refellipsoid[1] < 1 );
225 ARTS_ASSERT( r > 0 );
226 ARTS_ASSERT( za >= 90 );
227 // e=1e-7 corresponds to that polar radius
228 if( refellipsoid[1] < 1e-7 ) // less than 1 um smaller than equatorial
229 { // one for the Earth
230 r_tan = geometrical_ppc( r, za );
231 if( za > 0 )
232 { lat_tan = geompath_lat_at_za( za, lat, 90 ); }
233 else
234 { lat_tan = geompath_lat_at_za( za, lat, -90 ); }
235 }
236
237 else
238 {
239 ARTS_ASSERT( 0 ); // To be implemented
240 }
241}
242*/
243
245
261void line_circle_intersect(Numeric& x,
262 Numeric& z,
263 const Numeric& xl,
264 const Numeric& zl,
265 const Numeric& dx,
266 const Numeric& dz,
267 const Numeric& xc,
268 const Numeric& zc,
269 const Numeric& r) {
270 const Numeric a = dx * dx + dz * dz;
271 const Numeric b = 2 * (dx * (xl - xc) + dz * (zl - zc));
272 const Numeric c =
273 xc * xc + zc * zc + xl * xl + zl * zl - 2 * (xc * xl + zc * zl) - r * r;
274
275 Numeric d = b * b - 4 * a * c;
276 ARTS_ASSERT(d > 0);
277
278 const Numeric a2 = 2 * a;
279 const Numeric b2 = -b / a2;
280 const Numeric e = sqrt(d) / a2;
281
282 const Numeric l1 = b2 + e;
283 const Numeric l2 = b2 - e;
284
285 Numeric l;
286 if (l1 < 0) {
287 l = l2;
288 } else if (l2 < 0) {
289 l = l1;
290 } else {
291 l = min(l1, l2);
292 ARTS_ASSERT(l >= 0);
293 }
294
295 x = xl + l * dx;
296 z = zl + l * dz;
297}
298
300
314void pol2cart(Numeric& x, Numeric& z, const Numeric& r, const Numeric& lat) {
315 ARTS_ASSERT(r > 0);
316
317 const Numeric latrad = DEG2RAD * lat;
318
319 x = r * cos(latrad);
320 z = r * sin(latrad);
321}
322
324
338void poslos2cart(Numeric& x,
339 Numeric& z,
340 Numeric& dx,
341 Numeric& dz,
342 const Numeric& r,
343 const Numeric& lat,
344 const Numeric& za) {
345 ARTS_ASSERT(r > 0);
346 ARTS_ASSERT(za >= -180 && za <= 180);
347
348 const Numeric latrad = DEG2RAD * lat;
349 const Numeric zarad = DEG2RAD * za;
350
351 const Numeric coslat = cos(latrad);
352 const Numeric sinlat = sin(latrad);
353 const Numeric cosza = cos(zarad);
354 const Numeric sinza = sin(zarad);
355
356 // This part as pol2cart but uses local variables
357 x = r * coslat;
358 z = r * sinlat;
359
360 const Numeric dr = cosza;
361 const Numeric dlat = sinza; // r-term cancel out below
362
363 dx = coslat * dr - sinlat * dlat;
364 dz = sinlat * dr + coslat * dlat;
365}
366
367/*===========================================================================
368 === 3D functions
369 ===========================================================================*/
370
372
404void cart2poslos(Numeric& r,
405 Numeric& lat,
406 Numeric& lon,
407 Numeric& za,
408 Numeric& aa,
409 const Numeric& x,
410 const Numeric& y,
411 const Numeric& z,
412 const Numeric& dx,
413 const Numeric& dy,
414 const Numeric& dz,
415 const Numeric& ppc,
416 const Numeric& x0,
417 const Numeric& y0,
418 const Numeric& z0,
419 const Numeric& lat0,
420 const Numeric& lon0,
421 const Numeric& za0,
422 const Numeric& aa0) {
423 // Radius of new point
424 r = sqrt(x * x + y * y + z * z);
425
426 // Zenith and nadir cases
427 if (za0 < ANGTOL || za0 > 180 - ANGTOL) {
428 lat = lat0;
429 lon = lon0;
430 za = za0;
431 aa = aa0;
432 }
433
434 else {
435 lat = RAD2DEG * asin(z / r);
436 lon = RAD2DEG * atan2(y, x);
437
438 bool ns_case = false;
439 bool lon_flip = false;
440
441 // Make sure that lon is maintained for N-S cases (if not
442 // starting on a pole)
443 if ((abs(aa0) < ANGTOL || abs(180 - aa0) < ANGTOL) &&
444 abs(lat0) <= POLELAT) {
445 ns_case = true;
446 // Check that not lon changed with 180 deg
447 if (abs(abs(lon - lon0) - 180) < 5) {
448 lon_flip = true;
449 if (lon0 > 0) {
450 lon = lon0 - 180;
451 } else {
452 lon = lon0 + 180;
453 }
454 } else {
455 lon = lon0;
456 }
457 }
458
459 const Numeric latrad = DEG2RAD * lat;
460 const Numeric lonrad = DEG2RAD * lon;
461 const Numeric coslat = cos(latrad);
462 const Numeric sinlat = sin(latrad);
463 const Numeric coslon = cos(lonrad);
464 const Numeric sinlon = sin(lonrad);
465
466 // Set za by ppc for max accuracy, but this does not resolve
467 // za and 180-za. This was first resolved by dr, but using l and lmax was
468 // found to be more stable.
469 za = RAD2DEG * asin(ppc / r);
470
471 // Correct and check za
472 if (std::isnan(za)) {
473 za = 90;
474 }
475 // If za0 > 90, then correct za could be 180-za. Resolved by checking if
476 // the tangent point is passed or not
477 if (za0 > 90) {
478 const Numeric l =
479 sqrt(pow(x - x0, 2.0) + pow(y - y0, 2.0) + pow(z - z0, 2.0));
480 const Numeric r0 = sqrt(x0 * x0 + y0 * y0 + z0 * z0);
481 const Numeric ltan = geompath_l_at_r(ppc, r0);
482 if (l < ltan) {
483 za = 180.0 - za;
484 }
485 }
486
487 // For lat = +- 90 the azimuth angle gives the longitude along which
488 // the LOS goes
489 if (abs(lat) >= POLELAT) {
490 aa = RAD2DEG * atan2(dy, dx);
491 }
492
493 // N-S cases, not starting at a pole
494 else if (ns_case) {
495 if (!lon_flip) {
496 aa = aa0;
497 } else {
498 if (abs(aa0) < ANGTOL) {
499 aa = 180;
500 } else {
501 aa = 0;
502 }
503 }
504 }
505
506 else {
507 const Numeric dlat = -sinlat * coslon / r * dx -
508 sinlat * sinlon / r * dy + coslat / r * dz;
509 const Numeric dlon = -sinlon / coslat / r * dx + coslon / coslat / r * dy;
510
511 aa = RAD2DEG * acos(r * dlat / sin(DEG2RAD * za));
512
513 if (std::isnan(aa)) {
514 if (dlat >= 0) {
515 aa = 0;
516 } else {
517 aa = 180;
518 }
519 } else if (dlon < 0) {
520 aa = -aa;
521 }
522 }
523 }
524}
525void cart2poslos_plain(Numeric& r,
526 Numeric& lat,
527 Numeric& lon,
528 Numeric& za,
529 Numeric& aa,
530 const Numeric& x,
531 const Numeric& y,
532 const Numeric& z,
533 const Numeric& dx,
534 const Numeric& dy,
535 const Numeric& dz) {
536 cart2sph_plain(r, lat, lon, x, y, z);
537
538
539 const Numeric latrad = DEG2RAD * lat;
540 const Numeric lonrad = DEG2RAD * lon;
541 const Numeric coslat = cos(latrad);
542 const Numeric sinlat = sin(latrad);
543 const Numeric coslon = cos(lonrad);
544 const Numeric sinlon = sin(lonrad);
545
546 const Numeric dr = coslat*coslon*dx + sinlat*dz + coslat*sinlon*dy;
547 const Numeric dlat = -sinlat*coslon/r*dx + coslat/r*dz - sinlat*sinlon/r*dy;
548 const Numeric dlon = -sinlon/coslat/r*dx + coslon/coslat/r*dy;
549
550 za = acos( dr );
551 aa = RAD2DEG * acos( r * dlat / sin( za ) );
552 za *= RAD2DEG;
553
554 // Corrections of aa
555 if (std::isnan(aa)) {
556 if (dlat >= 0)
557 aa = 0;
558 else
559 aa = 180;
560 } else if (dlon < 0) {
561 aa = -aa;
562 }
563}
564
566
585void cart2sph(Numeric& r,
586 Numeric& lat,
587 Numeric& lon,
588 const Numeric& x,
589 const Numeric& y,
590 const Numeric& z,
591 const Numeric& lat0,
592 const Numeric& lon0,
593 const Numeric& za0,
594 const Numeric& aa0) {
595 r = sqrt(x * x + y * y + z * z);
596
597 // Zenith and nadir cases
598 if (za0 < ANGTOL || za0 > 180 - ANGTOL) {
599 lat = lat0;
600 lon = lon0;
601 }
602
603 else {
604 lat = RAD2DEG * asin(z / r);
605 lon = RAD2DEG * atan2(y, x);
606
607 // Make sure that lon is maintained for N-S cases (if not
608 // starting on a pole)
609 if ((abs(aa0) < ANGTOL || abs(180 - aa0) < ANGTOL) &&
610 abs(lat0) <= POLELAT) {
611 // Check that not lon changed with 180 deg
612 if (abs(lon - lon0) < 1) {
613 lon = lon0;
614 } else {
615 if (lon0 > 0) {
616 lon = lon0 - 180;
617 } else {
618 lon = lon0 + 180;
619 }
620 }
621 }
622 }
623}
624void cart2sph_plain(Numeric& r,
625 Numeric& lat,
626 Numeric& lon,
627 const Numeric& x,
628 const Numeric& y,
629 const Numeric& z) {
630 r = sqrt(x * x + y * y + z * z);
631 lat = RAD2DEG * asin(z / r);
632 lon = RAD2DEG * atan2(y, x);
633}
634
635
636
638
652void distance3D(Numeric& l,
653 const Numeric& r1,
654 const Numeric& lat1,
655 const Numeric& lon1,
656 const Numeric& r2,
657 const Numeric& lat2,
658 const Numeric& lon2) {
659 Numeric x1, y1, z1, x2, y2, z2;
660 sph2cart(x1, y1, z1, r1, lat1, lon1);
661 sph2cart(x2, y2, z2, r2, lat2, lon2);
662
663 const Numeric dx = x2 - x1;
664 const Numeric dy = y2 - y1;
665 const Numeric dz = z2 - z1;
666 l = sqrt(dx * dx + dy * dy + dz * dz);
667}
668
670
689void geompath_tanpos_3d(Numeric& r_tan,
690 Numeric& lat_tan,
691 Numeric& lon_tan,
692 Numeric& l_tan,
693 const Numeric& r,
694 const Numeric& lat,
695 const Numeric& lon,
696 const Numeric& za,
697 const Numeric& aa,
698 const Numeric& ppc) {
699 ARTS_ASSERT(za >= 90);
700 ARTS_ASSERT(r >= ppc);
701
702 Numeric x, y, z, dx, dy, dz;
703
704 poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
705
706 l_tan = sqrt(r * r - ppc * ppc);
707
708 cart2sph(r_tan,
709 lat_tan,
710 lon_tan,
711 x + dx * l_tan,
712 y + dy * l_tan,
713 z + dz * l_tan,
714 lat,
715 lon,
716 za,
717 aa);
718}
719
721
751/*
752void geomtanpoint(
753 Numeric& r_tan,
754 Numeric& lat_tan,
755 Numeric& lon_tan,
756 ConstVectorView refellipsoid,
757 const Numeric& r,
758 const Numeric& lat,
759 const Numeric& lon,
760 const Numeric& za,
761 const Numeric& aa )
762{
763 ARTS_ASSERT( refellipsoid.nelem() == 2 );
764 ARTS_ASSERT( refellipsoid[0] > 0 );
765 ARTS_ASSERT( refellipsoid[1] >= 0 );
766 ARTS_ASSERT( refellipsoid[1] < 1 );
767 ARTS_ASSERT( r > 0 );
768 ARTS_ASSERT( za >= 90 );
769
770 if( refellipsoid[1] < 1e-7 ) // e=1e-7 corresponds to that polar radius
771 { // less than 1 um smaller than equatorial
772 Numeric x, y, z, dx, dy, dz; // one for the Earth
773
774 poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
775
776 const Numeric ppc = r * sin( DEG2RAD * abs(za) );
777 const Numeric l_tan = sqrt( r*r - ppc*ppc );
778
779 cart2sph( r_tan, lat_tan, lon_tan, x+dx*l_tan, y+dy*l_tan, z+dz*l_tan );
780 }
781
782 else
783 {
784 // Equatorial and polar radii squared:
785 const Numeric a2 = refellipsoid[0]*refellipsoid[0];
786 const Numeric b2 = a2 * ( 1 - refellipsoid[1]*refellipsoid[1] );
787
788 Vector X(3), xunit(3), yunit(3), zunit(3);
789
790 poslos2cart( X[0], X[1], X[2], xunit[0], xunit[1], xunit[2],
791 r, lat, lon, za, aa );
792 cross( zunit, xunit, X );
793 unitl( zunit ); // Normalisation of length to 1
794
795 cross( yunit, zunit, xunit );
796 unitl( yunit ); // Normalisation of length to 1
797
798 Numeric x = X[0];
799 Numeric y = X[1];
800 const Numeric w11 = xunit[0];
801 const Numeric w12 = yunit[0];
802 const Numeric w21 = xunit[1];
803 const Numeric w22 = yunit[1];
804 const Numeric w31 = xunit[2];
805 const Numeric w32 = yunit[2];
806
807 const Numeric yr = X * yunit;
808 const Numeric xr = X * xunit;
809
810 const Numeric A = (w11*w11 + w21*w21)/a2 + w31*w31/b2;
811 const Numeric B = 2.0*((w11*w12 + w21*w22)/a2 + (w31*w32)/b2);
812 const Numeric C = (w12*w12 + w22*w22)/a2 + w32*w32/b2;
813
814 if( B == 0.0 )
815 { x = 0.0; }
816 else
817 {
818 const Numeric K = -2.0*A/B;
819 const Numeric factor = 1.0/(A+(B+C*K)*K);
820 x = sqrt(factor);
821 y = K*x;
822 }
823
824 const Numeric dist1 = (xr-X[0])*(xr-X[0]) + (yr-y)*(yr-y);
825 const Numeric dist2 = (xr+X[0])*(xr+X[0]) + (yr+y)*(yr+y);
826
827 if( dist1 > dist2 )
828 { x = -x; }
829
830 cart2sph( r_tan, lat_tan, lon_tan, w11*x + w12*yr, w21*x + w22*yr,
831 w31*x + w32*yr );
832 }
833}
834*/
835
837
857void line_refellipsoid_intersect(Numeric& l,
858 const Vector& refellipsoid,
859 const Numeric& x,
860 const Numeric& y,
861 const Numeric& z,
862 const Numeric& dx,
863 const Numeric& dy,
864 const Numeric& dz) {
865 // Code taken from Atmlab's ellipsoid_intersection
866
867 // Spherical case
868 if (refellipsoid[1] < 1e-7) {
869 const Numeric p = x*dx + y*dy + z*dz;
870 const Numeric pp = p*p;
871 const Numeric q = x*x + y*y + z*z - refellipsoid[0]*refellipsoid[0];
872 if (q>pp)
873 l = -1.0;
874 else {
875 const Numeric sq = sqrt(pp - q);
876 if (-p > sq)
877 l = -p - sq;
878 else
879 l = -p + sq;
880 }
881 }
882
883 // Ellipsoid case
884 else {
885 // Based on https://medium.com/@stephenhartzell/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6
886 const Numeric a = refellipsoid[0];
887 const Numeric b = refellipsoid[0];
888 const Numeric c = refellipsoid[0] * sqrt(1-refellipsoid[1]*refellipsoid[1]);
889 const Numeric a2 = a*a;
890 const Numeric b2 = b*b;
891 const Numeric c2 = c*c;
892 const Numeric x2 = x*x;
893 const Numeric y2 = y*y;
894 const Numeric z2 = z*z;
895 const Numeric dx2 = dx*dx;
896 const Numeric dy2 = dy*dy;
897 const Numeric dz2 = dz*dz;
898 const Numeric rad = a2*b2*dz2 + a2*c2*dy2 - a2*dy2*z2 + 2*a2*dy*dz*y*z -
899 a2*dz2*y2 + b2*c2*dx2 - b2*dx2*z2 + 2*b2*dx*dz*x*z -
900 b2*dz2*x2 - c2*dx2*y2 + 2*c2*dx*dy*x*y - c2*dy2*x2;
901 if (rad<0)
902 l = -1.0;
903 else {
904 const Numeric val = -a2*b2*dz*z - a2*c2*dy*y - b2*c2*dx*x;
905 const Numeric mag = a2*b2*dz2 + a2*c2*dy2 + b2*c2*dx2;
906 const Numeric abc = a*b*c*sqrt(rad);
907 if (val > abc)
908 l = (val - abc) / mag;
909 else
910 l = (val + abc) / mag;
911 }
912 }
913}
914
916
936void line_sphere_intersect(Numeric& x,
937 Numeric& y,
938 Numeric& z,
939 const Numeric& xl,
940 const Numeric& yl,
941 const Numeric& zl,
942 const Numeric& dx,
943 const Numeric& dy,
944 const Numeric& dz,
945 const Numeric& xc,
946 const Numeric& yc,
947 const Numeric& zc,
948 const Numeric& r) {
949 const Numeric a = dx * dx + dy * dy + dz * dz;
950 const Numeric b = 2 * (dx * (xl - xc) + dy * (yl - yc) + dz * (zl - zc));
951 const Numeric c = xc * xc + yc * yc + zc * zc + xl * xl + yl * yl + zl * zl -
952 2 * (xc * xl + yc * yl + zc * zl) - r * r;
953
954 Numeric d = b * b - 4 * a * c;
955 ARTS_ASSERT(d > 0);
956
957 const Numeric a2 = 2 * a;
958 const Numeric b2 = -b / a2;
959 const Numeric e = sqrt(d) / a2;
960
961 const Numeric l1 = b2 + e;
962 const Numeric l2 = b2 - e;
963
964 Numeric l;
965 if (l1 < 0) {
966 l = l2;
967 } else if (l2 < 0) {
968 l = l1;
969 } else {
970 l = min(l1, l2);
971 ARTS_ASSERT(l >= 0);
972 }
973
974 x = xl + l * dx;
975 y = yl + l * dy;
976 z = zl + l * dz;
977}
978
980
996void latlon_at_aa(Numeric& lat2,
997 Numeric& lon2,
998 const Numeric& lat1,
999 const Numeric& lon1,
1000 const Numeric& aa,
1001 const Numeric& ddeg) {
1002 // Code from http://www.movable-type.co.uk/scripts/latlong.html
1003 // (but with short-cuts, such as asin(sin(lat2)) = lat2)
1004 // Note that lat1 here is another latitude
1005
1006 const Numeric dang = DEG2RAD * ddeg;
1007 const Numeric cosdang = cos(dang);
1008 const Numeric sindang = sin(dang);
1009 const Numeric latrad = DEG2RAD * lat1;
1010 const Numeric coslat = cos(latrad);
1011 const Numeric sinlat = sin(latrad);
1012 const Numeric aarad = DEG2RAD * aa;
1013
1014 lat2 = sinlat * cosdang + coslat * sindang * cos(aarad);
1015 lon2 = lon1 + RAD2DEG * (atan2(sin(aarad) * sindang * coslat,
1016 cosdang - sinlat * lat2));
1017 lat2 = RAD2DEG * asin(lat2);
1018}
1019
1021
1042void los2xyz(Numeric& za,
1043 Numeric& aa,
1044 const Numeric& r1,
1045 const Numeric& lat1,
1046 const Numeric& lon1,
1047 const Numeric& x1,
1048 const Numeric& y1,
1049 const Numeric& z1,
1050 const Numeric& x2,
1051 const Numeric& y2,
1052 const Numeric& z2) {
1053 Numeric dx = x2 - x1, dy = y2 - y1, dz = z2 - z1;
1054 const Numeric ldxyz = sqrt(dx * dx + dy * dy + dz * dz);
1055 dx /= ldxyz;
1056 dy /= ldxyz;
1057 dz /= ldxyz;
1058
1059 // All below extracted from 3D version of cart2poslos:
1060 const Numeric latrad = DEG2RAD * lat1;
1061 const Numeric lonrad = DEG2RAD * lon1;
1062 const Numeric coslat = cos(latrad);
1063 const Numeric sinlat = sin(latrad);
1064 const Numeric coslon = cos(lonrad);
1065 const Numeric sinlon = sin(lonrad);
1066
1067 const Numeric dr = coslat * coslon * dx + coslat * sinlon * dy + sinlat * dz;
1068 const Numeric dlat =
1069 -sinlat * coslon / r1 * dx - sinlat * sinlon / r1 * dy + coslat / r1 * dz;
1070 const Numeric dlon = -sinlon / coslat / r1 * dx + coslon / coslat / r1 * dy;
1071
1072 za = RAD2DEG * acos(dr);
1073 aa = RAD2DEG * acos(r1 * dlat / sin(DEG2RAD * za));
1074 if (std::isnan(aa)) {
1075 if (dlat >= 0) {
1076 aa = 0;
1077 } else {
1078 aa = 180;
1079 }
1080 } else if (dlon < 0) {
1081 aa = -aa;
1082 }
1083}
1084
1086
1110void poslos2cart(Numeric& x,
1111 Numeric& y,
1112 Numeric& z,
1113 Numeric& dx,
1114 Numeric& dy,
1115 Numeric& dz,
1116 const Numeric& r,
1117 const Numeric& lat,
1118 const Numeric& lon,
1119 const Numeric& za,
1120 const Numeric& aa) {
1121 ARTS_ASSERT(r > 0);
1122 ARTS_ASSERT(abs(lat) <= 90);
1123 //ARTS_ASSERT( abs( lon ) <= 360 );
1124 ARTS_ASSERT(za >= 0 && za <= 180);
1125
1126 // lat = +-90
1127 // For lat = +- 90 the azimuth angle gives the longitude along which the
1128 // LOS goes
1129 if (abs(lat) > POLELAT) {
1130 const Numeric s = sign(lat);
1131
1132 x = 0;
1133 y = 0;
1134 z = s * r;
1135
1136 dz = s * cos(DEG2RAD * za);
1137 dx = sin(DEG2RAD * za);
1138 dy = dx * sin(DEG2RAD * aa);
1139 dx = dx * cos(DEG2RAD * aa);
1140 }
1141
1142 else {
1143 const Numeric latrad = DEG2RAD * lat;
1144 const Numeric lonrad = DEG2RAD * lon;
1145 const Numeric zarad = DEG2RAD * za;
1146 const Numeric aarad = DEG2RAD * aa;
1147
1148 const Numeric coslat = cos(latrad);
1149 const Numeric sinlat = sin(latrad);
1150 const Numeric coslon = cos(lonrad);
1151 const Numeric sinlon = sin(lonrad);
1152 const Numeric cosza = cos(zarad);
1153 const Numeric sinza = sin(zarad);
1154 const Numeric cosaa = cos(aarad);
1155 const Numeric sinaa = sin(aarad);
1156
1157 // This part as sph2cart but uses local variables
1158 x = r * coslat; // Common term for x and y
1159 y = x * sinlon;
1160 x = x * coslon;
1161 z = r * sinlat;
1162
1163 const Numeric dr = cosza;
1164 const Numeric dlat = sinza * cosaa; // r-term cancel out below
1165 const Numeric dlon = sinza * sinaa / coslat;
1166
1167 dx = coslat * coslon * dr - sinlat * coslon * dlat - coslat * sinlon * dlon;
1168 dz = sinlat * dr + coslat * dlat;
1169 dy = coslat * sinlon * dr - sinlat * sinlon * dlat + coslat * coslon * dlon;
1170 }
1171}
1172
1174
1194Numeric pos2refell_r(const Index& atmosphere_dim,
1195 ConstVectorView refellipsoid,
1196 ConstVectorView lat_grid,
1197 ConstVectorView lon_grid,
1198 ConstVectorView rte_pos) {
1199 if (atmosphere_dim == 1) {
1200 return refellipsoid[0];
1201 } else {
1202 ARTS_ASSERT(rte_pos.nelem() > 1);
1203
1204 bool inside = true;
1205
1206 if (rte_pos[1] < lat_grid[0] || rte_pos[1] > last(lat_grid)) {
1207 inside = false;
1208 } else if (atmosphere_dim == 3) {
1209 ARTS_ASSERT(rte_pos.nelem() == 3);
1210 if (rte_pos[2] < lon_grid[0] || rte_pos[2] > last(lon_grid)) {
1211 inside = false;
1212 }
1213 }
1214
1215 if (inside) {
1216 GridPos gp_lat;
1217 gridpos(gp_lat, lat_grid, rte_pos[1]);
1218 return refell2d(refellipsoid, lat_grid, gp_lat);
1219 } else {
1220 return refell2r(refellipsoid, rte_pos[1]);
1221 }
1222 }
1223}
1224
1226
1248Numeric refell2r(ConstVectorView refellipsoid, const Numeric& lat) {
1249 ARTS_ASSERT(refellipsoid.nelem() == 2);
1250 ARTS_ASSERT(refellipsoid[0] > 0);
1251 ARTS_ASSERT(refellipsoid[1] >= 0);
1252 ARTS_ASSERT(refellipsoid[1] < 1);
1253
1254 if (refellipsoid[1] < 1e-7) // e=1e-7 corresponds to that polar radius
1255 { // less than 1 um smaller than equatorial
1256 return refellipsoid[0]; // one for the Earth
1257 }
1258
1259 else {
1260 const Numeric c = 1 - refellipsoid[1] * refellipsoid[1];
1261 const Numeric b = refellipsoid[0] * sqrt(c);
1262 const Numeric v = DEG2RAD * lat;
1263 const Numeric ct = cos(v);
1264 const Numeric st = sin(v);
1265
1266 return b / sqrt(c * ct * ct + st * st);
1267 }
1268}
1269
1271
1287Numeric refell2d(ConstVectorView refellipsoid,
1288 ConstVectorView lat_grid,
1289 const GridPos gp) {
1290 if (gp.fd[0] == 0)
1291 return refell2r(refellipsoid, lat_grid[gp.idx]);
1292 else if (gp.fd[0] == 1)
1293 return refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1294 else
1295 return gp.fd[1] * refell2r(refellipsoid, lat_grid[gp.idx]) +
1296 gp.fd[0] * refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1297}
1298
1300
1318Numeric sphdist(const Numeric& lat1,
1319 const Numeric& lon1,
1320 const Numeric& lat2,
1321 const Numeric& lon2) {
1322 // Equations taken from http://www.movable-type.co.uk/scripts/latlong.html
1323 const Numeric slat = sin(DEG2RAD * (lat2 - lat1) / 2.0);
1324 const Numeric slon = sin(DEG2RAD * (lon2 - lon1) / 2.0);
1325 const Numeric a =
1326 slat * slat + cos(DEG2RAD * lat1) * cos(DEG2RAD * lat2) * slon * slon;
1327
1328 return RAD2DEG * 2 * atan2(sqrt(a), sqrt(1 - a));
1329}
1330
1332
1349void sph2cart(Numeric& x,
1350 Numeric& y,
1351 Numeric& z,
1352 const Numeric& r,
1353 const Numeric& lat,
1354 const Numeric& lon) {
1355 ARTS_ASSERT(r > 0);
1356 ARTS_ASSERT(abs(lat) <= 90);
1357 ARTS_ASSERT(abs(lon) <= 360);
1358
1359 const Numeric latrad = DEG2RAD * lat;
1360 const Numeric lonrad = DEG2RAD * lon;
1361
1362 x = r * cos(latrad); // Common term for x and z
1363 y = x * sin(lonrad);
1364 x = x * cos(lonrad);
1365 z = r * sin(latrad);
1366}
1367
1368
1369
1370/*===========================================================================
1371 === Fixes for longitudes
1372 ===========================================================================*/
1373
1375
1391void lon_shiftgrid(Vector& longrid_out,
1392 ConstVectorView longrid_in,
1393 const Numeric lon) {
1394 longrid_out = longrid_in;
1395 if (longrid_in[longrid_in.nelem() - 1] >= lon + 360.)
1396 longrid_out += -360.;
1397 else if (longrid_in[0] <= lon - 360.)
1398 longrid_out += 360.;
1399}
1400
1402/* This function wraps around the given latitude and longitude coordinates
1403 * to the ranges
1404 * - lat in [-90.0, 90.0]
1405 * - lon in [0.0, 360.0]
1406 * Only in the range [-180.0, 180.0] are accepted otherwise an error will
1407 * be thrown
1408 *
1409 * \param[in, out] lat The latitude coordinate.
1410 * \param[in, out] lon The longitude coordiante.
1411 *
1412 * \author Simon Pfreundschuh
1413 * \date 2018-05-22
1414*/
1415void cycle_lat_lon(Numeric& lat, Numeric& lon) {
1416 ARTS_USER_ERROR_IF (lat < -180.0,
1417 "Latitude values < -180.0 are not supported.");
1418 ARTS_USER_ERROR_IF (lat > 180.0,
1419 "Latitude values > 180.0 are not supported.");
1420
1421 if (lat < -90.0) {
1422 lat = -180.0 - lat;
1423 lon += 180.0;
1424 }
1425 if (lat > 90.0) {
1426 lat = 180.0 - lat;
1427 lon += 180.0;
1428 }
1429
1430 while (lon < 0.0) {
1431 lon += 360.0;
1432 }
1433 while (lon > 360.0) {
1434 lon -= 360.0;
1435 }
1436}
1437
1438
1439
1440/*===========================================================================
1441 === Functions involving geodetic latitudes (based on functions in Atmlab)
1442 ===========================================================================*/
1443
1445/*
1446 * \param[out] h Geodetic altitude
1447 * \param[out] la Geodetic latitude
1448 * \param[out] lon Longitude
1449 * \param[in] x x-coordinate (ECEF)
1450 * \param[in] y y-coordinate (ECEF)
1451 * \param[in] z z-coordinate (ECEF)
1452 * \param[in] refellipsoid As the WSV with the same name.
1453 *
1454 * \author Patrick Eriksson
1455 * \date 2020-09-17
1456*/
1457void cart2geodetic(Numeric& h,
1458 Numeric& lat,
1459 Numeric& lon,
1460 const Numeric& x,
1461 const Numeric& y,
1462 const Numeric& z,
1463 const Vector& refellipsoid ) {
1464 // Use geocentric function if geoid is spherical
1465 if (refellipsoid[1] < 1e-7)
1466 { Numeric r;
1467 cart2sph_plain(r, lat, lon, x, y, z);
1468 h = r - refellipsoid[0];
1469 }
1470 else
1471 {
1472 lon = RAD2DEG * atan2(y,x);
1473
1474 const Numeric sq = sqrt(x*x+y*y);
1475 Numeric B0 = atan2(z,sq);
1476 Numeric B = B0-1, N;
1477 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1478
1479 while (abs(B-B0)>1e-10) {
1480 N = refellipsoid[0] / sqrt(1-e2*sin(B0)*sin(B0));
1481 h = sq / cos(B0) - N;
1482 B = B0;
1483 B0 = atan((z/sq) * 1/(1-e2*N/(N+h)));
1484 }
1485 lat = RAD2DEG * B;
1486 }
1487}
1488
1489
1490
1492/*
1493 * \param[out] x x-coordinate (ECEF)
1494 * \param[out] y y-coordinate (ECEF)
1495 * \param[out] z z-coordinate (ECEF)
1496 * \param[in] h Geodetic altitude
1497 * \param[in] lat Geodetic latitude
1498 * \param[in] lon Longitude
1499 * \param[in] refellipsoid As the WSV with the same name.
1500 *
1501 * \author Patrick Eriksson
1502 * \date 2020-09-17
1503*/
1504void geodetic2cart(Numeric& x,
1505 Numeric& y,
1506 Numeric& z,
1507 const Numeric& h,
1508 const Numeric& lat,
1509 const Numeric& lon,
1510 const Vector& refellipsoid ) {
1511 // Use geocentric function if geoid is spherical
1512 if (refellipsoid[1] < 1e-7)
1513 { sph2cart(x, y, z, h+refellipsoid[0], lat, lon); }
1514 else
1515 {
1516 const Numeric a = refellipsoid[0];
1517 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1518 const Numeric sinlat = sin(DEG2RAD*lat);
1519 const Numeric coslat = cos(DEG2RAD*lat);
1520 const Numeric N = a / sqrt(1 - e2*sinlat*sinlat);
1521
1522 x = (N + h) * coslat * cos(DEG2RAD*lon);
1523 y = (N + h) * coslat * sin(DEG2RAD*lon);
1524 z = (N*(1 - e2) + h) * sinlat;
1525 }
1526}
1527
1528
1529
1531
1549void geodeticposlos2cart(Numeric& x,
1550 Numeric& y,
1551 Numeric& z,
1552 Numeric& dx,
1553 Numeric& dy,
1554 Numeric& dz,
1555 const Numeric& h,
1556 const Numeric& lat,
1557 const Numeric& lon,
1558 const Numeric& za,
1559 const Numeric& aa,
1560 const Vector& refellipsoid ) {
1561 ARTS_ASSERT(abs(lat) <= 90);
1562 ARTS_ASSERT(za >= 0 && za <= 180);
1563
1564 // lat = +-90
1565 // For lat = +- 90 the azimuth angle gives the longitude along which the
1566 // LOS goes
1567 // At the poles, no difference between geocentric and geodetic zenith
1568 if (abs(lat) > POLELAT) {
1569 const Numeric s = sign(lat);
1570
1571 x = 0;
1572 y = 0;
1573 z = s * (h + refellipsoid[0]*sqrt(1-refellipsoid[1]*refellipsoid[1]));
1574
1575 dz = s * cos(DEG2RAD * za);
1576 dx = sin(DEG2RAD * za);
1577 dy = dx * sin(DEG2RAD * aa);
1578 dx = dx * cos(DEG2RAD * aa);
1579 }
1580
1581 else {
1582 const Numeric coslat = cos(DEG2RAD * lat);
1583 const Numeric sinlat = sin(DEG2RAD * lat);
1584 const Numeric coslon = cos(DEG2RAD * lon);
1585 const Numeric sinlon = sin(DEG2RAD * lon);
1586
1587 geodetic2cart(x, y, z, h, lat, lon, refellipsoid);
1588
1589 Numeric de, dn, du;
1590 zaaa2enu(de, dn, du, za, aa);
1591
1592 // See
1593 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU
1594 dx = -sinlon*de - sinlat*coslon*dn + coslat*coslon*du;
1595 dy = coslon*de - sinlat*sinlon*dn + coslat*sinlon*du;
1596 dz = coslat* dn + sinlat* du;
1597 }
1598}
1599
1600
1601
1603
1621void cart2geodeticposlos(Numeric& h,
1622 Numeric& lat,
1623 Numeric& lon,
1624 Numeric& za,
1625 Numeric& aa,
1626 const Numeric& x,
1627 const Numeric& y,
1628 const Numeric& z,
1629 const Numeric& dx,
1630 const Numeric& dy,
1631 const Numeric& dz,
1632 const Vector& refellipsoid ) {
1633
1634 cart2geodetic(h, lat, lon, x, y, z, refellipsoid );
1635
1636 const Numeric latrad = DEG2RAD * lat;
1637 const Numeric lonrad = DEG2RAD * lon;
1638 const Numeric coslat = cos(latrad);
1639 const Numeric sinlat = sin(latrad);
1640 const Numeric coslon = cos(lonrad);
1641 const Numeric sinlon = sin(lonrad);
1642
1643 // See
1644 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU
1645 const Numeric de = -sinlon*dx + coslon*dy;
1646 const Numeric dn = -sinlat*coslon*dx - sinlat*sinlon*dy + coslat*dz;
1647 const Numeric du = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
1648
1649 enu2zaaa(za, aa, de, dn, du);
1650}
min
base min(const Array< base > &x)
Min function.
Definition array.h:144
arts_conversions.h
Common ARTS conversions.
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
distance3D
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition geodetic.cc:652
distance2D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition geodetic.cc:165
cart2geodeticposlos
void cart2geodeticposlos(Numeric &h, Numeric &lat, Numeric &lon, Numeric &za, Numeric &aa, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Vector &refellipsoid)
cart2geodeticposlos
Definition geodetic.cc:1621
sph2cart
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition geodetic.cc:1349
los2xyz
void los2xyz(Numeric &za, Numeric &aa, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &x1, const Numeric &y1, const Numeric &z1, const Numeric &x2, const Numeric &y2, const Numeric &z2)
los2xyz
Definition geodetic.cc:1042
cart2poslos_plain
void cart2poslos_plain(Numeric &r, Numeric &lat, Numeric &lon, Numeric &za, Numeric &aa, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Definition geodetic.cc:525
cart2sph
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition geodetic.cc:585
DEG2RAD
constexpr Numeric DEG2RAD
Definition geodetic.cc:25
pol2cart
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition geodetic.cc:314
refell2r
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition geodetic.cc:1248
geodeticposlos2cart
void geodeticposlos2cart(Numeric &x, Numeric &y, Numeric &z, Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &h, const Numeric &lat, const Numeric &lon, const Numeric &za, const Numeric &aa, const Vector &refellipsoid)
geodeticposlos2cart
Definition geodetic.cc:1549
cart2geodetic
void cart2geodetic(Numeric &h, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Vector &refellipsoid)
Conversion from cartesian to geodetic coordinates.
Definition geodetic.cc:1457
cycle_lat_lon
void cycle_lat_lon(Numeric &lat, Numeric &lon)
Cyclic latitude longitude coordinates.
Definition geodetic.cc:1415
pos2refell_r
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition geodetic.cc:1194
line_circle_intersect
void line_circle_intersect(Numeric &x, Numeric &z, const Numeric &xl, const Numeric &zl, const Numeric &dx, const Numeric &dz, const Numeric &xc, const Numeric &zc, const Numeric &r)
geomtanpoint2d
Definition geodetic.cc:261
line_refellipsoid_intersect
void line_refellipsoid_intersect(Numeric &l, const Vector &refellipsoid, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
geomtanpoint
Definition geodetic.cc:857
RAD2DEG
constexpr Numeric RAD2DEG
Definition geodetic.cc:26
sphdist
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition geodetic.cc:1318
latlon_at_aa
void latlon_at_aa(Numeric &lat2, Numeric &lon2, const Numeric &lat1, const Numeric &lon1, const Numeric &aa, const Numeric &ddeg)
latlon_at_aa
Definition geodetic.cc:996
poslos2cart
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition geodetic.cc:338
lon_shiftgrid
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition geodetic.cc:1391
geodetic2cart
void geodetic2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &h, const Numeric &lat, const Numeric &lon, const Vector &refellipsoid)
Conversion from geodetic to cartesian coordinates.
Definition geodetic.cc:1504
refell2d
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition geodetic.cc:1287
cart2sph_plain
void cart2sph_plain(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z)
Definition geodetic.cc:624
cart2poslos
void cart2poslos(Numeric &r, Numeric &lat, Numeric &za, const Numeric &x, const Numeric &z, const Numeric &dx, const Numeric &dz, const Numeric &ppc, const Numeric &lat0, const Numeric &za0)
cart2poslos
Definition geodetic.cc:96
cart2pol
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition geodetic.cc:57
geompath_tanpos_3d
void geompath_tanpos_3d(Numeric &r_tan, Numeric &lat_tan, Numeric &lon_tan, Numeric &l_tan, const Numeric &r, const Numeric &lat, const Numeric &lon, const Numeric &za, const Numeric &aa, const Numeric &ppc)
geompath_tanpos_3d
Definition geodetic.cc:689
line_sphere_intersect
void line_sphere_intersect(Numeric &x, Numeric &y, Numeric &z, const Numeric &xl, const Numeric &yl, const Numeric &zl, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &xc, const Numeric &yc, const Numeric &zc, const Numeric &r)
line_sphere_intersect
Definition geodetic.cc:936
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition interpolation.cc:141
last
Numeric last(ConstVectorView x)
last
Definition math_funcs.cc:142
sign
Numeric sign(const Numeric &x)
sign
Definition math_funcs.cc:445
math_funcs.h
Conversion::deg2rad
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
Definition arts_conversions.h:23
Conversion::rad2deg
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
Definition arts_conversions.h:26
zaaa2enu
void zaaa2enu(Numeric &de, Numeric &dn, Numeric &du, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to ENU unit vector.
Definition ppath.cc:319
enu2zaaa
void enu2zaaa(Numeric &za, Numeric &aa, const Numeric &de, const Numeric &dn, const Numeric &du)
Converts ENU unit vector vector to zenith and azimuth.
Definition ppath.cc:310
geompath_l_at_r
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
Calculates the length from the tangent point for the given radius.
Definition ppath.cc:142
ppath.h
Propagation path structure and functions.
POLELAT
constexpr Numeric POLELAT
Size of north and south poles.
Definition ppath.h:28
ANGTOL
constexpr Numeric ANGTOL
Width of zenith and nadir directions.
Definition ppath.h:39
GridPos
Structure to store a grid position.
Definition interpolation.h:56
GridPos::fd
std::array< Numeric, 2 > fd
Definition interpolation.h:58
GridPos::idx
Index idx
Definition interpolation.h:57
d
#define d
v
#define v
a
#define a
c
#define c
b
#define b