84 lines
1.7 KiB
Bash
84 lines
1.7 KiB
Bash
#!/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
|