From 0ed3296563bf8b1cac6d555248b667392a10adcc Mon Sep 17 00:00:00 2001 From: Wetitpig Date: Mon, 24 Aug 2020 18:18:14 +0800 Subject: [PATCH] Integrate --- .gitignore | 3 +- Makefile | 14 ++-- include/functions.h | 7 +- include/haversine.h | 10 +++ include/vincenty.h | 15 ++++ src/haversine.c | 51 +----------- src/main.c | 133 ++++++++++++++++++++++++++++++++ src/vincenty.c | 183 +++++++++++++++++++------------------------- 8 files changed, 253 insertions(+), 163 deletions(-) create mode 100644 include/haversine.h create mode 100644 include/vincenty.h create mode 100644 src/main.c diff --git a/.gitignore b/.gitignore index 86476ef..e660fd9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1 @@ -haversine -vincenty +bin/ diff --git a/Makefile b/Makefile index 1519c09..f54d1fe 100644 --- a/Makefile +++ b/Makefile @@ -3,15 +3,15 @@ CFLAGS=-O3 -Iinclude LDFLAGS=-lm -all: haversine vincenty +HEADERS = constants.h functions.h haversine.h vincenty.h -haversine: - @mkdir -p bin - $(CC) $(CFLAGS) -o bin/haversine src/haversine.c src/math.c src/io.c $(LDFLAGS) +%.o: src/%.c $(HEADERS) + $(CC) $(CFLAGS) -c -o $@ $< -vincenty: +all: src/haversine.o src/vincenty.o src/main.o src/io.o src/math.o @mkdir -p bin - $(CC) $(CFLAGS) -o bin/vincenty src/vincenty.c src/math.c src/io.c $(LDFLAGS) + $(CC) $(CFLAGS) -o bin/geodesic src/haversine.o src/vincenty.o src/main.o src/io.o src/math.o -lm + @strip bin/geodesic clean: - rm -rf bin + rm -rf src/*.o diff --git a/include/functions.h b/include/functions.h index f5caf74..6835c86 100644 --- a/include/functions.h +++ b/include/functions.h @@ -1,7 +1,12 @@ #include "constants.h" +#ifndef __HAVE_FUNCTIONS_H__ + +#define __HAVE_FUNCTIONS_H__ + int scan(int argc, char *argv, struct Coordinates *loc); -void print(int order, long double s, long double start, long double end); long double sqr(long double operand); long double atan2_modified(long double y, long double x); + +#endif diff --git a/include/haversine.h b/include/haversine.h new file mode 100644 index 0000000..d062017 --- /dev/null +++ b/include/haversine.h @@ -0,0 +1,10 @@ +#include "constants.h" + +#ifndef __HAVE_HAVERSINE_H__ + +#define __HAVE_HAVERSINE_H__ + +long double haversine_distance(struct Coordinates *location); +long double haversine_bearing(struct Coordinates *start, struct Coordinates *end); + +#endif diff --git a/include/vincenty.h b/include/vincenty.h new file mode 100644 index 0000000..841f466 --- /dev/null +++ b/include/vincenty.h @@ -0,0 +1,15 @@ +#include "constants.h" + +#ifndef __HAVE_VINCENTY_H__ + +#define __HAVE_VINCENTY_H__ + +struct vincenty_result { + long double distance; + long double start; + long double end; +}; + +void vincenty(struct vincenty_result *result, struct Coordinates *location); + +#endif diff --git a/src/haversine.c b/src/haversine.c index 361f4b0..cd438c7 100644 --- a/src/haversine.c +++ b/src/haversine.c @@ -1,11 +1,9 @@ #include -#include -#include #include "constants.h" - #include "functions.h" +#include "haversine.h" -long double distance(struct Coordinates *location) +long double haversine_distance(struct Coordinates *location) { long double londiff, latdiff; long double a, b; @@ -20,7 +18,7 @@ long double distance(struct Coordinates *location) return 2 * atan2_modified(sqrtl(a), sqrtl(1-a)) * RADIUS; } -long double bearing(struct Coordinates *start, struct Coordinates *end) +long double haversine_bearing(struct Coordinates *start, struct Coordinates *end) { long double londiff, y, x; @@ -29,46 +27,3 @@ long double bearing(struct Coordinates *start, struct Coordinates *end) x = cosl(start->lat) * sinl(end->lat) - sinl(start->lat) * cosl(end->lat) * cosl(londiff); return atan2_modified(y, x) / (RAD); } - -int main(int argc, char **argv) -{ - if (argc == 2) { - printf("Usage:\n\t%s [coordinate 1] [coordinate 2] ...\n", argv[0]); - return 1; - } - - struct Coordinates *location = malloc(sizeof(struct Coordinates) * 2); - - long double c, total = 0, start, end; - int i, j; - - scan(argc, argv[1], location); - - location->lat = location->lat * RAD; - location->lon = location->lon * RAD; - - puts("{"); - - for (i = 2; i < argc || argc == 1; ++i) { - if (scan(argc, argv[i], (location + 1)) == -1) - break; - - (location + 1)->lat = (location + 1)->lat * RAD; - (location + 1)->lon = (location + 1)->lon * RAD; - - c = distance(location); - start = NORMALISE(bearing(location, location + 1)); - end = NORMALISE(bearing(location + 1, location) - 180); - - total += c; - memcpy(location, location + 1, sizeof(struct Coordinates)); - - print(i - 2, c, start, end); - } - - free(location); - - printf(" \"total_distance\": %Lf\n}\n", total); - - return 0; -} diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..742580b --- /dev/null +++ b/src/main.c @@ -0,0 +1,133 @@ +#include +#include +#include +#include +#include + +#include "constants.h" +#include "functions.h" +#include "haversine.h" +#include "vincenty.h" + +void help(char *name) +{ + puts("Usage:"); + printf("\n\t%s [options] [coordinate 1] [coordinate 2] ...\n", name); + printf("\t-h Show usage.\n\t-f [haversine|vincenty] Set formula to haversine or vincenty.\n\t-s Compute distance.\n\t-o Compute azimuths.\n\n"); + puts("More info in README.md"); + return; +} + + +int main(int argc, char **argv) +{ + int i, j; + int distance = 0, azimuth = 0; + + while ((i = getopt(argc, argv, "f:soh")) != -1) { + switch (i) + { + case 'f': + if (strcmp(optarg, "haversine") == 0) + j = 1; + else if (strcmp(optarg, "vincenty") == 0) + j = 2; + else { + puts("Invalid formula. Abort."); + exit(1); + } + break; + case 's': + distance = 1; + break; + case 'o': + azimuth = 1; + break; + case 'h': + help(argv[0]); + exit(1); + break; + default: + exit(1); + break; + } + } + + if (distance == 0 && azimuth == 0) { + puts("Nothing to be shown. Abort."); + exit(1); + } + + struct Coordinates *location = malloc(sizeof(struct Coordinates) * 2); + + long double c, total = 0, start, end; + struct vincenty_result *res; + + scan(optind == argc, argv[optind], location); + + location->lat = location->lat * RAD; + location->lon = location->lon * RAD; + + putchar('{'); + + for (i = 1; optind == argc || (optind + i) < argc; ++i) { + if (optind != argc) + scan(argc, argv[optind + i], location + 1); + else if (scan(1, NULL, location + 1) == -1) + break; + + (location + 1)->lat = (location + 1)->lat * RAD; + (location + 1)->lon = (location + 1)->lon * RAD; + + if (azimuth == 1 && i != 1) + printf(","); + + printf("\n \"%d\": {\n", i - 1); + + switch (j) + { + case 1: + if (distance == 1) + c = haversine_distance(location); + if (azimuth == 1) { + start = NORMALISE(haversine_bearing(location, location + 1)); + end = NORMALISE(haversine_bearing(location + 1, location) - 180); + } + break; + case 2: + res = malloc(sizeof(struct vincenty_result)); + vincenty(res, location); + c = res->distance; + start = res->start; + end = res->end; + free(res); + break; + } + + if (distance == 1) { + total += c; + printf(" \"distance\": %Lf", c); + if (azimuth == 1) + printf(",\n"); + else + printf("\n },"); + } + + if (azimuth == 1) + printf(" \"start_azimuth\": %Lf,\n \"end_azimuth\": %Lf\n }", start, end); + + memcpy(location, location + 1, sizeof(struct Coordinates)); + } + + free(location); + + if (distance == 1) { + if (azimuth == 1 && total != 0) + printf(","); + printf("\n \"total_distance\": %Lf\n}\n", total); + } + else + printf("\n}\n"); + + exit(0); +} diff --git a/src/vincenty.c b/src/vincenty.c index 3bece35..530487d 100644 --- a/src/vincenty.c +++ b/src/vincenty.c @@ -1,142 +1,115 @@ -#include -#include #include #include "constants.h" - #include "functions.h" +#include "vincenty.h" -int main(int argc, char **argv) +void vincenty(struct vincenty_result *result, struct Coordinates *location) { - if (argc == 2) { - printf("Usage:\n\t%s [coordinate 1] [coordinate 2] ...\n", argv[0]); - return 1; - } - - struct Coordinates *location = malloc(sizeof(struct Coordinates) * 2); long double londiff, lambda, oldvalue, u1, u2; long double ssig, csig, sig, salp, calp, cos2, C; long double usq, k1, a, b, dsig, s; long double start, end; long double d = 0; - long double total = 0; - scan(argc, argv[1], location); + londiff = (location + 1)->lon - location->lon; + lambda = londiff; - location->lat = location->lat * RAD; - location->lon = location->lon * RAD; - puts("{"); + if (fabsl(location->lat) == M_PI_2) + u1 = location->lat > 0 ? M_PI_2 : -M_PI_2; + else + u1 = atanl((1 - FLAT) * tanl(location->lat)); - for (int i = 2; i < argc || argc == 1; ++i) { - if (scan(argc, argv[i], (location + 1)) == -1) - break; + if (fabsl((location + 1)->lat) == M_PI_2) + u2 = (location + 1)->lat > 0 ? M_PI_2 : -M_PI_2; + else + u2 = atanl((1 - FLAT) * tanl((location + 1)->lat)); - (location + 1)->lat = (location + 1)->lat * RAD; - (location + 1)->lon = (location + 1)->lon * RAD; - londiff = (location + 1)->lon - location->lon; - lambda = londiff; + do { + oldvalue = lambda; - if (fabsl(location->lat) == M_PI_2) - u1 = location->lat > 0 ? M_PI_2 : -M_PI_2; - else - u1 = atanl((1 - FLAT) * tanl(location->lat)); - - if (fabsl((location + 1)->lat) == M_PI_2) - u2 = (location + 1)->lat > 0 ? M_PI_2 : -M_PI_2; - else - u2 = atanl((1 - FLAT) * tanl((location + 1)->lat)); - - do { - oldvalue = lambda; + ssig = sqrtl(sqr(cosl(u2) * sinl(lambda)) + sqr(cosl(u1) * sinl(u2) - sinl(u1) * cosl(u2) * cosl(lambda))); + csig = sinl(u1) * sinl(u2) + cosl(u1) * cosl(u2) * cosl(lambda); - ssig = sqrtl(sqr(cosl(u2) * sinl(lambda)) + sqr(cosl(u1) * sinl(u2) - sinl(u1) * cosl(u2) * cosl(lambda))); - csig = sinl(u1) * sinl(u2) + cosl(u1) * cosl(u2) * cosl(lambda); - - sig = atan2_modified(ssig, csig); - - salp = cosl(u1) * cosl(u2) * sinl(lambda) / ssig; - calp = 1 - sqr(salp); - cos2 = csig - 2 * sinl(u1) * sinl(u2) / calp; - if (cos2 < -1 || isnan(cos2)) - cos2 = -1; - - C = FLAT / 16 * calp * (4 + FLAT * (4 - 3 * calp)); + sig = atan2_modified(ssig, csig); - lambda = londiff + (1 - C) * FLAT * salp * (sig + C * ssig * (cos2 + C * csig * (2 * sqr(cos2) - 1))); + salp = cosl(u1) * cosl(u2) * sinl(lambda) / ssig; + calp = 1 - sqr(salp); + cos2 = csig - 2 * sinl(u1) * sinl(u2) / calp; + if (cos2 < -1 || isnan(cos2)) + cos2 = -1; - if (fabsl(lambda) > M_PI) { - d = 1; - break; - } - } while (fabsl(oldvalue - lambda) >= powl(10,-15)); + C = FLAT / 16 * calp * (4 + FLAT * (4 - 3 * calp)); - if (d == 1) { - londiff = (londiff > 0 ? M_PI : -M_PI) - londiff; - lambda = 0; + lambda = londiff + (1 - C) * FLAT * salp * (sig + C * ssig * (cos2 + C * csig * (2 * sqr(cos2) - 1))); - calp = 0.5; - cos2 = 0; + if (fabsl(lambda) > M_PI) { + d = 1; + break; + } + } while (fabsl(oldvalue - lambda) >= powl(10,-15)); - sig = M_PI - fabsl(u1 + u2); - ssig = sinl(sig); - csig = cosl(sig); + if (d == 1) { + londiff = (londiff > 0 ? M_PI : -M_PI) - londiff; + lambda = 0; - do { - oldvalue = salp; + calp = 0.5; + cos2 = 0; - C = FLAT / 16 * calp * (4 + FLAT * (4 - 3 * calp)); - cos2 = csig - 2 * sinl(u1) * sinl(u2) / calp; + sig = M_PI - fabsl(u1 + u2); + ssig = sinl(sig); + csig = cosl(sig); - d = (1 - C) * FLAT * (sig + C * ssig * (cos2 + C * csig * (2 * sqr(cos2) - 1))); - salp = (londiff - asinl(lambda)) / d; - calp = 1 - sqr(salp); + do { + oldvalue = salp; + C = FLAT / 16 * calp * (4 + FLAT * (4 - 3 * calp)); + cos2 = csig - 2 * sinl(u1) * sinl(u2) / calp; - lambda = salp * ssig / (cosl(u1) * cosl(u2)); + d = (1 - C) * FLAT * (sig + C * ssig * (cos2 + C * csig * (2 * sqr(cos2) - 1))); + salp = (londiff - asinl(lambda)) / d; + calp = 1 - sqr(salp); - ssig = sqr(cosl(u2) * lambda) + sqr(cosl(u1) * sinl(u2) + sinl(u1) * cosl(u2) * sqrtl(1 - sqr(lambda))); - csig = -sqrtl(1 - ssig); - ssig = sqrtl(ssig); - sig = M_PI - asinl(ssig); + lambda = salp * ssig / (cosl(u1) * cosl(u2)); - } while (fabsl(oldvalue - salp) >= powl(10,-15)); - } + ssig = sqr(cosl(u2) * lambda) + sqr(cosl(u1) * sinl(u2) + sinl(u1) * cosl(u2) * sqrtl(1 - sqr(lambda))); + csig = -sqrtl(1 - ssig); + ssig = sqrtl(ssig); + sig = M_PI - asinl(ssig); - usq = calp * (sqr(RAD_MAJ) / sqr(RAD_MIN) - 1); - k1 = sqrtl(1 + usq); - k1 = (k1 - 1) / (k1 + 1); - a = (1 + sqr(k1) / 4) / (1 - k1); - b = k1 * (1 - 3 / 8 * sqr(k1)); + } while (fabsl(oldvalue - salp) >= powl(10,-15)); + } - dsig = b * ssig * (cos2 + b / 4 * (csig * (2 * sqr(cos2) - 1) - b / 6 * cos2 * (4 * sqr(ssig) - 3) * (4 * sqr(cos2) - 3))); - s = RAD_MIN * a * (sig - dsig); + usq = calp * (sqr(RAD_MAJ) / sqr(RAD_MIN) - 1); + k1 = sqrtl(1 + usq); + k1 = (k1 - 1) / (k1 + 1); + a = (1 + sqr(k1) / 4) / (1 - k1); + b = k1 * (1 - 3 / 8 * sqr(k1)); - if (fabsl(oldvalue - salp) > fabsl(oldvalue - lambda)) { - a = cosl(u2) * sinl(lambda); - b = cosl(u1) * sinl(u2) - sinl(u1) * cosl(u2) * cosl(lambda); - start = NORMALISE(atan2_modified(a, b) / (RAD)); + dsig = b * ssig * (cos2 + b / 4 * (csig * (2 * sqr(cos2) - 1) - b / 6 * cos2 * (4 * sqr(ssig) - 3) * (4 * sqr(cos2) - 3))); + s = RAD_MIN * a * (sig - dsig); - a = cosl(u1) * sinl(lambda); - b = cosl(u1) * sinl(u2) * cosl(lambda) - sinl(u1) * cosl(u2); - end = NORMALISE(atan2_modified(a, b) / (RAD)); - } - else { - a = salp / cosl(u1); - b = sqrtl(1 - sqr(a)); - if (cosl(u1) * sinl(u2) + sinl(u1) * cosl(u2) * cosl(lambda) < 0) - b = b * -1; - - start = fmodl(atan2_modified(a, b) / (RAD) + 360, 360); - end = fmodl(atan2_modified(salp, -sinl(u1) * ssig + cosl(u1) * csig * b) / (RAD) + 360, 360); - } + if (fabsl(oldvalue - salp) > fabsl(oldvalue - lambda)) { + a = cosl(u2) * sinl(lambda); + b = cosl(u1) * sinl(u2) - sinl(u1) * cosl(u2) * cosl(lambda); + start = NORMALISE(atan2_modified(a, b) / (RAD)); - total += s; - memcpy(location, location + 1, sizeof(struct Coordinates)); - print(i - 2, s, start, end); + a = cosl(u1) * sinl(lambda); + b = cosl(u1) * sinl(u2) * cosl(lambda) - sinl(u1) * cosl(u2); + end = NORMALISE(atan2_modified(a, b) / (RAD)); + } + else { + a = salp / cosl(u1); + b = sqrtl(1 - sqr(a)); + if (cosl(u1) * sinl(u2) + sinl(u1) * cosl(u2) * cosl(lambda) < 0) + b = b * -1; + + start = fmodl(atan2_modified(a, b) / (RAD) + 360, 360); + end = fmodl(atan2_modified(salp, -sinl(u1) * ssig + cosl(u1) * csig * b) / (RAD) + 360, 360); } - free(location); - printf(" \"total_distance\": %Lf\n}\n", total); - - return 0; + result->distance = s; + result->start = start; + result->end = end; + return; }