/**
 * \file geodesic.h
 * \brief Header for the geodesic routines in C
 *
 * This an implementation in C of the geodesic algorithms described in
 * - C. F. F. Karney,
 *   <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
 *   Algorithms for geodesics</a>,
 *   J. Geodesy <b>87</b>, 43--55 (2013);
 *   DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
 *   10.1007/s00190-012-0578-z</a>;
 *   addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
 *   geod-addenda.html</a>.
 * .
 * The principal advantages of these algorithms over previous ones (e.g.,
 * Vincenty, 1975) are
 * - accurate to round off for |<i>f</i>| &lt; 1/50;
 * - the solution of the inverse problem is always found;
 * - differential and integral properties of geodesics are computed.
 *
 * The shortest path between two points on the ellipsoid at (\e lat1, \e
 * lon1) and (\e lat2, \e lon2) is called the geodesic.  Its length is
 * \e s12 and the geodesic from point 1 to point 2 has forward azimuths
 * \e azi1 and \e azi2 at the two end points.
 *
 * Traditionally two geodesic problems are considered:
 * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
 *   determine \e lat2, \e lon2, and \e azi2.  This is solved by the function
 *   geod_direct().
 * - the inverse problem -- given \e lat1, \e lon1, and \e lat2, \e lon2,
 *   determine \e s12, \e azi1, and \e azi2.  This is solved by the function
 *   geod_inverse().
 *
 * The ellipsoid is specified by its equatorial radius \e a (typically in
 * meters) and flattening \e f.  The routines are accurate to round off with
 * double precision arithmetic provided that |<i>f</i>| &lt; 1/50; for the
 * WGS84 ellipsoid, the errors are less than 15 nanometers.  (Reasonably
 * accurate results are obtained for |<i>f</i>| &lt; 1/5.)  For a prolate
 * ellipsoid, specify \e f &lt; 0.
 *
 * The routines also calculate several other quantities of interest
 * - \e S12 is the area between the geodesic from point 1 to point 2 and the
 *   equator; i.e., it is the area, measured counter-clockwise, of the
 *   quadrilateral with corners (\e lat1,\e lon1), (0,\e lon1), (0,\e lon2),
 *   and (\e lat2,\e lon2).
 * - \e m12, the reduced length of the geodesic is defined such that if
 *   the initial azimuth is perturbed by \e dazi1 (radians) then the
 *   second point is displaced by \e m12 \e dazi1 in the direction
 *   perpendicular to the geodesic.  On a curved surface the reduced
 *   length obeys a symmetry relation, \e m12 + \e m21 = 0.  On a flat
 *   surface, we have \e m12 = \e s12.
 * - \e M12 and \e M21 are geodesic scales.  If two geodesics are
 *   parallel at point 1 and separated by a small distance \e dt, then
 *   they are separated by a distance \e M12 \e dt at point 2.  \e M21
 *   is defined similarly (with the geodesics being parallel to one
 *   another at point 2).  On a flat surface, we have \e M12 = \e M21
 *   = 1.
 * - \e a12 is the arc length on the auxiliary sphere.  This is a
 *   construct for converting the problem to one in spherical
 *   trigonometry.  \e a12 is measured in degrees.  The spherical arc
 *   length from one equator crossing to the next is always 180&deg;.
 *
 * If points 1, 2, and 3 lie on a single geodesic, then the following
 * addition rules hold:
 * - \e s13 = \e s12 + \e s23
 * - \e a13 = \e a12 + \e a23
 * - \e S13 = \e S12 + \e S23
 * - \e m13 = \e m12 \e M23 + \e m23 \e M21
 * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e
 *   m23 / \e m12
 * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e
 *   m12 / \e m23
 *
 * The shortest distance returned by the solution of the inverse problem is
 * (obviously) uniquely defined.  However, in a few special cases there are
 * multiple azimuths which yield the same shortest distance.  Here is a
 * catalog of those cases:
 * - \e lat1 = &minus;\e lat2 (with neither at a pole).  If \e azi1 = \e
 *   azi2, the geodesic is unique.  Otherwise there are two geodesics
 *   and the second one is obtained by setting [\e azi1, \e azi2] = [\e
 *   azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 =
 *   &minus;\e S12.  (This occurs when the longitude difference is near
 *   &plusmn;180&deg; for oblate ellipsoids.)
 * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither at a pole).  If
 *   \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
 *   Otherwise there are two geodesics and the second one is obtained by
 *   setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e
 *   S12 = &minus;\e S12.  (This occurs when the \e lat2 is near
 *   &minus;\e lat1 for prolate ellipsoids.)
 * - Points 1 and 2 at opposite poles.  There are infinitely many
 *   geodesics which can be generated by setting [\e azi1, \e azi2] =
 *   [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e d.  (For
 *   spheres, this prescription applies when points 1 and 2 are
 *   antipodal.)
 * - \e s12 = 0 (coincident points).  There are infinitely many geodesics
 *   which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e
 *   azi2] + [\e d, \e d], for arbitrary \e d.
 *
 * These routines are a simple transcription of the corresponding C++ classes
 * in <a href="http://geographiclib.sf.net"> GeographicLib</a>.  The "class
 * data" is represented by the structs geod_geodesic, geod_geodesicline,
 * geod_polygon and pointers to these objects are passed as initial arguments
 * to the member functions.  Most of the internal comments have been retained.
 * However, in the process of transcription some documentation has been lost
 * and the documentation for the C++ classes, GeographicLib::Geodesic,
 * GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be
 * consulted.  The C++ code remains the "reference implementation".  Think
 * twice about restructuring the internals of the C code since this may make
 * porting fixes from the C++ code more difficult.
 *
 * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
 * under the MIT/X11 License.  For more information, see
 * http://geographiclib.sourceforge.net/
 *
 * This library was distributed with
 * <a href="../index.html">GeographicLib</a> 1.32.
 **********************************************************************/

