aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--HOWTO-RELEASE2
-rw-r--r--man/man1/geod.126
-rw-r--r--man/man3/Makefile.am2
-rw-r--r--man/man3/Makefile.in2
-rw-r--r--src/Makefile.am8
-rw-r--r--src/Makefile.in14
-rw-r--r--src/geod_interface.c8
-rw-r--r--src/geod_interface.h4
-rw-r--r--src/geodesic.c346
-rw-r--r--src/geodesic.h659
-rw-r--r--src/makefile.vc8
-rw-r--r--src/proj.def117
12 files changed, 799 insertions, 397 deletions
diff --git a/HOWTO-RELEASE b/HOWTO-RELEASE
index 8ae52131..41bc3090 100644
--- a/HOWTO-RELEASE
+++ b/HOWTO-RELEASE
@@ -22,7 +22,7 @@
- If any interfaces have been removed since the last public release, then
set age to 0.
-4.5) Run "autgen.sh" (hopefully on the same machine it was last run on)
+4.5) Run "autogen.sh" (hopefully on the same machine it was last run on)
5) Add a note to the ChangeLog that a new release is being issued, and what
the release number is.
diff --git a/man/man1/geod.1 b/man/man1/geod.1
index ab6c74f7..0a2522ea 100644
--- a/man/man1/geod.1
+++ b/man/man1/geod.1
@@ -1,8 +1,8 @@
-.\" @(#)geod.1 - 1.1
-.nr LL 5.5i
+.\" @(#)geod.1
+.nr LL 6.5i
.ad b
.hy 1
-.TH GEOD 1 "2012/11/29 Rel. 4.8"
+.TH GEOD 1 "2013/02/27 Rel. 4.8"
.SH NAME
geod \- direct geodesic computations
.br
@@ -171,7 +171,7 @@ and number of points to be determined.
.RE
.SH EXAMPLE
The following script determines the geodesic azimuths and distance in
-U.S. stature miles from Boston, MA, to Portland, OR:
+U.S. statute miles from Boston, MA, to Portland, OR:
.RS 5
\f(CWgeod +ellps=clrk66 <<EOF -I +units=us-mi
42d15'N 71d07'W 45d31'N 123d41'W
@@ -199,9 +199,19 @@ which gives:
Note: lack of precision in the distance value compromises
the precision of the Portland location.
.SH SEE ALSO
-C. F. F. Karney,
-.I "Algorithms for Geodesics," J. Geodesy (2012);
-DOI: 10.1007/s00190-012-0578-z,
+.B geodesic(3)
+.br
+The \fBGeodSolve\fR utility in \fBGeographicLib\fR,
+http://geographiclib.sf.net. With the \fB-E\fR option, this solves the
+geodesic problems in terms of elliptic integrals; the results are
+accurate for arbitrary \fIf\fR.
+.br
+C. F. F. Karney, \fIAlgorithms for Geodesics\fR,
+.br
+J. Geodesy \fB87\fR, 43-55 (2013);
+.br
+DOI: http://dx.doi.org/10.1007/s00190-012-0578-z,
+.br
http://geographiclib.sf.net/geod-addenda.html.
.SH HOME PAGE
-http://www.remotesensing.org/proj
+http://proj.osgeo.org
diff --git a/man/man3/Makefile.am b/man/man3/Makefile.am
index 2a41471a..4a3128a8 100644
--- a/man/man3/Makefile.am
+++ b/man/man3/Makefile.am
@@ -1,3 +1,3 @@
-man_MANS = pj_init.3
+man_MANS = pj_init.3 geodesic.3
EXTRA_DIST = $(man_MANS)
diff --git a/man/man3/Makefile.in b/man/man3/Makefile.in
index faa81dea..c9e3da99 100644
--- a/man/man3/Makefile.in
+++ b/man/man3/Makefile.in
@@ -187,7 +187,7 @@ target_alias = @target_alias@
top_build_prefix = @top_build_prefix@
top_builddir = @top_builddir@
top_srcdir = @top_srcdir@
-man_MANS = pj_init.3
+man_MANS = pj_init.3 geodesic.3
EXTRA_DIST = $(man_MANS)
all: all-am
diff --git a/src/Makefile.am b/src/Makefile.am
index bbf28ba6..c1d1998a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -4,15 +4,14 @@ EXTRA_PROGRAMS = multistresstest
INCLUDES = -DPROJ_LIB=\"$(pkgdatadir)\" \
-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
-include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h
+include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h geodesic.h
EXTRA_DIST = makefile.vc proj.def
proj_SOURCES = proj.c gen_cheb.c p_series.c
cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h \
- geodesic.c geodesic.h
+geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
multistresstest_SOURCES = multistresstest.c
proj_LDADD = libproj.la
@@ -64,8 +63,7 @@ libproj_la_SOURCES = \
nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
- jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c
-
+ jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c
install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/Makefile.in b/src/Makefile.in
index eae1c90c..dd32d94b 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -107,7 +107,7 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \
- pj_apply_vgridshift.lo
+ pj_apply_vgridshift.lo geodesic.lo
libproj_la_OBJECTS = $(am_libproj_la_OBJECTS)
libproj_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
@@ -118,7 +118,7 @@ am_cs2cs_OBJECTS = cs2cs.$(OBJEXT) gen_cheb.$(OBJEXT) \
cs2cs_OBJECTS = $(am_cs2cs_OBJECTS)
cs2cs_DEPENDENCIES = libproj.la
am_geod_OBJECTS = geod.$(OBJEXT) geod_set.$(OBJEXT) \
- geod_interface.$(OBJEXT) geodesic.$(OBJEXT)
+ geod_interface.$(OBJEXT)
geod_OBJECTS = $(am_geod_OBJECTS)
geod_DEPENDENCIES = libproj.la
am_multistresstest_OBJECTS = multistresstest.$(OBJEXT)
@@ -270,14 +270,12 @@ top_srcdir = @top_srcdir@
INCLUDES = -DPROJ_LIB=\"$(pkgdatadir)\" \
-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
-include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h
+include_HEADERS = proj_api.h org_proj4_Projections.h org_proj4_PJ.h geodesic.h
EXTRA_DIST = makefile.vc proj.def
proj_SOURCES = proj.c gen_cheb.c p_series.c
cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h \
- geodesic.c geodesic.h
-
+geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
multistresstest_SOURCES = multistresstest.c
proj_LDADD = libproj.la
cs2cs_LDADD = libproj.la
@@ -325,7 +323,7 @@ libproj_la_SOURCES = \
nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
- jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c
+ jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c
all: proj_config.h
$(MAKE) $(AM_MAKEFLAGS) all-am
@@ -583,7 +581,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod_interface.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geod_set.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geodesic.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geodesic.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/jniproj.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mk_cheby.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/multistresstest.Po@am__quote@
diff --git a/src/geod_interface.c b/src/geod_interface.c
index c144a35d..dc757282 100644
--- a/src/geod_interface.c
+++ b/src/geod_interface.c
@@ -2,20 +2,20 @@
#include "geod_interface.h"
void geod_ini(void) {
- GeodesicInit(&GlobalGeodesic, geod_a, geod_f);
+ geod_init(&GlobalGeodesic, geod_a, geod_f);
}
void geod_pre(void) {
double
degree = PI/180,
lat1 = phi1 / degree, lon1 = lam1 /degree, azi1 = al12 / degree;
- GeodesicLineInit(&GlobalGeodesicLine, &GlobalGeodesic,
+ geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic,
lat1, lon1, azi1, 0U);
}
void geod_for(void) {
double degree = PI/180, s12 = geod_S, lat2, lon2, azi2;
- Position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
+ geod_position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
phi2 = lat2 * degree;
lam2 = lon2 * degree;
@@ -28,7 +28,7 @@ void geod_inv(void) {
lat1 = phi1 / degree, lon1 = lam1 / degree,
lat2 = phi2 / degree, lon2 = lam2 / degree,
azi1, azi2, s12;
- Inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
+ geod_inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
al12 = azi1 * degree; al21 = azi2 * degree; geod_S = s12;
}
diff --git a/src/geod_interface.h b/src/geod_interface.h
index 161c7572..255d505a 100644
--- a/src/geod_interface.h
+++ b/src/geod_interface.h
@@ -27,8 +27,8 @@ GEOD_EXTERN struct geodesic {
# define al21 GEODESIC.ALPHA21
# define geod_S GEODESIC.DIST
-GEOD_EXTERN struct Geodesic GlobalGeodesic;
-GEOD_EXTERN struct GeodesicLine GlobalGeodesicLine;
+GEOD_EXTERN struct geod_geodesic GlobalGeodesic;
+GEOD_EXTERN struct geod_geodesicline GlobalGeodesicLine;
GEOD_EXTERN int n_alpha, n_S;
GEOD_EXTERN double to_meter, fr_meter, del_alpha;
diff --git a/src/geodesic.c b/src/geodesic.c
index 794439a6..c5d481ee 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -1,22 +1,27 @@
+/**
+ * \file geodesic.c
+ * \brief Implementation of the geodesic routines in C
+ **********************************************************************/
+
/*
* This is a C implementation of the geodesic algorithms described in
*
- * C. F. F. Karney
- * Algorithms for geodesics
- * J. Geodesy (2012)
+ * C. F. F. Karney,
+ * Algorithms for geodesics,
+ * J. Geodesy <b>87</b>, 43--55 (2013);
* http://dx.doi.org/10.1007/s00190-012-0578-z
* Addenda: http://geographiclib.sf.net/geod-addenda.html
*
* See the comments in geodesic.h for documentation.
*
- * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed
+ * 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 file was distributed with GeographicLib 1.27.
*/
#include "geodesic.h"
+
+/** @cond SKIP */
#include <math.h>
#define GEOGRAPHICLIB_GEODESIC_ORDER 6
@@ -117,6 +122,19 @@ static real cbrtx(real x) {
return x < 0 ? -y : y;
}
+static real sumx(real u, real v, real* t) {
+ volatile real s = u + v;
+ volatile real up = s - v;
+ volatile real vpp = s - up;
+ up -= u;
+ vpp -= v;
+ *t = -(up + vpp);
+ /* error-free sum:
+ * u + v = s + t
+ * = round(u + v) + t */
+ return s;
+}
+
static real minx(real x, real y)
{ return x < y ? x : y; }
@@ -137,21 +155,30 @@ static real AngNormalize(real x)
static real AngNormalize2(real x)
{ return AngNormalize(fmod(x, (real)(360))); }
+static real AngDiff(real x, real y) {
+ real t, d = sumx(-x, y, &t);
+ if ((d - (real)(180)) + t > (real)(0)) /* y - x > 180 */
+ d -= (real)(360); /* exact */
+ else if ((d + (real)(180)) + t <= (real)(0)) /* y - x <= -180 */
+ d += (real)(360); /* exact */
+ return d + t;
+}
+
static real AngRound(real x) {
- const real z = (real)(0.0625); /* 1/16 */
+ const real z = 1/(real)(16);
volatile real y = fabs(x);
/* The compiler mustn't "simplify" z - (z - y) to y */
y = y < z ? z - (z - y) : y;
return x < 0 ? -y : y;
}
-static void A3coeff(struct Geodesic* g);
-static void C3coeff(struct Geodesic* g);
-static void C4coeff(struct Geodesic* g);
+static void A3coeff(struct geod_geodesic* g);
+static void C3coeff(struct geod_geodesic* g);
+static void C4coeff(struct geod_geodesic* g);
static real SinCosSeries(boolx sinp,
real sinx, real cosx,
const real c[], int n);
-static void Lengths(const struct Geodesic* g,
+static void Lengths(const struct geod_geodesic* g,
real eps, real sig12,
real ssig1, real csig1, real dn1,
real ssig2, real csig2, real dn2,
@@ -161,7 +188,7 @@ static void Lengths(const struct Geodesic* g,
/* Scratch areas of the right size */
real C1a[], real C2a[]);
static real Astroid(real x, real y);
-static real InverseStart(const struct Geodesic* g,
+static real InverseStart(const struct geod_geodesic* g,
real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real lam12,
@@ -170,7 +197,7 @@ static real InverseStart(const struct Geodesic* g,
real* psalp2, real* pcalp2,
/* Scratch areas of the right size */
real C1a[], real C2a[]);
-static real Lambda12(const struct Geodesic* g,
+static real Lambda12(const struct geod_geodesic* g,
real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real salp1, real calp1,
@@ -182,16 +209,19 @@ static real Lambda12(const struct Geodesic* g,
boolx diffp, real* pdlam12,
/* Scratch areas of the right size */
real C1a[], real C2a[], real C3a[]);
-static real A3f(const struct Geodesic* g, real eps);
-static void C3f(const struct Geodesic* g, real eps, real c[]);
-static void C4f(const struct Geodesic* g, real eps, real c[]);
+static real A3f(const struct geod_geodesic* g, real eps);
+static void C3f(const struct geod_geodesic* g, real eps, real c[]);
+static void C4f(const struct geod_geodesic* g, real eps, real c[]);
static real A1m1f(real eps);
static void C1f(real eps, real c[]);
static void C1pf(real eps, real c[]);
static real A2m1f(real eps);
static void C2f(real eps, real c[]);
+static int transit(real lon1, real lon2);
+
+/** @endcond */
-void GeodesicInit(struct Geodesic* g, real a, real f) {
+void geod_init(struct geod_geodesic* g, real a, real f) {
if (!init) Init();
g->a = a;
g->f = f <= 1 ? f : 1/f;
@@ -211,9 +241,9 @@ void GeodesicInit(struct Geodesic* g, real a, real f) {
C4coeff(g);
}
-void GeodesicLineInit(struct GeodesicLine* l,
- const struct Geodesic* g,
- real lat1, real lon1, real azi1, unsigned caps) {
+void geod_lineinit(struct geod_geodesicline* l,
+ const struct geod_geodesic* g,
+ real lat1, real lon1, real azi1, unsigned caps) {
real alp1, cbet1, sbet1, phi, eps;
l->a = g->a;
l->f = g->f;
@@ -221,23 +251,12 @@ void GeodesicLineInit(struct GeodesicLine* l,
l->c2 = g->c2;
l->f1 = g->f1;
/* If caps is 0 assume the standard direct calculation */
- l->caps = (caps ? caps : DISTANCE_IN | LONGITUDE) |
- LATITUDE | AZIMUTH; /* Always allow latitude and azimuth */
+ l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
+ GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
- azi1 = AngNormalize(azi1);
+ /* Guard against underflow in salp0 */
+ azi1 = AngRound(AngNormalize(azi1));
lon1 = AngNormalize(lon1);
- if (lat1 == 90) {
- lon1 += lon1 < 0 ? 180 : -180;
- lon1 = AngNormalize(lon1 - azi1);
- azi1 = -180;
- } else if (lat1 == -90) {
- lon1 = AngNormalize(lon1 + azi1);
- azi1 = 0;
- } else {
- /* Guard against underflow in salp0 */
- azi1 = AngRound(azi1);
- }
-
l->lat1 = lat1;
l->lon1 = lon1;
l->azi1 = azi1;
@@ -245,7 +264,7 @@ void GeodesicLineInit(struct GeodesicLine* l,
alp1 = azi1 * degree;
/* Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing
* problems directly than to skirt them. */
- l->salp1 = azi1 == -180 ? 0 : sin(alp1);
+ l->salp1 = azi1 == -180 ? 0 : sin(alp1);
l->calp1 = fabs(azi1) == 90 ? 0 : cos(alp1);
phi = lat1 * degree;
/* Ensure cbet1 = +epsilon at poles */
@@ -312,12 +331,12 @@ void GeodesicLineInit(struct GeodesicLine* l,
}
}
-real GenPosition(const struct GeodesicLine* l,
- boolx arcmode, real s12_a12,
- real* plat2, real* plon2, real* pazi2,
- real* ps12, real* pm12,
- real* pM12, real* pM21,
- real* pS12) {
+real geod_genposition(const struct geod_geodesicline* l,
+ boolx arcmode, real s12_a12,
+ real* plat2, real* plon2, real* pazi2,
+ real* ps12, real* pm12,
+ real* pM12, real* pM21,
+ real* pS12) {
real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
m12 = 0, M12 = 0, M21 = 0, S12 = 0;
/* Avoid warning about uninitialized B12. */
@@ -325,16 +344,17 @@ real GenPosition(const struct GeodesicLine* l,
real omg12, lam12, lon12;
real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
unsigned outmask =
- (plat2 ? LATITUDE : 0U) |
- (plon2 ? LONGITUDE : 0U) |
- (pazi2 ? AZIMUTH : 0U) |
- (ps12 ? DISTANCE : 0U) |
- (pm12 ? REDUCEDLENGTH : 0U) |
- (pM12 || pM21 ? GEODESICSCALE : 0U) |
- (pS12 ? AREA : 0U);
+ (plat2 ? GEOD_LATITUDE : 0U) |
+ (plon2 ? GEOD_LONGITUDE : 0U) |
+ (pazi2 ? GEOD_AZIMUTH : 0U) |
+ (ps12 ? GEOD_DISTANCE : 0U) |
+ (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+ (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+ (pS12 ? GEOD_AREA : 0U);
outmask &= l->caps & OUT_ALL;
- if (!( TRUE /*Init()*/ && (arcmode || (l->caps & DISTANCE_IN & OUT_ALL)) ))
+ if (!( TRUE /*Init()*/ &&
+ (arcmode || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
/* Uninitialized or impossible distance calculation requested */
return NaN;
@@ -397,7 +417,7 @@ real GenPosition(const struct GeodesicLine* l,
ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
dn2 = sqrt(1 + l->k2 * sq(ssig2));
- if (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) {
+ if (outmask & (GEOD_DISTANCE | GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
if (arcmode || fabs(l->f) > 0.01)
B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
AB1 = (1 + l->A1m1) * (B12 - l->B11);
@@ -417,10 +437,10 @@ real GenPosition(const struct GeodesicLine* l,
omg12 = atan2(somg2 * l->comg1 - comg2 * l->somg1,
comg2 * l->comg1 + somg2 * l->somg1);
- if (outmask & DISTANCE)
+ if (outmask & GEOD_DISTANCE)
s12 = arcmode ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
- if (outmask & LONGITUDE) {
+ if (outmask & GEOD_LONGITUDE) {
lam12 = omg12 + l->A3c *
( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
- l->B31));
@@ -431,31 +451,31 @@ real GenPosition(const struct GeodesicLine* l,
lon2 = AngNormalize(l->lon1 + lon12);
}
- if (outmask & LATITUDE)
+ if (outmask & GEOD_LATITUDE)
lat2 = atan2(sbet2, l->f1 * cbet2) / degree;
- if (outmask & AZIMUTH)
+ if (outmask & GEOD_AZIMUTH)
/* minus signs give range [-180, 180). 0- converts -0 to +0. */
azi2 = 0 - atan2(-salp2, calp2) / degree;
- if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) {
+ if (outmask & (GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
real
B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
AB2 = (1 + l->A2m1) * (B22 - l->B21),
J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
- if (outmask & REDUCEDLENGTH)
+ if (outmask & GEOD_REDUCEDLENGTH)
/* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
* accurate cancellation in the case of coincident points. */
m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
- l->csig1 * csig2 * J12);
- if (outmask & GEODESICSCALE) {
+ if (outmask & GEOD_GEODESICSCALE) {
real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) / (l->dn1 + dn2);
M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
}
}
- if (outmask & AREA) {
+ if (outmask & GEOD_AREA) {
real
B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
real salp12, calp12;
@@ -488,66 +508,66 @@ real GenPosition(const struct GeodesicLine* l,
S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
}
- if (outmask & LATITUDE)
+ if (outmask & GEOD_LATITUDE)
*plat2 = lat2;
- if (outmask & LONGITUDE)
+ if (outmask & GEOD_LONGITUDE)
*plon2 = lon2;
- if (outmask & AZIMUTH)
+ if (outmask & GEOD_AZIMUTH)
*pazi2 = azi2;
- if (outmask & DISTANCE)
+ if (outmask & GEOD_DISTANCE)
*ps12 = s12;
- if (outmask & REDUCEDLENGTH)
+ if (outmask & GEOD_REDUCEDLENGTH)
*pm12 = m12;
- if (outmask & GEODESICSCALE) {
+ if (outmask & GEOD_GEODESICSCALE) {
if (pM12) *pM12 = M12;
if (pM21) *pM21 = M21;
}
- if (outmask & AREA)
+ if (outmask & GEOD_AREA)
*pS12 = S12;
return arcmode ? s12_a12 : sig12 / degree;
}
-void Position(const struct GeodesicLine* l, real s12,
- real* plat2, real* plon2, real* pazi2) {
- GenPosition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
+void geod_position(const struct geod_geodesicline* l, real s12,
+ real* plat2, real* plon2, real* pazi2) {
+ geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
}
-real GenDirect(const struct Geodesic* g,
- real lat1, real lon1, real azi1,
- boolx arcmode, real s12_a12,
- real* plat2, real* plon2, real* pazi2,
- real* ps12, real* pm12, real* pM12, real* pM21,
- real* pS12) {
- struct GeodesicLine l;
+real geod_gendirect(const struct geod_geodesic* g,
+ real lat1, real lon1, real azi1,
+ boolx arcmode, real s12_a12,
+ real* plat2, real* plon2, real* pazi2,
+ real* ps12, real* pm12, real* pM12, real* pM21,
+ real* pS12) {
+ struct geod_geodesicline l;
unsigned outmask =
- (plat2 ? LATITUDE : 0U) |
- (plon2 ? LONGITUDE : 0U) |
- (pazi2 ? AZIMUTH : 0U) |
- (ps12 ? DISTANCE : 0U) |
- (pm12 ? REDUCEDLENGTH : 0U) |
- (pM12 || pM21 ? GEODESICSCALE : 0U) |
- (pS12 ? AREA : 0U);
-
- GeodesicLineInit(&l, g, lat1, lon1, azi1,
- /* Automatically supply DISTANCE_IN if necessary */
- outmask | (arcmode ? NONE : DISTANCE_IN));
- return GenPosition(&l, arcmode, s12_a12,
- plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
+ (plat2 ? GEOD_LATITUDE : 0U) |
+ (plon2 ? GEOD_LONGITUDE : 0U) |
+ (pazi2 ? GEOD_AZIMUTH : 0U) |
+ (ps12 ? GEOD_DISTANCE : 0U) |
+ (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+ (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+ (pS12 ? GEOD_AREA : 0U);
+
+ geod_lineinit(&l, g, lat1, lon1, azi1,
+ /* Automatically supply GEOD_DISTANCE_IN if necessary */
+ outmask | (arcmode ? GEOD_NONE : GEOD_DISTANCE_IN));
+ return geod_genposition(&l, arcmode, s12_a12,
+ plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
}
-void Direct(const struct Geodesic* g,
- real lat1, real lon1, real azi1,
- real s12,
- real* plat2, real* plon2, real* pazi2) {
- GenDirect(g, lat1, lon1, azi1, FALSE, s12, plat2, plon2, pazi2,
- 0, 0, 0, 0, 0);
+void geod_direct(const struct geod_geodesic* g,
+ real lat1, real lon1, real azi1,
+ real s12,
+ real* plat2, real* plon2, real* pazi2) {
+ geod_gendirect(g, lat1, lon1, azi1, FALSE, s12, plat2, plon2, pazi2,
+ 0, 0, 0, 0, 0);
}
-real GenInverse(const struct Geodesic* g,
- real lat1, real lon1, real lat2, real lon2,
- real* ps12, real* pazi1, real* pazi2,
- real* pm12, real* pM12, real* pM21, real* pS12) {
+real geod_geninverse(const struct geod_geodesic* g,
+ real lat1, real lon1, real lat2, real lon2,
+ real* ps12, real* pazi1, real* pazi2,
+ real* pm12, real* pM12, real* pM21, real* pS12) {
real s12 = 0, azi1 = 0, azi2 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
real lon12;
int latsign, lonsign, swapp;
@@ -560,23 +580,22 @@ real GenInverse(const struct Geodesic* g,
real omg12 = 0;
unsigned outmask =
- (ps12 ? DISTANCE : 0U) |
- (pazi1 || pazi2 ? AZIMUTH : 0U) |
- (pm12 ? REDUCEDLENGTH : 0U) |
- (pM12 || pM21 ? GEODESICSCALE : 0U) |
- (pS12 ? AREA : 0U);
+ (ps12 ? GEOD_DISTANCE : 0U) |
+ (pazi1 || pazi2 ? GEOD_AZIMUTH : 0U) |
+ (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
+ (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
+ (pS12 ? GEOD_AREA : 0U);
outmask &= OUT_ALL;
- lon12 = AngNormalize(AngNormalize(lon2) -
- AngNormalize(lon1));
- /* If very close to being on the same meridian, then make it so.
- * Not sure this is necessary... */
+ /* Compute longitude difference (AngDiff does this carefully). Result is
+ * in [-180, 180] but -180 is only for west-going geodesics. 180 is for
+ * east-going and meridional geodesics. */
+ lon12 = AngDiff(AngNormalize(lon1), AngNormalize(lon2));
+ /* If very close to being on the same half-meridian, then make it so. */
lon12 = AngRound(lon12);
/* Make longitude difference positive. */
lonsign = lon12 >= 0 ? 1 : -1;
lon12 *= lonsign;
- if (lon12 == 180)
- lonsign = 1;
/* If really close to the equator, treat as on equator. */
lat1 = AngRound(lat1);
lat2 = AngRound(lat2);
@@ -637,7 +656,6 @@ real GenInverse(const struct Geodesic* g,
slam12 = lon12 == 180 ? 0 : sin(lam12);
clam12 = cos(lam12); /* lon12 == 90 isn't interesting */
-
meridian = lat1 == -90 || slam12 == 0;
if (meridian) {
@@ -660,7 +678,7 @@ real GenInverse(const struct Geodesic* g,
real dummy;
Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
cbet1, cbet2, &s12x, &m12x, &dummy,
- (outmask & GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
+ (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
}
/* Add the check for sig12 since zero length geodesics might yield m12 <
* 0. Test case was
@@ -688,7 +706,7 @@ real GenInverse(const struct Geodesic* g,
s12x = g->a * lam12;
sig12 = omg12 = lam12 / g->f1;
m12x = g->b * sin(sig12);
- if (outmask & GEODESICSCALE)
+ if (outmask & GEOD_GEODESICSCALE)
M12 = M21 = cos(sig12);
a12 = lon12 / g->f1;
@@ -708,7 +726,7 @@ real GenInverse(const struct Geodesic* g,
real dnm = (dn1 + dn2) / 2;
s12x = sig12 * g->b * dnm;
m12x = sq(dnm) * g->b * sin(sig12 / dnm);
- if (outmask & GEODESICSCALE)
+ if (outmask & GEOD_GEODESICSCALE)
M12 = M21 = cos(sig12 / dnm);
a12 = sig12 / degree;
omg12 = lam12 / (g->f1 * dnm);
@@ -739,13 +757,13 @@ real GenInverse(const struct Geodesic* g,
&eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a)
- lam12);
/* 2 * tol0 is approximately 1 ulp for a number in [0, pi]. */
- if (tripb || fabs(v) < (tripn ? 8 : 2) * tol0) break;
+ /* Reversed test to allow escape with NaNs */
+ if (tripb || !(fabs(v) >= (tripn ? 8 : 2) * tol0)) break;
/* Update bracketing values */
- if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b)) {
- salp1b = salp1; calp1b = calp1;
- } else if (numit > maxit1 || calp1/salp1 < calp1a/salp1a) {
- salp1a = salp1; calp1a = calp1;
- }
+ if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
+ { salp1b = salp1; calp1b = calp1; }
+ else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))
+ { salp1a = salp1; calp1a = calp1; }
if (numit < maxit1 && dv > 0) {
real
dalp1 = -v/dv;
@@ -782,7 +800,7 @@ real GenInverse(const struct Geodesic* g,
real dummy;
Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
cbet1, cbet2, &s12x, &m12x, &dummy,
- (outmask & GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
+ (outmask & GEOD_GEODESICSCALE) != 0U, &M12, &M21, C1a, C2a);
}
m12x *= g->b;
s12x *= g->b;
@@ -791,13 +809,13 @@ real GenInverse(const struct Geodesic* g,
}
}
- if (outmask & DISTANCE)
+ if (outmask & GEOD_DISTANCE)
s12 = 0 + s12x; /* Convert -0 to 0 */
- if (outmask & REDUCEDLENGTH)
+ if (outmask & GEOD_REDUCEDLENGTH)
m12 = 0 + m12x; /* Convert -0 to 0 */
- if (outmask & AREA) {
+ if (outmask & GEOD_AREA) {
real
/* From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) */
salp0 = salp1 * cbet1,
@@ -860,44 +878,46 @@ real GenInverse(const struct Geodesic* g,
if (swapp < 0) {
swapx(&salp1, &salp2);
swapx(&calp1, &calp2);
- if (outmask & GEODESICSCALE)
+ if (outmask & GEOD_GEODESICSCALE)
swapx(&M12, &M21);
}
salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
- if (outmask & AZIMUTH) {
+ if (outmask & GEOD_AZIMUTH) {
/* minus signs give range [-180, 180). 0- converts -0 to +0. */
azi1 = 0 - atan2(-salp1, calp1) / degree;
azi2 = 0 - atan2(-salp2, calp2) / degree;
}
- if (outmask & DISTANCE)
+ if (outmask & GEOD_DISTANCE)
*ps12 = s12;
- if (outmask & AZIMUTH) {
+ if (outmask & GEOD_AZIMUTH) {
if (pazi1) *pazi1 = azi1;
if (pazi2) *pazi2 = azi2;
}
- if (outmask & REDUCEDLENGTH)
+ if (outmask & GEOD_REDUCEDLENGTH)
*pm12 = m12;
- if (outmask & GEODESICSCALE) {
+ if (outmask & GEOD_GEODESICSCALE) {
if (pM12) *pM12 = M12;
if (pM21) *pM21 = M21;
}
- if (outmask & AREA)
+ if (outmask & GEOD_AREA)
*pS12 = S12;
/* Returned value in [0, 180] */
return a12;
}
-void Inverse(const struct Geodesic* g,
- double lat1, double lon1, double lat2, double lon2,
- double* ps12, double* pazi1, double* pazi2) {
- GenInverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
+void geod_inverse(const struct geod_geodesic* g,
+ double lat1, double lon1, double lat2, double lon2,
+ double* ps12, double* pazi1, double* pazi2) {
+ geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
}
+/** @cond SKIP */
+
real SinCosSeries(boolx sinp,
real sinx, real cosx,
const real c[], int n) {
@@ -922,7 +942,7 @@ real SinCosSeries(boolx sinp,
: cosx * (y0 - y1); /* cos(x) * (y0 - y1) */
}
-void Lengths(const struct Geodesic* g,
+void Lengths(const struct geod_geodesic* g,
real eps, real sig12,
real ssig1, real csig1, real dn1,
real ssig2, real csig2, real dn2,
@@ -1019,7 +1039,7 @@ real Astroid(real x, real y) {
return k;
}
-real InverseStart(const struct Geodesic* g,
+real InverseStart(const struct geod_geodesic* g,
real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real lam12,
@@ -1187,7 +1207,7 @@ real InverseStart(const struct Geodesic* g,
return sig12;
}
-real Lambda12(const struct Geodesic* g,
+real Lambda12(const struct geod_geodesic* g,
real sbet1, real cbet1, real dn1,
real sbet2, real cbet2, real dn2,
real salp1, real calp1,
@@ -1286,7 +1306,7 @@ real Lambda12(const struct Geodesic* g,
return lam12;
}
-real A3f(const struct Geodesic* g, real eps) {
+real A3f(const struct geod_geodesic* g, real eps) {
/* Evaluate sum(A3x[k] * eps^k, k, 0, nA3x-1) by Horner's method */
real v = 0;
int i;
@@ -1295,7 +1315,7 @@ real A3f(const struct Geodesic* g, real eps) {
return v;
}
-void C3f(const struct Geodesic* g, real eps, real c[]) {
+void C3f(const struct geod_geodesic* g, real eps, real c[]) {
/* Evaluate C3 coeffs by Horner's method
* Elements c[1] thru c[nC3 - 1] are set */
int i, j, k;
@@ -1313,7 +1333,7 @@ void C3f(const struct Geodesic* g, real eps, real c[]) {
}
}
-void C4f(const struct Geodesic* g, real eps, real c[]) {
+void C4f(const struct geod_geodesic* g, real eps, real c[]) {
/* Evaluate C4 coeffs by Horner's method
* Elements c[0] thru c[nC4 - 1] are set */
int i, j, k;
@@ -1404,7 +1424,7 @@ void C2f(real eps, real c[]) {
}
/* The scale factor A3 = mean value of (d/dsigma)I3 */
-void A3coeff(struct Geodesic* g) {
+void A3coeff(struct geod_geodesic* g) {
g->A3x[0] = 1;
g->A3x[1] = (g->n-1)/2;
g->A3x[2] = (g->n*(3*g->n-1)-2)/8;
@@ -1414,7 +1434,7 @@ void A3coeff(struct Geodesic* g) {
}
/* The coefficients C3[l] in the Fourier expansion of B3 */
-void C3coeff(struct Geodesic* g) {
+void C3coeff(struct geod_geodesic* g) {
g->C3x[0] = (1-g->n)/4;
g->C3x[1] = (1-g->n*g->n)/8;
g->C3x[2] = ((3-g->n)*g->n+3)/64;
@@ -1435,7 +1455,7 @@ void C3coeff(struct Geodesic* g) {
/* Generated by Maxima on 2012-10-19 08:02:34-04:00 */
/* The coefficients C4[l] in the Fourier expansion of I4 */
-void C4coeff(struct Geodesic* g) {
+void C4coeff(struct geod_geodesic* g) {
g->C4x[0] = (g->n*(g->n*(g->n*(g->n*(100*g->n+208)+572)+3432)-12012)+30030)/
45045;
g->C4x[1] = (g->n*(g->n*(g->n*(64*g->n+624)-4576)+6864)-3003)/15015;
@@ -1459,3 +1479,41 @@ void C4coeff(struct Geodesic* g) {
g->C4x[19] = -128/(real)(135135);
g->C4x[20] = 128/(real)(99099);
}
+
+int transit(real lon1, real lon2) {
+ real lon12;
+ /* Return 1 or -1 if crossing prime meridian in east or west direction.
+ * Otherwise return zero. */
+ /* Compute lon12 the same way as Geodesic::Inverse. */
+ lon1 = AngNormalize(lon1);
+ lon2 = AngNormalize(lon2);
+ lon12 = AngDiff(lon1, lon2);
+ return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
+ (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
+}
+
+/** @endcond */
+
+void geod_polygonarea(const struct geod_geodesic* g,
+ real lats[], real lons[], int n,
+ real* pA, real* pP) {
+ int i, crossings = 0;
+ real area0 = 4 * pi * g->c2, A = 0, P = 0;
+ for (i = 0; i < n; ++i) {
+ real s12, S12;
+ geod_geninverse(g, lats[i], lons[i], lats[(i + 1) % n], lons[(i + 1) % n],
+ &s12, 0, 0, 0, 0, 0, &S12);
+ P += s12;
+ A -= S12; /* The minus sign is due to the counter-clockwise convention */
+ crossings += transit(lons[i], lons[(i + 1) % n]);
+ }
+ if (crossings & 1)
+ A += (A < 0 ? 1 : -1) * area0/2;
+ /* Put area in (-area0/2, area0/2] */
+ if (A > area0/2)
+ A -= area0;
+ else if (A <= -area0/2)
+ A += area0;
+ if (pA) *pA = A;
+ if (pP) *pP = P;
+}
diff --git a/src/geodesic.h b/src/geodesic.h
index acf5d875..2c49d0dc 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -1,138 +1,123 @@
-/*
- * This is a C implementation of the geodesic algorithms described in
- *
- * C. F. F. Karney,
- * Algorithms for geodesics,
- * J. Geodesy (2012);
- * http://dx.doi.org/10.1007/s00190-012-0578-z
- * Addenda: http://geographiclib.sf.net/geod-addenda.html
+/**
+ * \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 abs(f) < 1/50;
- * * the solution of the inverse problem is always found;
- * * differential and integral properties of geodesics are computed.
+ * - 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 (lat1, lon1) and
- * (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic
- * from point 1 to point 2 has forward azimuths azi1 and azi2 at the two end
- * points.
+ * 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 lat1, lon1, s12, and azi1, determine
- * lat2, lon2, and azi2.
- * * the inverse problem -- given lat1, lon1, lat2, lon2, determine
- * s12, azi1, and azi2.
+ * - 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 a (typically in meters)
- * and flattening f. The routines are accurate to round off with double
- * precision arithmetic provided that abs(f) < 1/50; for the WGS84 ellipsoid,
- * the errors are less than 15 nanometers. (Reasonably accurate results are
- * obtained for abs(f) < 1/5.) Latitudes, longitudes, and azimuths are in
- * degrees. Latitudes must lie in [-90,90] and longitudes and azimuths must
- * lie in [-540,540). The returned values for longitude and azimuths are in
- * [-180,180). The distance s12 is measured in meters (more precisely the same
- * units as a).
+ * 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
- * * SS12 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 (lat1,lon1), (0,lon1), (0,lon2), and
- * (lat2,lon2). It is given in meters^2.
- * * m12, the reduced length of the geodesic is defined such that if the
- * initial azimuth is perturbed by dazi1 (radians) then the second point is
- * displaced by m12 dazi1 in the direction perpendicular to the geodesic.
- * m12 is given in meters. On a curved surface the reduced length obeys a
- * symmetry relation, m12 + m21 = 0. On a flat surface, we have m12 = s12.
- * * MM12 and MM21 are geodesic scales. If two geodesics are parallel at
- * point 1 and separated by a small distance dt, then they are separated by
- * a distance MM12 dt at point 2. MM21 is defined similarly (with the
- * geodesics being parallel to one another at point 2). MM12 and MM21 are
- * dimensionless quantities. On a flat surface, we have MM12 = MM21 = 1.
- * * a12 is the arc length on the auxiliary sphere. This is a construct for
- * converting the problem to one in spherical trigonometry. a12 is
- * measured in degrees. The spherical arc length from one equator crossing
- * to the next is always 180 degrees.
- *
- * Simple interface
- *
- * #include "geodesic.h"
- *
- * double a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12;
- * struct Geodesic g;
- *
- * GeodesicInit(&g, a, f);
- * Direct(&g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2);
- * Inverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
- *
- * GeodesicInit initalizes g for the ellipsoid. Subsequent calls to Direct and
- * Inverse solve the direct and inverse geodesic problems.
- *
- * Returning auxiliary quantities (continued from the previous example):
- *
- * double a12, s12_a12, m12, M12, M21, S12;
- * a12 = GenDirect(&g, lat1, lon1, azi1, 0, s12,
- * &lat1, &lat2, &azi2, 0, &m12, &M12, &M21, &S21);
- * GenDirect(&g, lat1, lon1, azi1, 1, a12,
- * &lat1, &lat2, &azi2, &s12, &m12, &M12, &M21, &S21);
- * a12 = GenInverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2,
- * &m12, &M12, &M21, &S12);
- *
- * GenDirect is a generalized version of Direct allowing the return of the
- * auxiliary quantities. With the first variant (arcmode = 0), the length of
- * the geodesic is specified by the true length s12 and the arc length a12 is
- * returned as the function value. In the second variant (arcmode = 1), the
- * length is specified by the arc length a12 (in degrees), and the true length
- * is returned via &s12.
- *
- * a12 = GenInverse(&g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2,
- * &m12, &M12, &M21, &S12);
- *
- * GenInverse is a generalized version of Inverse allowing the return of the
- * auxiliary quantities.
- *
- * Any of the "return" arguments &s12, etc. in these routines may be replaced
- * by 0, if you do not need some quantities computed.
- *
- * Computing multiple points on a geodesic. This may be accomplished by
- * repeated invocations of Direct varying s12. However, it is more efficient
- * to create a GeodesicLine object, as follows.
- *
- * struct GeodesicLine l;
- * int caps = 0;
+ * - \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;.
*
- * GeodesicLineInit(&l, &g, a, f, lat1, lon1, azi1, caps);
- * Position(l, s12, &lat2, &lon2, &azi2)
+ * 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
*
- * caps is a bit mask specifying the capabilities of the GeodesicLine object.
- * It should be an or'ed combination of
+ * 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.
*
- * LATITUDE compute lat2 (in effect this is always set)
- * LONGITUDE compute lon2
- * AZIMUTH compute azi2 (in effect this is always set)
- * DISTANCE compute s12
- * DISTANCE_IN allow the length to be specified in terms of distance
- * REDUCEDLENGTH compute m12
- * GEODESICSCALE compute M12 and M21
- * AREA compute S12
- * ALL all of the above
+ * 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 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.
*
- * caps = 0 is treated as LATITUDE | LONGITUDE | AZIMUTH | DISTANCE_IN (to
- * support the solution "standard" direct problem).
- *
- * There's also a generalized version of Position
- *
- * a12 = GenPosition(&l, arcmode, s12_a12,
- * &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
- *
- * See the documentation on GenDirect for the meaning of arcmode.
- *
- * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed
+ * 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 file was distributed with GeographicLib 1.27.
- */
+ * This library was distributed with
+ * <a href="../index.html">GeographicLib</a> 1.30.
+ **********************************************************************/
#if !defined(GEODESIC_H)
#define GEODESIC_H 1
@@ -141,62 +126,404 @@
extern "C" {
#endif
- struct Geodesic {
- double a, f, f1, e2, ep2, n, b, c2, etol2;
+ /**
+ * 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 */
};
- struct GeodesicLine {
- double lat1, lon1, azi1;
- double a, f, b, c2, f1, salp0, calp0, k2,
+ /**
+ * 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];
- unsigned caps;
+ /**< @endcond */
+ unsigned caps; /**< the capabilities */
};
- void GeodesicInit(struct Geodesic* g, double a, double f);
- void GeodesicLineInit(struct GeodesicLine* l,
- const struct Geodesic* g,
- double lat1, double lon1, double azi1, unsigned caps);
-
- void Direct(const struct Geodesic* g,
- double lat1, double lon1, double azi1, double s12,
- double* plat2, double* plon2, double* pazi2);
- void Inverse(const struct Geodesic* g,
- double lat1, double lon1, double lat2, double lon2,
- double* ps12, double* pazi1, double* pazi2);
- void Position(const struct GeodesicLine* l, double s12,
- double* plat2, double* plon2, double* pazi2);
-
- double GenDirect(const struct 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);
- double GenInverse(const struct Geodesic* g,
+ /**
+ * 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 [&minus;90&deg;, 90&deg;]; \e azi1 should be in the
+ * range [&minus;540&deg;, 540&deg;).
+ *
+ * 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 [&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 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 = &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 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 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 = &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,
- double* pm12, double* pM12, double* pM21, double* pS12);
- double GenPosition(const struct GeodesicLine* l,
- int arcmode, double s12_a12,
- double* plat2, double* plon2, double* pazi2,
- double* ps12, double* pm12,
- double* pM12, double* pM21,
- double* pS12);
-
- enum mask {
- NONE = 0U,
- LATITUDE = 1U<<7 | 0U,
- LONGITUDE = 1U<<8 | 1U<<3,
- AZIMUTH = 1U<<9 | 0U,
- DISTANCE = 1U<<10 | 1U<<0,
- DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1,
- REDUCEDLENGTH = 1U<<12 | 1U<<0 | 1U<<2,
- GEODESICSCALE = 1U<<13 | 1U<<0 | 1U<<2,
- AREA = 1U<<14 | 1U<<4,
- ALL = 0x7F80U| 0x1FU
+ 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 [&minus;180&deg;, 180&deg;). 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
+ * (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 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
+ * (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 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
+ * (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
+ * 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 (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 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)
diff --git a/src/makefile.vc b/src/makefile.vc
index c9b68cd5..bd042332 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -55,11 +55,12 @@ support = \
pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \
proj_mdist.obj pj_mutex.obj pj_initcache.obj \
pj_ctx.obj pj_log.obj pj_apply_vgridshift.obj
-
-LIBOBJ = $(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc)
+geodesic = geodesic.obj
+LIBOBJ = $(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \
+ $(geodesic)
PROJEXE_OBJ = proj.obj gen_cheb.obj p_series.obj emess.obj
CS2CSEXE_OBJ = cs2cs.obj gen_cheb.obj p_series.obj emess.obj
-GEODEXE_OBJ = geod.obj geod_set.obj geod_interface.obj geodesic.obj emess.obj
+GEODEXE_OBJ = geod.obj geod_set.obj geod_interface.obj emess.obj
PROJ_DLL = proj$(VERSION).dll
PROJ_EXE = proj.exe
CS2CS_EXE = cs2cs.exe
@@ -134,4 +135,5 @@ install: all
copy *.dll $(INSTDIR)\bin
copy *.lib $(INSTDIR)\lib
copy proj_api.h $(INSTDIR)\include
+ copy geodesic.h $(INSTDIR)\include
diff --git a/src/proj.def b/src/proj.def
index f4d60cab..7ae3daaa 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -1,56 +1,65 @@
-VERSION 1.2
+VERSION 1.2
EXPORTS
- pj_init @1
- pj_fwd @2
- pj_inv @3
- pj_free @4
- pj_transform @5
- pj_geocentric_to_geodetic @6
- pj_geodetic_to_geocentric @7
- pj_deallocate_grids @8
- pj_init_plus @9
- pj_latlong_from_proj @10
- pj_is_latlong @11
- pj_get_errno_ref @12
- pj_set_finder @13
- pj_strerrno @14
- pj_errno @15
- pj_get_def @16
- pj_dalloc @17
- pj_is_geocent @18
- pj_get_release @19
- pj_malloc @20
- pj_pr_list @21
- pj_compare_datums @22
- pj_apply_gridshift @23
- pj_datum_transform @24
- pj_set_searchpath @25
- dmstor @26
- pj_get_ellps_ref @27
- pj_get_datums_ref @28
- pj_get_units_ref @29
- pj_get_list_ref @30
+ pj_init @1
+ pj_fwd @2
+ pj_inv @3
+ pj_free @4
+ pj_transform @5
+ pj_geocentric_to_geodetic @6
+ pj_geodetic_to_geocentric @7
+ pj_deallocate_grids @8
+ pj_init_plus @9
+ pj_latlong_from_proj @10
+ pj_is_latlong @11
+ pj_get_errno_ref @12
+ pj_set_finder @13
+ pj_strerrno @14
+ pj_errno @15
+ pj_get_def @16
+ pj_dalloc @17
+ pj_is_geocent @18
+ pj_get_release @19
+ pj_malloc @20
+ pj_pr_list @21
+ pj_compare_datums @22
+ pj_apply_gridshift @23
+ pj_datum_transform @24
+ pj_set_searchpath @25
+ dmstor @26
+ pj_get_ellps_ref @27
+ pj_get_datums_ref @28
+ pj_get_units_ref @29
+ pj_get_list_ref @30
pj_get_prime_meridians_ref @31
- rtodms @32
- set_rtodms @33
- pj_factors @34
- mk_cheby @35
- adjlon @36
- pj_param @37
- pj_ell_set @38
- pj_mkparam @39
- pj_init_ctx @40
- pj_init_plus_ctx @41
- pj_get_default_ctx @42
- pj_get_ctx @43
- pj_set_ctx @44
- pj_ctx_alloc @45
- pj_ctx_free @46
- pj_ctx_get_errno @47
- pj_ctx_set_errno @48
- pj_ctx_set_debug @49
- pj_ctx_set_logger @50
- pj_ctx_set_app_data @51
- pj_ctx_get_app_data @52
- pj_log @53
- pj_clear_initcache @54
+ rtodms @32
+ set_rtodms @33
+ pj_factors @34
+ mk_cheby @35
+ adjlon @36
+ pj_param @37
+ pj_ell_set @38
+ pj_mkparam @39
+ pj_init_ctx @40
+ pj_init_plus_ctx @41
+ pj_get_default_ctx @42
+ pj_get_ctx @43
+ pj_set_ctx @44
+ pj_ctx_alloc @45
+ pj_ctx_free @46
+ pj_ctx_get_errno @47
+ pj_ctx_set_errno @48
+ pj_ctx_set_debug @49
+ pj_ctx_set_logger @50
+ pj_ctx_set_app_data @51
+ pj_ctx_get_app_data @52
+ pj_log @53
+ pj_clear_initcache @54
+ geod_init @55
+ geod_lineinit @56
+ geod_genposition @57
+ geod_position @58
+ geod_gendirect @59
+ geod_direct @60
+ geod_geninverse @61
+ geod_inverse @62
+ geod_polygonarea @63
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 17:10:24 +0000