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