/**
* \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,
*
* Algorithms for geodesics,
* J. Geodesy 87, 43--55 (2013);
* DOI:
* 10.1007/s00190-012-0578-z;
* addenda:
* geod-addenda.html.
* .
* The principal advantages of these algorithms over previous ones (e.g.,
* Vincenty, 1975) are
* - accurate to round off for |f| < 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, \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 |f| < 1/50; for the
* WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably
* accurate results are obtained for |f| < 1/5.) For a prolate
* ellipsoid, specify \e f < 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°.
*
* 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 − (1 − \e M12 \e M21) \e
* m23 / \e m12
* - \e M31 = \e M32 \e M21 − (1 − \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 = −\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 =
* −\e S12. (This occurs when the longitude difference is near
* ±180° for oblate ellipsoids.)
* - \e lon2 = \e lon1 ± 180° (with neither at a pole). If
* \e azi1 = 0° or ±180°, the geodesic is unique.
* Otherwise there are two geodesics and the second one is obtained by
* setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e
* S12 = −\e S12. (This occurs when the \e lat2 is near
* −\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, −\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 GeographicLib. The
* "class data" is represented by the structs geod_geodesic and
* geod_geodesicline and these are passed as initial arguments to the member
* functions. (But note that the PolygonArea class has been replaced by a
* simple function interface and the area computation does use the accurate
* summing providing by the Accumulator class.) 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) and licensed
* under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/
*
* This library was distributed with
* GeographicLib 1.30.
**********************************************************************/
#if !defined(GEODESIC_H)
#define GEODESIC_H 1
#if defined(__cplusplus)
extern "C" {
#endif
/**
* The struct containing information about the ellipsoid.
**********************************************************************/
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.
**********************************************************************/
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 */
};
/**
* Initialize a geod_geodesic object
*
* @param[out] g 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 the object to be initialized.
* @param[in] g 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 [−90°, 90°]; \e azi1 should be in the
* range [−540°, 540°).
*
* The geod_mask values are
* - \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 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 [−90°, 90°]; \e lon1 and \e azi1
* should be in the range [−540°, 540°). The values of \e lon2
* and \e azi2 returned are in the range [−180°, 180°). Any of
* the "return" arguments 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 and writing \e lat = ±(90° − ε)
* and taking the limit ε → 0+. An arc length greater that
* 180° 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°.)
*
* 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 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 [−90°, 90°]; \e lon1 and
* \e lon2 should be in the range [−540°, 540°). The values of
* \e azi1 and \e azi2 returned are in the range [−180°, 180°).
* Any of the "return" arguments 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 and writing \e lat = ±(90° − ε)
* and taking the limit ε → 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 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 the \e 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 [−180°, 180°). Any of the "return" arguments
* 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 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
* (meters2).
* @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 [−90°, 90°]; \e lon1 and \e azi1
* should be in the range [−540°, 540°). The function value \e
* a12 equals \e s12_a12 is \e arcmode is non-zero. Any of the "return"
* arguments 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 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
* (meters2).
* @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 [−90°, 90°]; \e lon1 and
* \e lon2 should be in the range [−540°, 540°). Any of the
* "return" arguments 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 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
* (meters2); 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 [−180°, 180°). Any of the "return" arguments
* 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, computer 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);
/**
* The area of a geodesic polygon.
*
* @param[in] g 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 (meters2).
* @param[out] pP pointer to the perimeter of the polygon (meters).
*
* \e lats should be in the range [−90°, 90°]; \e lons should
* be in the range [−540°, 540°).
*
* 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 geod_geodesicline capabilities.
**********************************************************************/
enum geod_mask {
GEOD_NONE = 0U,
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