aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--.gitignore14
-rw-r--r--LICENSE15
-rw-r--r--Makefile96
-rw-r--r--README23
-rw-r--r--TODO3
-rw-r--r--common.h19
-rw-r--r--config.mk8
-rw-r--r--libtellurian.h440
-rw-r--r--libtellurian_azimuth.c19
-rw-r--r--libtellurian_azimuth_radians.c13
-rw-r--r--libtellurian_azimuthal_radius.c11
-rw-r--r--libtellurian_azimuthal_radius_radians.c15
-rw-r--r--libtellurian_coarse_distance.c14
-rw-r--r--libtellurian_coarse_distance_radians.c13
-rw-r--r--libtellurian_distance.c21
-rw-r--r--libtellurian_distance_radians.c15
-rw-r--r--libtellurian_elevated_gravity.c10
-rw-r--r--libtellurian_elevated_gravity_radians.c14
-rw-r--r--libtellurian_end_point.c20
-rw-r--r--libtellurian_end_point_radians.c79
-rw-r--r--libtellurian_gaussian_radius.c10
-rw-r--r--libtellurian_gaussian_radius_radians.c16
-rw-r--r--libtellurian_meridan_radius.c10
-rw-r--r--libtellurian_meridan_radius_radians.c14
-rw-r--r--libtellurian_normal_gravity.c10
-rw-r--r--libtellurian_normal_gravity_radians.c18
-rw-r--r--libtellurian_sea_level.c10
-rw-r--r--libtellurian_sea_level_radians.c17
-rw-r--r--libtellurian_transverse_radius.c10
-rw-r--r--libtellurian_transverse_radius_radians.c14
-rw-r--r--libtellurian_vincenty_inverse__.c104
-rw-r--r--mk/linux.mk6
-rw-r--r--mk/macos.mk6
-rw-r--r--mk/windows.mk6
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
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..f930ea6
--- /dev/null
+++ b/LICENSE
@@ -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
diff --git a/README b/README
new file mode 100644
index 0000000..7dc7dae
--- /dev/null
+++ b/README
@@ -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.
diff --git a/TODO b/TODO
new file mode 100644
index 0000000..09145cf
--- /dev/null
+++ b/TODO
@@ -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 = :