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