ARTS  2.4.0(git:4fb77825)
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 
42 extern const Numeric DEG2RAD;
43 extern 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 
74 void cart2pol(Numeric& r,
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  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 /*
230 void 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  assert( refellipsoid.nelem() == 2 );
239  assert( refellipsoid[0] > 0 );
240  assert( refellipsoid[1] >= 0 );
241  assert( refellipsoid[1] < 1 );
242  assert( r > 0 );
243  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  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  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  assert(l >= 0);
310  }
311 
312  x = xl + l * dx;
313  z = zl + l * dz;
314 }
315 
317 
331 void pol2cart(Numeric& x, Numeric& z, const Numeric& r, const Numeric& lat) {
332  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  assert(r > 0);
363  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 }
542 
544 
564  Numeric& lat,
565  Numeric& lon,
566  const Numeric& x,
567  const Numeric& y,
568  const Numeric& z,
569  const Numeric& lat0,
570  const Numeric& lon0,
571  const Numeric& za0,
572  const Numeric& aa0) {
573  r = sqrt(x * x + y * y + z * z);
574 
575  // Zenith and nadir cases
576  if (za0 < ANGTOL || za0 > 180 - ANGTOL) {
577  lat = lat0;
578  lon = lon0;
579  }
580 
581  else {
582  lat = RAD2DEG * asin(z / r);
583  lon = RAD2DEG * atan2(y, x);
584 
585  // Make sure that lon is maintained for N-S cases (if not
586  // starting on a pole)
587  if ((abs(aa0) < ANGTOL || abs(180 - aa0) < ANGTOL) &&
588  abs(lat0) <= POLELAT) {
589  // Check that not lon changed with 180 deg
590  if (abs(lon - lon0) < 1) {
591  lon = lon0;
592  } else {
593  if (lon0 > 0) {
594  lon = lon0 - 180;
595  } else {
596  lon = lon0 + 180;
597  }
598  }
599  }
600  }
601 }
602 
604 
619  const Numeric& r1,
620  const Numeric& lat1,
621  const Numeric& lon1,
622  const Numeric& r2,
623  const Numeric& lat2,
624  const Numeric& lon2) {
625  Numeric x1, y1, z1, x2, y2, z2;
626  sph2cart(x1, y1, z1, r1, lat1, lon1);
627  sph2cart(x2, y2, z2, r2, lat2, lon2);
628 
629  const Numeric dx = x2 - x1;
630  const Numeric dy = y2 - y1;
631  const Numeric dz = z2 - z1;
632  l = sqrt(dx * dx + dy * dy + dz * dz);
633 }
634 
636 
656  Numeric& lat_tan,
657  Numeric& lon_tan,
658  Numeric& l_tan,
659  const Numeric& r,
660  const Numeric& lat,
661  const Numeric& lon,
662  const Numeric& za,
663  const Numeric& aa,
664  const Numeric& ppc) {
665  assert(za >= 90);
666  assert(r >= ppc);
667 
668  Numeric x, y, z, dx, dy, dz;
669 
670  poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
671 
672  l_tan = sqrt(r * r - ppc * ppc);
673 
674  cart2sph(r_tan,
675  lat_tan,
676  lon_tan,
677  x + dx * l_tan,
678  y + dy * l_tan,
679  z + dz * l_tan,
680  lat,
681  lon,
682  za,
683  aa);
684 }
685 
687 
717 /*
718 void geomtanpoint(
719  Numeric& r_tan,
720  Numeric& lat_tan,
721  Numeric& lon_tan,
722  ConstVectorView refellipsoid,
723  const Numeric& r,
724  const Numeric& lat,
725  const Numeric& lon,
726  const Numeric& za,
727  const Numeric& aa )
728 {
729  assert( refellipsoid.nelem() == 2 );
730  assert( refellipsoid[0] > 0 );
731  assert( refellipsoid[1] >= 0 );
732  assert( refellipsoid[1] < 1 );
733  assert( r > 0 );
734  assert( za >= 90 );
735 
736  if( refellipsoid[1] < 1e-7 ) // e=1e-7 corresponds to that polar radius
737  { // less than 1 um smaller than equatorial
738  Numeric x, y, z, dx, dy, dz; // one for the Earth
739 
740  poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
741 
742  const Numeric ppc = r * sin( DEG2RAD * abs(za) );
743  const Numeric l_tan = sqrt( r*r - ppc*ppc );
744 
745  cart2sph( r_tan, lat_tan, lon_tan, x+dx*l_tan, y+dy*l_tan, z+dz*l_tan );
746  }
747 
748  else
749  {
750  // Equatorial and polar radii squared:
751  const Numeric a2 = refellipsoid[0]*refellipsoid[0];
752  const Numeric b2 = a2 * ( 1 - refellipsoid[1]*refellipsoid[1] );
753 
754  Vector X(3), xunit(3), yunit(3), zunit(3);
755 
756  poslos2cart( X[0], X[1], X[2], xunit[0], xunit[1], xunit[2],
757  r, lat, lon, za, aa );
758  cross( zunit, xunit, X );
759  unitl( zunit ); // Normalisation of length to 1
760 
761  cross( yunit, zunit, xunit );
762  unitl( yunit ); // Normalisation of length to 1
763 
764  Numeric x = X[0];
765  Numeric y = X[1];
766  const Numeric w11 = xunit[0];
767  const Numeric w12 = yunit[0];
768  const Numeric w21 = xunit[1];
769  const Numeric w22 = yunit[1];
770  const Numeric w31 = xunit[2];
771  const Numeric w32 = yunit[2];
772 
773  const Numeric yr = X * yunit;
774  const Numeric xr = X * xunit;
775 
776  const Numeric A = (w11*w11 + w21*w21)/a2 + w31*w31/b2;
777  const Numeric B = 2.0*((w11*w12 + w21*w22)/a2 + (w31*w32)/b2);
778  const Numeric C = (w12*w12 + w22*w22)/a2 + w32*w32/b2;
779 
780  if( B == 0.0 )
781  { x = 0.0; }
782  else
783  {
784  const Numeric K = -2.0*A/B;
785  const Numeric factor = 1.0/(A+(B+C*K)*K);
786  x = sqrt(factor);
787  y = K*x;
788  }
789 
790  const Numeric dist1 = (xr-X[0])*(xr-X[0]) + (yr-y)*(yr-y);
791  const Numeric dist2 = (xr+X[0])*(xr+X[0]) + (yr+y)*(yr+y);
792 
793  if( dist1 > dist2 )
794  { x = -x; }
795 
796  cart2sph( r_tan, lat_tan, lon_tan, w11*x + w12*yr, w21*x + w22*yr,
797  w31*x + w32*yr );
798  }
799 }
800 */
801 
803 
824  Numeric& y,
825  Numeric& z,
826  const Numeric& xl,
827  const Numeric& yl,
828  const Numeric& zl,
829  const Numeric& dx,
830  const Numeric& dy,
831  const Numeric& dz,
832  const Numeric& xc,
833  const Numeric& yc,
834  const Numeric& zc,
835  const Numeric& r) {
836  const Numeric a = dx * dx + dy * dy + dz * dz;
837  const Numeric b = 2 * (dx * (xl - xc) + dy * (yl - yc) + dz * (zl - zc));
838  const Numeric c = xc * xc + yc * yc + zc * zc + xl * xl + yl * yl + zl * zl -
839  2 * (xc * xl + yc * yl + zc * zl) - r * r;
840 
841  Numeric d = b * b - 4 * a * c;
842  assert(d > 0);
843 
844  const Numeric a2 = 2 * a;
845  const Numeric b2 = -b / a2;
846  const Numeric e = sqrt(d) / a2;
847 
848  const Numeric l1 = b2 + e;
849  const Numeric l2 = b2 - e;
850 
851  Numeric l;
852  if (l1 < 0) {
853  l = l2;
854  } else if (l2 < 0) {
855  l = l1;
856  } else {
857  l = min(l1, l2);
858  assert(l >= 0);
859  }
860 
861  x = xl + l * dx;
862  y = yl + l * dy;
863  z = zl + l * dz;
864 }
865 
867 
884  Numeric& lon2,
885  const Numeric& lat1,
886  const Numeric& lon1,
887  const Numeric& aa,
888  const Numeric& ddeg) {
889  // Code from http://www.movable-type.co.uk/scripts/latlong.html
890  // (but with short-cuts, such as asin(sin(lat2)) = lat2)
891  // Note that lat1 here is another latitude
892 
893  const Numeric dang = DEG2RAD * ddeg;
894  const Numeric cosdang = cos(dang);
895  const Numeric sindang = sin(dang);
896  const Numeric latrad = DEG2RAD * lat1;
897  const Numeric coslat = cos(latrad);
898  const Numeric sinlat = sin(latrad);
899  const Numeric aarad = DEG2RAD * aa;
900 
901  lat2 = sinlat * cosdang + coslat * sindang * cos(aarad);
902  lon2 = lon1 + RAD2DEG * (atan2(sin(aarad) * sindang * coslat,
903  cosdang - sinlat * lat2));
904  lat2 = RAD2DEG * asin(lat2);
905 }
906 
908 
929 void los2xyz(Numeric& za,
930  Numeric& aa,
931  const Numeric& r1,
932  const Numeric& lat1,
933  const Numeric& lon1,
934  const Numeric& x1,
935  const Numeric& y1,
936  const Numeric& z1,
937  const Numeric& x2,
938  const Numeric& y2,
939  const Numeric& z2) {
940  Numeric dx = x2 - x1, dy = y2 - y1, dz = z2 - z1;
941  const Numeric ldxyz = sqrt(dx * dx + dy * dy + dz * dz);
942  dx /= ldxyz;
943  dy /= ldxyz;
944  dz /= ldxyz;
945 
946  // All below extracted from 3D version of cart2poslos:
947  const Numeric latrad = DEG2RAD * lat1;
948  const Numeric lonrad = DEG2RAD * lon1;
949  const Numeric coslat = cos(latrad);
950  const Numeric sinlat = sin(latrad);
951  const Numeric coslon = cos(lonrad);
952  const Numeric sinlon = sin(lonrad);
953 
954  const Numeric dr = coslat * coslon * dx + coslat * sinlon * dy + sinlat * dz;
955  const Numeric dlat =
956  -sinlat * coslon / r1 * dx - sinlat * sinlon / r1 * dy + coslat / r1 * dz;
957  const Numeric dlon = -sinlon / coslat / r1 * dx + coslon / coslat / r1 * dy;
958 
959  za = RAD2DEG * acos(dr);
960  aa = RAD2DEG * acos(r1 * dlat / sin(DEG2RAD * za));
961  if (std::isnan(aa)) {
962  if (dlat >= 0) {
963  aa = 0;
964  } else {
965  aa = 180;
966  }
967  } else if (dlon < 0) {
968  aa = -aa;
969  }
970 }
971 
973 
998  Numeric& y,
999  Numeric& z,
1000  Numeric& dx,
1001  Numeric& dy,
1002  Numeric& dz,
1003  const Numeric& r,
1004  const Numeric& lat,
1005  const Numeric& lon,
1006  const Numeric& za,
1007  const Numeric& aa) {
1008  assert(r > 0);
1009  assert(abs(lat) <= 90);
1010  //assert( abs( lon ) <= 360 );
1011  assert(za >= 0 && za <= 180);
1012 
1013  // lat = +-90
1014  // For lat = +- 90 the azimuth angle gives the longitude along which the
1015  // LOS goes
1016  if (abs(lat) > POLELAT) {
1017  const Numeric s = sign(lat);
1018 
1019  x = 0;
1020  y = 0;
1021  z = s * r;
1022 
1023  dz = s * cos(DEG2RAD * za);
1024  dx = sin(DEG2RAD * za);
1025  dy = dx * sin(DEG2RAD * aa);
1026  dx = dx * cos(DEG2RAD * aa);
1027  }
1028 
1029  else {
1030  const Numeric latrad = DEG2RAD * lat;
1031  const Numeric lonrad = DEG2RAD * lon;
1032  const Numeric zarad = DEG2RAD * za;
1033  const Numeric aarad = DEG2RAD * aa;
1034 
1035  const Numeric coslat = cos(latrad);
1036  const Numeric sinlat = sin(latrad);
1037  const Numeric coslon = cos(lonrad);
1038  const Numeric sinlon = sin(lonrad);
1039  const Numeric cosza = cos(zarad);
1040  const Numeric sinza = sin(zarad);
1041  const Numeric cosaa = cos(aarad);
1042  const Numeric sinaa = sin(aarad);
1043 
1044  // This part as sph2cart but uses local variables
1045  x = r * coslat; // Common term for x and y
1046  y = x * sinlon;
1047  x = x * coslon;
1048  z = r * sinlat;
1049 
1050  const Numeric dr = cosza;
1051  const Numeric dlat = sinza * cosaa; // r-term cancel out below
1052  const Numeric dlon = sinza * sinaa / coslat;
1053 
1054  dx = coslat * coslon * dr - sinlat * coslon * dlat - coslat * sinlon * dlon;
1055  dz = sinlat * dr + coslat * dlat;
1056  dy = coslat * sinlon * dr - sinlat * sinlon * dlat + coslat * coslon * dlon;
1057  }
1058 }
1059 
1061 
1086  if (atmosphere_dim == 1) {
1087  return refellipsoid[0];
1088  } else {
1089  assert(rte_pos.nelem() > 1);
1090 
1091  bool inside = true;
1092 
1093  if (rte_pos[1] < lat_grid[0] || rte_pos[1] > last(lat_grid)) {
1094  inside = false;
1095  } else if (atmosphere_dim == 3) {
1096  assert(rte_pos.nelem() == 3);
1097  if (rte_pos[2] < lon_grid[0] || rte_pos[2] > last(lon_grid)) {
1098  inside = false;
1099  }
1100  }
1101 
1102  if (inside) {
1103  GridPos gp_lat;
1104  gridpos(gp_lat, lat_grid, rte_pos[1]);
1105  return refell2d(refellipsoid, lat_grid, gp_lat);
1106  } else {
1107  return refell2r(refellipsoid, rte_pos[1]);
1108  }
1109  }
1110 }
1111 
1113 
1136  assert(refellipsoid.nelem() == 2);
1137  assert(refellipsoid[0] > 0);
1138  assert(refellipsoid[1] >= 0);
1139  assert(refellipsoid[1] < 1);
1140 
1141  if (refellipsoid[1] < 1e-7) // e=1e-7 corresponds to that polar radius
1142  { // less than 1 um smaller than equatorial
1143  return refellipsoid[0]; // one for the Earth
1144  }
1145 
1146  else {
1147  const Numeric c = 1 - refellipsoid[1] * refellipsoid[1];
1148  const Numeric b = refellipsoid[0] * sqrt(c);
1149  const Numeric v = DEG2RAD * lat;
1150  const Numeric ct = cos(v);
1151  const Numeric st = sin(v);
1152 
1153  return b / sqrt(c * ct * ct + st * st);
1154  }
1155 }
1156 
1158 
1176  const GridPos gp) {
1177  if (gp.fd[0] == 0)
1178  return refell2r(refellipsoid, lat_grid[gp.idx]);
1179  else if (gp.fd[0] == 1)
1180  return refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1181  else
1182  return gp.fd[1] * refell2r(refellipsoid, lat_grid[gp.idx]) +
1183  gp.fd[0] * refell2r(refellipsoid, lat_grid[gp.idx + 1]);
1184 }
1185 
1187 
1206  const Numeric& lon1,
1207  const Numeric& lat2,
1208  const Numeric& lon2) {
1209  // Equations taken from http://www.movable-type.co.uk/scripts/latlong.html
1210  const Numeric slat = sin(DEG2RAD * (lat2 - lat1) / 2.0);
1211  const Numeric slon = sin(DEG2RAD * (lon2 - lon1) / 2.0);
1212  const Numeric a =
1213  slat * slat + cos(DEG2RAD * lat1) * cos(DEG2RAD * lat2) * slon * slon;
1214 
1215  return RAD2DEG * 2 * atan2(sqrt(a), sqrt(1 - a));
1216 }
1217 
1219 
1237  Numeric& y,
1238  Numeric& z,
1239  const Numeric& r,
1240  const Numeric& lat,
1241  const Numeric& lon) {
1242  assert(r > 0);
1243  assert(abs(lat) <= 90);
1244  assert(abs(lon) <= 360);
1245 
1246  const Numeric latrad = DEG2RAD * lat;
1247  const Numeric lonrad = DEG2RAD * lon;
1248 
1249  x = r * cos(latrad); // Common term for x and z
1250  y = x * sin(lonrad);
1251  x = x * cos(lonrad);
1252  z = r * sin(latrad);
1253 }
1254 
1255 /*===========================================================================
1256  === coordinate transformations
1257  ===========================================================================*/
1258 
1260 
1276 void lon_shiftgrid(Vector& longrid_out,
1277  ConstVectorView longrid_in,
1278  const Numeric lon) {
1279  longrid_out = longrid_in;
1280  if (longrid_in[longrid_in.nelem() - 1] >= lon + 360.)
1281  longrid_out += -360.;
1282  else if (longrid_in[0] <= lon - 360.)
1283  longrid_out += 360.;
1284 }
1285 
1287 /* This function wraps around the given latitude and longitude coordinates
1288  * to the ranges
1289  * - lat in [-90.0, 90.0]
1290  * - lon in [0.0, 360.0]
1291  * Only in the range [-180.0, 180.0] are accepted otherwise an error will
1292  * be thrown
1293  *
1294  * \param[in, out] lat The latitude coordinate.
1295  * \param[in, out] lon The longitude coordiante.
1296  *
1297  * \author Simon Pfreundschuh
1298  * \date 2018-05-22
1299 */
1301  if (lat < -180.0) {
1302  throw std::runtime_error("Latitude values < -180.0 are not supported.");
1303  }
1304  if (lat > 180.0) {
1305  throw std::runtime_error("Latitude values > 180.0 are not supported.");
1306  }
1307 
1308  if (lat < -90.0) {
1309  lat = -180.0 - lat;
1310  lon += 180.0;
1311  }
1312  if (lat > 90.0) {
1313  lat = 180.0 - lat;
1314  lon += 180.0;
1315  }
1316 
1317  while (lon < 0.0) {
1318  lon += 360.0;
1319  }
1320  while (lon > 360.0) {
1321  lon -= 360.0;
1322  }
1323 }
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
sph2cart
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition: geodetic.cc:1236
ANGTOL
const Numeric ANGTOL
Width of zenith and nadir directions.
Definition: ppath.h:108
cycle_lat_lon
void cycle_lat_lon(Numeric &lat, Numeric &lon)
Cyclic latitude longitude coordinates.
Definition: geodetic.cc:1300
ARTS::Var::lat
Numeric lat(Workspace &ws) noexcept
Definition: autoarts.h:3933
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
sign
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:423
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
refell2d
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1174
last
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
abs
#define abs(x)
Definition: legacy_continua.cc:20626
sphdist
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1205
DEG2RAD
const Numeric DEG2RAD
POLELAT
const Numeric POLELAT
Size of north and south poles.
Definition: ppath.h:97
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
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:655
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)
geomtanpoint
Definition: geodetic.cc:823
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:563
x0
#define x0
Definition: lineshapemodel.h:413
b2
#define b2
Definition: complex.h:59
GridPos::fd
Numeric fd[2]
Definition: interpolation.h:75
ARTS::Var::lon
Numeric lon(Workspace &ws) noexcept
Definition: autoarts.h:4064
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
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:113
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
cart2pol
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:74
min
#define min(a, b)
Definition: legacy_continua.cc:20628
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
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:1081
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:278
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
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:158
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:929
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
RAD2DEG
const Numeric RAD2DEG
math_funcs.h
lon_shiftgrid
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition: geodetic.cc:1276
geodetic.h
ARTS::Var::refellipsoid
Vector refellipsoid(Workspace &ws) noexcept
Definition: autoarts.h:5519
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ppath.h
Propagation path structure and functions.
pol2cart
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition: geodetic.cc:331
x1
#define x1
Definition: lineshapemodel.h:414
dx
#define dx
Definition: legacy_continua.cc:22180
poslos2cart
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:355
GridPos::idx
Index idx
Definition: interpolation.h:74
a2
#define a2
Definition: complex.h:58
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
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:883
Vector
The Vector class.
Definition: matpackI.h:860
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:618
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
x2
#define x2
Definition: lineshapemodel.h:415
refell2r
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1135
distance2D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:182