Geodesic.hpp

Go to the documentation of this file.
00001 /**
00002  * \file Geodesic.hpp
00003  * \brief Header for GeographicLib::Geodesic class
00004  *
00005  * Copyright (c) Charles Karney (2009, 2010, 2011) <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 6950 2011-02-11 04:09:24Z 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 
00030    * The shortest path between two points on a ellipsoid at (\e lat1, \e lon1)
00031    * and (\e lat2, \e lon2) is called the geodesic.  Its length is \e s12 and
00032    * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
00033    * the two end points.  (The azimuth is the heading measured clockwise from
00034    * north.  \e azi2 is the "forward" azimuth, i.e., the heading that takes you
00035    * beyond point 2 not back to point 1.)
00036    *
00037    * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
00038    * lon2, and \e azi2.  This is the \e direct geodesic problem and its
00039    * solution is given by the function Geodesic::Direct.  (If \e s12 is
00040    * sufficiently large that the geodesic wraps more than halfway around the
00041    * earth, there will be another geodesic between the points with a smaller \e
00042    * s12.)
00043    *
00044    * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e
00045    * azi2, and \e s12.  This is the \e inverse geodesic problem, whose solution
00046    * is given by Geodesic::Inverse.  Usually, the solution to the inverse
00047    * problem is unique.  In cases where there are muliple solutions (all with
00048    * the same \e s12, of course), all the solutions can be easily generated
00049    * once a particular solution is provided.
00050    *
00051    * The standard way of specifying the direct problem is the specify the
00052    * distance \e s12 to the second point.  However it is sometimes useful
00053    * instead to specify the the arc length \e a12 (in degrees) on the auxiliary
00054    * sphere.  This is a mathematical construct used in solving the geodesic
00055    * problems.  The solution of the direct problem in this form is provide by
00056    * Geodesic::ArcDirect.  An arc length in excess of 180<sup>o</sup> indicates
00057    * that the geodesic is not a shortest path.  In addition, the arc length
00058    * between an equatorial crossing and the next extremum of latitude for a
00059    * geodesic is 90<sup>o</sup>.
00060    *
00061    * This class can also calculate several other quantities related to
00062    * geodesics.  These are:
00063    * - <i>reduced length</i>.  If we fix the first point and increase \e azi1
00064    *   by \e dazi1 (radians), the the second point is displaced \e m12 \e dazi1
00065    *   in the direction \e azi2 + 90<sup>o</sup>.  The quantity \e m12 is
00066    *   called the "reduced length" and is symmetric under interchange of the
00067    *   two points.  On a flat surface, we have \e m12 = \e s12.  The ratio \e
00068    *   s12/\e m12 gives the azimuthal scale for an azimuthal equidistant
00069    *   projection.
00070    * - <i>geodesic scale</i>.  Consider a reference geodesic and a second
00071    *   geodesic parallel to this one at point 1 and separated by a small
00072    *   distance \e dt.  The separation of the two geodesics at point 2 is \e
00073    *   M12 \e dt where \e M12 is called the "geodesic scale".  \e M21 is
00074    *   defined similarly (with the geodesics being parallel at point 2).  On a
00075    *   flat surface, we have \e M12 = \e M21 = 1.  The quantity 1/\e M12 gives
00076    *   the scale of the Cassini-Soldner projection.
00077    * - <i>area</i>.  Consider the quadrilateral bounded by the following lines:
00078    *   the geodesic from point 1 to point 2, the meridian from point 2 to the
00079    *   equator, the equator from \e lon2 to \e lon1, the meridian from the
00080    *   equator to point 1.  The area of this quadrilateral is represented by \e
00081    *   S12 with a clockwise traversal of the perimeter counting as a positive
00082    *   area and it can be used to compute the area of any simple geodesic
00083    *   polygon.
00084    *
00085    * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and
00086    * Geodesic::Inverse allow these quantities to be returned.  In addition
00087    * there are general functions Geodesic::GenDirect, and Geodesic::GenInverse
00088    * which allow an arbitrary set of results to be computed.
00089    *
00090    * Additional functionality if provided by the GeodesicLine class, which
00091    * allows a sequence of points along a geodesic to be computed.
00092    *
00093    * The calculations are accurate to better than 15 nm.  See Sec. 9 of
00094    * <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a> for details.
00095    *
00096    * The algorithms are described in
00097    * - C. F. F. Karney,
00098    *   <a href="http://arxiv.org/abs/1102.1215">Geodesics
00099    *   on an ellipsoid of revolution</a>,
00100    *   Feb. 2011;
00101    *   preprint
00102    *   <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a>.
00103    * .
00104    * For more information on geodesics see \ref geodesic.
00105    **********************************************************************/
00106 
00107   class Geodesic {
00108   private:
00109     typedef Math::real real;
00110     friend class GeodesicLine;
00111     static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD,
00112       nA2 = GEOD_ORD, nC2 = GEOD_ORD,
00113       nA3 = GEOD_ORD, nA3x = nA3,
00114       nC3 = GEOD_ORD, nC3x = (nC3 * (nC3 - 1)) / 2,
00115       nC4 = GEOD_ORD, nC4x = (nC4 * (nC4 + 1)) / 2;
00116     static const unsigned maxit = 50;
00117 
00118     static inline real sq(real x) throw() { return x * x; }
00119     void Lengths(real eps, real sig12,
00120                  real ssig1, real csig1, real ssig2, real csig2,
00121                  real cbet1, real cbet2,
00122                  real& s12s, real& m12a, real& m0,
00123                  bool scalep, real& M12, real& M21,
00124                  real tc[], real zc[]) const throw();
00125     static real Astroid(real R, real z) throw();
00126     real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2,
00127                       real lam12,
00128                       real& salp1, real& calp1,
00129                       real& salp2, real& calp2,
00130                       real C1a[], real C2a[]) const throw();
00131     real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2,
00132                   real salp1, real calp1,
00133                   real& salp2, real& calp2, real& sig12,
00134                   real& ssig1, real& csig1, real& ssig2, real& csig2,
00135                   real& eps, real& domg12, bool diffp, real& dlam12,
00136                   real C1a[], real C2a[], real C3a[])
00137       const throw();
00138 
00139     static const real eps2, tol0, tol1, tol2, xthresh;
00140     const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
00141     real _A3x[nA3x], _C3x[nC3x], _C4x[nC4x];
00142     static real SinCosSeries(bool sinp,
00143                              real sinx, real cosx, const real c[], int n)
00144       throw();
00145     static inline real AngNormalize(real x) throw() {
00146       // Place angle in [-180, 180).  Assumes x is in [-540, 540).
00147       return x >= 180 ? x - 360 : x < -180 ? x + 360 : x;
00148     }
00149     static inline real AngRound(real x) throw() {
00150       // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
00151       // for reals = 0.7 pm on the earth if x is an angle in degrees.  (This
00152       // is about 1000 times more resolution than we get with angles around 90
00153       // degrees.)  We use this to avoid having to deal with near singular
00154       // cases when x is non-zero but tiny (e.g., 1.0e-200).
00155       const real z = real(0.0625); // 1/16
00156       volatile real y = std::abs(x);
00157       // The compiler mustn't "simplify" z - (z - y) to y
00158       y = y < z ? z - (z - y) : y;
00159       return x < 0 ? -y : y;
00160     }
00161     static inline void SinCosNorm(real& sinx, real& cosx) throw() {
00162       real r = Math::hypot(sinx, cosx);
00163       sinx /= r;
00164       cosx /= r;
00165     }
00166     // These are Maxima generated functions to provide series approximations to
00167     // the integrals for the ellipsoidal geodesic.
00168     static real A1m1f(real eps) throw();
00169     static void C1f(real eps, real c[]) throw();
00170     static void C1pf(real eps, real c[]) throw();
00171     static real A2m1f(real eps) throw();
00172     static void C2f(real eps, real c[]) throw();
00173     void A3coeff() throw();
00174     real A3f(real eps) const throw();
00175     void C3coeff() throw();
00176     void C3f(real eps, real c[]) const throw();
00177     void C4coeff() throw();
00178     void C4f(real k2, real c[]) const throw();
00179 
00180     enum captype {
00181       CAP_NONE = 0U,
00182       CAP_C1   = 1U<<0,
00183       CAP_C1p  = 1U<<1,
00184       CAP_C2   = 1U<<2,
00185       CAP_C3   = 1U<<3,
00186       CAP_C4   = 1U<<4,
00187       CAP_ALL  = 0x1FU,
00188       OUT_ALL  = 0x7F80U,
00189     };
00190   public:
00191 
00192     /**
00193      * Bit masks for what calculations to do.  These masks do double duty.
00194      * They signify to the GeodesicLine::GeodesicLine constructor and to
00195      * Geodesic::Line what capabilities should be included in the GeodesicLine
00196      * object.  They also specify which results to return in the general
00197      * routines Geodesic::GenDirect and Geodesic::GenInverse routines.
00198      * GeodesicLine::mask is a duplication of this enum.
00199      **********************************************************************/
00200     enum mask {
00201       /**
00202        * No capabilities, no output.
00203        * @hideinitializer
00204        **********************************************************************/
00205       NONE          = 0U,
00206       /**
00207        * Calculate latitude \e lat2.  (It's not necessary to include this as a
00208        * capability to GeodesicLine because this is included by default.)
00209        * @hideinitializer
00210        **********************************************************************/
00211       LATITUDE      = 1U<<7  | CAP_NONE,
00212       /**
00213        * Calculate longitude \e lon2.
00214        * @hideinitializer
00215        **********************************************************************/
00216       LONGITUDE     = 1U<<8  | CAP_C3,
00217       /**
00218        * Calculate azimuths \e azi1 and \e azi2.  (It's not necessary to
00219        * include this as a capability to GeodesicLine because this is included
00220        * by default.)
00221        * @hideinitializer
00222        **********************************************************************/
00223       AZIMUTH       = 1U<<9  | CAP_NONE,
00224       /**
00225        * Calculate distance \e s12.
00226        * @hideinitializer
00227        **********************************************************************/
00228       DISTANCE      = 1U<<10 | CAP_C1,
00229       /**
00230        * Allow distance \e s12 to be used as input in the direct geodesic
00231        * problem.
00232        * @hideinitializer
00233        **********************************************************************/
00234       DISTANCE_IN   = 1U<<11 | CAP_C1 | CAP_C1p,
00235       /**
00236        * Calculate reduced length \e m12.
00237        * @hideinitializer
00238        **********************************************************************/
00239       REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
00240       /**
00241        * Calculate geodesic scales \e M12 and \e M21.
00242        * @hideinitializer
00243        **********************************************************************/
00244       GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
00245       /**
00246        * Calculate area \e S12.
00247        * @hideinitializer
00248        **********************************************************************/
00249       AREA          = 1U<<14 | CAP_C4,
00250       /**
00251        * All capabilities.  Calculate everything.
00252        * @hideinitializer
00253        **********************************************************************/
00254       ALL           = OUT_ALL| CAP_ALL,
00255     };
00256 
00257     /** \name Constructor
00258      **********************************************************************/
00259     ///@{
00260     /**
00261      * Constructor for a ellipsoid with
00262      *
00263      * @param[in] a equatorial radius (meters)
00264      * @param[in] r reciprocal flattening.  Setting \e r = 0 implies \e r = inf
00265      *   or flattening = 0 (i.e., a sphere).  Negative \e r indicates a prolate
00266      *   ellipsoid.
00267      *
00268      * An exception is thrown if either of the axes of the ellipsoid is
00269      * non-positive.
00270      **********************************************************************/
00271     Geodesic(real a, real r);
00272     ///@}
00273 
00274     /** \name Direct geodesic problem specified in terms of distance.
00275      **********************************************************************/
00276     ///@{
00277     /**
00278      * Perform the direct geodesic calculation where the length of the geodesic
00279      * is specify in terms of distance.
00280      *
00281      * @param[in] lat1 latitude of point 1 (degrees).
00282      * @param[in] lon1 longitude of point 1 (degrees).
00283      * @param[in] azi1 azimuth at point 1 (degrees).
00284      * @param[in] s12 distance between point 1 and point 2 (meters); it can be
00285      *   signed.
00286      * @param[out] lat2 latitude of point 2 (degrees).
00287      * @param[out] lon2 longitude of point 2 (degrees).
00288      * @param[out] azi2 (forward) azimuth at point 2 (degrees).
00289      * @param[out] m12 reduced length of geodesic (meters).
00290      * @param[out] M12 geodesic scale of point 2 relative to point 1
00291      *   (dimensionless).
00292      * @param[out] M21 geodesic scale of point 1 relative to point 2
00293      *   (dimensionless).
00294      * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
00295      * @return \e a12 arc length of between point 1 and point 2 (degrees).
00296      *
00297      * If either point is at a pole, the azimuth is defined by keeping the
00298      * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
00299      * taking the limit \e eps -> 0 from above.  An arc length greater that 180
00300      * degrees signifies a geodesic which is not a shortest path.  (For a
00301      * prolate ellipsoid, an additional condition is necessary for a shortest
00302      * path: the longitudinal extent must not exceed of 180 degrees.)
00303      *
00304      * The following functions are overloaded versions of Geodesic::Direct
00305      * which omit some of the output parameters.  Note, however, that the arc
00306      * length is always computed and returned as the function value.
00307      **********************************************************************/
00308     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00309                       real& lat2, real& lon2, real& azi2,
00310                       real& m12, real& M12, real& M21, real& S12)
00311       const throw() {
00312       real t;
00313       return GenDirect(lat1, lon1, azi1, false, s12,
00314                        LATITUDE | LONGITUDE | AZIMUTH |
00315                        REDUCEDLENGTH | GEODESICSCALE | AREA,
00316                        lat2, lon2, azi2, t, m12, M12, M21, S12);
00317     }
00318 
00319     /**
00320      * See the documentation for Geodesic::Direct.
00321      **********************************************************************/
00322     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00323                       real& lat2, real& lon2)
00324       const throw() {
00325       real t;
00326       return GenDirect(lat1, lon1, azi1, false, s12,
00327                        LATITUDE | LONGITUDE,
00328                        lat2, lon2, t, t, t, t, t, t);
00329     }
00330 
00331     /**
00332      * See the documentation for Geodesic::Direct.
00333      **********************************************************************/
00334     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00335                       real& lat2, real& lon2, real& azi2)
00336       const throw() {
00337       real t;
00338       return GenDirect(lat1, lon1, azi1, false, s12,
00339                        LATITUDE | LONGITUDE | AZIMUTH,
00340                        lat2, lon2, azi2, t, t, t, t, t);
00341     }
00342 
00343     /**
00344      * See the documentation for Geodesic::Direct.
00345      **********************************************************************/
00346     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00347                       real& lat2, real& lon2, real& azi2, real& m12)
00348       const throw() {
00349       real t;
00350       return GenDirect(lat1, lon1, azi1, false, s12,
00351                        LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
00352                        lat2, lon2, azi2, t, m12, t, t, t);
00353     }
00354 
00355     /**
00356      * See the documentation for Geodesic::Direct.
00357      **********************************************************************/
00358     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00359                       real& lat2, real& lon2, real& azi2,
00360                       real& M12, real& M21)
00361       const throw() {
00362       real t;
00363       return GenDirect(lat1, lon1, azi1, false, s12,
00364                        LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
00365                        lat2, lon2, azi2, t, t, M12, M21, t);
00366     }
00367 
00368     /**
00369      * See the documentation for Geodesic::Direct.
00370      **********************************************************************/
00371     Math::real Direct(real lat1, real lon1, real azi1, real s12,
00372                       real& lat2, real& lon2, real& azi2,
00373                       real& m12, real& M12, real& M21)
00374       const throw() {
00375       real t;
00376       return GenDirect(lat1, lon1, azi1, false, s12,
00377                        LATITUDE | LONGITUDE | AZIMUTH |
00378                        REDUCEDLENGTH | GEODESICSCALE,
00379                        lat2, lon2, azi2, t, m12, M12, M21, t);
00380     }
00381     ///@}
00382 
00383     /** \name Direct geodesic problem specified in terms of arc length.
00384      **********************************************************************/
00385     ///@{
00386     /**
00387      * Perform the direct geodesic calculation where the length of the geodesic
00388      * is specify in terms of arc length.
00389      *
00390      * @param[in] lat1 latitude of point 1 (degrees).
00391      * @param[in] lon1 longitude of point 1 (degrees).
00392      * @param[in] azi1 azimuth at point 1 (degrees).
00393      * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
00394      *   be signed.
00395      * @param[out] lat2 latitude of point 2 (degrees).
00396      * @param[out] lon2 longitude of point 2 (degrees).
00397      * @param[out] azi2 (forward) azimuth at point 2 (degrees).
00398      * @param[out] s12 distance between point 1 and point 2 (meters).
00399      * @param[out] m12 reduced length of geodesic (meters).
00400      * @param[out] M12 geodesic scale of point 2 relative to point 1
00401      *   (dimensionless).
00402      * @param[out] M21 geodesic scale of point 1 relative to point 2
00403      *   (dimensionless).
00404      * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
00405      *
00406      * If either point is at a pole, the azimuth is defined by keeping the
00407      * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
00408      * taking the limit \e eps -> 0 from above.  An arc length greater that 180
00409      * degrees signifies a geodesic which is not a shortest path.  (For a
00410      * prolate ellipsoid, an additional condition is necessary for a shortest
00411      * path: the longitudinal extent must not exceed of 180 degrees.)
00412      *
00413      * The following functions are overloaded versions of Geodesic::Direct
00414      * which omit some of the output parameters.
00415      **********************************************************************/
00416     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00417                    real& lat2, real& lon2, real& azi2, real& s12,
00418                    real& m12, real& M12, real& M21, real& S12)
00419       const throw() {
00420       GenDirect(lat1, lon1, azi1, true, a12,
00421                 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
00422                 REDUCEDLENGTH | GEODESICSCALE | AREA,
00423                 lat2, lon2, azi2, s12, m12, M12, M21, S12);
00424     }
00425 
00426     /**
00427      * See the documentation for Geodesic::ArcDirect.
00428      **********************************************************************/
00429     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00430                    real& lat2, real& lon2) const throw() {
00431       real t;
00432       GenDirect(lat1, lon1, azi1, true, a12,
00433                 LATITUDE | LONGITUDE,
00434                 lat2, lon2, t, t, t, t, t, t);
00435     }
00436 
00437     /**
00438      * See the documentation for Geodesic::ArcDirect.
00439      **********************************************************************/
00440     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00441                    real& lat2, real& lon2, real& azi2) const throw() {
00442       real t;
00443       GenDirect(lat1, lon1, azi1, true, a12,
00444                 LATITUDE | LONGITUDE | AZIMUTH,
00445                 lat2, lon2, azi2, t, t, t, t, t);
00446     }
00447 
00448     /**
00449      * See the documentation for Geodesic::ArcDirect.
00450      **********************************************************************/
00451     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00452                    real& lat2, real& lon2, real& azi2, real& s12)
00453       const throw() {
00454       real t;
00455       GenDirect(lat1, lon1, azi1, true, a12,
00456                 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
00457                 lat2, lon2, azi2, s12, t, t, t, t);
00458     }
00459 
00460     /**
00461      * See the documentation for Geodesic::ArcDirect.
00462      **********************************************************************/
00463     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00464                    real& lat2, real& lon2, real& azi2,
00465                    real& s12, real& m12) const throw() {
00466       real t;
00467       GenDirect(lat1, lon1, azi1, true, a12,
00468                 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
00469                 REDUCEDLENGTH,
00470                 lat2, lon2, azi2, s12, m12, t, t, t);
00471     }
00472 
00473     /**
00474      * See the documentation for Geodesic::ArcDirect.
00475      **********************************************************************/
00476     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00477                    real& lat2, real& lon2, real& azi2, real& s12,
00478                    real& M12, real& M21) const throw() {
00479       real t;
00480       GenDirect(lat1, lon1, azi1, true, a12,
00481                 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
00482                 GEODESICSCALE,
00483                 lat2, lon2, azi2, s12, t, M12, M21, t);
00484     }
00485 
00486     /**
00487      * See the documentation for Geodesic::ArcDirect.
00488      **********************************************************************/
00489     void ArcDirect(real lat1, real lon1, real azi1, real a12,
00490                    real& lat2, real& lon2, real& azi2, real& s12,
00491                    real& m12, real& M12, real& M21) const throw() {
00492       real t;
00493       GenDirect(lat1, lon1, azi1, true, a12,
00494                 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
00495                 REDUCEDLENGTH | GEODESICSCALE,
00496                 lat2, lon2, azi2, s12, m12, M12, M21, t);
00497     }
00498     ///@}
00499 
00500     /** \name General version of the direct geodesic solution.
00501      **********************************************************************/
00502     ///@{
00503 
00504     /**
00505      * The general direct geodesic calculation.  Geodesic::Direct and
00506      * Geodesic::ArcDirect are defined in terms of this function.
00507      *
00508      * @param[in] lat1 latitude of point 1 (degrees).
00509      * @param[in] lon1 longitude of point 1 (degrees).
00510      * @param[in] azi1 azimuth at point 1 (degrees).
00511      * @param[in] arcmode boolean flag determining the meaning of the second
00512      *   parameter.
00513      * @param[in] s12_a12 if \e arcmode is false, this is the distance between
00514      *   point 1 and point 2 (meters); otherwise it is the arc length between
00515      *   point 1 and point 2 (degrees); it can be signed.
00516      * @param[in] outmask a bitor'ed combination of Geodesic::mask values
00517      *   specifying which of the following parameters should be set.
00518      * @param[out] lat2 latitude of point 2 (degrees).
00519      * @param[out] lon2 longitude of point 2 (degrees).
00520      * @param[out] azi2 (forward) azimuth at point 2 (degrees).
00521      * @param[out] s12 distance between point 1 and point 2 (meters).
00522      * @param[out] m12 reduced length of geodesic (meters).
00523      * @param[out] M12 geodesic scale of point 2 relative to point 1
00524      *   (dimensionless).
00525      * @param[out] M21 geodesic scale of point 1 relative to point 2
00526      *   (dimensionless).
00527      * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
00528      * @return \e a12 arc length of between point 1 and point 2 (degrees).
00529      *
00530      * The Geodesic::mask values possible for \e outmask are
00531      * - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2.
00532      * - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2.
00533      * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2.
00534      * - \e outmask |= Geodesic::DISTANCE for the distance \e s12.
00535      * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
00536      *   m12.
00537      * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
00538      *   M12 and \e M21.
00539      * - \e outmask |= Geodesic::AREA for the area \e S12.
00540      * .
00541      * The function value \e a12 is always computed and returned and this
00542      * equals \e s12_a12 is \e arcmode is true.  If \e outmask includes
00543      * Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12.
00544      * It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this
00545      * is automatically included is \e arcmode is false.
00546      **********************************************************************/
00547     Math::real GenDirect(real lat1, real lon1, real azi1,
00548                          bool arcmode, real s12_a12, unsigned outmask,
00549                          real& lat2, real& lon2, real& azi2,
00550                          real& s12, real& m12, real& M12, real& M21,
00551                          real& S12) const throw();
00552     ///@}
00553 
00554     /** \name Inverse geodesic problem.
00555      **********************************************************************/
00556     ///@{
00557     /**
00558      * Perform the inverse geodesic calculation.
00559      *
00560      * @param[in] lat1 latitude of point 1 (degrees).
00561      * @param[in] lon1 longitude of point 1 (degrees).
00562      * @param[in] lat2 latitude of point 2 (degrees).
00563      * @param[in] lon2 longitude of point 2 (degrees).
00564      * @param[out] s12 distance between point 1 and point 2 (meters).
00565      * @param[out] azi1 azimuth at point 1 (degrees).
00566      * @param[out] azi2 (forward) azimuth at point 1 (degrees).
00567      * @param[out] m12 reduced length of geodesic (meters).
00568      * @param[out] M12 geodesic scale of point 2 relative to point 1
00569      *   (dimensionless).
00570      * @param[out] M21 geodesic scale of point 1 relative to point 2
00571      *   (dimensionless).
00572      * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
00573      * @return \e a12 arc length of between point 1 and point 2 (degrees).
00574      *
00575      * If either point is at a pole, the azimuth is defined by keeping the
00576      * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
00577      * taking the limit \e eps -> 0 from above.  If the routine fails to
00578      * converge, then all the requested outputs are set to Math::NaN().  This
00579      * is not expected to happen with ellipsoidal models of the earth; please
00580      * report all cases where this occurs.
00581      *
00582      * The following functions are overloaded versions of Geodesic::Inverse
00583      * which omit some of the output parameters.  Note, however, that the arc
00584      * length is always computed and returned as the function value.
00585      **********************************************************************/
00586     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00587                        real& s12, real& azi1, real& azi2, real& m12,
00588                        real& M12, real& M21, real& S12) const throw() {
00589       return GenInverse(lat1, lon1, lat2, lon2,
00590                         DISTANCE | AZIMUTH |
00591                         REDUCEDLENGTH | GEODESICSCALE | AREA,
00592                         s12, azi1, azi2, m12, M12, M21, S12);
00593     }
00594 
00595     /**
00596      * See the documentation for Geodesic::Inverse.
00597      **********************************************************************/
00598     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00599                        real& s12) const throw() {
00600       real t;
00601       return GenInverse(lat1, lon1, lat2, lon2,
00602                         DISTANCE,
00603                         s12, t, t, t, t, t, t);
00604     }
00605 
00606     /**
00607      * See the documentation for Geodesic::Inverse.
00608      **********************************************************************/
00609     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00610                        real& azi1, real& azi2) const throw() {
00611       real t;
00612       return GenInverse(lat1, lon1, lat2, lon2,
00613                         AZIMUTH,
00614                         t, azi1, azi2, t, t, t, t);
00615     }
00616 
00617     /**
00618      * See the documentation for Geodesic::Inverse.
00619      **********************************************************************/
00620     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00621                        real& s12, real& azi1, real& azi2)
00622       const throw() {
00623       real t;
00624       return GenInverse(lat1, lon1, lat2, lon2,
00625                         DISTANCE | AZIMUTH,
00626                         s12, azi1, azi2, t, t, t, t);
00627     }
00628 
00629     /**
00630      * See the documentation for Geodesic::Inverse.
00631      **********************************************************************/
00632     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00633                        real& s12, real& azi1, real& azi2, real& m12)
00634       const throw() {
00635       real t;
00636       return GenInverse(lat1, lon1, lat2, lon2,
00637                         DISTANCE | AZIMUTH | REDUCEDLENGTH,
00638                         s12, azi1, azi2, m12, t, t, t);
00639     }
00640 
00641     /**
00642      * See the documentation for Geodesic::Inverse.
00643      **********************************************************************/
00644     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00645                        real& s12, real& azi1, real& azi2,
00646                        real& M12, real& M21) const throw() {
00647       real t;
00648       return GenInverse(lat1, lon1, lat2, lon2,
00649                         DISTANCE | AZIMUTH | GEODESICSCALE,
00650                         s12, azi1, azi2, t, M12, M21, t);
00651     }
00652 
00653     /**
00654      * See the documentation for Geodesic::Inverse.
00655      **********************************************************************/
00656     Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
00657                        real& s12, real& azi1, real& azi2, real& m12,
00658                        real& M12, real& M21) const throw() {
00659       real t;
00660       return GenInverse(lat1, lon1, lat2, lon2,
00661                         DISTANCE | AZIMUTH |
00662                         REDUCEDLENGTH | GEODESICSCALE,
00663                         s12, azi1, azi2, m12, M12, M21, t);
00664     }
00665     ///@}
00666 
00667     /** \name General version of inverse geodesic solution.
00668      **********************************************************************/
00669     ///@{
00670     /**
00671      * The general inverse geodesic calculation.  Geodesic::Inverse is defined
00672      * in terms of this function.
00673      *
00674      * @param[in] lat1 latitude of point 1 (degrees).
00675      * @param[in] lon1 longitude of point 1 (degrees).
00676      * @param[in] lat2 latitude of point 2 (degrees).
00677      * @param[in] lon2 longitude of point 2 (degrees).
00678      * @param[in] outmask a bitor'ed combination of Geodesic::mask values
00679      *   specifying which of the following parameters should be set.
00680      * @param[out] s12 distance between point 1 and point 2 (meters).
00681      * @param[out] azi1 azimuth at point 1 (degrees).
00682      * @param[out] azi2 (forward) azimuth at point 1 (degrees).
00683      * @param[out] m12 reduced length of geodesic (meters).
00684      * @param[out] M12 geodesic scale of point 2 relative to point 1
00685      *   (dimensionless).
00686      * @param[out] M21 geodesic scale of point 1 relative to point 2
00687      *   (dimensionless).
00688      * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
00689      * @return \e a12 arc length of between point 1 and point 2 (degrees).
00690      *
00691      * The Geodesic::mask values possible for \e outmask are
00692      * - \e outmask |= Geodesic::DISTANCE for the distance \e s12.
00693      * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2.
00694      * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
00695      *   m12.
00696      * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
00697      *   M12 and \e M21.
00698      * - \e outmask |= Geodesic::AREA for the area \e S12.
00699      * .
00700      * The arc length is always computed and returned as the function value.
00701      **********************************************************************/
00702     Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
00703                        unsigned outmask,
00704                        real& s12, real& azi1, real& azi2,
00705                        real& m12, real& M12, real& M21, real& S12)
00706       const throw();
00707     ///@}
00708 
00709     /** \name Interface to GeodesicLine.
00710      **********************************************************************/
00711     ///@{
00712 
00713     /**
00714      * Set up to do a series of ranges.
00715      *
00716      * @param[in] lat1 latitude of point 1 (degrees).
00717      * @param[in] lon1 longitude of point 1 (degrees).
00718      * @param[in] azi1 azimuth at point 1 (degrees).
00719      * @param[in] caps bitor'ed combination of Geodesic::mask values
00720      *   specifying the capabilities the GeodesicLine object should possess,
00721      *   i.e., which quantities can be returned in calls to
00722      *   GeodesicLib::Position.
00723      *
00724      * The Geodesic::mask values are
00725      * - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
00726      *   added automatically
00727      * - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2
00728      * - \e caps |= Geodesic::AZIMUTH for the latitude \e azi2; this is
00729      *   added automatically
00730      * - \e caps |= Geodesic::DISTANCE for the distance \e s12
00731      * - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12
00732      * - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
00733      *   and \e M21
00734      * - \e caps |= Geodesic::AREA for the area \e S12
00735      * - \e caps |= Geodesic::DISTANCE_IN permits the length of the
00736      *   geodesic to be given in terms of \e s12; without this capability the
00737      *   length can only be specified in terms of arc length.
00738      * .
00739      * The default value of \e caps is Geodesic::ALL which turns on all the
00740      * capabilities.
00741      *
00742      * If the point is at a pole, the azimuth is defined by keeping the \e lon1
00743      * fixed and writing \e lat1 = 90 - \e eps or -90 + \e eps and taking the
00744      * limit \e eps -> 0 from above.
00745      **********************************************************************/
00746     GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
00747       const throw();
00748 
00749     ///@}
00750 
00751     /** \name Inspector functions.
00752      **********************************************************************/
00753     ///@{
00754 
00755     /**
00756      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00757      *   the value used in the constructor.
00758      **********************************************************************/
00759     Math::real MajorRadius() const throw() { return _a; }
00760 
00761     /**
00762      * @return \e r the inverse flattening of the ellipsoid.  This is the
00763      *   value used in the constructor.  A value of 0 is returned for a sphere
00764      *   (infinite inverse flattening).
00765      **********************************************************************/
00766     Math::real InverseFlattening() const throw() { return _r; }
00767 
00768     /**
00769      * @return total area of ellipsoid in meters<sup>2</sup>.  The area of a
00770      *   polygon encircling a pole can be found by adding
00771      *   Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
00772      *   polygon.
00773      **********************************************************************/
00774     Math::real EllipsoidArea() const throw()
00775     { return 4 * Math::pi<real>() * _c2; }
00776     ///@}
00777 
00778     /**
00779      * A global instantiation of Geodesic with the parameters for the WGS84
00780      * ellipsoid.
00781      **********************************************************************/
00782     static const Geodesic WGS84;
00783 
00784 
00785     /** \name Deprecated function.
00786      **********************************************************************/
00787     ///@{
00788     /**
00789      * <b>DEPRECATED</b> Perform the direct geodesic calculation.  Given a
00790      * latitude, \e lat1, longitude, \e lon1, and azimuth \e azi1 (degrees) for
00791      * point 1 and a range, \e s12 (meters) from point 1 to point 2, return the
00792      * latitude, \e lat2, longitude, \e lon2, and forward azimuth, \e azi2
00793      * (degrees) for point 2 and the reduced length \e m12 (meters).  If either
00794      * point is at a pole, the azimuth is defined by keeping the longitude
00795      * fixed and writing \e lat = 90 - \e eps or -90 + \e eps and taking the
00796      * limit \e eps -> 0 from above.  If \e arcmode (default false) is set to
00797      * true, \e s12 is interpreted as the arc length \e a12 (degrees) on the
00798      * auxiliary sphere.  An arc length greater that 180 degrees results in a
00799      * geodesic which is not a shortest path.  For a prolate ellipsoid, an
00800      * additional condition is necessary for a shortest path: the longitudinal
00801      * extent must not exceed of 180 degrees.  Returned value is the arc length
00802      * \e a12 (degrees) if \e arcmode is false, otherwise it is the distance \e
00803      * s12 (meters).
00804      **********************************************************************/
00805     Math::real Direct(real lat1, real lon1, real azi1, real s12_a12,
00806                       real& lat2, real& lon2, real& azi2, real& m12,
00807                       bool arcmode) const throw() {
00808       if (arcmode) {
00809         real a12 = s12_a12, s12;
00810         ArcDirect(lat1, lon1, azi1, a12, lat2, lon2, azi2, s12, m12);
00811         return s12;
00812       } else {
00813         real s12 = s12_a12;
00814         return Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2, m12);
00815       }
00816     }
00817     ///@}
00818   };
00819 
00820 } // namespace GeographicLib
00821 #endif