Haversine formula: Difference between revisions

=={{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.