Separate meridian-parallel polygons

Makefile

@@ -11,7 +11,7 @@ HEADERS = constants.h functions.h haversine.h vincenty.h greatcircle.h karney.h
 all: debug
 	@strip bin/geodesic
 
-debug: src/haversine.o src/vincenty.o src/main.o src/io.o src/math.o src/greatcircle.o src/karney.o
+debug: src/haversine.o src/vincenty.o src/main.o src/io.o src/math.o src/greatcircle.o src/karney.o src/mpblock.o
 	@mkdir -p bin
 	$(CC) $(CFLAGS) -o bin/geodesic $^ $(LDFLAGS)
 

include/geodesic.h

@@ -13,6 +13,7 @@
 #define FLAT (1.0l/298.257223563l)
 
 #define RAD_MIN (RAD_MAJ * (1.0l - FLAT))
+#define FLAT_3 ((RAD_MAJ - RAD_MIN) / (RAD_MAJ + RAD_MIN))
 #define ECC (FLAT * (2.0l - FLAT))
 #define ECC2 (ECC / (1.0l - ECC))
 

include/mpblock.h

@@ -0,0 +1,9 @@
+#include "geodesic.h"
+
+#ifndef __HAVE_MPBLOCK_H__
+
+#define __HAVE_MPBLOCK_H__
+
+void mpblock(struct Coordinates *vertex, int i, long double s, long double a, long double *res);
+
+#endif

src/karney.c

@@ -83,84 +83,16 @@ long double I4(long double csig, long double salp)
 	return i;
 }
 
