diff options
34 files changed, 1113 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a071ed4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +*\#* +*~ +*.o +*.a +*.lo +*.su +*.so +*.so.* +*.dll +*.dylib +*.gch +*.gcov +*.gcno +*.gcda @@ -0,0 +1,15 @@ +ISC License + +© 2024 Mattias Andrée <m@maandree.se> + +Permission to use, copy, modify, and/or distribute this software for any +purpose with or without fee is hereby granted, provided that the above +copyright notice and this permission notice appear in all copies. + +THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES +WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR +ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES +WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN +ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF +OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..d79bd06 --- /dev/null +++ b/Makefile @@ -0,0 +1,96 @@ +.POSIX: + +CONFIGFILE = config.mk +include $(CONFIGFILE) + +OS = linux +# Linux: linux +# Mac OS: macos +# Windows: windows +include mk/$(OS).mk + + +LIB_MAJOR = 1 +LIB_MINOR = 0 +LIB_VERSION = $(LIB_MAJOR).$(LIB_MINOR) +LIB_NAME = tellurian + + +OBJ_PUBLIC =\ + libtellurian_sea_level.o\ + libtellurian_sea_level_radians.o\ + libtellurian_coarse_distance.o\ + libtellurian_coarse_distance_radians.o\ + libtellurian_distance.o\ + libtellurian_distance_radians.o\ + libtellurian_azimuth.o\ + libtellurian_azimuth_radians.o\ + libtellurian_end_point.o\ + libtellurian_end_point_radians.o\ + libtellurian_normal_gravity.o\ + libtellurian_normal_gravity_radians.o\ + libtellurian_elevated_gravity.o\ + libtellurian_elevated_gravity_radians.o\ + libtellurian_meridan_radius.o\ + libtellurian_meridan_radius_radians.o\ + libtellurian_transverse_radius.o\ + libtellurian_transverse_radius_radians.o\ + libtellurian_azimuthal_radius.o\ + libtellurian_azimuthal_radius_radians.o\ + libtellurian_gaussian_radius.o\ + libtellurian_gaussian_radius_radians.o + +OBJ =\ + $(OBJ_PUBLIC)\ + libtellurian_vincenty_inverse__.o + +HDR =\ + libtellurian.h\ + common.h + +LOBJ = $(OBJ:.o=.lo) + + +all: libtellurian.a libtellurian.$(LIBEXT) +$(OBJ): $(HDR) +$(LOBJ): $(HDR) + +.c.o: + $(CC) -c -o $@ $< $(CFLAGS) $(CPPFLAGS) + +.c.lo: + $(CC) -fPIC -c -o $@ $< $(CFLAGS) $(CPPFLAGS) + +libtellurian.a: $(OBJ) + @rm -f -- $@ + $(AR) rc $@ $(OBJ) + $(AR) ts $@ > /dev/null + +libtellurian.$(LIBEXT): $(LOBJ) + $(CC) $(LIBFLAGS) -o $@ $(LOBJ) $(LDFLAGS) + +install: libtellurian.a libtellurian.$(LIBEXT) + mkdir -p -- "$(DESTDIR)$(PREFIX)/lib" + mkdir -p -- "$(DESTDIR)$(PREFIX)/include" + cp -- libtellurian.a "$(DESTDIR)$(PREFIX)/lib/" + cp -- libtellurian.$(LIBEXT) "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBMINOREXT)" + $(FIX_INSTALL_NAME) "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBMINOREXT)" + ln -sf -- libtellurian.$(LIBMINOREXT) "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBMAJOREXT)" + ln -sf -- libtellurian.$(LIBMAJOREXT) "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBEXT)" + cp -- libtellurian.h "$(DESTDIR)$(PREFIX)/include/" + +uninstall: + -rm -f -- "$(DESTDIR)$(PREFIX)/lib/libtellurian.a" + -rm -f -- "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBMAJOREXT)" + -rm -f -- "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBMINOREXT)" + -rm -f -- "$(DESTDIR)$(PREFIX)/lib/libtellurian.$(LIBEXT)" + -rm -f -- "$(DESTDIR)$(PREFIX)/include/libtellurian.h" + +clean: + -rm -f -- *.o *.a *.lo *.su *.so *.so.* *.dll *.dylib + -rm -f -- *.gch *.gcov *.gcno *.gcda *.$(LIBEXT) + +.SUFFIXES: +.SUFFIXES: .lo .o .c + +.PHONY: all install uninstall clean @@ -0,0 +1,23 @@ +NAME + libtellurian - Geodesy library + +SYNPOSIS + #include <libtellurian.h> + + Link with -ltellurian -lm. + +DESCRIPTION + libtellurian provides a collection of geodesy functions + and constants. + + libtellurian always uses meters as the length unit, + meters per square-second as the acceleration unit, + and degrees for the angle unit, except when the + function name uses the suffix "_radians", in which + case the function uses radians instead of degrees. + All functions that input or output angles have a + version that uses degress and a version that uses + radians. Unless otherwise specified, functions that + use degress convert the input and output to and from + radians and call the version of the function that + uses radians. @@ -0,0 +1,3 @@ +Add man pages +Finish libtellurian_end_point_radians +Finish libtellurian_vincenty_inverse__ diff --git a/common.h b/common.h new file mode 100644 index 0000000..36ca06f --- /dev/null +++ b/common.h @@ -0,0 +1,19 @@ +/* See LICENSE file for copyright and license details. */ +#include "libtellurian.h" + +#include <stddef.h> +#include <math.h> + + +#if __GNUC__ +# pragma GCC diagnostic ignored "-Wunsuffixed-float-constants" +#endif + + +#define radians(DEG) ((double)M_PI / 180.0 * (DEG)) +#define degrees(RAD) (180.0 / (double)M_PI * (RAD)) +#define haversin(X) (fma(-0.5, cos(X), 0.5)) + + +void libtellurian_vincenty_inverse__(double latitude1, double longitude1, double latitude2, double longitude2, + double *distance_out, double *azimuth1_out, double *azimuth2_out); diff --git a/config.mk b/config.mk new file mode 100644 index 0000000..037fe35 --- /dev/null +++ b/config.mk @@ -0,0 +1,8 @@ +PREFIX = /usr +MANPREFIX = $(PREFIX)/share/man + +CC = c99 + +CPPFLAGS = -D_DEFAULT_SOURCE -D_BSD_SOURCE -D_XOPEN_SOURCE=700 -D_GNU_SOURCE +CFLAGS = +LDFLAGS = -lm diff --git a/libtellurian.h b/libtellurian.h new file mode 100644 index 0000000..f700691 --- /dev/null +++ b/libtellurian.h @@ -0,0 +1,440 @@ +/* See LICENSE file for copyright and license details. */ +#ifndef LIBTELLURIAN_H +#define LIBTELLURIAN_H + + +#if defined(__GNUC__) +# define LIBTELLURIAN_CONST__ __attribute__((__const__, __warn_unused_result__)) +# define LIBTELLURIAN_WUR__ __attribute__((__warn_unused_result__)) +#else +# define LIBTELLURIAN_CONST__ +# define LIBTELLURIAN_WUR__ +#endif + + +/** + * The Earth's equatorial radius, in meters + */ +#define LIBTELLURIAN_EQUATORIAL_RADIUS 6378137.0 /* a */ + +/** + * The Earth's polar radius, in meters + */ +#define LIBTELLURIAN_POLAR_RADIUS 6356752.314245 /* b */ + +/** + * The Earth's mean radius (arithmetic mean), in meters + */ +#define LIBTELLURIAN_MEAN_RADIUS 6371008.771415 /* (2a + b) / 3 */ + +/** + * The Earth's volumetric radius (geometric mean), in meters + */ +#define LIBTELLURIAN_VOLUMETRIC_RADIUS 6371000.7900090935 /* ∛(a² * b) */ + +/** + * The Earth's authalic radius (equal-area mean), in meters + */ +#define LIBTELLURIAN_AUTHALIC_RADIUS 6371007.180918414 /* √(a²/2 + (b² artanh e) / 2e) */ + +/** + * The Earth's rectifying radius, in meters + */ +#define LIBTELLURIAN_RECTIFYING_RADIUS 6367449.14582 /* 2π⁻¹∫{0 → ½π} √(a² cos² φ + b² sin² φ) dφ ≅ ∛(½(√a³+√b³))² */ + +/** + * The Earth's nominal equatorial radius, in meters + */ +#define LIBTELLURIAN_NOMINAL_EQUATORIAL_RADIUS 6378100. + +/** + * The Earth's nominal polar radius, in meters + */ +#define LIBTELLURIAN_NOMINAL_POLAR_RADIUS 6356800. + +/** + * The Earth's nominal radius, in meters + */ +#define LIBTELLURIAN_NOMINAL_RADIUS LIBTELLURIAN_NOMINAL_EQUATORIAL_RADIUS + +/** + * The Earth's equatorial circumference, in meters + */ +#define LIBTELLURIAN_EQUATORIAL_CIRCUMFERENCE 40075016.68557849 /* 2bπ */ + +/** + * The Earth's meridional circumference, in meters + */ +#define LIBTELLURIAN_MERIDIONAL_CIRCUMFERENCE 39940652.74224401 /* 2bπ */ + +/** + * The Earth's mean circumference (arithmetic mean), in meters + */ +#define LIBTELLURIAN_MEAN_CIRCUMFERENCE 40030228.70446699 /* 2π(2a + b) / 3 */ + +/** + * The Earth's volumetric circumference (geometric mean), in meters + */ +#define LIBTELLURIAN_VOLUMETRIC_CIRCUMFERENCE 40030178.555814676 /* 2π∛(a² * b) */ + +/** + * The Earth's authalic circumference (equal-area mean), in meters + */ +#define LIBTELLURIAN_AUTHALIC_CIRCUMFERENCE 40030218.71108221 /* 2π√(a²/2 + (b² artanh e) / 2e) */ + +/** + * The Earth's rectifying circumference, in meters + */ +#define LIBTELLURIAN_RECTIFYING_CIRCUMFERENCE 40007862.91722943 /* 4 ∫{0 → ½π} √(a² cos² φ + b² sin² φ) dφ */ + +/** + * The Earth's gravity at the equator, in meters per square-second + */ +#define LIBTELLURIAN_EQUATORIAL_GRAVITY 9.7803253359 + +/** + * The Earth's gravity at the poles, in meters per square-second + */ +#define LIBTELLURIAN_POLAR_GRAVITY 9.8321849378 + +/** + * The Earth's normal gravity at the equator, in meters per square-second + */ +#define LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY 9.7803267715 + +/** + * The Earth's normal gravity at the poles, in meters per square-second + */ +#define LIBTELLURIAN_NORMAL_POLAR_GRAVITY 9.8321863685 + +/** + * The Earth's mass, in kilograms + */ +#define LIBTELLURIAN_MASS_OF_EARTH 5.972168e24 + + +/** + * Calculate the distance of the sea level (geocentric radius), + * for some point on the Earth's surface, from the centre of + * the Earth + * + * @param latitude GPS latitude coordinate, in degrees + * @return The geocentric altitude of the sea level, in meters + */ +LIBTELLURIAN_CONST__ +double libtellurian_sea_level(double latitude); + +/** + * Calculate the distance of the sea level (geocentric radius), + * for some point on the Earth's surface, from the centre of + * the Earth + * + * @param latitude GPS latitude coordinate, in radians + * @return The geocentric altitude of the sea level, in meters + */ +LIBTELLURIAN_CONST__ +double libtellurian_sea_level_radians(double latitude); + +/** + * Calculate the distance between two points on the Earth's surface + * + * This function is gives an approximate value using + * a sphere as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in degrees + * @param longitude1 GPS longitude coordinate for the first point, in degrees + * @param latitude2 GPS latitude coordinate for the second point, in degrees + * @param longitude2 GPS longitude coordinate for the second point, in degrees + * @return Approximate distance between the two points + */ +LIBTELLURIAN_CONST__ +double libtellurian_coarse_distance(double latitude1, double longitude1, + double latitude2, double longitude2); + +/** + * Calculate the distance between two points on the Earth's surface + * + * This function is gives an approximate value using + * a sphere as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in radians + * @param longitude1 GPS longitude coordinate for the first point, in radians + * @param latitude2 GPS latitude coordinate for the second point, in radians + * @param longitude2 GPS longitude coordinate for the second point, in radians + * @return Approximate distance between the two points + */ +LIBTELLURIAN_CONST__ +double libtellurian_coarse_distance_radians(double latitude1, double longitude1, + double latitude2, double longitude2); + +/** + * Calculate the distance and azimuths between two points on the Earth's surface + * + * This function is gives good approximate values + * using an oblate spheroid as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in degrees + * @param longitude1 GPS longitude coordinate for the first point, in degrees + * @param latitude2 GPS latitude coordinate for the second point, in degrees + * @param longitude2 GPS longitude coordinate for the second point, in degrees + * @param azimuth1_out Output parameter for the forward azimuth, in degrees, + * at the first point; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth, in degrees, + * at the second point; or `NULL` + * @return Approximate distance between the two points + * + * Calculating the the forward azimuths is not free, but it + * is cheap to compute them (especially the first one) when + * most of the computations for the distance have been + * performed. If you have no need for an azimuth you can set + * the corresponding output parameter to `NULL`, and the + * function will not compute it. + */ +LIBTELLURIAN_WUR__ +double libtellurian_distance(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out); + +/** + * Calculate the distance and azimuths between two points on the Earth's surface + * + * This function is gives good approximate values + * using an oblate spheroid as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in radians + * @param longitude1 GPS longitude coordinate for the first point, in radians + * @param latitude2 GPS latitude coordinate for the second point, in radians + * @param longitude2 GPS longitude coordinate for the second point, in radians + * @param azimuth1_out Output parameter for the forward azimuth, in radians, + * at the first point; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth, in radians, + * at the second point; or `NULL` + * @return Approximate distance between the two points + * + * Calculating the the forward azimuths is not free, but it + * is cheap to compute them (especially the first one) when + * most of the computations for the distance have been + * performed. If you have no need for an azimuth you can set + * the corresponding output parameter to `NULL`, and the + * function will not compute it. + */ +LIBTELLURIAN_WUR__ +double libtellurian_distance_radians(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out); + +/** + * Calculate the azimuths between two points on the Earth's surface + * + * This function is gives good approximate values + * using an oblate spheroid as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in degrees + * @param longitude1 GPS longitude coordinate for the first point, in degrees + * @param latitude2 GPS latitude coordinate for the second point, in degrees + * @param longitude2 GPS longitude coordinate for the second point, in degrees + * @param azimuth1_out Output parameter for the forward azimuth, in degrees, + * at the first point; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth, in degrees, + * at the second point; or `NULL` + * + * This function is identical to libtellurian_distance` + * except it does not compute the distance between the + * points + */ +void libtellurian_azimuth(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out); + +/** + * Calculate the azimuths between two points on the Earth's surface + * + * This function is gives good approximate values + * using an oblate spheroid as a model for the Earth + * + * @param latitude1 GPS latitude coordinate for the first point, in radians + * @param longitude1 GPS longitude coordinate for the first point, in radians + * @param latitude2 GPS latitude coordinate for the second point, in radians + * @param longitude2 GPS longitude coordinate for the second point, in radians + * @param azimuth1_out Output parameter for the forward azimuth, in radians, + * at the first point; or `NULL` + * @param azimuth2_out Output parameter for the forward azimuth, in radians, + * at the second point; or `NULL` + * + * + * This function is identical to libtellurian_distance_radians` + * except it does not compute the distance between the points + */ +void libtellurian_azimuth_radians(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out); + +/** + * Calculate the location that is some distance away, + * in some direction, from some point + * + * @param latitude1 GPS latitude coordinate for the starting point, in degrees + * @param longitude1 GPS longitude coordinate for the starting point, in degrees + * @param azimuth1 The direction to travel, in degrees + * @param distance The distance to travel, in meters + * @param latitude2_out Output parameter for the end point's GPS latitude + * coordinate, in degrees; or `NULL` + * @param longitude2_out Output parameter for the end point's GPS longitude + * coordinate, in degrees; or `NULL` + * @param azimuth2_out Output parameter for the direction from the end point + * to the starting point, in degrees; or `NULL` + */ +void libtellurian_end_point(double latitude1, double longitude1, double azimuth1, double distance, + double *latitude2_out, double *longitude2_out, double *azimuth2_out); + +/** + * Calculate the location that is some distance away, + * in some direction, from some point + * + * @param latitude1 GPS latitude coordinate for the starting point, in radians + * @param longitude1 GPS longitude coordinate for the starting point, in radians + * @param azimuth1 The direction to travel, in radians + * @param distance The distance to travel, in meters + * @param latitude2_out Output parameter for the end point's GPS latitude + * coordinate, in radians; or `NULL` + * @param longitude2_out Output parameter for the end point's GPS longitude + * coordinate, in radians; or `NULL` + * @param azimuth2_out Output parameter for the direction from the end point + * to the starting point, in radians; or `NULL` + */ +void libtellurian_end_point_radians(double latitude1, double longitude1, double azimuth1, double distance, + double *latitude2_out, double *longitude2_out, double *azimuth2_out); + +/** + * Calculate the normal gravity for some point on + * the Earth's surface, where the Earth's is model + * as an oblate spheroid + * + * @param latitude GPS latitude coordinate, in degrees + * @return The normal gravity, in meters per square-second + */ +LIBTELLURIAN_CONST__ +double libtellurian_normal_gravity(double latitude); + +/** + * Calculate the normal gravity for some point on + * the Earth's surface, where the Earth's is model + * as an oblate spheroid + * + * @param latitude GPS latitude coordinate, in radians + * @return The normal gravity, in meters per square-second + */ +LIBTELLURIAN_CONST__ +double libtellurian_normal_gravity_radians(double latitude); + +/** + * Calculate the gravity adjusted for the elevation + * above the altitude where the gravity is measure + * + * Altitudes above circa 100000 meters is out of range + * for this function (that would be in outer space) + * + * @param gravity The gravity at sea level, in meters per square-second + * @param latitude GPS latitude coordinate, in degrees + * @param altitude Elevation above the gravity's measurement point, in meters + * @return The height-adjusted gravity, in meters per square-second + */ +LIBTELLURIAN_CONST__ +double libtellurian_elevated_gravity(double gravity, double latitude, double altitude); + +/** + * Calculate the gravity adjusted for the elevation + * above the altitude where the gravity is measure + * + * Altitudes above circa 100000 meters is out of range + * for this function (that would be in outer space) + * + * @param gravity The gravity at sea level, in meters per square-second + * @param latitude GPS latitude coordinate, in radians + * @param altitude Elevation above the gravity's measurement point, in meters + * @return The height-adjusted gravity, in meters per square-second + */ +LIBTELLURIAN_CONST__ +double libtellurian_elevated_gravity_radians(double gravity, double latitude, double altitude); + +/** + * Calculate the Earth's meridan radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in degrees + * @return The meridan radius of curvature, in meters (for radians!) + */ +LIBTELLURIAN_CONST__ +double libtellurian_meridan_radius(double latitude); + +/** + * Calculate the Earth's meridan radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in radians + * @return The meridan radius of curvature, in meters (for radians) + */ +LIBTELLURIAN_CONST__ +double libtellurian_meridan_radius_radians(double latitude); + +/** + * Calculate the Earth's transverse radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in degrees + * @return The transverse radius of curvature, in meters (for radians!) + */ +LIBTELLURIAN_CONST__ +double libtellurian_transverse_radius(double latitude); + +/** + * Calculate the Earth's transverse radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in radians + * @return The transverse radius of curvature, in meters (for radians) + */ +LIBTELLURIAN_CONST__ +double libtellurian_transverse_radius_radians(double latitude); + +/** + * Calculate the Earth's azimuthal radius of curvature + * at some latitude for some azimuth + * + * @param latitude GPS latitude coordinate, in degrees + * @param azimuth The azimuth, in degrees + * @return The transverse radius of curvature, in meters (for radians!) + */ +LIBTELLURIAN_CONST__ +double libtellurian_azimuthal_radius(double latitude, double azimuth); + +/** + * Calculate the Earth's azimuthal radius of curvature + * at some latitude for some azimuth + * + * @param latitude GPS latitude coordinate, in radians + * @param azimuth The azimuth, in radians + * @return The transverse radius of curvature, in meters (for radians) + */ +LIBTELLURIAN_CONST__ +double libtellurian_azimuthal_radius_radians(double latitude, double azimuth); + +/** + * Calculate the Earth's Gaussian radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in degrees + * @return The Gaussian radius of curvature, in meters (for radians!) + */ +LIBTELLURIAN_CONST__ +double libtellurian_gaussian_radius(double latitude); + +/** + * Calculate the Earth's Gaussian radius of curvature at some latitude + * + * @param latitude GPS latitude coordinate, in radians + * @return The Gaussian radius of curvature, in meters (for radians) + */ +LIBTELLURIAN_CONST__ +double libtellurian_gaussian_radius_radians(double latitude); + + +#undef LIBTELLURIAN_CONST__ +#undef LIBTELLURIAN_WUR__ + +#endif diff --git a/libtellurian_azimuth.c b/libtellurian_azimuth.c new file mode 100644 index 0000000..cc08037 --- /dev/null +++ b/libtellurian_azimuth.c @@ -0,0 +1,19 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +void +libtellurian_azimuth(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out) +{ + latitude1 = radians(latitude1); + longitude1 = radians(longitude1); + latitude2 = radians(latitude2); + longitude2 = radians(longitude2); + libtellurian_azimuth_radians(latitude1, longitude1, latitude2, longitude2, azimuth1_out, azimuth2_out); + if (azimuth1_out) + *azimuth1_out = degrees(*azimuth1_out); + if (azimuth2_out) + *azimuth2_out = degrees(*azimuth2_out); +} diff --git a/libtellurian_azimuth_radians.c b/libtellurian_azimuth_radians.c new file mode 100644 index 0000000..fa729d6 --- /dev/null +++ b/libtellurian_azimuth_radians.c @@ -0,0 +1,13 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +void +libtellurian_azimuth_radians(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out) +{ + libtellurian_vincenty_inverse__(latitude1, longitude1, + latitude2, longitude2, + NULL, azimuth1_out, azimuth2_out); +} diff --git a/libtellurian_azimuthal_radius.c b/libtellurian_azimuthal_radius.c new file mode 100644 index 0000000..3ad5841 --- /dev/null +++ b/libtellurian_azimuthal_radius.c @@ -0,0 +1,11 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_azimuthal_radius(double latitude, double azimuth) +{ + latitude = radians(latitude); + azimuth = radians(azimuth); + return libtellurian_azimuthal_radius_radians(latitude, azimuth); +} diff --git a/libtellurian_azimuthal_radius_radians.c b/libtellurian_azimuthal_radius_radians.c new file mode 100644 index 0000000..00671af --- /dev/null +++ b/libtellurian_azimuthal_radius_radians.c @@ -0,0 +1,15 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_azimuthal_radius_radians(double latitude, double azimuth) +{ + double m = libtellurian_meridan_radius_radians(latitude); + double n = libtellurian_transverse_radius_radians(latitude); + double c2 = cos(azimuth) * cos(azimuth); + double s2 = sin(azimuth) * sin(azimuth); + double x = c2 / m; + double y = s2 / n; + return 1.0 / (x + y); +} diff --git a/libtellurian_coarse_distance.c b/libtellurian_coarse_distance.c new file mode 100644 index 0000000..8f056d1 --- /dev/null +++ b/libtellurian_coarse_distance.c @@ -0,0 +1,14 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_coarse_distance(double latitude1, double longitude1, + double latitude2, double longitude2) +{ + latitude1 = radians(latitude1); + longitude1 = radians(longitude1); + latitude2 = radians(latitude2); + longitude2 = radians(longitude2); + return libtellurian_coarse_distance_radians(latitude1, longitude1, latitude2, longitude2); +} diff --git a/libtellurian_coarse_distance_radians.c b/libtellurian_coarse_distance_radians.c new file mode 100644 index 0000000..eb9cac0 --- /dev/null +++ b/libtellurian_coarse_distance_radians.c @@ -0,0 +1,13 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_coarse_distance_radians(double latitude1, double longitude1, + double latitude2, double longitude2) +{ + double h = fma(haversin(longitude2 - longitude1), + cos(latitude1) * cos(latitude2), + haversin(latitude2 - latitude1)); + return 2.0 * LIBTELLURIAN_AUTHALIC_RADIUS * asin(sqrt(h)); +} diff --git a/libtellurian_distance.c b/libtellurian_distance.c new file mode 100644 index 0000000..0a290d1 --- /dev/null +++ b/libtellurian_distance.c @@ -0,0 +1,21 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_distance(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out) +{ + double r; + latitude1 = radians(latitude1); + longitude1 = radians(longitude1); + latitude2 = radians(latitude2); + longitude2 = radians(longitude2); + r = libtellurian_distance_radians(latitude1, longitude1, latitude2, longitude2, azimuth1_out, azimuth2_out); + if (azimuth1_out) + *azimuth1_out = degrees(*azimuth1_out); + if (azimuth2_out) + *azimuth2_out = degrees(*azimuth2_out); + return r; +} diff --git a/libtellurian_distance_radians.c b/libtellurian_distance_radians.c new file mode 100644 index 0000000..b2b2134 --- /dev/null +++ b/libtellurian_distance_radians.c @@ -0,0 +1,15 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_distance_radians(double latitude1, double longitude1, + double latitude2, double longitude2, + double *azimuth1_out, double *azimuth2_out) +{ + double s; + libtellurian_vincenty_inverse__(latitude1, longitude1, + latitude2, longitude2, + &s, azimuth1_out, azimuth2_out); + return s; +} diff --git a/libtellurian_elevated_gravity.c b/libtellurian_elevated_gravity.c new file mode 100644 index 0000000..575c0e8 --- /dev/null +++ b/libtellurian_elevated_gravity.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_elevated_gravity(double gravity, double latitude, double altitude) +{ + latitude = radians(latitude); + return libtellurian_elevated_gravity_radians(gravity, latitude, altitude); +} diff --git a/libtellurian_elevated_gravity_radians.c b/libtellurian_elevated_gravity_radians.c new file mode 100644 index 0000000..ffd1a5b --- /dev/null +++ b/libtellurian_elevated_gravity_radians.c @@ -0,0 +1,14 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_elevated_gravity_radians(double gravity, double latitude, double altitude) +{ + double k1 = -3.1570429877205807e-07; + double k2 = 2.1026896504579084e-09; + double k3 = -7.374516772941995e-14; + double sin2_phi = sin(latitude) * sin(latitude); + double f = fma(k3, altitude, fma(k2, sin2_phi, k1)); + return fma(f * altitude, gravity, gravity); +} diff --git a/libtellurian_end_point.c b/libtellurian_end_point.c new file mode 100644 index 0000000..fb2c8be --- /dev/null +++ b/libtellurian_end_point.c @@ -0,0 +1,20 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +void +libtellurian_end_point(double latitude1, double longitude1, double azimuth1, double distance, + double *latitude2_out, double *longitude2_out, double *azimuth2_out) +{ + latitude1 = radians(latitude1); + longitude1 = radians(longitude1); + azimuth1 = radians(azimuth1); + libtellurian_end_point_radians(latitude1, longitude1, azimuth1, distance, + latitude2_out, longitude2_out, azimuth2_out); + if (latitude2_out) + *latitude2_out = degrees(*latitude2_out); + if (longitude2_out) + *longitude2_out = degrees(*longitude2_out); + if (azimuth2_out) + *azimuth2_out = degrees(*azimuth2_out); +} diff --git a/libtellurian_end_point_radians.c b/libtellurian_end_point_radians.c new file mode 100644 index 0000000..79e759f --- /dev/null +++ b/libtellurian_end_point_radians.c @@ -0,0 +1,79 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +/* TODO complete the implementation */ + +void +libtellurian_end_point_radians(double latitude1, double longitude1, double azimuth1, double distance, + double *latitude2_out, double *longitude2_out, double *azimuth2_out) +{ + /* + * Vincenty's formula for the "direct problem" + */ + double t, sigma, cos_2sigma_m, sin_sigma, cos_sigma; + double x, x0, y, C, L, v, lambda; + + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double c = b / a; + double f = 1.0 - c; + + double u1 = atan(c * tan(latitude1)); + + double sin_u1 = sin(u1); + double cos_u1 = cos(u1); + double sin_alpha1 = sin(azimuth1); + double cos_alpha1 = cos(azimuth1); + + double sigma1 = atan2(tan(u1), cos_alpha1); + + double sin_alpha = cos_u1 * sin_alpha1; + double cos2_alpha = fma(-sin_alpha, sin_alpha, 1.0); + + double uu = cos2_alpha * fma(a / b, a / b, -1.0); + + double A = fma(fma(fma(fma(-175.0, uu, 320.0), uu, -768.0), uu, 4096.0), uu / 16384.0, 1.0); + double B = fma(fma(fma(-47.0, uu, 74.0), uu, -128.0), uu, 256.0) * (uu / 1024.0); + + double sigma_0 = distance / (b * A); + + sigma = sigma_0; + + /* { */ + + cos_2sigma_m = cos(fma(2.0, sigma1, sigma)); + sin_sigma = sin(sigma); + cos_sigma = cos(sigma); + + v = cos_sigma * fma(2 * cos_2sigma_m, cos_2sigma_m, -1.0); + t = fma(2.0 * cos_2sigma_m, 2.0 * cos_2sigma_m, -3.0); + t *= fma(2.0 * sin_sigma, 2.0 * sin_sigma, -3.0); + t = fma(B * cos_2sigma_m / -6.0, t, v); + t = fma(0.25 * B, t, cos_2sigma_m); + sigma = fma(B * sin_sigma, t, sigma_0); + + /* } repeat until sigma converges */ + + if (latitude2_out || azimuth2_out) + x0 = fma(cos_sigma * cos_alpha1, cos_u1, -sin_u1 * sin_alpha); + + if (latitude2_out) { + y = fma(sin_sigma * cos_alpha1, cos_u1, sin_u1 * cos_sigma); + x = sqrt(fma(sin_alpha, sin_alpha, x0 * x0)) * c; + *latitude2_out = atan2(y, x); + } + + if (azimuth2_out) + *azimuth2_out = atan2(sin_alpha, x0); + + if (longitude2_out) { + C = fma(fma(-0.75, cos2_alpha, 1.0), f, 1.0) * cos2_alpha * f / 4.0; + y = sin_sigma * sin_alpha1; + x = fma(sin_sigma * cos_alpha1, sin_u1, cos_u1 * cos_sigma); + lambda = atan2(y, x); + t = fma(C * sin_sigma, fma(C, v, cos_2sigma_m), sigma); + L = fma(fma(C, f, -f) * sin_alpha, t, lambda); + *longitude2_out = longitude1 + L; + } +} diff --git a/libtellurian_gaussian_radius.c b/libtellurian_gaussian_radius.c new file mode 100644 index 0000000..24c9b54 --- /dev/null +++ b/libtellurian_gaussian_radius.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_gaussian_radius(double latitude) +{ + latitude = radians(latitude); + return libtellurian_gaussian_radius_radians(latitude); +} diff --git a/libtellurian_gaussian_radius_radians.c b/libtellurian_gaussian_radius_radians.c new file mode 100644 index 0000000..d61319a --- /dev/null +++ b/libtellurian_gaussian_radius_radians.c @@ -0,0 +1,16 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_gaussian_radius_radians(double latitude) +{ + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double f = 1.0 - b / a; + double neg_e2 = (f - 2.0) * f; + double neg_e2_plus_1 = fma(f - 2.0, f, 1.0); + double s = sin(latitude); + double denom = fma(neg_e2, s * s, 1.0); + return a * sqrt(neg_e2_plus_1) / denom; +} diff --git a/libtellurian_meridan_radius.c b/libtellurian_meridan_radius.c new file mode 100644 index 0000000..01c597f --- /dev/null +++ b/libtellurian_meridan_radius.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_meridan_radius(double latitude) +{ + latitude = radians(latitude); + return libtellurian_meridan_radius_radians(latitude); +} diff --git a/libtellurian_meridan_radius_radians.c b/libtellurian_meridan_radius_radians.c new file mode 100644 index 0000000..0a832e5 --- /dev/null +++ b/libtellurian_meridan_radius_radians.c @@ -0,0 +1,14 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_meridan_radius_radians(double latitude) +{ + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double f = 1.0 - b / a; + double neg_e2 = (f - 2.0) * f; + double s = sin(latitude); + return fma(a, neg_e2, a) / pow(fma(neg_e2, s * s, 1.0), 1.5); +} diff --git a/libtellurian_normal_gravity.c b/libtellurian_normal_gravity.c new file mode 100644 index 0000000..27722b1 --- /dev/null +++ b/libtellurian_normal_gravity.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_normal_gravity(double latitude) +{ + latitude = radians(latitude); + return libtellurian_normal_gravity_radians(latitude); +} diff --git a/libtellurian_normal_gravity_radians.c b/libtellurian_normal_gravity_radians.c new file mode 100644 index 0000000..3fdb037 --- /dev/null +++ b/libtellurian_normal_gravity_radians.c @@ -0,0 +1,18 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_normal_gravity_radians(double latitude) +{ + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double neg_e2 = fma(b / a, b / a, -1.0); + double ag = a * LIBTELLURIAN_EQUATORIAL_GRAVITY; + double bg = b * LIBTELLURIAN_POLAR_GRAVITY; + double k = bg / ag - 1.0; + double sin2_phi = sin(latitude) * sin(latitude); + double num = fma(k, sin2_phi, 1.0); + double denom2 = fma(neg_e2, sin2_phi, 1.0); + return LIBTELLURIAN_EQUATORIAL_GRAVITY * num / sqrt(denom2); +} diff --git a/libtellurian_sea_level.c b/libtellurian_sea_level.c new file mode 100644 index 0000000..76f817f --- /dev/null +++ b/libtellurian_sea_level.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_sea_level(double latitude) +{ + latitude = radians(latitude); + return libtellurian_sea_level_radians(latitude); +} diff --git a/libtellurian_sea_level_radians.c b/libtellurian_sea_level_radians.c new file mode 100644 index 0000000..a23b222 --- /dev/null +++ b/libtellurian_sea_level_radians.c @@ -0,0 +1,17 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_sea_level_radians(double latitude) +{ + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double c = cos(latitude); + double s = sin(latitude); + double x = a * c * a; + double y = b * s * b; + double num = fma(x, x, y * y); + double denom = fma(x, c, y * s); + return sqrt(num / denom); +} diff --git a/libtellurian_transverse_radius.c b/libtellurian_transverse_radius.c new file mode 100644 index 0000000..d94ca86 --- /dev/null +++ b/libtellurian_transverse_radius.c @@ -0,0 +1,10 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_transverse_radius(double latitude) +{ + latitude = radians(latitude); + return libtellurian_transverse_radius_radians(latitude); +} diff --git a/libtellurian_transverse_radius_radians.c b/libtellurian_transverse_radius_radians.c new file mode 100644 index 0000000..dae8f04 --- /dev/null +++ b/libtellurian_transverse_radius_radians.c @@ -0,0 +1,14 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +double +libtellurian_transverse_radius_radians(double latitude) +{ + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double f = 1.0 - b / a; + double neg_e2 = (f - 2.0) * f; + double s = sin(latitude); + return pow(fma(neg_e2, s * s, 1.0), -0.5); +} diff --git a/libtellurian_vincenty_inverse__.c b/libtellurian_vincenty_inverse__.c new file mode 100644 index 0000000..746011f --- /dev/null +++ b/libtellurian_vincenty_inverse__.c @@ -0,0 +1,104 @@ +/* See LICENSE file for copyright and license details. */ +#include "common.h" + + +/* TODO complete the implementation */ + +void +libtellurian_vincenty_inverse__(double latitude1, double longitude1, double latitude2, double longitude2, + double *distance_out, double *azimuth1_out, double *azimuth2_out) +{ + /* + * Vincenty's formula for the "inverse problem" + */ + + double lambda, t, uu, A, B, sigma_minus_delta_sigma; + double x, y, cos_lambda, sin_lambda, sin_sigma, cos_sigma, sigma; + double sin_alpha, cos2_alpha, cos_2sigma_m, C, cos2_2sigma_m; + + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double c = b / a; + double f = 1.0 - c; + + double u1 = atan(c * tan(latitude1)); + double u2 = atan(c * tan(latitude2)); + + double cos_u1 = cos(u1), sin_u1 = sin(u1); + double cos_u2 = cos(u2), sin_u2 = sin(u2); + + double cos_u1_sin_u2 = cos_u1 * sin_u2; + double sin_u1_cos_u2 = sin_u1 * cos_u2; + double cos_u1_cos_u2 = cos_u1 * cos_u2; + double sin_u1_sin_u2 = sin_u1 * sin_u2; + + double L = longitude2 - longitude1; + + lambda = L; + + /* { */ + + cos_lambda = cos(lambda); + sin_lambda = sin(lambda); + + y = cos_u2 * sin_lambda; + x = fma(-sin_u1_cos_u2, cos_lambda, cos_u1_sin_u2); + + sin_sigma = sqrt(fma(y, y, x * x)); + cos_sigma = fma(cos_u1_cos_u2, cos_lambda, sin_u1_sin_u2); + sigma = atan2(sin_sigma, cos_sigma); + + sin_alpha = (cos_u1_cos_u2 * sin_lambda) / sin_sigma; + cos2_alpha = fma(-sin_alpha, sin_alpha, 1.0); + /* + * If sin σ = 0 the value of sin α is indeterminate. + * It represents an end point coincident with, or + * diametrically opposed to, the start point. + */ + + cos_2sigma_m = fma(-2.0 / cos2_alpha, sin_u1_sin_u2, cos_sigma); + C = 0.25 * f * cos2_alpha * fma(f, fma(-0.75, cos2_alpha, 1.0), 1.0); + + cos2_2sigma_m = cos_2sigma_m * cos_2sigma_m; + t = fma(2.0, cos2_2sigma_m, -1.0); + t = fma(C * cos_sigma, t, cos_2sigma_m); + t = fma(C * sin_sigma, t, sigma); + lambda = fma(fma(C, f, -f) * sin_alpha, t, L); + + /* } repeat until lambda converges */ + + if (distance_out) { + uu = cos2_alpha * fma(a / b, a / b, -1.0); + A = fma(fma(fma(fma(-175.0, uu, 320.0), uu, -768.0), uu, 4096.0), uu / 16384.0, 1.0); + B = fma(fma(fma(-47.0, uu, 74.0), uu, -128.0), uu, 256.0) * (uu / 1024.0); + + t = fma(4.0, cos2_2sigma_m, -3.0); + t *= fma(2.0 * sin_sigma, 2.0 * sin_sigma, -3.0); + t *= cos_2sigma_m * B / -6.0; + t = fma(cos_sigma, fma(2.0, cos2_2sigma_m, -1.0), t); + t = fma(t, B / 4.0, cos_2sigma_m); + sigma_minus_delta_sigma = fma(t, -B * sin_sigma, sigma); + + *distance_out = b * A * sigma_minus_delta_sigma; + } + + if (azimuth1_out) + *azimuth1_out = atan2(y, x); + + if (azimuth2_out) { + y = cos_u1 * sin_lambda; + x = fma(cos_u1_sin_u2, cos_lambda, -sin_u1_cos_u2); + *azimuth1_out = atan2(y, x); + } + + /* + * Between two nearly antipodal points, the iterative formula may fail to converge; + * this will occur when the first guess at λ as computed by the equation above is + * greater than π in absolute value. + */ + + /* + * https://en.wikipedia.org/wiki/Vincenty's_formulae#Vincenty's_modification + * https://en.wikipedia.org/wiki/Vincenty's_formulae#Nearly_antipodal_points + */ +} diff --git a/mk/linux.mk b/mk/linux.mk new file mode 100644 index 0000000..ad58f69 --- /dev/null +++ b/mk/linux.mk @@ -0,0 +1,6 @@ +LIBEXT = so +LIBFLAGS = -shared -Wl,-soname,lib$(LIB_NAME).$(LIBEXT).$(LIB_MAJOR) +LIBMAJOREXT = $(LIBEXT).$(LIB_MAJOR) +LIBMINOREXT = $(LIBEXT).$(LIB_VERSION) + +FIX_INSTALL_NAME = : diff --git a/mk/macos.mk b/mk/macos.mk new file mode 100644 index 0000000..72be640 --- /dev/null +++ b/mk/macos.mk @@ -0,0 +1,6 @@ +LIBEXT = dylib +LIBFLAGS = -dynamiclib -Wl,-compatibility_version,$(LIB_MAJOR) -Wl,-current_version,$(LIB_VERSION) +LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT) +LIBMINOREXT = $(LIB_VERSION).$(LIBEXT) + +FIX_INSTALL_NAME = install_name_tool -id "$(PREFIX)/lib/libtellurian.$(LIBMAJOREXT)" diff --git a/mk/windows.mk b/mk/windows.mk new file mode 100644 index 0000000..ed5ec8d --- /dev/null +++ b/mk/windows.mk @@ -0,0 +1,6 @@ +LIBEXT = dll +LIBFLAGS = -shared +LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT) +LIBMINOREXT = $(LIB_VERSION).$(LIBEXT) + +FIX_INSTALL_NAME = : |