00001 /** 00002 * \file Geodesic.hpp 00003 * \brief Header for GeographicLib::Geodesic class 00004 * 00005 * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> 00006 * and licensed under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_GEODESIC_HPP) 00011 #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6827 2010-05-20 19:56:18Z karney $" 00012 00013 #include "GeographicLib/Constants.hpp" 00014 00015 #if !defined(GEOD_ORD) 00016 /** 00017 * The order of the expansions used by Geodesic. 00018 **********************************************************************/ 00019 #define GEOD_ORD (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 3 : 7) 00020 #endif 00021 00022 namespace GeographicLib { 00023 00024 class GeodesicLine; 00025 00026 /** 00027 * \brief %Geodesic calculations 00028 * 00029 * The shortest path between two points on a ellipsoid at (\e lat1, \e lon1) 00030 * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and 00031 * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at 00032 * the two end points. (The azimuth is the heading measured clockwise from 00033 * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you 00034 * beyond point 2 not back to point 1.) 00035 * 00036 * If we fix the first point and increase \e s12 by \e ds12, then the second 00037 * point is displaced \e ds12 in the direction \e azi2. Similarly we 00038 * increase \e azi1 by \e dazi1 (radians), the the second point is displaced 00039 * \e m12 \e dazi1 in the direction \e azi2 + 90<sup>o</sup>. The quantity 00040 * \e m12 is called the "reduced length" and is symmetric under interchange 00041 * of the two points. On a flat surface, he have \e m12 = \e s12. The ratio 00042 * \e s12/\e m12 gives the azimuthal scale for an azimuthal equidistant 00043 * projection. 00044 * 00045 * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e 00046 * lon2, \e azi2, \e m12. This is the \e direct geodesic problem. (If \e 00047 * s12 is sufficiently large that the geodesic wraps more than halfway around 00048 * the earth, there will be another geodesic between the points with a 00049 * smaller \e s12.) 00050 * 00051 * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e 00052 * azi2, \e s12, \e m12. This is the \e inverse geodesic problem. Usually, 00053 * the solution to the inverse problem is unique. In cases where there are 00054 * muliple solutions (all with the same \e s12, of course), all the solutions 00055 * can be easily generated once a particular solution is provided. 00056 * 00057 * As an alternative to using distance to measure \e s12, the class can also 00058 * use the arc length \e a12 (in degrees) on the auxiliary sphere. This is a 00059 * mathematical construct used in solving the geodesic problems. However, an 00060 * arc length in excess of 180<sup>o</sup> indicates that the geodesic is not 00061 * a shortest path. In addition, the arc length between an equatorial 00062 * crossing and the next extremum of latitude for a geodesic is 00063 * 90<sup>o</sup>. 00064 * 00065 * The calculations are accurate to better than 15 nm. (See \ref geoderrors 00066 * for details.) 00067 * 00068 * For more information on geodesics see \ref geodesic. 00069 **********************************************************************/ 00070 00071 class Geodesic { 00072 private: 00073 typedef Math::real real; 00074 friend class GeodesicLine; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-54; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-54; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-54; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-54; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-54; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087