#!/bin/sh
#-
# © 2021 mirabilos Ⓕ CC0; implementation of Haversine GCD from public sources
#
# now developed online:
# https://evolvis.org/plugins/scmgit/cgi-bin/gitweb.cgi?p=useful-scripts/mirkarte.git;a=blob;f=geo.sh;hb=HEAD

if test "$#" -ne 4; then
	echo >&2 "E: syntax: $0 lat1 lon1 lat2 lon2"
	exit 1
fi

set -e

# make GNU bc use POSIX mode and shut up
BC_ENV_ARGS=-qs
export BC_ENV_ARGS

# assignment of constants, variables and functions
# p: multiply with to convert from degrees to radians (π/180)
# r: earth radius in metres
# d: distance
# h: haversine intermediate
# i,j: (lat,lon) point 1
# x,y: (lat,lon) point 2
# k: delta lat
# l: delta lon
# m: sin(k/2) (square root of hav(k))
# n: sin(l/2) (  partial haversine  )
# n(x): arcsin(x)
# r(x,n): round x to n decimal digits
# v(x): sign (Vorzeichen)
# w(x): min(1, sqrt(x)) (Wurzel)

bc -l <<-EOF
scale=64
define n(x) {
	if (x == -1) return (-2 * a(1))
	if (x == 1) return (2 * a(1))
	return (a(x / sqrt(1 - x*x)))
}
define v(x) {
	if (x < 0) return (-1)
	if (x > 0) return (1)
	return (0)
}
define r(x, n) {
	auto o
	o = scale
	if (scale < (n + 1)) scale = (n + 1)
	x += v(x) * 0.5 * A^-n
	scale = n
	x /= 1
	scale = o
	return (x)
}
define w(x) {
	if (x >= 1) return (1)
	return (sqrt(x))
}
/* WGS84 reference ellipsoid: große Halbachse (metres), Abplattung */
i = 6378137.000
x = 1/298.257223563
/* other axis */
j = i * (1 - x)
/* mean radius resulting */
r = (2 * i + j) / 3
/* coordinates */
p = (4 * a(1) / 180)
i = (p * $1)
j = (p * $2)
x = (p * $3)
y = (p * $4)
/* calculation */
k = (x - i)
l = (y - j)
m = s(k / 2)
n = s(l / 2)
h = ((m * m) + (c(i) * c(x) * n * n))
d = 2 * r * n(w(h))
r(d, 3)
EOF

# output is in metres, rounded to millimetres, error < ¼% in WGS84