diff options
Diffstat (limited to 'src/geodesic.h')
| -rw-r--r-- | src/geodesic.h | 481 |
1 files changed, 302 insertions, 179 deletions
diff --git a/src/geodesic.h b/src/geodesic.h index 7bd8270f..c3f28c79 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -1,6 +1,6 @@ /** * \file geodesic.h - * \brief Header for the geodesic routines in C + * \brief API for the geodesic routines in C * * This an implementation in C of the geodesic algorithms described in * - C. F. F. Karney, @@ -9,7 +9,7 @@ * J. Geodesy <b>87</b>, 43--55 (2013); * DOI: <a href="https://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"> + * addenda: <a href="http://geographiclib.sourceforge.net/geod-addenda.html"> * geod-addenda.html</a>. * . * The principal advantages of these algorithms over previous ones (e.g., @@ -75,30 +75,29 @@ * (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 point 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 point 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 \e lat2 is near −\e lat1 for + * - \e lat1 = −\e lat2 (with neither point 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 point 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 \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. + * - 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 <a href="http://geographiclib.sf.net"> GeographicLib</a>. The "class - * data" is represented by the structs geod_geodesic, geod_geodesicline, + * in <a href="http://geographiclib.sourceforge.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 @@ -108,12 +107,12 @@ * 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-2015) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2016) <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.44. + * <a href="../index.html">GeographicLib</a> 1.46. **********************************************************************/ #if !defined(GEODESIC_H) @@ -128,12 +127,12 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 44 +#define GEODESIC_VERSION_MINOR 46 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_PATCH 0 +#define GEODESIC_VERSION_PATCH 1 /** * Pack the version components into a single integer. Users should not rely on @@ -177,7 +176,8 @@ extern "C" { /** * The struct containing information about a single geodesic. This must be - * initialized by geod_lineinit() before use. + * initialized by geod_lineinit(), geod_directline(), geod_gendirectline(), + * or geod_inverseline() before use. **********************************************************************/ struct geod_geodesicline { double lat1; /**< the starting latitude */ @@ -185,9 +185,13 @@ extern "C" { double azi1; /**< the starting azimuth */ double a; /**< the equatorial radius */ double f; /**< the flattening */ + double salp1; /**< sine of \e azi1 */ + double calp1; /**< cosine of \e azi1 */ + double a13; /**< arc length to reference point */ + double s13; /**< distance to reference point */ /**< @cond SKIP */ double b, c2, f1, salp0, calp0, k2, - salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1, + 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 */ @@ -223,46 +227,6 @@ extern "C" { 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 [−90°, 90°]. - * - * 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 @@ -270,7 +234,7 @@ extern "C" { * @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 + * @param[in] s12 distance from point 1 to 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). @@ -303,6 +267,51 @@ extern "C" { 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] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_UNROLL "unrolls" \e lon2. + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance + * from point 1 to point 2 (meters); otherwise it is the arc length + * from point 1 to 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 from point 1 to 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 from point 1 to point 2 (degrees). + * + * \e g must have been initialized with a call to geod_init(). \e lat1 + * should be in the range [−90°, 90°]. The function value \e + * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return" + * arguments, \e plat2, etc., may be replaced by 0, if you do not need some + * quantities computed. + * + * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so + * that the quantity \e lon2 − \e lon1 indicates how many times and in + * what sense the geodesic encircles the ellipsoid. + **********************************************************************/ + double geod_gendirect(const struct geod_geodesic* g, + double lat1, double lon1, double azi1, + unsigned flags, double s12_a12, + double* plat2, double* plon2, double* pazi2, + double* ps12, double* pm12, double* pM12, double* pM21, + double* pS12); + + /** * Solve the inverse geodesic problem. * * @param[in] g a pointer to the geod_geodesic object specifying the @@ -311,7 +320,7 @@ extern "C" { * @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 + * @param[out] ps12 pointer to the distance from point 1 to 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). @@ -345,22 +354,180 @@ extern "C" { double* ps12, double* pazi1, double* pazi2); /** + * 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 from point 1 to 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 from point 1 to 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°]. 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); + + /** + * 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 [−90°, 90°]. + * + * 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). + * + * When initialized by this function, point 3 is undefined (l->s13 = l->a13 = + * NaN). + **********************************************************************/ + void geod_lineinit(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, unsigned caps); + + /** + * Initialize a geod_geodesicline object in terms of the direct geodesic + * problem. + * + * @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] s12 distance from point 1 to point 2 (meters); it can be + * negative. + * @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(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the direct geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_directline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, double s12, + unsigned caps); + + /** + * Initialize a geod_geodesicline object in terms of the direct geodesic + * problem spacified in terms of either distance or arc length. + * + * @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] flags either GEOD_NOFLAGS or GEOD_ARCMODE to determining the + * meaning of the \e s12_a12. + * @param[in] s12_a12 if \e flags = GEOD_NOFLAGS, this is the distance + * from point 1 to point 2 (meters); if \e flags = GEOD_ARCMODE, it is + * the arc length from point 1 to point 2 (degrees); it can be + * negative. + * @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(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the direct geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_gendirectline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double azi1, + unsigned flags, double s12_a12, + unsigned caps); + + /** + * Initialize a geod_geodesicline object in terms of the inverse geodesic + * problem. + * + * @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] lat2 latitude of point 2 (degrees). + * @param[in] lon2 longitude of point 2 (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(). + * + * This function sets point 3 of the geod_geodesicline to correspond to point + * 2 of the inverse geodesic problem. See geod_lineinit() for more + * informaion. + **********************************************************************/ + void geod_inverseline(struct geod_geodesicline* l, + const struct geod_geodesic* g, + double lat1, double lon1, double lat2, double lon2, + unsigned caps); + + /** * 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 + * @param[in] s12 distance from point 1 to 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 [−180°, 180°). Any of the "return" arguments - * \e plat2, etc., may be replaced by 0, if you do not need some quantities - * computed. + * \e l must have been initialized with a call, e.g., 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 \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(): @@ -379,13 +546,12 @@ extern "C" { @code{.c} struct geod_geodesic g; struct geod_geodesicline l; - double s12, azi1, lat[101],lon[101]; + double 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); + geod_inverseline(&l, &g, 40.64, -73.78, 1.36, 103.99, 0); + for (i = 0; i <= 100; ++i) { + geod_position(&l, i * l.s13 * 0.01, lat + i, lon + i, 0); printf("%.5f %.5f\n", lat[i], lon[i]); } @endcode @@ -394,84 +560,6 @@ extern "C" { 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] flags bitor'ed combination of geod_flags(); \e flags & - * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & - * GEOD_LONG_UNROLL "unrolls" \e lon2. - * @param[in] s12_a12 if \e flags & GEOD_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 [−90°, 90°]. The function value \e - * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return" - * arguments, \e plat2, etc., may be replaced by 0, if you do not need some - * quantities computed. - * - * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so - * that the quantity \e lon2 − \e lon1 indicates how many times and in - * what sense the geodesic encircles the ellipsoid. - **********************************************************************/ - double geod_gendirect(const struct geod_geodesic* g, - double lat1, double lon1, double azi1, - unsigned flags, 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 [−90°, 90°]. 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 @@ -481,14 +569,14 @@ extern "C" { * GEOD_LONG_UNROLL "unrolls" \e lon2; if \e flags & GEOD_ARCMODE is 0, * then \e l must have been initialized with \e caps |= GEOD_DISTANCE_IN. * @param[in] s12_a12 if \e flags & GEOD_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 + * distance from point 1 to point 2 (meters); otherwise it is the + * arc length from point 1 to 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 + * @param[out] ps12 pointer to the distance from point 1 to 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); @@ -502,7 +590,7 @@ extern "C" { * @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). + * @return \e a12 arc length from point 1 to point 2 (degrees). * * \e l must have been initialized with a call to geod_lineinit() with \e * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range @@ -515,22 +603,21 @@ extern "C" { * that the quantity \e lon2 − \e lon1 indicates how many times and in * what sense the geodesic encircles the ellipsoid. * - * 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. + * 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 spaced in distance). This is + * faster than using geod_position() and would be appropriate if drawing the + * path on a map. @code{.c} struct geod_geodesic g; struct geod_geodesicline l; - double a12, azi1, lat[101], lon[101]; + double 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, + geod_inverseline(&l, &g, 40.64, -73.78, 1.36, 103.99, + GEOD_LATITUDE | GEOD_LONGITUDE); + for (i = 0; i <= 100; ++i) { + geod_genposition(&l, GEOD_ARCMODE, i * l.a13 * 0.01, lat + i, lon + i, 0, 0, 0, 0, 0, 0); printf("%.5f %.5f\n", lat[i], lon[i]); } @@ -544,6 +631,36 @@ extern "C" { double* pS12); /** + * Specify position of point 3 in terms of distance. + * + * @param[inout] l a pointer to the geod_geodesicline object. + * @param[in] s13 the distance from point 1 to point 3 (meters); it + * can be negative. + * + * This is only useful if the geod_geodesicline object has been constructed + * with \e caps |= GEOD_DISTANCE_IN. + **********************************************************************/ + void geod_setdistance(struct geod_geodesicline* l, double s13); + + /** + * Specify position of point 3 in terms of either distance or arc length. + * + * @param[inout] l a pointer to the geod_geodesicline object. + * @param[in] flags either GEOD_NOFLAGS or GEOD_ARCMODE to determining the + * meaning of the \e s13_a13. + * @param[in] s13_a13 if \e flags = GEOD_NOFLAGS, this is the distance + * from point 1 to point 3 (meters); if \e flags = GEOD_ARCMODE, it is + * the arc length from point 1 to point 3 (degrees); it can be + * negative. + * + * If flags = GEOD_NOFLAGS, this calls geod_setdistance(). If flags = + * GEOD_ARCMODE, the \e s13 is only set if the geod_geodesicline object has + * been constructed with \e caps |= GEOD_DISTANCE. + **********************************************************************/ + void geod_gensetdistance(struct geod_geodesicline* l, + unsigned flags, double s13_a13); + + /** * Initialize a geod_polygon object. * * @param[out] p a pointer to the object to be initialized. @@ -565,6 +682,13 @@ extern "C" { void geod_polygon_init(struct geod_polygon* p, int polylinep); /** + * Clear the polygon, allowing a new polygon to be started. + * + * @param[in,out] p a pointer to the object to be cleared. + **********************************************************************/ + void geod_polygon_clear(struct geod_polygon* p); + + /** * Add a point to the polygon or polyline. * * @param[in] g a pointer to the geod_geodesic object specifying the @@ -630,6 +754,8 @@ extern "C" { * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding * quantity returned. * + * More points can be added to the polygon after this call. + * * Example, compute the perimeter and area of the geodesic triangle with * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). @code{.c} @@ -774,10 +900,7 @@ extern "C" { enum geod_flags { GEOD_NOFLAGS = 0U, /**< No flags */ GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */ - GEOD_LONG_UNROLL = 1U<<15, /**< Unroll the longitude */ - /**< @cond SKIP */ - GEOD_LONG_NOWRAP = GEOD_LONG_UNROLL /* For backward compatibility only */ - /**< @endcond */ + GEOD_LONG_UNROLL = 1U<<15 /**< Unroll the longitude */ }; #if defined(__cplusplus) |