ARTS  2.2.66
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 
21 
22 /*===========================================================================
23  === File description
24  ===========================================================================*/
25 
37 /*===========================================================================
38  === External declarations
39  ===========================================================================*/
40 
41 #include <cmath>
42 #include <stdexcept>
43 #include "geodetic.h"
44 #include "math_funcs.h"
45 #include "ppath.h"
46 
47 extern const Numeric DEG2RAD;
48 extern const Numeric RAD2DEG;
49 
50 
51 
52 /*===========================================================================
53  === 2D functions
54  ===========================================================================*/
55 
56 // The 2D case is treated as being the 3D x/z-plane. That is, y-coordinate is
57 // skipped. For simplicity, the angle coordinate is denoted as latitude.
58 // However, the latitude is here not limited to [-90,90]. It is cyclic and can
59 // have any value. The input *lat0* is used to shift the output from atan2 with
60 // n*360 to return what should be the expected latitude. That is, it is assumed
61 // that no operation moves the latitude more than 180 degrees from the initial
62 // value *lat0*. Negative zeniath angles are handled, following ARTS definition
63 // of 2D geometry.
64 
65 
67 
82 void cart2pol(
83  Numeric& r,
84  Numeric& lat,
85  const Numeric& x,
86  const Numeric& z,
87  const Numeric& lat0,
88  const Numeric& za0 )
89 {
90  r = sqrt( x*x + z*z );
91 
92  // Zenith and nadir cases
93  const Numeric absza = abs( za0 );
94  if( absza < ANGTOL || absza > 180-ANGTOL )
95  { lat = lat0; }
96 
97  else
98  { // Latitude inside [0,360]
99  lat = RAD2DEG * atan2( z, x );
100  // Shift with n*360 to get as close as possible to lat0
101  lat = lat - 360.0 * Numeric( round( ( lat - lat0 ) / 360.0 ) );
102  }
103 }
104 
105 
106 
108 
126  Numeric& r,
127  Numeric& lat,
128  Numeric& za,
129  const Numeric& x,
130  const Numeric& z,
131  const Numeric& dx,
132  const Numeric& dz,
133  const Numeric& ppc,
134  const Numeric& lat0,
135  const Numeric& za0 )
136 {
137  r = sqrt( x*x + z*z );
138 
139  // Zenith and nadir cases
140  const Numeric absza = abs( za0 );
141  if( absza < ANGTOL || absza > 180-ANGTOL )
142  {
143  lat = lat0;
144  za = za0;
145  }
146 
147  else
148  {
149  lat = RAD2DEG * atan2( z, x );
150 
151  const Numeric latrad = DEG2RAD * lat;
152  const Numeric coslat = cos( latrad );
153  const Numeric sinlat = sin( latrad );
154  const Numeric dr = coslat*dx + sinlat*dz;
155 
156  // Use ppc for max accuracy, but dr required to resolve if up-
157  // and downward cases.
158 
159  // Another possibility to obtain (absolute value of) za is
160  // RAD2DEG*acos(dr).
161  // It is checked that the two ways give consistent results, but
162  // occasionally deviate with 1e-4 deg (due to numerical issues).
163 
164  za = RAD2DEG * asin( ppc / r );
165  if( za0 > 0 )
166  {
167  if( isnan( za ) )
168  { za = 90; }
169  else if( dr < 0 )
170  { za = 180.0 - za; }
171  }
172  else
173  {
174  if( isnan( za ) )
175  { za = -90; }
176  else if( dr < 0 )
177  { za = -180.0 + za; }
178  else
179  { za = -za; }
180  }
181  }
182 }
183 
184 
185 
187 
202  Numeric& l,
203  const Numeric& r1,
204  const Numeric& lat1,
205  const Numeric& r2,
206  const Numeric& lat2 )
207 {
208  assert( abs( lat2 - lat1 ) <= 180 );
209 
210  Numeric x1, z1, x2, z2;
211  pol2cart( x1, z1, r1, lat1 );
212  pol2cart( x2, z2, r2, lat2 );
213 
214  const Numeric dx = x2 - x1;
215  const Numeric dz = z2 - z1;
216  l = sqrt( dx*dx + dz*dz );
217 }
218 
219 
220 
222 
252 /*
253 void geomtanpoint2d(
254  Numeric& r_tan,
255  Numeric& lat_tan,
256  ConstVectorView refellipsoid,
257  const Numeric& r,
258  const Numeric& lat,
259  const Numeric& za )
260 {
261  assert( refellipsoid.nelem() == 2 );
262  assert( refellipsoid[0] > 0 );
263  assert( refellipsoid[1] >= 0 );
264  assert( refellipsoid[1] < 1 );
265  assert( r > 0 );
266  assert( za >= 90 );
267  // e=1e-7 corresponds to that polar radius
268  if( refellipsoid[1] < 1e-7 ) // less than 1 um smaller than equatorial
269  { // one for the Earth
270  r_tan = geometrical_ppc( r, za );
271  if( za > 0 )
272  { lat_tan = geompath_lat_at_za( za, lat, 90 ); }
273  else
274  { lat_tan = geompath_lat_at_za( za, lat, -90 ); }
275  }
276 
277  else
278  {
279  assert( 0 ); // To be implemented
280  }
281 }
282 */
283 
284 
286 
303  Numeric& x,
304  Numeric& z,
305  const Numeric& xl,
306  const Numeric& zl,
307  const Numeric& dx,
308  const Numeric& dz,
309  const Numeric& xc,
310  const Numeric& zc,
311  const Numeric& r )
312 {
313  const Numeric a = dx*dx + dz*dz;
314  const Numeric b = 2*( dx*(xl-xc) + dz*(zl-zc) );
315  const Numeric c = xc*xc + zc*zc + xl*xl + zl*zl - 2*(xc*xl + zc*zl) - r*r;
316 
317  Numeric d = b*b - 4*a*c;
318  assert( d > 0 );
319 
320  const Numeric a2 = 2*a;
321  const Numeric b2 = -b / a2;
322  const Numeric e = sqrt( d ) / a2;
323 
324  const Numeric l1 = b2 + e;
325  const Numeric l2 = b2 - e;
326 
327  Numeric l;
328  if( l1 < 0 )
329  { l = l2; }
330  else if( l2 < 0 )
331  { l = l1; }
332  else
333  { l = min(l1,l2); assert( l>=0 ); }
334 
335  x = xl + l*dx;
336  z = zl + l*dz;
337 }
338 
339 
340 
342 
356 void pol2cart(
357  Numeric& x,
358  Numeric& z,
359  const Numeric& r,
360  const Numeric& lat )
361 {
362  assert( r > 0 );
363 
364  const Numeric latrad = DEG2RAD * lat;
365 
366  x = r * cos( latrad );
367  z = r * sin( latrad );
368 }
369 
370 
371 
373 
388  Numeric& x,
389  Numeric& z,
390  Numeric& dx,
391  Numeric& dz,
392  const Numeric& r,
393  const Numeric& lat,
394  const Numeric& za )
395 {
396  assert( r > 0 );
397  assert( za >= -180 && za<=180 );
398 
399  const Numeric latrad = DEG2RAD * lat;
400  const Numeric zarad = DEG2RAD * za;
401 
402  const Numeric coslat = cos( latrad );
403  const Numeric sinlat = sin( latrad );
404  const Numeric cosza = cos( zarad );
405  const Numeric sinza = sin( zarad );
406 
407  // This part as pol2cart but uses local variables
408  x = r * coslat;
409  z = r * sinlat;
410 
411  const Numeric dr = cosza;
412  const Numeric dlat = sinza; // r-term cancel out below
413 
414  dx = coslat * dr - sinlat * dlat;
415  dz = sinlat * dr + coslat * dlat;
416 }
417 
418 
419 
420 
421 
422 /*===========================================================================
423  === 3D functions
424  ===========================================================================*/
425 
427 
460  Numeric& r,
461  Numeric& lat,
462  Numeric& lon,
463  Numeric& za,
464  Numeric& aa,
465  const Numeric& x,
466  const Numeric& y,
467  const Numeric& z,
468  const Numeric& dx,
469  const Numeric& dy,
470  const Numeric& dz,
471  const Numeric& ppc,
472  const Numeric& lat0,
473  const Numeric& lon0,
474  const Numeric& za0,
475  const Numeric& aa0 )
476 {
477  r = sqrt( x*x + y*y + z*z );
478 
479  // Zenith and nadir cases
480  if( za0 < ANGTOL || za0 > 180-ANGTOL )
481  {
482  lat = lat0;
483  lon = lon0;
484  za = za0;
485  aa = aa0;
486  }
487 
488  else
489  {
490  lat = RAD2DEG * asin( z / r );
491  lon = RAD2DEG * atan2( y, x );
492 
493  bool ns_case = false;
494  bool lon_flip = false;
495 
496  // Make sure that lon is maintained for N-S cases (if not
497  // starting on a pole)
498  if( ( abs(aa0) < ANGTOL || abs(180-aa0) < ANGTOL ) &&
499  abs( lat0 ) <= POLELAT )
500  {
501  ns_case = true;
502  // Check that not lon changed with 180 deg
503  if( abs(lon-lon0) < 1 )
504  { lon = lon0; }
505  else
506  {
507  lon_flip = true;
508  if( lon0 > 0 )
509  { lon = lon0 - 180; }
510  else
511  { lon = lon0 + 180; }
512  }
513  }
514 
515  const Numeric latrad = DEG2RAD * lat;
516  const Numeric lonrad = DEG2RAD * lon;
517  const Numeric coslat = cos( latrad );
518  const Numeric sinlat = sin( latrad );
519  const Numeric coslon = cos( lonrad );
520  const Numeric sinlon = sin( lonrad );
521  const Numeric dr = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
522 
523  // Use ppc for max accuracy, but dr required to resolve if up-
524  // and downward cases
525  za = RAD2DEG * asin( ppc / r );
526  if( isnan( za ) )
527  { za = 90; }
528  if( dr < 0 )
529  { za = 180.0 - za; }
530 
531  // The difference below can at least be 3e-6 for tangent points
532  assert( abs( za - RAD2DEG*acos(dr) ) < 1e-4 );
533 
534  // For lat = +- 90 the azimuth angle gives the longitude along which
535  // the LOS goes
536  if( abs( lat ) >= POLELAT )
537  { aa = RAD2DEG * atan2( dy, dx ); }
538 
539  // N-S cases, not starting at a pole
540  else if( ns_case )
541  {
542  if( !lon_flip )
543  { aa = aa0; }
544  else
545  {
546  if( abs(aa0) < ANGTOL )
547  { aa = 180; }
548  else
549  { aa = 0; }
550  }
551  }
552 
553  else
554  {
555  const Numeric dlat = -sinlat*coslon/r*dx - sinlat*sinlon/r*dy +
556  coslat/r*dz;
557  const Numeric dlon = -sinlon/coslat/r*dx + coslon/coslat/r*dy;
558 
559  aa = RAD2DEG * acos( r * dlat / sin( DEG2RAD * za ) );
560 
561  if( isnan( aa ) )
562  {
563  if( dlat >= 0 )
564  { aa = 0; }
565  else
566  { aa = 180; }
567  }
568  else if( dlon < 0 )
569  { aa = -aa; }
570  }
571  }
572 }
573 
574 
575 
577 
596 void cart2sph(
597  Numeric& r,
598  Numeric& lat,
599  Numeric& lon,
600  const Numeric& x,
601  const Numeric& y,
602  const Numeric& z,
603  const Numeric& lat0,
604  const Numeric& lon0,
605  const Numeric& za0,
606  const Numeric& aa0 )
607 {
608  r = sqrt( x*x + y*y + z*z );
609 
610  // Zenith and nadir cases
611  if( za0 < ANGTOL || za0 > 180-ANGTOL )
612  {
613  lat = lat0;
614  lon = lon0;
615  }
616 
617  else
618  {
619  lat = RAD2DEG * asin( z / r );
620  lon = RAD2DEG * atan2( y, x );
621 
622  // Make sure that lon is maintained for N-S cases (if not
623  // starting on a pole)
624  if( ( abs(aa0) < ANGTOL || abs(180-aa0) < ANGTOL ) &&
625  abs( lat0 ) <= POLELAT )
626  {
627  // Check that not lon changed with 180 deg
628  if( abs(lon-lon0) < 1 )
629  { lon = lon0; }
630  else
631  {
632  if( lon0 > 0 )
633  { lon = lon0 - 180; }
634  else
635  { lon = lon0 + 180; }
636  }
637  }
638  }
639 }
640 
641 
642 
644 
659  Numeric& l,
660  const Numeric& r1,
661  const Numeric& lat1,
662  const Numeric& lon1,
663  const Numeric& r2,
664  const Numeric& lat2,
665  const Numeric& lon2 )
666 {
667  Numeric x1, y1, z1, x2, y2, z2;
668  sph2cart( x1, y1, z1, r1, lat1, lon1 );
669  sph2cart( x2, y2, z2, r2, lat2, lon2 );
670 
671  const Numeric dx = x2 - x1;
672  const Numeric dy = y2 - y1;
673  const Numeric dz = z2 - z1;
674  l = sqrt( dx*dx + dy*dy + dz*dz );
675 }
676 
677 
678 
680 
700  Numeric& r_tan,
701  Numeric& lat_tan,
702  Numeric& lon_tan,
703  Numeric& l_tan,
704  const Numeric& r,
705  const Numeric& lat,
706  const Numeric& lon,
707  const Numeric& za,
708  const Numeric& aa,
709  const Numeric& ppc )
710 {
711  assert( za >= 90 );
712  assert( r >= ppc );
713 
714  Numeric x, y, z, dx, dy, dz;
715 
716  poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
717 
718  l_tan = sqrt( r*r - ppc*ppc );
719 
720  cart2sph( r_tan, lat_tan, lon_tan, x+dx*l_tan, y+dy*l_tan, z+dz*l_tan,
721  lat, lon, za, aa );
722 }
723 
724 
725 
727 
757 /*
758 void geomtanpoint(
759  Numeric& r_tan,
760  Numeric& lat_tan,
761  Numeric& lon_tan,
762  ConstVectorView refellipsoid,
763  const Numeric& r,
764  const Numeric& lat,
765  const Numeric& lon,
766  const Numeric& za,
767  const Numeric& aa )
768 {
769  assert( refellipsoid.nelem() == 2 );
770  assert( refellipsoid[0] > 0 );
771  assert( refellipsoid[1] >= 0 );
772  assert( refellipsoid[1] < 1 );
773  assert( r > 0 );
774  assert( za >= 90 );
775 
776  if( refellipsoid[1] < 1e-7 ) // e=1e-7 corresponds to that polar radius
777  { // less than 1 um smaller than equatorial
778  Numeric x, y, z, dx, dy, dz; // one for the Earth
779 
780  poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
781 
782  const Numeric ppc = r * sin( DEG2RAD * abs(za) );
783  const Numeric l_tan = sqrt( r*r - ppc*ppc );
784 
785  cart2sph( r_tan, lat_tan, lon_tan, x+dx*l_tan, y+dy*l_tan, z+dz*l_tan );
786  }
787 
788  else
789  {
790  // Equatorial and polar radii squared:
791  const Numeric a2 = refellipsoid[0]*refellipsoid[0];
792  const Numeric b2 = a2 * ( 1 - refellipsoid[1]*refellipsoid[1] );
793 
794  Vector X(3), xunit(3), yunit(3), zunit(3);
795 
796  poslos2cart( X[0], X[1], X[2], xunit[0], xunit[1], xunit[2],
797  r, lat, lon, za, aa );
798  cross( zunit, xunit, X );
799  unitl( zunit ); // Normalisation of length to 1
800 
801  cross( yunit, zunit, xunit );
802  unitl( yunit ); // Normalisation of length to 1
803 
804  Numeric x = X[0];
805  Numeric y = X[1];
806  const Numeric w11 = xunit[0];
807  const Numeric w12 = yunit[0];
808  const Numeric w21 = xunit[1];
809  const Numeric w22 = yunit[1];
810  const Numeric w31 = xunit[2];
811  const Numeric w32 = yunit[2];
812 
813  const Numeric yr = X * yunit;
814  const Numeric xr = X * xunit;
815 
816  const Numeric A = (w11*w11 + w21*w21)/a2 + w31*w31/b2;
817  const Numeric B = 2.0*((w11*w12 + w21*w22)/a2 + (w31*w32)/b2);
818  const Numeric C = (w12*w12 + w22*w22)/a2 + w32*w32/b2;
819 
820  if( B == 0.0 )
821  { x = 0.0; }
822  else
823  {
824  const Numeric K = -2.0*A/B;
825  const Numeric factor = 1.0/(A+(B+C*K)*K);
826  x = sqrt(factor);
827  y = K*x;
828  }
829 
830  const Numeric dist1 = (xr-X[0])*(xr-X[0]) + (yr-y)*(yr-y);
831  const Numeric dist2 = (xr+X[0])*(xr+X[0]) + (yr+y)*(yr+y);
832 
833  if( dist1 > dist2 )
834  { x = -x; }
835 
836  cart2sph( r_tan, lat_tan, lon_tan, w11*x + w12*yr, w21*x + w22*yr,
837  w31*x + w32*yr );
838  }
839 }
840 */
841 
842 
843 
845 
866  Numeric& x,
867  Numeric& y,
868  Numeric& z,
869  const Numeric& xl,
870  const Numeric& yl,
871  const Numeric& zl,
872  const Numeric& dx,
873  const Numeric& dy,
874  const Numeric& dz,
875  const Numeric& xc,
876  const Numeric& yc,
877  const Numeric& zc,
878  const Numeric& r )
879 {
880  const Numeric a = dx*dx + dy*dy + dz*dz;
881  const Numeric b = 2*( dx*(xl-xc) + dy*(yl-yc) + dz*(zl-zc) );
882  const Numeric c = xc*xc + yc*yc + zc*zc +
883  xl*xl + yl*yl + zl*zl - 2*(xc*xl + yc*yl + zc*zl) - r*r;
884 
885  Numeric d = b*b - 4*a*c;
886  assert( d > 0 );
887 
888  const Numeric a2 = 2*a;
889  const Numeric b2 = -b / a2;
890  const Numeric e = sqrt( d ) / a2;
891 
892  const Numeric l1 = b2 + e;
893  const Numeric l2 = b2 - e;
894 
895  Numeric l;
896  if( l1 < 0 )
897  { l = l2; }
898  else if( l2 < 0 )
899  { l = l1; }
900  else
901  { l = min(l1,l2); assert( l>=0 ); }
902 
903  x = xl + l*dx;
904  y = yl + l*dy;
905  z = zl + l*dz;
906 }
907 
908 
909 
911 
928  Numeric& lat2,
929  Numeric& lon2,
930  const Numeric& lat1,
931  const Numeric& lon1,
932  const Numeric& aa,
933  const Numeric& ddeg )
934 {
935  // Code from http://www.movable-type.co.uk/scripts/latlong.html
936  // (but with short-cuts, such as asin(sin(lat2)) = lat2)
937  // Note that lat1 here is another latitude
938 
939  const Numeric dang = DEG2RAD * ddeg;
940  const Numeric cosdang= cos( dang );
941  const Numeric sindang= sin( dang );
942  const Numeric latrad = DEG2RAD * lat1;
943  const Numeric coslat = cos( latrad );
944  const Numeric sinlat = sin( latrad );
945  const Numeric aarad = DEG2RAD * aa;
946 
947  lat2 = sinlat*cosdang + coslat*sindang*cos(aarad);
948  lon2 = lon1 + RAD2DEG*( atan2( sin(aarad)*sindang*coslat,
949  cosdang-sinlat*lat2 ) );
950  lat2 = RAD2DEG * asin( lat2 );
951 }
952 
953 
954 
956 
977 void los2xyz(
978  Numeric& za,
979  Numeric& aa,
980  const Numeric& r1,
981  const Numeric& lat1,
982  const Numeric& lon1,
983  const Numeric& x1,
984  const Numeric& y1,
985  const Numeric& z1,
986  const Numeric& x2,
987  const Numeric& y2,
988  const Numeric& z2 )
989 {
990  Numeric dx = x2-x1, dy = y2-y1, dz = z2-z1;
991  const Numeric ldxyz = sqrt( dx*dx + dy*dy + dz*dz );
992  dx /= ldxyz;
993  dy /= ldxyz;
994  dz /= ldxyz;
995 
996  // All below extracted from 3D version of cart2poslos:
997  const Numeric latrad = DEG2RAD * lat1;
998  const Numeric lonrad = DEG2RAD * lon1;
999  const Numeric coslat = cos( latrad );
1000  const Numeric sinlat = sin( latrad );
1001  const Numeric coslon = cos( lonrad );
1002  const Numeric sinlon = sin( lonrad );
1003 
1004  const Numeric dr = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
1005  const Numeric dlat = -sinlat*coslon/r1*dx - sinlat*sinlon/r1*dy +
1006  coslat/r1*dz;
1007  const Numeric dlon = -sinlon/coslat/r1*dx + coslon/coslat/r1*dy;
1008 
1009  za = RAD2DEG * acos( dr );
1010  aa = RAD2DEG * acos( r1 * dlat / sin( DEG2RAD * za ) );
1011  if( isnan( aa ) )
1012  {
1013  if( dlat >= 0 )
1014  { aa = 0; }
1015  else
1016  { aa = 180; }
1017  }
1018  else if( dlon < 0 )
1019  { aa = -aa; }
1020 }
1021 
1022 
1023 
1025 
1050  Numeric& x,
1051  Numeric& y,
1052  Numeric& z,
1053  Numeric& dx,
1054  Numeric& dy,
1055  Numeric& dz,
1056  const Numeric& r,
1057  const Numeric& lat,
1058  const Numeric& lon,
1059  const Numeric& za,
1060  const Numeric& aa )
1061 {
1062  assert( r > 0 );
1063  assert( abs( lat ) <= 90 );
1064  assert( abs( lon ) <= 360 );
1065  assert( za >= 0 && za<=180 );
1066 
1067  // lat = +-90
1068  // For lat = +- 90 the azimuth angle gives the longitude along which the
1069  // LOS goes
1070  if( abs( lat ) > POLELAT )
1071  {
1072  const Numeric s = sign( lat );
1073 
1074  x = 0;
1075  y = 0;
1076  z = s * r;
1077 
1078  dz = s * cos( DEG2RAD * za );
1079  dx = sin( DEG2RAD * za );
1080  dy = dx * sin( DEG2RAD * aa );
1081  dx = dx * cos( DEG2RAD * aa );
1082  }
1083 
1084  else
1085  {
1086  const Numeric latrad = DEG2RAD * lat;
1087  const Numeric lonrad = DEG2RAD * lon;
1088  const Numeric zarad = DEG2RAD * za;
1089  const Numeric aarad = DEG2RAD * aa;
1090 
1091  const Numeric coslat = cos( latrad );
1092  const Numeric sinlat = sin( latrad );
1093  const Numeric coslon = cos( lonrad );
1094  const Numeric sinlon = sin( lonrad );
1095  const Numeric cosza = cos( zarad );
1096  const Numeric sinza = sin( zarad );
1097  const Numeric cosaa = cos( aarad );
1098  const Numeric sinaa = sin( aarad );
1099 
1100  // This part as sph2cart but uses local variables
1101  x = r * coslat; // Common term for x and y
1102  y = x * sinlon;
1103  x = x * coslon;
1104  z = r * sinlat;
1105 
1106  const Numeric dr = cosza;
1107  const Numeric dlat = sinza * cosaa; // r-term cancel out below
1108  const Numeric dlon = sinza * sinaa / coslat;
1109 
1110  dx = coslat*coslon * dr - sinlat*coslon * dlat - coslat*sinlon * dlon;
1111  dz = sinlat * dr + coslat * dlat;
1112  dy = coslat*sinlon * dr - sinlat*sinlon * dlat + coslat*coslon * dlon;
1113  }
1114 }
1115 
1116 
1117 
1119 
1140  const Index& atmosphere_dim,
1141  ConstVectorView refellipsoid,
1142  ConstVectorView lat_grid,
1143  ConstVectorView lon_grid,
1144  ConstVectorView rte_pos )
1145 {
1146  if( atmosphere_dim == 1 )
1147  { return refellipsoid[0]; }
1148  else
1149  {
1150  assert( rte_pos.nelem() > 1 );
1151 
1152  bool inside = true;
1153 
1154  if( rte_pos[1] < lat_grid[0] || rte_pos[1] > last(lat_grid) )
1155  { inside = false; }
1156  else if( atmosphere_dim == 3 )
1157  {
1158  assert( rte_pos.nelem() == 3 );
1159  if( rte_pos[2] < lon_grid[0] || rte_pos[2] > last(lon_grid) )
1160  { inside = false; }
1161  }
1162 
1163  if( inside )
1164  {
1165  GridPos gp_lat;
1166  gridpos( gp_lat, lat_grid, rte_pos[1] );
1167  return refell2d( refellipsoid, lat_grid, gp_lat );
1168  }
1169  else
1170  { return refell2r( refellipsoid, rte_pos[1] ); }
1171  }
1172 }
1173 
1174 
1175 
1177 
1200  ConstVectorView refellipsoid,
1201  const Numeric& lat )
1202 {
1203  assert( refellipsoid.nelem() == 2 );
1204  assert( refellipsoid[0] > 0 );
1205  assert( refellipsoid[1] >= 0 );
1206  assert( refellipsoid[1] < 1 );
1207 
1208  if( refellipsoid[1] < 1e-7 ) // e=1e-7 corresponds to that polar radius
1209  { // less than 1 um smaller than equatorial
1210  return refellipsoid[0]; // one for the Earth
1211  }
1212 
1213  else
1214  {
1215  const Numeric c = 1 - refellipsoid[1]*refellipsoid[1];
1216  const Numeric b = refellipsoid[0] * sqrt( c );
1217  const Numeric v = DEG2RAD * lat;
1218  const Numeric ct = cos( v );
1219  const Numeric st = sin( v );
1220 
1221  return b / sqrt( c*ct*ct + st*st );
1222  }
1223 }
1224 
1225 
1226 
1228 
1245  ConstVectorView refellipsoid,
1246  ConstVectorView lat_grid,
1247  const GridPos gp )
1248 {
1249  if( gp.fd[0] == 0 )
1250  return refell2r(refellipsoid,lat_grid[gp.idx]);
1251  else if( gp.fd[0] == 1 )
1252  return refell2r(refellipsoid,lat_grid[gp.idx+1]);
1253  else
1254  return gp.fd[1] * refell2r(refellipsoid,lat_grid[gp.idx]) +
1255  gp.fd[0] * refell2r(refellipsoid,lat_grid[gp.idx+1]);
1256 }
1257 
1258 
1259 
1261 
1277  const Numeric& lat1,
1278  const Numeric& lon1,
1279  const Numeric& lat2,
1280  const Numeric& lon2 )
1281 {
1282  // Equations taken from http://www.movable-type.co.uk/scripts/latlong.html
1283  const Numeric slat = sin( DEG2RAD*(lat2-lat1) / 2.0 );
1284  const Numeric slon = sin( DEG2RAD*(lon2-lon1) / 2.0 );
1285  const Numeric a = slat*slat + cos(DEG2RAD*lat1)*cos(DEG2RAD*lat2)*slon*slon;
1286 
1287  return RAD2DEG * 2 * atan2( sqrt(a), sqrt(1-a) );
1288 }
1289 
1290 
1291 
1292 
1293 
1295 
1313  Numeric& x,
1314  Numeric& y,
1315  Numeric& z,
1316  const Numeric& r,
1317  const Numeric& lat,
1318  const Numeric& lon )
1319 {
1320  assert( r > 0 );
1321  assert( abs( lat ) <= 90 );
1322  assert( abs( lon ) <= 360 );
1323 
1324  const Numeric latrad = DEG2RAD * lat;
1325  const Numeric lonrad = DEG2RAD * lon;
1326 
1327  x = r * cos( latrad ); // Common term for x and z
1328  y = x * sin( lonrad );
1329  x = x * cos( lonrad );
1330  z = r * sin( latrad );
1331 }
1332 
1333 
1334 
1335 
1336 
1337 /*===========================================================================
1338  === coordinate transformations
1339  ===========================================================================*/
1340 
1342 
1359  Vector& longrid_out,
1360  ConstVectorView longrid_in,
1361  const Numeric lon )
1362 {
1363  longrid_out = longrid_in;
1364  if (longrid_in[longrid_in.nelem()-1] >= lon+360.)
1365  longrid_out += -360.;
1366  else if (longrid_in[0] <= lon-360.)
1367  longrid_out += 360.;
1368 }
1369 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
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:927
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:977
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:201
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition: geodetic.cc:356
The Vector class.
Definition: matpackI.h:556
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:302
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:183
const Numeric RAD2DEG
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1244
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Structure to store a grid position.
Definition: interpolation.h:74
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:473
const Numeric ANGTOL
Definition: ppath.h:103
Numeric fd[2]
Definition: interpolation.h:76
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
Definition: geodetic.cc:1276
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:125
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
#define abs(x)
Definition: continua.cc:20458
Index idx
Definition: interpolation.h:75
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define dx
Definition: continua.cc:21928
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:658
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
Propagation path structure and functions.
const Numeric DEG2RAD
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
Definition: geodetic.cc:1358
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:82
#define min(a, b)
Definition: continua.cc:20460
A constant view of a Vector.
Definition: matpackI.h:292
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:596
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition: geodetic.cc:1312
const Numeric POLELAT
Definition: ppath.h:94
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:699
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:387
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:865
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1139