aboutsummaryrefslogtreecommitdiff
path: root/src/geodesic.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/geodesic.c')
-rw-r--r--src/geodesic.c275
1 files changed, 242 insertions, 33 deletions
diff --git a/src/geodesic.c b/src/geodesic.c
index 957f2144..bd9fc960 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -1,8 +1,12 @@
/**
* \file geodesic.c
* \brief Implementation of the geodesic routines in C
+ *
+ * For the full documentation see geodesic.h.
**********************************************************************/
+/** @cond SKIP */
+
/*
* This is a C implementation of the geodesic algorithms described in
*
@@ -20,8 +24,6 @@
*/
#include "geodesic.h"
-
-/** @cond SKIP */
#include <math.h>
#define GEOGRAPHICLIB_GEODESIC_ORDER 6
@@ -221,9 +223,10 @@ static real A2m1f(real eps);
static void C2f(real eps, real c[]);
static int transit(real lon1, real lon2);
static void accini(real s[]);
+static void acccopy(const real s[], real t[]);
static void accadd(real s[], real y);
-
-/** @endcond */
+static real accsum(const real s[], real y);
+static void accneg(real s[]);
void geod_init(struct geod_geodesic* g, real a, real f) {
if (!init) Init();
@@ -925,16 +928,12 @@ real geod_geninverse(const struct geod_geodesic* g,
}
void geod_inverse(const struct geod_geodesic* g,
- double lat1, double lon1, double lat2, double lon2,
- double* ps12, double* pazi1, double* pazi2) {
+ real lat1, real lon1, real lat2, real lon2,
+ real* ps12, real* pazi1, real* 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) {
+real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
/* Evaluate
* y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
* sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
@@ -1522,6 +1521,11 @@ void accini(real s[]) {
s[0] = s[1] = 0;
}
+void acccopy(const real s[], real t[]) {
+ /* Copy an accumulator; t = s. */
+ t[0] = s[0]; t[1] = s[1];
+}
+
void accadd(real s[], real y) {
/* Add y to an accumulator. */
real u, z = sumx(y, s[1], &u);
@@ -1532,30 +1536,235 @@ void accadd(real s[], real y) {
s[1] = s[1] + u;
}
-/** @endcond */
+real accsum(const real s[], real y) {
+ /* Return accumulator + y (but don't add to accumulator). */
+ real t[2];
+ acccopy(s, t);
+ accadd(t, y);
+ return t[0];
+}
-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[2], P[2];
- accini(A); accini(P);
- for (i = 0; i < n; ++i) {
+void accneg(real s[]) {
+ /* Negate an accumulator. */
+ s[0] = -s[0]; s[1] = -s[1];
+}
+
+void geod_polygon_init(struct geod_polygon* p, boolx polylinep) {
+ p->lat0 = p->lon0 = p->lat = p->lon = NaN;
+ p->polyline = (polylinep != 0);
+ accini(p->P);
+ accini(p->A);
+ p->num = p->crossings = 0;
+}
+
+void geod_polygon_addpoint(const struct geod_geodesic* g,
+ struct geod_polygon* p,
+ real lat, real lon) {
+ lon = AngNormalize(lon);
+ if (p->num == 0) {
+ p->lat0 = p->lat = lat;
+ p->lon0 = p->lon = lon;
+ } else {
+ real s12, S12;
+ geod_geninverse(g, p->lat, p->lon, lat, lon,
+ &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+ accadd(p->P, s12);
+ if (!p->polyline) {
+ accadd(p->A, S12);
+ p->crossings += transit(p->lon, lon);
+ }
+ p->lat = lat; p->lon = lon;
+ }
+ ++p->num;
+}
+
+void geod_polygon_addedge(const struct geod_geodesic* g,
+ struct geod_polygon* p,
+ real azi, real s) {
+ if (p->num) { /* Do nothing is num is zero */
+ real lat, lon, S12;
+ geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+ &lat, &lon, 0,
+ 0, 0, 0, 0, p->polyline ? 0 : &S12);
+ accadd(p->P, s);
+ if (!p->polyline) {
+ accadd(p->A, S12);
+ p->crossings += transit(p->lon, lon);
+ }
+ p->lat = lat; p->lon = lon;
+ ++p->num;
+ }
+}
+
+unsigned geod_polygon_compute(const struct geod_geodesic* g,
+ const struct geod_polygon* p,
+ boolx reverse, boolx sign,
+ real* pA, real* pP) {
+ real s12, S12, t[2], area0;
+ int crossings;
+ if (p->num < 2) {
+ if (pP) *pP = 0;
+ if (!p->polyline && pA) *pA = 0;
+ return p->num;
+ }
+ if (p->polyline) {
+ if (pP) *pP = p->P[0];
+ return p->num;
+ }
+ geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0,
+ &s12, 0, 0, 0, 0, 0, &S12);
+ if (pP) *pP = accsum(p->P, s12);
+ acccopy(p->A, t);
+ accadd(t, S12);
+ crossings = p->crossings + transit(p->lon, p->lon0);
+ area0 = 4 * pi * g->c2;
+ if (crossings & 1)
+ accadd(t, (t[0] < 0 ? 1 : -1) * area0/2);
+ /* area is with the clockwise sense. If !reverse convert to
+ * counter-clockwise convention. */
+ if (!reverse)
+ accneg(t);
+ /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+ if (sign) {
+ if (t[0] > area0/2)
+ accadd(t, -area0);
+ else if (t[0] <= -area0/2)
+ accadd(t, +area0);
+ } else {
+ if (t[0] >= area0)
+ accadd(t, -area0);
+ else if (t[0] < 0)
+ accadd(t, +area0);
+ }
+ if (pA) *pA = 0 + t[0];
+ return p->num;
+}
+
+unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
+ const struct geod_polygon* p,
+ real lat, real lon,
+ boolx reverse, boolx sign,
+ real* pA, real* pP) {
+ real perimeter, tempsum, area0;
+ int crossings, i;
+ unsigned num = p->num + 1;
+ if (num == 1) {
+ if (pP) *pP = 0;
+ if (!p->polyline && pA) *pA = 0;
+ return num;
+ }
+ perimeter = p->P[0];
+ tempsum = p->polyline ? 0 : p->A[0];
+ crossings = p->crossings;
+ for (i = 0; i < (p->polyline ? 1 : 2); ++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);
- accadd(P, s12);
- /* The minus sign is due to the counter-clockwise convention */
- accadd(A, -S12);
- crossings += transit(lons[i], lons[(i + 1) % n]);
+ geod_geninverse(g,
+ i == 0 ? p->lat : lat, i == 0 ? p->lon : lon,
+ i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
+ &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+ perimeter += s12;
+ if (!p->polyline) {
+ tempsum += S12;
+ crossings += transit(i == 0 ? p->lon : lon,
+ i != 0 ? p->lon0 : lon);
+ }
}
+
+ if (pP) *pP = perimeter;
+ if (p->polyline)
+ return num;
+
+ area0 = 4 * pi * g->c2;
if (crossings & 1)
- accadd(A, (A[0] < 0 ? 1 : -1) * area0/2);
- /* Put area in (-area0/2, area0/2] */
- if (A[0] > area0/2)
- accadd(A, -area0);
- else if (A[0] <= -area0/2)
- accadd(A, +area0);
- if (pA) *pA = A[0];
- if (pP) *pP = P[0];
+ tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+ /* area is with the clockwise sense. If !reverse convert to
+ * counter-clockwise convention. */
+ if (!reverse)
+ tempsum *= -1;
+ /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+ if (sign) {
+ if (tempsum > area0/2)
+ tempsum -= area0;
+ else if (tempsum <= -area0/2)
+ tempsum += area0;
+ } else {
+ if (tempsum >= area0)
+ tempsum -= area0;
+ else if (tempsum < 0)
+ tempsum += area0;
+ }
+ if (pA) *pA = 0 + tempsum;
+ return num;
}
+
+unsigned geod_polygon_testedge(const struct geod_geodesic* g,
+ const struct geod_polygon* p,
+ real azi, real s,
+ boolx reverse, boolx sign,
+ real* pA, real* pP) {
+ real perimeter, tempsum, area0;
+ int crossings;
+ unsigned num = p->num + 1;
+ if (num == 1) { /* we don't have a starting point! */
+ if (pP) *pP = NaN;
+ if (!p->polyline && pA) *pA = NaN;
+ return 0;
+ }
+ perimeter = p->P[0] + s;
+ if (p->polyline) {
+ if (pP) *pP = perimeter;
+ return num;
+ }
+
+ tempsum = p->A[0];
+ crossings = p->crossings;
+ {
+ real lat, lon, s12, S12;
+ geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+ &lat, &lon, 0,
+ 0, 0, 0, 0, &S12);
+ tempsum += S12;
+ crossings += transit(p->lon, lon);
+ geod_geninverse(g, lat, lon, p->lat0, p->lon0,
+ &s12, 0, 0, 0, 0, 0, &S12);
+ perimeter += s12;
+ tempsum += S12;
+ crossings += transit(lon, p->lon0);
+ }
+
+ area0 = 4 * pi * g->c2;
+ if (crossings & 1)
+ tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+ /* area is with the clockwise sense. If !reverse convert to
+ * counter-clockwise convention. */
+ if (!reverse)
+ tempsum *= -1;
+ /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+ if (sign) {
+ if (tempsum > area0/2)
+ tempsum -= area0;
+ else if (tempsum <= -area0/2)
+ tempsum += area0;
+ } else {
+ if (tempsum >= area0)
+ tempsum -= area0;
+ else if (tempsum < 0)
+ tempsum += area0;
+ }
+ if (pP) *pP = perimeter;
+ if (pA) *pA = 0 + tempsum;
+ return num;
+}
+
+void geod_polygonarea(const struct geod_geodesic* g,
+ real lats[], real lons[], int n,
+ real* pA, real* pP) {
+ int i;
+ struct geod_polygon p;
+ geod_polygon_init(&p, FALSE);
+ for (i = 0; i < n; ++i)
+ geod_polygon_addpoint(g, &p, lats[i], lons[i]);
+ geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP);
+}
+
+/** @endcond */
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 20:50:39 +0000