ARTS 2.5.4 (git: bcd8c674)
geodetic.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
19/*===========================================================================
20 === File description
21 ===========================================================================*/
22
32/*===========================================================================
33 === External declarations
34 ===========================================================================*/
35
36#include "geodetic.h"
37#include <cmath>
38#include <stdexcept>
39#include "math_funcs.h"
40#include "ppath.h"
41
42extern const Numeric DEG2RAD;
43extern const Numeric RAD2DEG;
44
45/*===========================================================================
46 === 2D functions
47 ===========================================================================*/
48
49// The 2D case is treated as being the 3D x/z-plane. That is, y-coordinate is
50// skipped. For simplicity, the angle coordinate is denoted as latitude.
51// However, the latitude is here not limited to [-90,90]. It is cyclic and can
52// have any value. The input *lat0* is used to shift the output from atan2 with
53// n*360 to return what should be the expected latitude. That is, it is assumed
54// that no operation moves the latitude more than 180 degrees from the initial
55// value *lat0*. Negative zeniath angles are handled, following ARTS definition
56// of 2D geometry.
57
59
75 Numeric& lat,
76 const Numeric& x,
77 const Numeric& z,
78 const Numeric& lat0,
79 const Numeric& za0) {
80 r = sqrt(x * x + z * z);
81
82 // Zenith and nadir cases
83 const Numeric absza = abs(za0);
84 if (absza < ANGTOL || absza > 180 - ANGTOL) {
85 lat = lat0;
86 }
87
88 else { // Latitude inside [0,360]
89 lat = RAD2DEG * atan2(z, x);
90 // Shift with n*360 to get as close as possible to lat0
91 lat = lat - 360.0 * Numeric(round((lat - lat0) / 360.0));
92 }
93}
94
96
114 Numeric& lat,
115 Numeric& za,
116 const Numeric& x,
117 const Numeric& z,
118 const Numeric& dx,
119 const Numeric& dz,
120 const Numeric& ppc,
121 const Numeric& lat0,
122 const Numeric& za0) {
123 r = sqrt(x * x + z * z);
124
125 // Zenith and nadir cases
126 const Numeric absza = abs(za0);
127 if (absza < ANGTOL || absza > 180 - ANGTOL) {
128 lat = lat0;
129 za = za0;
130 }
131
132 else {
133 lat = RAD2DEG * atan2(z, x);
134
135 const Numeric latrad = DEG2RAD * lat;
136 const Numeric coslat = cos(latrad);
137 const Numeric sinlat = sin(latrad);
138 const Numeric dr = coslat * dx + sinlat * dz;
139
140 // Use ppc for max accuracy, but dr required to resolve if up-
141 // and downward cases.
142
143 // Another possibility to obtain (absolute value of) za is
144 // RAD2DEG*acos(dr).
145 // It is checked that the two ways give consistent results, but
146 // occasionally deviate with 1e-4 deg (due to numerical issues).
147
148 za = RAD2DEG * asin(ppc / r);
149 if (za0 > 0) {
150 if (std::isnan(za)) {
151 za = 90;
152 } else if (dr < 0) {
153 za = 180.0 - za;
154 }
155 } else {
156 if (std::isnan(za)) {
157 za = -90;
158 } else if (dr < 0) {
159 za = -180.0 + za;
160 } else {
161 za = -za;
162 }
163 }
164 }
165}
166
168
183 const Numeric& r1,
184 const Numeric& lat1,
185 const Numeric& r2,
186 const Numeric& lat2) {
187 ARTS_ASSERT(abs(lat2 - lat1) <= 180);
188
189 Numeric x1, z1, x2, z2;
190 pol2cart(x1, z1, r1, lat1);
191 pol2cart(x2, z2, r2, lat2);
192
193 const Numeric dx = x2 - x1;
194 const Numeric dz = z2 - z1;
195 l = sqrt(dx * dx + dz * dz);
196}
197
199
229/*
230void geomtanpoint2d(
231 Numeric& r_tan,
232 Numeric& lat_tan,
233 ConstVectorView refellipsoid,
234 const Numeric& r,
235 const Numeric& lat,
236 const Numeric& za )
237{
238 ARTS_ASSERT( refellipsoid.nelem() == 2 );
239 ARTS_ASSERT( refellipsoid[0] > 0 );
240 ARTS_ASSERT( refellipsoid[1] >= 0 );
241 ARTS_ASSERT( refellipsoid[1] < 1 );
242 ARTS_ASSERT( r > 0 );
243 ARTS_ASSERT( za >= 90 );
244 // e=1e-7 corresponds to that polar radius
245 if( refellipsoid[1] < 1e-7 ) // less than 1 um smaller than equatorial
246 { // one for the Earth
247 r_tan = geometrical_ppc( r, za );
248 if( za > 0 )
249 { lat_tan = geompath_lat_at_za( za, lat, 90 ); }
250 else
251 { lat_tan = geompath_lat_at_za( za, lat, -90 ); }
252 }
253
254 else
255 {
256 ARTS_ASSERT( 0 ); // To be implemented
257 }
258}
259*/
260
262
279 Numeric& z,
280 const Numeric& xl,
281 const Numeric& zl,
282 const Numeric& dx,
283 const Numeric& dz,
284 const Numeric& xc,
285 const Numeric& zc,
286 const Numeric& r) {
287 const Numeric a = dx * dx + dz * dz;
288 const Numeric b = 2 * (dx * (xl - xc) + dz * (zl - zc));
289 const Numeric c =
290 xc * xc + zc * zc + xl * xl + zl * zl - 2 * (xc * xl + zc * zl) - r * r;
291
292 Numeric d = b * b - 4 * a * c;
293 ARTS_ASSERT(d > 0);
294
295 const Numeric a2 = 2 * a;
296 const Numeric b2 = -b / a2;
297 const Numeric e = sqrt(d) / a2;
298
299 const Numeric l1 = b2 + e;
300 const Numeric l2 = b2 - e;
301
302 Numeric l;
303 if (l1 < 0) {
304 l = l2;
305 } else if (l2 < 0) {
306 l = l1;
307 } else {
308 l = min(l1, l2);
309 ARTS_ASSERT(l >= 0);
310 }
311
312 x = xl + l * dx;
313 z = zl + l * dz;
314}
315
317
331void pol2cart(Numeric& x, Numeric& z, const Numeric& r, const Numeric& lat) {
332 ARTS_ASSERT(r > 0);
333
334 const Numeric latrad = DEG2RAD * lat;
335
336 x = r * cos(latrad);
337 z = r * sin(latrad);
338}
339
341
356 Numeric& z,
357 Numeric& dx,
358 Numeric& dz,
359 const Numeric& r,
360 const Numeric& lat,
361 const Numeric& za) {
362 ARTS_ASSERT(r > 0);
363 ARTS_ASSERT(za >= -180 && za <= 180);
364
365 const Numeric latrad = DEG2RAD * lat;
366 const Numeric zarad = DEG2RAD * za;
367
368 const Numeric coslat = cos(latrad);
369 const Numeric sinlat = sin(latrad);
370 const Numeric cosza = cos(zarad);
371 const Numeric sinza = sin(zarad);
372
373 // This part as pol2cart but uses local variables
374 x = r * coslat;
375 z = r * sinlat;
376
377 const Numeric dr = cosza;
378 const Numeric dlat = sinza; // r-term cancel out below
379
380 dx = coslat * dr - sinlat * dlat;
381 dz = sinlat * dr + coslat * dlat;
382}
383
384/*===========================================================================
385 === 3D functions
386 ===========================================================================*/
387
389
422 Numeric& lat,
423 Numeric& lon,
424 Numeric& za,
425 Numeric& aa,
426 const Numeric& x,
427 const Numeric& y,
428 const Numeric& z,
429 const Numeric& dx,
430 const Numeric& dy,
431 const Numeric& dz,
432 const Numeric& ppc,
433 const Numeric& x0,
434 const Numeric& y0,
435 const Numeric& z0,
436 const Numeric& lat0,
437 const Numeric& lon0,
438 const Numeric& za0,
439 const Numeric& aa0) {
440 // Radius of new point
441 r = sqrt(x * x + y * y + z * z);
442
443 // Zenith and nadir cases
444 if (za0 < ANGTOL || za0 > 180 - ANGTOL) {
445 lat = lat0;
446 lon = lon0;
447 za = za0;
448 aa = aa0;
449 }
450
451 else {
452 lat = RAD2DEG * asin(z / r);
453 lon = RAD2DEG * atan2(y, x);
454
455 bool ns_case = false;
456 bool lon_flip = false;
457
458 // Make sure that lon is maintained for N-S cases (if not
459 // starting on a pole)
460 if ((abs(aa0) < ANGTOL || abs(180 - aa0) < ANGTOL) &&
461 abs(lat0) <= POLELAT) {
462 ns_case = true;
463 // Check that not lon changed with 180 deg
464 if (abs(abs(lon - lon0) - 180) < 5) {
465 lon_flip = true;
466 if (lon0 > 0) {
467 lon = lon0 - 180;
468 } else {
469 lon = lon0 + 180;
470 }
471 } else {
472 lon = lon0;
473 }
474 }
475
476 const Numeric latrad = DEG2RAD * lat;
477 const Numeric lonrad = DEG2RAD * lon;
478 const Numeric coslat = cos(latrad);
479 const Numeric sinlat = sin(latrad);
480 const Numeric coslon = cos(lonrad);
481 const Numeric sinlon = sin(lonrad);
482
483 // Set za by ppc for max accuracy, but this does not resolve
484 // za and 180-za. This was first resolved by dr, but using l and lmax was
485 // found to be more stable.
486 za = RAD2DEG * asin(ppc / r);
487
488 // Correct and check za
489 if (std::isnan(za)) {
490 za = 90;
491 }
492 // If za0 > 90, then correct za could be 180-za. Resolved by checking if
493 // the tangent point is passed or not
494 if (za0 > 90) {
495 const Numeric l =
496 sqrt(pow(x - x0, 2.0) + pow(y - y0, 2.0) + pow(z - z0, 2.0));
497 const Numeric r0 = sqrt(x0 * x0 + y0 * y0 + z0 * z0);
498 const Numeric ltan = geompath_l_at_r(ppc, r0);
499 if (l < ltan) {
500 za = 180.0 - za;
501 }
502 }
503
504 // For lat = +- 90 the azimuth angle gives the longitude along which
505 // the LOS goes
506 if (abs(lat) >= POLELAT) {
507 aa = RAD2DEG * atan2(dy, dx);
508 }
509
510 // N-S cases, not starting at a pole
511 else if (ns_case) {
512 if (!lon_flip) {
513 aa = aa0;
514 } else {
515 if (abs(aa0) < ANGTOL) {
516 aa = 180;
517 } else {
518 aa = 0;
519 }
520 }
521 }
522
523 else {
524 const Numeric dlat = -sinlat * coslon / r * dx -
525 sinlat * sinlon / r * dy + coslat / r * dz;
526 const Numeric dlon = -sinlon / coslat / r * dx + coslon / coslat / r * dy;
527
528 aa = RAD2DEG * acos(r * dlat / sin(DEG2RAD * za));
529
530 if (std::isnan(aa)) {
531 if (dlat >= 0) {
532 aa = 0;
533 } else {
534 aa = 180;
535 }
536 } else if (dlon < 0) {
537 aa = -aa;
538 }
539 }
540 }
541}
543 Numeric& lat,
544 Numeric& lon,
545 Numeric& za,
546 Numeric& aa,
547 const Numeric& x,
548 const Numeric& y,
549 const Numeric& z,
550 const Numeric& dx,
551 const Numeric& dy,
552 const Numeric& dz) {
553 cart2sph_plain(r, lat, lon, x, y, z);
554
555
556 const Numeric latrad = DEG2RAD * lat;
557 const Numeric lonrad = DEG2RAD * lon;
558 const Numeric coslat = cos(latrad);
559 const Numeric sinlat = sin(latrad);
560 const Numeric coslon = cos(lonrad);
561 const Numeric sinlon = sin(lonrad);
562
563 const Numeric dr = coslat*coslon*dx + sinlat*dz + coslat*sinlon*dy;
564 const Numeric dlat = -sinlat*coslon/r*dx + coslat/r*dz - sinlat*sinlon/r*dy;
565 const Numeric dlon = -sinlon/coslat/r*dx + coslon/coslat/r*dy;
566
567 za = acos( dr );
568 aa = RAD2DEG * acos( r * dlat / sin( za ) );
569 za *= RAD2DEG;
570
571 // Corrections of aa
572 if (std::isnan(aa)) {
573 if (dlat >= 0)
574 aa = 0;
575 else
576 aa = 180;
577 } else if (dlon < 0) {
578 aa = -aa;
579 }
580}
581
583
603 Numeric& lat,
604 Numeric& lon,
605 const Numeric& x,
606 const Numeric& y,
607 const Numeric& z,
608 const Numeric& lat0,
609 const Numeric& lon0,
610 const Numeric& za0,
611 const Numeric& aa0) {
612 r = sqrt(x * x + y * y + z * z);
613
614 // Zenith and nadir cases
615 if (za0 < ANGTOL || za0 > 180 - ANGTOL) {
616 lat = lat0;
617 lon = lon0;
618 }
619
620 else {
621 lat = RAD2DEG * asin(z / r);
622 lon = RAD2DEG * atan2(y, x);
623
624 // Make sure that lon is maintained for N-S cases (if not
625 // starting on a pole)
626 if ((abs(aa0) < ANGTOL || abs(180 - aa0) < ANGTOL) &&
627 abs(lat0) <= POLELAT) {
628 // Check that not lon changed with 180 deg
629 if (abs(lon - lon0) < 1) {
630 lon = lon0;
631 } else {
632 if (lon0 > 0) {
633 lon = lon0 - 180;
634 } else {
635 lon = lon0 + 180;
636 }
637 }
638 }
639 }
640}
642 Numeric& lat,
643 Numeric& lon,
644 const Numeric& x,
645 const Numeric& y,
646 const Numeric& z) {
647 r = sqrt(x * x + y * y + z * z);
648 lat = RAD2DEG * asin(z / r);
649 lon = RAD2DEG * atan2(y, x);
650}
651
652
653
655
670 const Numeric& r1,
671 const Numeric& lat1,
672 const Numeric& lon1,
673 const Numeric& r2,
674 const Numeric& lat2,
675 const Numeric& lon2) {
676 Numeric x1, y1, z1, x2, y2, z2;
677 sph2cart(x1, y1, z1, r1, lat1, lon1);
678 sph2cart(x2, y2, z2, r2, lat2, lon2);
679
680 const Numeric dx = x2 - x1;
681 const Numeric dy = y2 - y1;
682 const Numeric dz = z2 - z1;
683 l = sqrt(dx * dx + dy * dy + dz * dz);
684}
685
687
707 Numeric& lat_tan,
708 Numeric& lon_tan,
709 Numeric& l_tan,
710 const Numeric& r,
711 const Numeric& lat,
712 const Numeric& lon,
713 const Numeric& za,
714 const Numeric& aa,
715 const Numeric& ppc) {
716 ARTS_ASSERT(za >= 90);
717 ARTS_ASSERT(r >= ppc);
718
719 Numeric x, y, z, dx, dy, dz;
720
721 poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
722
723 l_tan = sqrt(r * r - ppc * ppc);
724
725 cart2sph(r_tan,
726 lat_tan,
727 lon_tan,
728 x + dx * l_tan,
729 y + dy * l_tan,
730 z + dz * l_tan,
731 lat,
732 lon,
733 za,
734 aa);
735}
736
738
768/*
769void geomtanpoint(
770 Numeric& r_tan,
771 Numeric& lat_tan,
772 Numeric& lon_tan,
773 ConstVectorView refellipsoid,
774 const Numeric& r,
775 const Numeric& lat,
776 const Numeric& lon,
777 const Numeric& za,
778 const Numeric& aa )
779{
780 ARTS_ASSERT( refellipsoid.nelem() == 2 );
781 ARTS_ASSERT( refellipsoid[0] > 0 );
782 ARTS_ASSERT( refellipsoid[1] >= 0 );
783 ARTS_ASSERT( refellipsoid[1] < 1 );
784 ARTS_ASSERT( r > 0 );
785 ARTS_ASSERT( za >= 90 );
786
787 if( refellipsoid[1] < 1e-7 ) // e=1e-7 corresponds to that polar radius
788 { // less than 1 um smaller than equatorial
789 Numeric x, y, z, dx, dy, dz; // one for the Earth
790
791 poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
792
793 const Numeric ppc = r * sin( DEG2RAD * abs(za) );
794 const Numeric l_tan = sqrt( r*r - ppc*ppc );
795
796 cart2sph( r_tan, lat_tan, lon_tan, x+dx*l_tan, y+dy*l_tan, z+dz*l_tan );
797 }
798
799 else
800 {
801 // Equatorial and polar radii squared:
802 const Numeric a2 = refellipsoid[0]*refellipsoid[0];
803 const Numeric b2 = a2 * ( 1 - refellipsoid[1]*refellipsoid[1] );
804
805 Vector X(3), xunit(3), yunit(3), zunit(3);
806
807 poslos2cart( X[0], X[1], X[2], xunit[0], xunit[1], xunit[2],
808 r, lat, lon, za, aa );
809 cross( zunit, xunit, X );
810 unitl( zunit ); // Normalisation of length to 1
811
812 cross( yunit, zunit, xunit );
813 unitl( yunit ); // Normalisation of length to 1
814
815 Numeric x = X[0];
816 Numeric y = X[1];
817 const Numeric w11 = xunit[0];
818 const Numeric w12 = yunit[0];
819 const Numeric w21 = xunit[1];
820 const Numeric w22 = yunit[1];
821 const Numeric w31 = xunit[2];
822 const Numeric w32 = yunit[2];
823
824 const Numeric yr = X * yunit;
825 const Numeric xr = X * xunit;
826
827 const Numeric A = (w11*w11 + w21*w21)/a2 + w31*w31/b2;
828 const Numeric B = 2.0*((w11*w12 + w21*w22)/a2 + (w31*w32)/b2);
829 const Numeric C = (w12*w12 + w22*w22)/a2 + w32*w32/b2;
830
831 if( B == 0.0 )
832 { x = 0.0; }
833 else
834 {
835 const Numeric K = -2.0*A/B;
836 const Numeric factor = 1.0/(A+(B+C*K)*K);
837 x = sqrt(factor);
838 y = K*x;
839 }
840
841 const Numeric dist1 = (xr-X[0])*(xr-X[0]) + (yr-y)*(yr-y);
842 const Numeric dist2 = (xr+X[0])*(xr+X[0]) + (yr+y)*(yr+y);
843
844 if( dist1 > dist2 )
845 { x = -x; }
846
847 cart2sph( r_tan, lat_tan, lon_tan, w11*x + w12*yr, w21*x + w22*yr,
848 w31*x + w32*yr );
849 }
850}
851*/
852
854
875 const Vector& refellipsoid,
876 const Numeric& x,
877 const Numeric& y,
878 const Numeric& z,
879 const Numeric& dx,
880 const Numeric& dy,
881 const Numeric& dz) {
882 // Code taken from Atmlab's ellipsoid_intersection
883
884 // Spherical case
885 if (refellipsoid[1] < 1e-7) {
886 const Numeric p = x*dx + y*dy + z*dz;
887 const Numeric pp = p*p;
888 const Numeric q = x*x + y*y + z*z - refellipsoid[0]*refellipsoid[0];
889 if (q>pp)
890 l = -1.0;
891 else {
892 const Numeric sq = sqrt(pp - q);
893 if (-p > sq)
894 l = -p - sq;
895 else
896 l = -p + sq;
897 }
898 }
899
900 // Ellipsoid case
901 else {
902 // Based on https://medium.com/@stephenhartzell/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6
903 const Numeric a = refellipsoid[0];
904 const Numeric b = refellipsoid[0];
905 const Numeric c = refellipsoid[0] * sqrt(1-refellipsoid[1]*refellipsoid[1]);
906 const Numeric a2 = a*a;
907 const Numeric b2 = b*b;
908 const Numeric c2 = c*c;
909 const Numeric x2 = x*x;
910 const Numeric y2 = y*y;
911 const Numeric z2 = z*z;
912 const Numeric dx2 = dx*dx;
913 const Numeric dy2 = dy*dy;
914 const Numeric dz2 = dz*dz;
915 const Numeric rad = a2*b2*dz2 + a2*c2*dy2 - a2*dy2*z2 + 2*a2*dy*dz*y*z -
916 a2*dz2*y2 + b2*c2*dx2 - b2*dx2*z2 + 2*b2*dx*dz*x*z -
917 b2*dz2*x2 - c2*dx2*y2 + 2*c2*dx*dy*x*y - c2*dy2*x2;
918 if (rad<0)
919 l = -1.0;
920 else {
921 const Numeric val = -a2*b2*dz*z - a2*c2*dy*y - b2*c2*dx*x;
922 const Numeric mag = a2*b2*dz2 + a2*c2*dy2 + b2*c2*dx2;
923 const Numeric abc = a*b*c*sqrt(rad);
924 if (val > abc)
925 l = (val - abc) / mag;
926 else
927 l = (val + abc) / mag;
928 }
929 }
930}
931
933
954 Numeric& y,
955 Numeric& z,
956 const Numeric& xl,
957 const Numeric& yl,
958 const Numeric& zl,
959 const Numeric& dx,
960 const Numeric& dy,
961 const Numeric& dz,
962 const Numeric& xc,
963 const Numeric& yc,
964 const Numeric& zc,
965 const Numeric& r) {
966 const Numeric a = dx * dx + dy * dy + dz * dz;
967 const Numeric b = 2 * (dx * (xl - xc) + dy * (yl - yc) + dz * (zl - zc));
968 const Numeric c = xc * xc + yc * yc + zc * zc + xl * xl + yl * yl + zl * zl -
969 2 * (xc * xl + yc * yl + zc * zl) - r * r;
970
971 Numeric d = b * b - 4 * a * c;
972 ARTS_ASSERT(d > 0);
973
974 const Numeric a2 = 2 * a;
975 const Numeric b2 = -b / a2;
976 const Numeric e = sqrt(d) / a2;
977
978 const Numeric l1 = b2 + e;
979 const Numeric l2 = b2 - e;
980
981 Numeric l;
982 if (l1 < 0) {
983 l = l2;
984 } else if (l2 < 0) {
985 l = l1;
986 } else {
987 l = min(l1, l2);
988 ARTS_ASSERT(l >= 0);
989 }
990
991 x = xl + l * dx;
992 y = yl + l * dy;
993 z = zl + l * dz;
994}
995
997
1014 Numeric& lon2,
1015 const Numeric& lat1,
1016 const Numeric& lon1,
1017 const Numeric& aa,
1018 const Numeric& ddeg) {
1019 // Code from http://www.movable-type.co.uk/scripts/latlong.html
1020 // (but with short-cuts, such as asin(sin(lat2)) = lat2)
1021 // Note that lat1 here is another latitude
1022
1023 const Numeric dang = DEG2RAD * ddeg;
1024 const Numeric cosdang = cos(dang);
1025 const Numeric sindang = sin(dang);
1026 const Numeric latrad = DEG2RAD * lat1;
1027 const Numeric coslat = cos(latrad);
1028 const Numeric sinlat = sin(latrad);
1029 const Numeric aarad = DEG2RAD * aa;
1030
1031 lat2 = sinlat * cosdang + coslat * sindang * cos(aarad);
1032 lon2 = lon1 + RAD2DEG * (atan2(sin(aarad) * sindang * coslat,
1033 cosdang - sinlat * lat2));
1034 lat2 = RAD2DEG * asin(lat2);
1035}
1036
1038
1060 Numeric& aa,
1061 const Numeric& r1,
1062 const Numeric& lat1,
1063 const Numeric& lon1,
1064 const Numeric& x1,
1065 const Numeric& y1,
1066 const Numeric& z1,
1067 const Numeric& x2,
1068 const Numeric& y2,
1069 const Numeric& z2) {
1070 Numeric dx = x2 - x1, dy = y2 - y1, dz = z2 - z1;
1071 const Numeric ldxyz = sqrt(dx * dx + dy * dy + dz * dz);
1072 dx /= ldxyz;
1073 dy /= ldxyz;
1074 dz /= ldxyz;
1075
1076 // All below extracted from 3D version of cart2poslos:
1077 const Numeric latrad = DEG2RAD * lat1;
1078 const Numeric lonrad = DEG2RAD * lon1;
1079 const Numeric coslat = cos(latrad);
1080 const Numeric sinlat = sin(latrad);
1081 const Numeric coslon = cos(lonrad);
1082 const Numeric sinlon = sin(lonrad);
1083
1084 const Numeric dr = coslat * coslon * dx + coslat * sinlon * dy + sinlat * dz;
1085 const Numeric dlat =
1086 -sinlat * coslon / r1 * dx - sinlat * sinlon / r1 * dy + coslat / r1 * dz;
1087 const Numeric dlon = -sinlon / coslat / r1 * dx + coslon / coslat / r1 * dy;
1088
1089 za = RAD2DEG * acos(dr);
1090 aa = RAD2DEG * acos(r1 * dlat / sin(DEG2RAD * za));
1091 if (std::isnan(aa)) {
1092 if (dlat >= 0) {
1093 aa = 0;
1094 } else {
1095 aa = 180;
1096 }
1097 } else if (dlon < 0) {
1098 aa = -aa;
1099 }
1100}
1101
1103
1128 Numeric& y,
1129 Numeric& z,
1130 Numeric& dx,
1131 Numeric& dy,
1132 Numeric& dz,
1133 const Numeric& r,
1134 const Numeric& lat,
1135 const Numeric& lon,
1136 const Numeric& za,
1137 const Numeric& aa) {
1138 ARTS_ASSERT(r > 0);
1139 ARTS_ASSERT(abs(lat) <= 90);
1140 //ARTS_ASSERT( abs( lon ) <= 360 );
1141 ARTS_ASSERT(za >= 0 && za <= 180);
1142
1143 // lat = +-90
1144 // For lat = +- 90 the azimuth angle gives the longitude along which the
1145 // LOS goes
1146 if (abs(lat) > POLELAT) {
1147 const Numeric s = sign(lat);
1148
1149 x = 0;
1150 y = 0;
1151 z = s * r;
1152
1153 dz = s * cos(DEG2RAD * za);
1154 dx = sin(DEG2RAD * za);
1155 dy = dx * sin(DEG2RAD * aa);
1156 dx = dx * cos(DEG2RAD * aa);
1157 }
1158
1159 else {
1160 const Numeric latrad = DEG2RAD * lat;
1161 const Numeric lonrad = DEG2RAD * lon;
1162 const Numeric zarad = DEG2RAD * za;
1163 const Numeric aarad = DEG2RAD * aa;
1164
1165 const Numeric coslat = cos(latrad);
1166 const Numeric sinlat = sin(latrad);
1167 const Numeric coslon = cos(lonrad);
1168 const Numeric sinlon = sin(lonrad);
1169 const Numeric cosza = cos(zarad);
1170 const Numeric sinza = sin(zarad);
1171 const Numeric cosaa = cos(aarad);
1172 const Numeric sinaa = sin(aarad);
1173
1174 // This part as sph2cart but uses local variables
1175 x = r * coslat; // Common term for x and y
1176 y = x * sinlon;
1177 x = x * coslon;
1178 z = r * sinlat;
1179
1180 const Numeric dr = cosza;
1181 const Numeric dlat = sinza * cosaa; // r-term cancel out below
1182 const Numeric dlon = sinza * sinaa / coslat;
1183
1184 dx = coslat * coslon * dr - sinlat * coslon * dlat - coslat * sinlon * dlon;
1185 dz = sinlat * dr + coslat * dlat;
1186 dy = coslat * sinlon * dr - sinlat * sinlon * dlat + coslat * coslon * dlon;
1187 }
1188}
1189
1191
1211Numeric pos2refell_r(const Index& atmosphere_dim,
1212 ConstVectorView refellipsoid,
1213 ConstVectorView lat_grid,
1214 ConstVectorView lon_grid,
1215 ConstVectorView rte_pos) {
1216 if (atmosphere_dim == 1) {
1217 return refellipsoid[0];
1218 } else {
1219 ARTS_ASSERT(rte_pos.nelem() > 1);
1220
1221 bool inside = true;
1222
1223 if (rte_pos[1] < lat_grid[0] || rte_pos[1] > last(lat_grid)) {
1224 inside = false;
1225 } else if (atmosphere_dim == 3) {
1226 ARTS_ASSERT(rte_pos.nelem() == 3);
1227 if (rte_pos[2] < lon_grid[0] || rte_pos[2] > last(lon_grid)) {
1228 inside = false;
1229 }
1230 }
1231
1232 if (inside) {
1233 GridPos gp_lat;
1234 gridpos(gp_lat, lat_grid, rte_pos[1]);
1235 return refell2d(refellipsoid, lat_grid, gp_lat);
1236 } else {
1237 return refell2r(refellipsoid, rte_pos[1]);
1238 }
1239 }
1240}
1241
1243
1265Numeric refell2r(ConstVectorView refellipsoid, const Numeric& lat) {
1266 ARTS_ASSERT(refellipsoid.nelem() == 2);
1267 ARTS_ASSERT(refellipsoid[0] > 0);
1268 ARTS_ASSERT(refellipsoid[1] >= 0);
1269 ARTS_ASSERT(refellipsoid[1] < 1);
1270
1271 if (refellipsoid[1] < 1e-7) // e=1e-7 corresponds to that polar radius
1272 { // less than 1 um smaller than equatorial
1273 return refellipsoid[0]; // one for the Earth
1274 }
1275
1276 else {
1277 const Numeric c = 1 - refellipsoid[1] * refellipsoid[1];
1278 const Numeric b = refellipsoid[0] * sqrt(c);
1279 const Numeric v = DEG2RAD * lat;
1280 const Numeric ct = cos(v);
1281 const Numeric st = sin(v);
1282
1283 return b / sqrt(c * ct * ct + st * st);
1284 }
1285}
1286
1288
1305 ConstVectorView lat_grid,
1306 const GridPos gp) {
1307 if (gp.fd[0] == 0)
1308 return refell2r(refellipsoid, lat_grid[gp.idx]);
1309 else if (gp.fd[0] == 1)
1310 return refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1311 else
1312 return gp.fd[1] * refell2r(refellipsoid, lat_grid[gp.idx]) +
1313 gp.fd[0] * refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1314}
1315
1317
1336 const Numeric& lon1,
1337 const Numeric& lat2,
1338 const Numeric& lon2) {
1339 // Equations taken from http://www.movable-type.co.uk/scripts/latlong.html
1340 const Numeric slat = sin(DEG2RAD * (lat2 - lat1) / 2.0);
1341 const Numeric slon = sin(DEG2RAD * (lon2 - lon1) / 2.0);
1342 const Numeric a =
1343 slat * slat + cos(DEG2RAD * lat1) * cos(DEG2RAD * lat2) * slon * slon;
1344
1345 return RAD2DEG * 2 * atan2(sqrt(a), sqrt(1 - a));
1346}
1347
1349
1367 Numeric& y,
1368 Numeric& z,
1369 const Numeric& r,
1370 const Numeric& lat,
1371 const Numeric& lon) {
1372 ARTS_ASSERT(r > 0);
1373 ARTS_ASSERT(abs(lat) <= 90);
1374 ARTS_ASSERT(abs(lon) <= 360);
1375
1376 const Numeric latrad = DEG2RAD * lat;
1377 const Numeric lonrad = DEG2RAD * lon;
1378
1379 x = r * cos(latrad); // Common term for x and z
1380 y = x * sin(lonrad);
1381 x = x * cos(lonrad);
1382 z = r * sin(latrad);
1383}
1384
1385
1386
1387/*===========================================================================
1388 === Fixes for longitudes
1389 ===========================================================================*/
1390
1392
1408void lon_shiftgrid(Vector& longrid_out,
1409 ConstVectorView longrid_in,
1410 const Numeric lon) {
1411 longrid_out = longrid_in;
1412 if (longrid_in[longrid_in.nelem() - 1] >= lon + 360.)
1413 longrid_out += -360.;
1414 else if (longrid_in[0] <= lon - 360.)
1415 longrid_out += 360.;
1416}
1417
1419/* This function wraps around the given latitude and longitude coordinates
1420 * to the ranges
1421 * - lat in [-90.0, 90.0]
1422 * - lon in [0.0, 360.0]
1423 * Only in the range [-180.0, 180.0] are accepted otherwise an error will
1424 * be thrown
1425 *
1426 * \param[in, out] lat The latitude coordinate.
1427 * \param[in, out] lon The longitude coordiante.
1428 *
1429 * \author Simon Pfreundschuh
1430 * \date 2018-05-22
1431*/
1433 ARTS_USER_ERROR_IF (lat < -180.0,
1434 "Latitude values < -180.0 are not supported.");
1435 ARTS_USER_ERROR_IF (lat > 180.0,
1436 "Latitude values > 180.0 are not supported.");
1437
1438 if (lat < -90.0) {
1439 lat = -180.0 - lat;
1440 lon += 180.0;
1441 }
1442 if (lat > 90.0) {
1443 lat = 180.0 - lat;
1444 lon += 180.0;
1445 }
1446
1447 while (lon < 0.0) {
1448 lon += 360.0;
1449 }
1450 while (lon > 360.0) {
1451 lon -= 360.0;
1452 }
1453}
1454
1455
1456
1457/*===========================================================================
1458 === Functions involving geodetic latitudes (based on functions in Atmlab)
1459 ===========================================================================*/
1460
1462/*
1463 * \param[out] h Geodetic altitude
1464 * \param[out] la Geodetic latitude
1465 * \param[out] lon Longitude
1466 * \param[in] x x-coordinate (ECEF)
1467 * \param[in] y y-coordinate (ECEF)
1468 * \param[in] z z-coordinate (ECEF)
1469 * \param[in] refellipsoid As the WSV with the same name.
1470 *
1471 * \author Patrick Eriksson
1472 * \date 2020-09-17
1473*/
1475 Numeric& lat,
1476 Numeric& lon,
1477 const Numeric& x,
1478 const Numeric& y,
1479 const Numeric& z,
1480 const Vector& refellipsoid ) {
1481 // Use geocentric function if geoid is spherical
1482 if (refellipsoid[1] < 1e-7)
1483 { Numeric r;
1484 cart2sph_plain(r, lat, lon, x, y, z);
1485 h = r - refellipsoid[0];
1486 }
1487 else
1488 {
1489 lon = RAD2DEG * atan2(y,x);
1490
1491 const Numeric sq = sqrt(x*x+y*y);
1492 Numeric B0 = atan2(z,sq);
1493 Numeric B = B0-1, N;
1494 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1495
1496 while (abs(B-B0)>1e-10) {
1497 N = refellipsoid[0] / sqrt(1-e2*sin(B0)*sin(B0));
1498 h = sq / cos(B0) - N;
1499 B = B0;
1500 B0 = atan((z/sq) * 1/(1-e2*N/(N+h)));
1501 }
1502 lat = RAD2DEG * B;
1503 }
1504}
1505
1506
1507
1509/*
1510 * \param[out] x x-coordinate (ECEF)
1511 * \param[out] y y-coordinate (ECEF)
1512 * \param[out] z z-coordinate (ECEF)
1513 * \param[in] h Geodetic altitude
1514 * \param[in] lat Geodetic latitude
1515 * \param[in] lon Longitude
1516 * \param[in] refellipsoid As the WSV with the same name.
1517 *
1518 * \author Patrick Eriksson
1519 * \date 2020-09-17
1520*/
1522 Numeric& y,
1523 Numeric& z,
1524 const Numeric& h,
1525 const Numeric& lat,
1526 const Numeric& lon,
1527 const Vector& refellipsoid ) {
1528 // Use geocentric function if geoid is spherical
1529 if (refellipsoid[1] < 1e-7)
1530 { sph2cart(x, y, z, h+refellipsoid[0], lat, lon); }
1531 else
1532 {
1533 const Numeric a = refellipsoid[0];
1534 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1535 const Numeric sinlat = sin(DEG2RAD*lat);
1536 const Numeric coslat = cos(DEG2RAD*lat);
1537 const Numeric N = a / sqrt(1 - e2*sinlat*sinlat);
1538
1539 x = (N + h) * coslat * cos(DEG2RAD*lon);
1540 y = (N + h) * coslat * sin(DEG2RAD*lon);
1541 z = (N*(1 - e2) + h) * sinlat;
1542 }
1543}
1544
1545
1546
1548
1567 Numeric& y,
1568 Numeric& z,
1569 Numeric& dx,
1570 Numeric& dy,
1571 Numeric& dz,
1572 const Numeric& h,
1573 const Numeric& lat,
1574 const Numeric& lon,
1575 const Numeric& za,
1576 const Numeric& aa,
1577 const Vector& refellipsoid ) {
1578 ARTS_ASSERT(abs(lat) <= 90);
1579 ARTS_ASSERT(za >= 0 && za <= 180);
1580
1581 // lat = +-90
1582 // For lat = +- 90 the azimuth angle gives the longitude along which the
1583 // LOS goes
1584 // At the poles, no difference between geocentric and geodetic zenith
1585 if (abs(lat) > POLELAT) {
1586 const Numeric s = sign(lat);
1587
1588 x = 0;
1589 y = 0;
1590 z = s * (h + refellipsoid[0]*sqrt(1-refellipsoid[1]*refellipsoid[1]));
1591
1592 dz = s * cos(DEG2RAD * za);
1593 dx = sin(DEG2RAD * za);
1594 dy = dx * sin(DEG2RAD * aa);
1595 dx = dx * cos(DEG2RAD * aa);
1596 }
1597
1598 else {
1599 const Numeric coslat = cos(DEG2RAD * lat);
1600 const Numeric sinlat = sin(DEG2RAD * lat);
1601 const Numeric coslon = cos(DEG2RAD * lon);
1602 const Numeric sinlon = sin(DEG2RAD * lon);
1603
1604 geodetic2cart(x, y, z, h, lat, lon, refellipsoid);
1605
1606 Numeric de, dn, du;
1607 zaaa2enu(de, dn, du, za, aa);
1608
1609 // See
1610 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU
1611 dx = -sinlon*de - sinlat*coslon*dn + coslat*coslon*du;
1612 dy = coslon*de - sinlat*sinlon*dn + coslat*sinlon*du;
1613 dz = coslat* dn + sinlat* du;
1614 }
1615}
1616
1617
1618
1620
1639 Numeric& lat,
1640 Numeric& lon,
1641 Numeric& za,
1642 Numeric& aa,
1643 const Numeric& x,
1644 const Numeric& y,
1645 const Numeric& z,
1646 const Numeric& dx,
1647 const Numeric& dy,
1648 const Numeric& dz,
1649 const Vector& refellipsoid ) {
1650
1651 cart2geodetic(h, lat, lon, x, y, z, refellipsoid );
1652
1653 const Numeric latrad = DEG2RAD * lat;
1654 const Numeric lonrad = DEG2RAD * lon;
1655 const Numeric coslat = cos(latrad);
1656 const Numeric sinlat = sin(latrad);
1657 const Numeric coslon = cos(lonrad);
1658 const Numeric sinlon = sin(lonrad);
1659
1660 // See
1661 // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU
1662 const Numeric de = -sinlon*dx + coslon*dy;
1663 const Numeric dn = -sinlat*coslon*dx - sinlat*sinlon*dy + coslat*dz;
1664 const Numeric du = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
1665
1666 enu2zaaa(za, aa, de, dn, du);
1667}
A constant view of a Vector.
Definition: matpackI.h:530
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
The Vector class.
Definition: matpackI.h:922
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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:669
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:182
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:1638
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition: geodetic.cc:1366
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:1059
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:542
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
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition: geodetic.cc:331
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1265
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:1566
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:1474
void cycle_lat_lon(Numeric &lat, Numeric &lon)
Cyclic latitude longitude coordinates.
Definition: geodetic.cc:1432
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1211
void line_circle_intersect(Numeric &x, Numeric &z, const Numeric &xl, const Numeric &zl, const Numeric &dx, const Numeric &dz, const Numeric &xc, const Numeric &zc, const Numeric &r)
geomtanpoint2d
Definition: geodetic.cc:278
const Numeric RAD2DEG
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:874
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1335
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
const Numeric DEG2RAD
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition: geodetic.cc:1408
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:1521
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1304
void cart2sph_plain(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z)
Definition: geodetic.cc:641
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
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:706
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:953
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
#define abs(x)
#define q
#define min(a, b)
#define dx
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:159
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:440
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
#define a2
#define b2
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
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:53
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 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
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
Propagation path structure and functions.
constexpr Numeric POLELAT
Size of north and south poles.
Definition: ppath.h:45
constexpr Numeric ANGTOL
Width of zenith and nadir directions.
Definition: ppath.h:56
#define d
#define v
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:736
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:728
#define N
Definition: rng.cc:164
Structure to store a grid position.
Definition: interpolation.h:73
std::array< Numeric, 2 > fd
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74