From aa024f39f281a6770538be8f7b5dcf79f6260e5a Mon Sep 17 00:00:00 2001 From: Wetitpig Date: Fri, 4 Sep 2020 23:48:11 +0800 Subject: [PATCH] Separate meridian-parallel polygons --- Makefile | 2 +- include/geodesic.h | 1 + include/mpblock.h | 9 ++++ src/karney.c | 68 ------------------------ src/main.c | 16 ++++++ src/mpblock.c | 126 +++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 153 insertions(+), 69 deletions(-) create mode 100644 include/mpblock.h create mode 100644 src/mpblock.c diff --git a/Makefile b/Makefile index 8646428..cba6de3 100644 --- a/Makefile +++ b/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) diff --git a/include/geodesic.h b/include/geodesic.h index 314ea4b..499eaa1 100644 --- a/include/geodesic.h +++ b/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)) diff --git a/include/mpblock.h b/include/mpblock.h new file mode 100644 index 0000000..e466117 --- /dev/null +++ b/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 diff --git a/src/karney.c b/src/karney.c index a45de65..e38011a 100644 --- a/src/karney.c +++ b/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; diff --git a/src/main.c b/src/main.c index e810f24..336b2f6 100644 --- a/src/main.c +++ b/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 { diff --git a/src/mpblock.c b/src/mpblock.c new file mode 100644 index 0000000..1877727 --- /dev/null +++ b/src/mpblock.c @@ -0,0 +1,126 @@ +#include +#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; +}