Haversine formula: Difference between revisions
Content added Content deleted
Drkameleon (talk | contribs) (Added Arturo implementation) |
(add bc(1) implementation) |
||
Line 420: | Line 420: | ||
Distance = 2887.25995 km |
Distance = 2887.25995 km |
||
</pre> |
</pre> |
||
=={{header|bc}}== |
|||
{{works with|GNU_bc|GNU bc(1)|1.07.1-2}} |
|||
{{works with|OpenBSD_bc|dc(1)-based MirBSD bc(1)|as of 2021-06-11}} |
|||
{{works with|Bc|bc(1)|POSIX}} |
|||
… supplied with a small POSIX shell wrapper to feed the arguments to <code>bc</code>. ([https://edugit.org/-/snippets/27 see also]) |
|||
<lang sh> |
|||
#!/bin/sh |
|||
#- |
|||
# © 2021 mirabilos Ⓕ CC0; implementation of Haversine GCD from public sources |
|||
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 |
|||
# 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=42 |
|||
define n(x) { |
|||
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)) |
|||
} |
|||
p = (3.1415926535897932 / 180) |
|||
r = 6371008.8 |
|||
i = (p * $1) |
|||
j = (p * $2) |
|||
x = (p * $3) |
|||
y = (p * $4) |
|||
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 |
|||
</lang> |
|||
{{out}} |
|||
<pre>$ sh dist.sh 36.12 -86.67 33.94 -118.4 |
|||
2886448.430</pre> |
|||
Note I used a more precise earth radius; this matches the other implementations when choosing the same. |
|||
=={{header|C}}== |
=={{header|C}}== |