#if !defined(GEODESIC_H)
#define GEODESIC_H 1

/**
 * The major version of the geodesic library.  (This tracks the version of
 * GeographicLib.)
 **********************************************************************/
#define GEODESIC_VERSION_MAJOR 1
/**
 * The minor version of the geodesic library.  (This tracks the version of
 * GeographicLib.)
 **********************************************************************/
#define GEODESIC_VERSION_MINOR 32
/**
 * The patch level of the geodesic library.  (This tracks the version of
 * GeographicLib.)
 **********************************************************************/
#define GEODESIC_VERSION_PATCH 0

#if defined(__cplusplus)
extern "C" {
#endif

//------------------
//WGS84 parameters
#define A_ 6378137
#define B_ 6356752.314245
#define F_ 0.00335281066
//-------------------
    /**
     * The struct containing information about the ellipsoid.  This must be
     * initialized by geod_init() before use.
     **********************************************************************/
    struct geod_geodesic {
        double a;                   /**< the equatorial radius */
        double f;                   /**< the flattening */
        /**< @cond SKIP */
        double f1, e2, ep2, n, b, c2, etol2;
        double A3x[6], C3x[15], C4x[21];
        /**< @endcond */
    };
    
    /**
     * The struct containing information about a single geodesic.  This must be
     * initialized by geod_lineinit() before use.
     **********************************************************************/
    struct geod_geodesicline {
        double lat1;                /**< the starting latitude */
        double lon1;                /**< the starting longitude */
        double azi1;                /**< the starting azimuth */
        double a;                   /**< the equatorial radius */
        double f;                   /**< the flattening */
        /**< @cond SKIP */
        double b, c2, f1, salp0, calp0, k2,
        salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1,
        A1m1, A2m1, A3c, B11, B21, B31, A4, B41;
        double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6];
        /**< @endcond */
        unsigned caps;              /**< the capabilities */
    };
    
    /**
     * The struct for accumulating information about a geodesic polygon.  This is
     * used for computing the perimeter and area of a polygon.  This must be
     * initialized by geod_polygon_init() before use.
     **********************************************************************/
    struct geod_polygon {
        double lat;                 /**< the current latitude */
        double lon;                 /**< the current longitude */
        /**< @cond SKIP */
        double lat0;
        double lon0;
        double A[2];
        double P[2];
        int polyline;
        int crossings;
        /**< @endcond */
        unsigned num;               /**< the number of points so far */
    };
    
    /**
     * Initialize a geod_geodesic object.
     *
     * @param[out] g a pointer to the object to be initialized.
     * @param[in] a the equatorial radius (meters).
     * @param[in] f the flattening.
     **********************************************************************/
    void geod_init(struct geod_geodesic* g, double a, double f);
    
    /**
     * Initialize a geod_geodesicline object.
     *
     * @param[out] l a pointer to the object to be initialized.
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lat1 latitude of point 1 (degrees).
     * @param[in] lon1 longitude of point 1 (degrees).
     * @param[in] azi1 azimuth at point 1 (degrees).
     * @param[in] caps bitor'ed combination of geod_mask() values specifying the
     *   capabilities the geod_geodesicline object should possess, i.e., which
     *   quantities can be returned in calls to geod_position() and
     *   geod_genposition().
     *
     * \e g must have been initialized with a call to geod_init().  \e lat1
     * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
     * should be in the range [&minus;540&deg;, 540&deg;).
     *
     * The geod_mask values are [see geod_mask()]:
     * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
     *   added automatically,
     * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2,
     * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is
     *   added automatically,
     * - \e caps |= GEOD_DISTANCE for the distance \e s12,
     * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12,
     * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12
     *   and \e M21,
     * - \e caps |= GEOD_AREA for the area \e S12,
     * - \e caps |= GEOD_DISTANCE_IN permits the length of the
     *   geodesic to be given in terms of \e s12; without this capability the
     *   length can only be specified in terms of arc length.
     * .
     * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE |
     * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard"
     * direct problem).
     **********************************************************************/
    void geod_lineinit(struct geod_geodesicline* l,
                       const struct geod_geodesic* g,
                       double lat1, double lon1, double azi1, unsigned caps);
    
    /**
     * Solve the direct geodesic problem.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lat1 latitude of point 1 (degrees).
     * @param[in] lon1 longitude of point 1 (degrees).
     * @param[in] azi1 azimuth at point 1 (degrees).
     * @param[in] s12 distance between point 1 and point 2 (meters); it can be
     *   negative.
     * @param[out] plat2 pointer to the latitude of point 2 (degrees).
     * @param[out] plon2 pointer to the longitude of point 2 (degrees).
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     *
     * \e g must have been initialized with a call to geod_init().  \e lat1
     * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
     * should be in the range [&minus;540&deg;, 540&deg;).  The values of \e lon2
     * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).  Any of
     * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
     * need some quantities computed.
     *
     * If either point is at a pole, the azimuth is defined by keeping the
     * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
     * taking the limit &epsilon; &rarr; 0+.  An arc length greater that 180&deg;
     * signifies a geodesic which is not a shortest path.  (For a prolate
     * ellipsoid, an additional condition is necessary for a shortest path: the
     * longitudinal extent must not exceed of 180&deg;.)
     *
     * Example, determine the point 10000 km NE of JFK:
     @code
     struct geod_geodesic g;
     double lat, lon;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0);
     printf("%.5f %.5f\n", lat, lon);
     @endcode
     **********************************************************************/
    void geod_direct(const struct geod_geodesic* g,
                     double lat1, double lon1, double azi1, double s12,
                     double* plat2, double* plon2, double* pazi2);
    
    /**
     * Solve the inverse geodesic problem.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lat1 latitude of point 1 (degrees).
     * @param[in] lon1 longitude of point 1 (degrees).
     * @param[in] lat2 latitude of point 2 (degrees).
     * @param[in] lon2 longitude of point 2 (degrees).
     * @param[out] ps12 pointer to the distance between point 1 and point 2
     *   (meters).
     * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     *
     * \e g must have been initialized with a call to geod_init().  \e lat1
     * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
     * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  The values of
     * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
     * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you
     * do not need some quantities computed.
     *
     * If either point is at a pole, the azimuth is defined by keeping the
     * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
     * taking the limit &epsilon; &rarr; 0+.
     *
     * The solution to the inverse problem is found using Newton's method.  If
     * this fails to converge (this is very unlikely in geodetic applications
     * but does occur for very eccentric ellipsoids), then the bisection method
     * is used to refine the solution.
     *
     * Example, determine the distance between JFK and Singapore Changi Airport:
     @code
     struct geod_geodesic g;
     double s12;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0);
     printf("%.3f\n", s12);
     @endcode
     **********************************************************************/
    void geod_inverse(const struct geod_geodesic* g,
                      double lat1, double lon1, double lat2, double lon2,
                      double* ps12, double* pazi1, double* pazi2);
    
    /**
     * Compute the position along a geod_geodesicline.
     *
     * @param[in] l a pointer to the geod_geodesicline object specifying the
     *   geodesic line.
     * @param[in] s12 distance between point 1 and point 2 (meters); it can be
     *   negative.
     * @param[out] plat2 pointer to the latitude of point 2 (degrees).
     * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
     *   that \e l was initialized with \e caps |= GEOD_LONGITUDE.
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     *
     * \e l must have been initialized with a call to geod_lineinit() with \e
     * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
     * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
     * \e plat2, etc., may be replaced by 0, if you do not need some quantities
     * computed.
     *
     * Example, compute way points between JFK and Singapore Changi Airport
     * the "obvious" way using geod_direct():
     @code
     struct geod_geodesic g;
     double s12, azi1, lat[101],lon[101];
     int i;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
     for (i = 0; i < 101; ++i) {
     geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0);
     printf("%.5f %.5f\n", lat[i], lon[i]);
     }
     @endcode
     * A faster way using geod_position():
     @code
     struct geod_geodesic g;
     struct geod_geodesicline l;
     double s12, azi1, lat[101],lon[101];
     int i;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
     geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0);
     for (i = 0; i < 101; ++i) {
     geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0);
     printf("%.5f %.5f\n", lat[i], lon[i]);
     }
     @endcode
     **********************************************************************/
    void geod_position(const struct geod_geodesicline* l, double s12,
                       double* plat2, double* plon2, double* pazi2);
    
    /**
     * The general direct geodesic problem.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lat1 latitude of point 1 (degrees).
     * @param[in] lon1 longitude of point 1 (degrees).
     * @param[in] azi1 azimuth at point 1 (degrees).
     * @param[in] arcmode flag determining the meaning of the \e
     *   s12_a12.
     * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
     *   point 1 and point 2 (meters); otherwise it is the arc length between
     *   point 1 and point 2 (degrees); it can be negative.
     * @param[out] plat2 pointer to the latitude of point 2 (degrees).
     * @param[out] plon2 pointer to the longitude of point 2 (degrees).
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     * @param[out] ps12 pointer to the distance between point 1 and point 2
     *   (meters).
     * @param[out] pm12 pointer to the reduced length of geodesic (meters).
     * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
     *   point 1 (dimensionless).
     * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
     *   point 2 (dimensionless).
     * @param[out] pS12 pointer to the area under the geodesic
     *   (meters<sup>2</sup>).
     * @return \e a12 arc length of between point 1 and point 2 (degrees).
     *
     * \e g must have been initialized with a call to geod_init().  \e lat1
     * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
     * should be in the range [&minus;540&deg;, 540&deg;).  The function value \e
     * a12 equals \e s12_a12 is \e arcmode is non-zero.  Any of the "return"
     * arguments \e plat2, etc., may be replaced by 0, if you do not need some
     * quantities computed.
     **********************************************************************/
    double geod_gendirect(const struct geod_geodesic* g,
                          double lat1, double lon1, double azi1,
                          int arcmode, double s12_a12,
                          double* plat2, double* plon2, double* pazi2,
                          double* ps12, double* pm12, double* pM12, double* pM21,
                          double* pS12);
    
    /**
     * The general inverse geodesic calculation.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lat1 latitude of point 1 (degrees).
     * @param[in] lon1 longitude of point 1 (degrees).
     * @param[in] lat2 latitude of point 2 (degrees).
     * @param[in] lon2 longitude of point 2 (degrees).
     * @param[out] ps12 pointer to the distance between point 1 and point 2
     *  (meters).
     * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     * @param[out] pm12 pointer to the reduced length of geodesic (meters).
     * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
     *   point 1 (dimensionless).
     * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
     *   point 2 (dimensionless).
     * @param[out] pS12 pointer to the area under the geodesic
     *   (meters<sup>2</sup>).
     * @return \e a12 arc length of between point 1 and point 2 (degrees).
     *
     * \e g must have been initialized with a call to geod_init().  \e lat1
     * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
     * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  Any of the
     * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
     * some quantities computed.
     **********************************************************************/
    double geod_geninverse(const struct geod_geodesic* g,
                           double lat1, double lon1, double lat2, double lon2,
                           double* ps12, double* pazi1, double* pazi2,
                           double* pm12, double* pM12, double* pM21,
                           double* pS12);
    
    /**
     * The general position function.
     *
     * @param[in] l a pointer to the geod_geodesicline object specifying the
     *   geodesic line.
     * @param[in] arcmode flag determining the meaning of the second parameter;
     *   if arcmode is 0, then \e l must have been initialized with \e caps |=
     *   GEOD_DISTANCE_IN.
     * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
     *   point 1 and point 2 (meters); otherwise it is the arc length between
     *   point 1 and point 2 (degrees); it can be negative.
     * @param[out] plat2 pointer to the latitude of point 2 (degrees).
     * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
     *   that \e l was initialized with \e caps |= GEOD_LONGITUDE.
     * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
     * @param[out] ps12 pointer to the distance between point 1 and point 2
     *   (meters); requires that \e l was initialized with \e caps |=
     *   GEOD_DISTANCE.
     * @param[out] pm12 pointer to the reduced length of geodesic (meters);
     *   requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH.
     * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
     *   point 1 (dimensionless); requires that \e l was initialized with \e caps
     *   |= GEOD_GEODESICSCALE.
     * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
     *   point 2 (dimensionless); requires that \e l was initialized with \e caps
     *   |= GEOD_GEODESICSCALE.
     * @param[out] pS12 pointer to the area under the geodesic
     *   (meters<sup>2</sup>); requires that \e l was initialized with \e caps |=
     *   GEOD_AREA.
     * @return \e a12 arc length of between point 1 and point 2 (degrees).
     *
     * \e l must have been initialized with a call to geod_lineinit() with \e
     * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
     * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
     * \e plat2, etc., may be replaced by 0, if you do not need some quantities
     * computed.  Requesting a value which \e l is not capable of computing is
     * not an error; the corresponding argument will not be altered.
     *
     * Example, compute way points between JFK and Singapore Changi Airport
     * using geod_genposition().  In this example, the points are evenly space in
     * arc length (and so only approximately equally space in distance).  This is
     * faster than using geod_position() would be appropriate if drawing the path
     * on a map.
     @code
     struct geod_geodesic g;
     struct geod_geodesicline l;
     double a12, azi1, lat[101],lon[101];
     int i;
     geod_init(&g, 6378137, 1/298.257223563);
     a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
     0, &azi1, 0, 0, 0, 0, 0);
     geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
     for (i = 0; i < 101; ++i) {
     geod_genposition(&l, 1, i * a12 * 0.01,
     lat + i, lon + i, 0, 0, 0, 0, 0, 0);
     printf("%.5f %.5f\n", lat[i], lon[i]);
     }
     @endcode
     **********************************************************************/
    double geod_genposition(const struct geod_geodesicline* l,
                            int arcmode, double s12_a12,
                            double* plat2, double* plon2, double* pazi2,
                            double* ps12, double* pm12,
                            double* pM12, double* pM21,
                            double* pS12);
    
    /**
     * Initialize a geod_polygon object.
     *
     * @param[out] p a pointer to the object to be initialized.
     * @param[in] polylinep non-zero if a polyline instead of a polygon.
     *
     * If \e polylinep is zero, then the sequence of vertices and edges added by
     * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
     * the perimeter and area are returned by geod_polygon_compute().  If \e
     * polylinep is non-zero, then the vertices and edges define a polyline and
     * only the perimeter is returned by geod_polygon_compute().
     *
     * An example of the use of this function is given in the documentation for
     * geod_polygon_compute().
     **********************************************************************/
    void geod_polygon_init(struct geod_polygon* p, int polylinep);
    
    /**
     * Add a point to the polygon or polyline.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in,out] p a pointer to the geod_polygon object specifying the
     *   polygon.
     * @param[in] lat the latitude of the point (degrees).
     * @param[in] lon the longitude of the point (degrees).
     *
     * \e g and \e p must have been initialized with calls to geod_init() and
     * geod_polygon_init(), respectively.  The same \e g must be used for all the
     * points and edges in a polygon.  \e lat should be in the range
     * [&minus;90&deg;, 90&deg;] and \e lon should be in the range
     * [&minus;540&deg;, 540&deg;).
     *
     * An example of the use of this function is given in the documentation for
     * geod_polygon_compute().
     **********************************************************************/
    void geod_polygon_addpoint(const struct geod_geodesic* g,
                               struct geod_polygon* p,
                               double lat, double lon);
    
    /**
     * Add an edge to the polygon or polyline.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in,out] p a pointer to the geod_polygon object specifying the
     *   polygon.
     * @param[in] azi azimuth at current point (degrees).
     * @param[in] s distance from current point to next point (meters).
     *
     * \e g and \e p must have been initialized with calls to geod_init() and
     * geod_polygon_init(), respectively.  The same \e g must be used for all the
     * points and edges in a polygon.  \e azi should be in the range
     * [&minus;540&deg;, 540&deg;).  This does nothing if no points have been
     * added yet.  The \e lat and \e lon fields of \e p give the location of
     * the new vertex.
     **********************************************************************/
    void geod_polygon_addedge(const struct geod_geodesic* g,
                              struct geod_polygon* p,
                              double azi, double s);
    
    /**
     * Return the results for a polygon.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] p a pointer to the geod_polygon object specifying the polygon.
     * @param[in] reverse if non-zero then clockwise (instead of
     *   counter-clockwise) traversal counts as a positive area.
     * @param[in] sign if non-zero then return a signed result for the area if
     *   the polygon is traversed in the "wrong" direction instead of returning
     *   the area for the rest of the earth.
     * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
     *   only set if \e polyline is non-zero in the call to geod_polygon_init().
     * @param[out] pP pointer to the perimeter of the polygon or length of the
     *   polyline (meters).
     * @return the number of points.
     *
     * Only simple polygons (which are not self-intersecting) are allowed.
     * There's no need to "close" the polygon by repeating the first vertex.  Set
     * \e pA or \e pP to zero, if you do not want the corresponding quantity
     * returned.
     *
     * Example, compute the perimeter and area of the geodesic triangle with
     * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
     @code
     double A, P;
     int n;
     struct geod_geodesic g;
     struct geod_polygon p;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_polygon_init(&p, 0);
     
     geod_polygon_addpoint(&g, &p,  0,  0);
     geod_polygon_addpoint(&g, &p,  0, 90);
     geod_polygon_addpoint(&g, &p, 90,  0);
     n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
     printf("%d %.8f %.3f\n", n, P, A);
     @endcode
     **********************************************************************/
    unsigned geod_polygon_compute(const struct geod_geodesic* g,
                                  const struct geod_polygon* p,
                                  int reverse, int sign,
                                  double* pA, double* pP);
    
    /**
     * Return the results assuming a tentative final test point is added;
     * however, the data for the test point is not saved.  This lets you report a
     * running result for the perimeter and area as the user moves the mouse
     * cursor.  Ordinary floating point arithmetic is used to accumulate the data
     * for the test point; thus the area and perimeter returned are less accurate
     * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] p a pointer to the geod_polygon object specifying the polygon.
     * @param[in] lat the latitude of the test point (degrees).
     * @param[in] lon the longitude of the test point (degrees).
     * @param[in] reverse if non-zero then clockwise (instead of
     *   counter-clockwise) traversal counts as a positive area.
     * @param[in] sign if non-zero then return a signed result for the area if
     *   the polygon is traversed in the "wrong" direction instead of returning
     *   the area for the rest of the earth.
     * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
     *   only set if \e polyline is non-zero in the call to geod_polygon_init().
     * @param[out] pP pointer to the perimeter of the polygon or length of the
     *   polyline (meters).
     * @return the number of points.
     *
     * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
     * lon should be in the range [&minus;540&deg;, 540&deg;).
     **********************************************************************/
    unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
                                    const struct geod_polygon* p,
                                    double lat, double lon,
                                    int reverse, int sign,
                                    double* pA, double* pP);
    
    /**
     * Return the results assuming a tentative final test point is added via an
     * azimuth and distance; however, the data for the test point is not saved.
     * This lets you report a running result for the perimeter and area as the
     * user moves the mouse cursor.  Ordinary floating point arithmetic is used
     * to accumulate the data for the test point; thus the area and perimeter
     * returned are less accurate than if geod_polygon_addedge() and
     * geod_polygon_compute() are used.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] p a pointer to the geod_polygon object specifying the polygon.
     * @param[in] azi azimuth at current point (degrees).
     * @param[in] s distance from current point to final test point (meters).
     * @param[in] reverse if non-zero then clockwise (instead of
     *   counter-clockwise) traversal counts as a positive area.
     * @param[in] sign if non-zero then return a signed result for the area if
     *   the polygon is traversed in the "wrong" direction instead of returning
     *   the area for the rest of the earth.
     * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
     *   only set if \e polyline is non-zero in the call to geod_polygon_init().
     * @param[out] pP pointer to the perimeter of the polygon or length of the
     *   polyline (meters).
     * @return the number of points.
     *
     * \e azi should be in the range [&minus;540&deg;, 540&deg;).
     **********************************************************************/
    unsigned geod_polygon_testedge(const struct geod_geodesic* g,
                                   const struct geod_polygon* p,
                                   double azi, double s,
                                   int reverse, int sign,
                                   double* pA, double* pP);
    
    /**
     * A simple interface for computing the area of a geodesic polygon.
     *
     * @param[in] g a pointer to the geod_geodesic object specifying the
     *   ellipsoid.
     * @param[in] lats an array of latitudes of the polygon vertices (degrees).
     * @param[in] lons an array of longitudes of the polygon vertices (degrees).
     * @param[in] n the number of vertices.
     * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>).
     * @param[out] pP pointer to the perimeter of the polygon (meters).
     *
     * \e lats should be in the range [&minus;90&deg;, 90&deg;]; \e lons should
     * be in the range [&minus;540&deg;, 540&deg;).
     *
     * Only simple polygons (which are not self-intersecting) are allowed.
     * There's no need to "close" the polygon by repeating the first vertex.  The
     * area returned is signed with counter-clockwise traversal being treated as
     * positive.
     *
     * Example, compute the area of Antarctic:
     @code
     double
     lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
     -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
     lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
     88, 59, 25, -4, -14, -33, -46, -61};
     struct geod_geodesic g;
     double A, P;
     geod_init(&g, 6378137, 1/298.257223563);
     geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
     printf("%.0f %.2f\n", A, P);
     @endcode
     **********************************************************************/
    void geod_polygonarea(const struct geod_geodesic* g,
                          double lats[], double lons[], int n,
                          double* pA, double* pP);
    
    /**
     * mask values for the the \e caps argument to geod_lineinit().
     **********************************************************************/
    enum geod_mask {
        GEOD_NONE         = 0U,                     /**< Calculate nothing */
        GEOD_LATITUDE     = 1U<<7  | 0U,            /**< Calculate latitude */
        GEOD_LONGITUDE    = 1U<<8  | 1U<<3,         /**< Calculate longitude */
        GEOD_AZIMUTH      = 1U<<9  | 0U,            /**< Calculate azimuth */
        GEOD_DISTANCE     = 1U<<10 | 1U<<0,         /**< Calculate distance */
        GEOD_DISTANCE_IN  = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input  */
        GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */
        GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */
        GEOD_AREA         = 1U<<14 | 1U<<4,         /**< Calculate reduced length */
        GEOD_ALL          = 0x7F80U| 0x1FU          /**< Calculate everything */
    };
    
#if defined(__cplusplus)
}
#endif

#endif