-int isblock(struct Coordinates *vertex)
-{
-	int block = 1;
-	block &= (vertex->lat == (vertex + 1)->lat && (vertex + 2)->lat == (vertex + 3)->lat && (vertex + 1)->lon == (vertex + 2)->lon && vertex->lon == (vertex + 3)->lon);
-	block |= (vertex->lon == (vertex + 1)->lon && (vertex + 2)->lon == (vertex + 3)->lon && (vertex + 1)->lat == (vertex + 2)->lat && vertex->lat == (vertex + 3)->lat);
-	return block;
-}
-
-int ispolariso(struct Coordinates *vertex)
-{
-	int k;
-	for (k = 0; k < 3; k++) {
-		if (fabsl((vertex + k)->lat) / RAD == 90 && (vertex + (k + 1) % 3)->lat == (vertex + (k + 2) % 3)->lat)
-			break;
-	}
-	return k;
-}
-
-long double ellipblock(struct Coordinates *vertex)
-{
-	long double lat[2], lon[2];
-	if (vertex->lat < (vertex + 2)->lat) {
-		lat[0] = vertex->lat;
-		lat[1] = (vertex + 2)->lat;
-	}
-	else {
-		lat[0] = (vertex + 2)->lat;
-		lat[1] = vertex->lat;
-	}
-
-	if (vertex->lon < (vertex + 2)->lon) {
-		lon[0] = vertex->lon;
-		lon[1] = (vertex + 2)->lon;
-	}
-	else {
-		lon[0] = (vertex + 2)->lon;
-		lon[1] = vertex->lon;
-	}
-
-	long double londiff = fabsl(normalise_c(lon[1] - lon[0]));
-	long double latdiff = sinl(lat[1]) / (1.0l - ECC * sqr(sinl(lat[1]))) - sinl(lat[0]) / (1.0l - ECC * sqr(sinl(lat[0])));
-	latdiff += logl((1.0l + sqrtl(ECC) * sinl(lat[1])) * (1.0l - sqrtl(ECC) * sinl(lat[0])) / (1.0l - sqrtl(ECC) * sinl(lat[1])) / (1.0l + sqrtl(ECC) * sinl(lat[0]))) / (2.0l * sqrtl(ECC));
-
-	return sqr(RAD_MIN) * londiff * latdiff / 2.0l;
-}
-
 void karney(struct Coordinates *vertex, int i, int s, int a, long double *res)
 {
 	long double *inter = malloc(sizeof(long double) * 8);
 	long double prev, next, excess = 0;
 	long double area, darea = 0, interarea;
 	long double sig0, sig1;
-
 	long double perimeter = 0;
 
 	int h, k;
 
-	if (a == 1) {
-		if (i == 4 && isblock(vertex)) {
-			area = ellipblock(vertex);
-			a = 0;
-		}
-		else if (i == 3 && (k = ispolariso(vertex)) < 3) {
-			struct Coordinates *triangle = malloc(sizeof(struct Coordinates) * 4);
-
-			triangle->lat = (vertex + k)->lat;
-			triangle->lon = (vertex + ((k + 2) % 3))->lon;
-			memcpy(triangle + 1, vertex + ((k + 1) % 3), sizeof(struct Coordinates));
-			for (h = 2; h < 4; h++) {
-				(triangle + h)->lat = (vertex + ((k + h) % 3))->lat;
-				(triangle + h)->lon = (triangle + (3 - h))->lon;
-			}
-			area = ellipblock(triangle);
-			free(triangle);
-			a = 0;
-		}
-	}
-
 	struct Coordinates *coor0, *coor1;
 	for (h = 0; h < i; h++) {
 		coor0 = vertex + h;

src/main.c

@@ -10,6 +10,7 @@
 #include "vincenty.h"
 #include "greatcircle.h"
 #include "karney.h"
+#include "mpblock.h"
 
 void help(char *name)
 {
@@ -83,6 +84,8 @@ int main(int argc, char **argv)
 				j = 1;
 			else if (strcmp(optarg, "ellipsoid") == 0)
 				j = 2;
+			else if (strcmp(optarg, "block") == 0)
+				j = 3;
 			else
 				error("Invalid formula.");
 			break;
@@ -144,6 +147,9 @@ int main(int argc, char **argv)
 				start = *(res + 1);
 				end = *(res + 2);
 				break;
+				default:
+				error("Formula not available for inverse problem.");
+				break;
 			}
 
 			if (isnan(c) && memcmp(location, location + 1, sizeof(struct Coordinates)) == 0)
@@ -200,6 +206,9 @@ int main(int argc, char **argv)
 				(point + 1)->lon = *(res + 1);
 				end = *(res + 2);
 				break;
+				default:
+				error("Formula not available for direct problem.");
+				break;
 			}
 
 			if ((isnan((point + 1)->lat) || isnan((point + 1)->lon)) && add->s == 0)
@@ -264,6 +273,13 @@ int main(int argc, char **argv)
 				perimeter = *res;
 				area = *(res + 1);
 				break;
+				case 3:
+				mpblock(vertex, i, distance, azimuth, res);
+				perimeter = *res;
+				area = *(res + 1);
+				if (isnan(area) && isnan(perimeter))
+					error("Not a meridian-parallel block.");
+				break;
 			}
 		}
 		else {

src/mpblock.c

@@ -0,0 +1,126 @@
+#include <math.h>
+#include "mpblock.h"
+
+long double double_fac(int x)
+{
+	if (x == -3)
+		return -1;
+	else {
+		long double y = 1;
+		while (x > 0) {
+			y *= x;
+			x -= 2;
+		}
+		return y;
+	}
+}
+
+int isblock(struct Coordinates *vertex)
+{
+	int block = 1;
+	block &= (vertex->lat == (vertex + 1)->lat && (vertex + 2)->lat == (vertex + 3)->lat && (vertex + 1)->lon == (vertex + 2)->lon && vertex->lon == (vertex + 3)->lon);
+	block |= (vertex->lon == (vertex + 1)->lon && (vertex + 2)->lon == (vertex + 3)->lon && (vertex + 1)->lat == (vertex + 2)->lat && vertex->lat == (vertex + 3)->lat);
+	return block;
+}
+
+int ispolariso(struct Coordinates *vertex)
+{
+	int k;
+	for (k = 0; k < 3; k++) {
+		if (fabsl((vertex + k)->lat) / RAD == 90 && (vertex + (k + 1) % 3)->lat == (vertex + (k + 2) % 3)->lat)
+			break;
+	}
+	return k;
+}
+
+long double ellipblock(long double lat0, long double lat1, long double lon0, long double lon1)
+{
+	long double londiff = fabsl(normalise_c(lon1 - lon0));
+	long double latdiff = sinl(lat1) / (1.0l - ECC * sqr(sinl(lat1))) - sinl(lat0) / (1.0l - ECC * sqr(sinl(lat0)));
+	latdiff += logl((1.0l + sqrtl(ECC) * sinl(lat1)) * (1.0l - sqrtl(ECC) * sinl(lat0)) / (1.0l - sqrtl(ECC) * sinl(lat1)) / (1.0l + sqrtl(ECC) * sinl(lat0))) / (2.0l * sqrtl(ECC));
+
+	return sqr(RAD_MIN) * londiff * latdiff / 2.0l;
+}
+
+long double parallel_length(long double lon0, long double lon1, long double lat)
+{
+	return cosl(lat) * fabsl(normalise_c(lon1 - lon0)) / sqrtl(1.0l - ECC * sqr(sinl(lat)));
+}
+
+long double meridian_arc(long double lat0, long double lat1)
+{
+	int k, j;
+	long double c = 0, m = 0;
+
+	for (j = 0; j < 11; j++)
+		c += sqr(double_fac(2 * j - 3) / double_fac(2 * j)) * powl(FLAT_3, 2 * j);
+	m += c * (lat1 - lat0);
+
+	for (k = 1; k < 6; k++) {
+		c = 0;
+		for (j = 0; j < 11; j++)
+			c += double_fac(2 * j - 3) / double_fac(2 * j) * double_fac(2 * j + 2 * k - 3) / double_fac(2 * j + 2 * k) * powl(FLAT_3, k + 2 * j);
+		c /= k;
+		m += powl(-1.0l, k) * c * (sin(2.0l * lat1) - sin(2.0l * lat0));
+	}
+
+	return (RAD_MAJ + RAD_MIN) / 2 * m;
+}
+
+void mpblock(struct Coordinates *vertex, int i, long double s, long double a, long double *res)
+{
+	int h, k;
+
+	if (i == 4 && isblock(vertex)) {
+		long double lat[2], lon[2];
+		if (vertex->lat < (vertex + 2)->lat) {
+			lat[0] = vertex->lat;
+			lat[1] = (vertex + 2)->lat;
+		}
+		else {
+			lat[0] = (vertex + 2)->lat;
+			lat[1] = vertex->lat;
+		}
+
+		if (vertex->lon < (vertex + 2)->lon) {
+			lon[0] = vertex->lon;
+			lon[1] = (vertex + 2)->lon;
+		}
+		else {
+			lon[0] = (vertex + 2)->lon;
+			lon[1] = vertex->lon;
+		}
+
+		if (a == 1)
+			*(res + 1) = ellipblock(lat[0], lat[1], lon[0], lon[1]);
+		if (s == 1) {
+			*res = parallel_length(lat[0], lat[1], lon[0]);
+			*res += parallel_length(lat[0], lat[1], lon[1]);
+			*res += meridian_arc(lat[0], lat[1]) * 2.0l;
+		}
+	}
+	else if (i == 3 && (k = ispolariso(vertex)) < 3) {
+		long double lon[2];
+
+		if ((vertex + k + 1)->lon < (vertex + ((k + 2) % 3))->lon) {
+			lon[0] = (vertex + k + 1)->lon;
+			lon[1] = (vertex + ((k + 2) % 3))->lon;
+		}
+		else {
+			lon[1] = (vertex + k + 1)->lon;
+			lon[0] = (vertex + ((k + 2) % 3))->lon;
+		}
+
+		if (a == 1)
+			*(res + 1) = ellipblock((vertex + k + 1)->lat, (vertex + k)->lat, lon[0], lon[1]);
+		if (s == 1) {
+			*res = parallel_length(lon[0], lon[1], (vertex + k + 1)->lat);
+			*res += meridian_arc((vertex + k + 1)->lat, (vertex + k)->lat) * 2.0l;
+		}
+	}
+	else {
+		*res = NAN;
+		*(res + 1) = NAN;
+	}
+	return;
+}