Examining the Haversine formula

[Mike Sanders]

Examining the Haversine formula - Michael Sanders 2016


SYNOPSIS


Looking for a way to calculate distance, I discovered the Haversine formula
which as it turns out, is indeed a good way to measure 'straight line'
distances. And because awk is an interpreted rather than compiled language,
I was able to test this formula quickly before deploying its use in other
languages, reducing write/compile/debug times...

So far, I've measured the accuracy of the equation against a single
source[i], but the results are close enough that I consider the awk
function as 'ready to use'. Most of my reading on the subject is based on
an article at NYU titled 'Computing Distances'[ii].


BACKGROUND


The Haversine formula is an equation important in navigation, giving great-
circle distances between two points on a sphere from their longitudes and
latitudes. It is a special case of a more general formula in spherical
trigonometry, the law of haversines, relating the sides and angles of
spherical triangles. The first table of haversines in English was published
by James Andrew in 1805[iii].


HAVERSINE FORMULA


Here's the formula as published by R.W. Sinnott in 1984[iv]...

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat / 2)) ^ 2 + cos(lat1) * cos(lat2) * (sin(dlon / 2)) ^ 2
c = 2 * arcsin(min(1, sqrt(a)))
d = R * c

where R is the radius of the Earth


CAVEAT


What the equation can not account for is variation in terrain such as a
mountain or valley...

For instance: Driving a car from point A to point B along the planet's
surface, the distance traveled would be greater than flying a plane between
the same two points. This straight line method of reckoning, 'as the crow
flies' or 'a beeline', means the distance between two points will tend to
be somewhat shorter. In any event, the equation calculates approximate
values only.

Another item affecting accuracy is the input radius. Below, I've used
earth's radius at the equator for no other reason than being easy to
visualize.


AWK SCRIPT


# Haversine formula as implemented in awk
# v1.00 - Michael Sanders 2016 <www.peanut-software.com>
# based on the equation by R.W. Sinnott
# invocation; awk -f haversine.awk

BEGIN{

# city | latitude | longitude
# nyc | 40.71278 | -74.00594
# la | 34.05223 | -118.24368

rkm = 6371; # avg. radius of earth at equator in kilometers
rmi = 3959; # avg. radius of earth at equator in miles

dkm = haversine(rkm, 34.05223, -118.24368, 40.71278, -74.00594);
dmi = haversine(rmi, 34.05223, -118.24368, 40.71278, -74.00594);

printf("%s%s\n", "distance in kilometers from la to nyc: ", dkm);
printf("%s%s\n", "distance in miles from la to nyc: ", dmi);

}

function haversine(radius, lat1, lon1, lat2, lon2) {

deglon = lon2 - lon1;
deglat = lat2 - lat1;

a = (sin(deglat/2))^2 + cos(lat1) * cos(lat2) * (sin(deglat/2))^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));

return (radius * c);

}

# eof


NOTES


i. <http://www.distancefromto.net/distance-from-new-york-to-los-angeles-us>

ii. <http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html>

iii. <https://en.wikipedia.org/wiki/Haversine_formula>

iv. Virtues of the Haversine: Sky and Telescope, vol. 68, no. 2,
1984 - R.W. Sinnott


LICENSE


[c]2016 Michael Sanders. Reprint freely so long as nothing is modified.


--
Mike Sanders
www.peanut-software.com

[Mike Sanders]

Hope posting these papers is not somehow out of the
norm for this group, just sharing the brain food...

--
Mike Sanders
www.peanut-software.com

[Kenny McCormack]

In article <slrnnip6b...@porkchop.nix>,

Mike Sanders <mi...@porkchop.nix> wrote:
>Hope posting these papers is not somehow out of the
>norm for this group, just sharing the brain food...

I think it is fine. In particular, this specific problem was one I've long
wanted a solution to - and having it in AWK is perfect for my use.

Just don't try this sort of thing in comp.lang.c, though... (heh heh)

--
Modern Christian: Someone who can take time out from
complaining about "welfare mothers popping out babies we
have to feed" to complain about welfare mothers getting
abortions that PREVENT more babies to be raised at public
expense.

[Marc de Bourget]

Hello Mike,

thank you a lot for sharing your ideas/scripts.
Please continue to do so, I enjoy reading them!

Marc de Bourget

[Mike Sanders]

Kenny McCormack <gaz...@shell.xmission.com> wrote:

> I think it is fine. In particular, this specific problem was one
> I've long wanted a solution to - and having it in AWK is perfect
> for my use.

Great, hope you can put that pup to work Kenny. Heavens knows I've
nabbed alot from this group myself.

> Just don't try this sort of thing in comp.lang.c, though... (heh heh)

[shudder] Not me, foul dragons lurk therein...

--
Mike Sanders
www.peanut-software.com

[Mike Sanders]

Marc de Bourget <marcde...@gmail.com> wrote:

> Hello Mike,
>
> thank you a lot for sharing your ideas/scripts.
> Please continue to do so, I enjoy reading them!

Hey Marc.

Thank you, appreciate your kind words. =)

--
Mike Sanders
www.peanut-software.com

[Anton Treuenfels]

"Mike Sanders" <mi...@porkchop.nix> wrote in message
news:slrnnip68...@porkchop.nix...

>
>
> Here's the formula as published by R.W. Sinnott in 1984[iv]...
>
> dlon = lon2 - lon1
> dlat = lat2 - lat1

Nitpicker's observation:

> a = (sin(dlat / 2)) ^ 2 + cos(lat1) * cos(lat2) * (sin(dlon / 2)) ^ 2


> c = 2 * arcsin(min(1, sqrt(a)))
> d = R * c
>
> where R is the radius of the Earth

[...]

> function haversine(radius, lat1, lon1, lat2, lon2) {
>
> deglon = lon2 - lon1;
> deglat = lat2 - lat1;
>

a difference here:

[Mike Sanders]

Anton Treuenfels <teamt...@yahoo.com> wrote:

> Nitpicker's observation...

Keen observations Anton.

Truth is, you're reading my attempt to understand the formula.
But since I've cited my sources, everyone *should* judge its
fitness accordingly because I'm no expert...

--
Mike Sanders
www.peanut-software.com

[Mike Sanders]

Anton Treuenfels (thank you Anton) noted this error/typo:

> a = (sin(deglat/2))^2 + cos(lat1) * cos(lat2) * (sin(deglat/2))^2;

It should instead read:

a = (sin(deglat/2))^2 + cos(lat1) * cos(lat2) * (sin(deglon/2))^2;

The updated paper is here:

http://www.peanut-software.com/haversine-formula.pdf

--
Mike Sanders
www.peanut-software.com

[Paulino Huerta]

Hi Mike,

It seems to forget conversion to radians.

To convert decimal degrees to radians, multiply the number
of degrees by pi/180 = 0.017453293 radians/degree.

Example:

rad = 0.017453293
...
...
deglon =( lon2 - lon1 ) * rad
...
cos(lat1 * rad)
...

Regards,

Paulino

[Mike Sanders]

Paulino Huerta <paulin...@gmail.com> wrote:

> Hi Mike,
>

> It seems to forget conversion to radians...

Hello Paulino.

ahh yes - I understand the *entire* idea now!!!

*decimals to radians*

***

result = deg * (pi / 180)

***

And thank you all. I will update the paper yet again
(I'm determined to work this out properly).

--
Mike Sanders
www.peanut-software.com

[Mike Sanders]

Okay! Heres my definitive version:

http://www.peanut-software.com/haversine-formula.pdf

Changes...

. pulled my head out (pi/180)

. added thanks section to pdf

. updated all foot notes

. just found another source (of course!) to check my work against:

https://rosettacode.org/wiki/Haversine_formula#AWK


And here's the updated script...

# Haversine formula as implemented in awk

# v1.02 - Michael Sanders 2016 <www.peanut-software.com>

# based on the equation by R.W. Sinnott
# invocation: awk -f haversine.awk

BEGIN {

# city | latitude | longitude

# la | 34.05223 | -118.24368
# ny | 40.71278 | -74.00594

rkm = 6371; # earth's radius at equator in kilometers
rmi = 3959; # earth's radius at equator in miles

dkm = haversine(rkm, 34.05223, -118.24368, 40.71278, -74.00594);
dmi = haversine(rmi, 34.05223, -118.24368, 40.71278, -74.00594);

printf("%s%s\n", "distance in kilometers from la to ny: ", dkm);
printf("%s%s\n", "distance in miles from la to ny: ", dmi);

}

function haversine(radius, lat1, lon1, lat2, lon2) {

rad = 0.017453293; # pi/180
deglon = (lon2 - lon1) * rad;
deglat = (lat2 - lat1) * rad;
lat1 = (lat1 * rad);
lat2 = (lat2 * rad);

a = (sin(deglat/2))^2 + cos(lat1) * cos(lat2) * (sin(deglon/2))^2;

c = 2 * atan2(sqrt(a), sqrt(1-a));

return (radius * c);

}

# eof

And my brain is now mush... =)

--
Mike Sanders
www.peanut-software.com

[Luuk]

On 07-05-16 16:09, Mike Sanders wrote:
> . just found another source (of course!) to check my work against:
>
> https://rosettacode.org/wiki/Haversine_formula#AWK

for a visual look:
http://www.distance.to/New-York/Los-Angeles

[Janis Papanagnou]

On 07.05.2016 16:09, Mike Sanders wrote:
> Okay! Heres my definitive version:

Famous words!

That's why version numbers have been introduced and used in
software development (e.g. major/minor/patch).

Janis

> [...]

[Mike Sanders]

Luuk <lu...@invalid.lan> wrote:

> for a visual look:
> http://www.distance.to/New-York/Los-Angeles

Thanks Luuk (a good site to bookmark in any event).

--
Mike Sanders
www.peanut-software.com

[Mike Sanders]

Janis Papanagnou <janis_pa...@hotmail.com> wrote:

>> Okay! Heres my definitive version:
>
> Famous words!
>
> That's why version numbers have been introduced and used in
> software development (e.g. major/minor/patch).

[In]Famous in my case Janis. Embarrassed to say I became too excited :/
Fortunately, the script in the article does have a version number.

--
Mike Sanders
www.peanut-software.com

[Janis Papanagnou]

On 07.05.2016 17:24, Mike Sanders wrote:
> Janis Papanagnou <janis_pa...@hotmail.com> wrote:
>
>>> Okay! Heres my definitive version:
>>
>> Famous words!
>>
>> That's why version numbers have been introduced and used in
>> software development (e.g. major/minor/patch).
>
> [In]Famous in my case Janis. Embarrassed to say I became too excited :/
> Fortunately, the script in the article does have a version number.

Yeah, I saw it; but don't count on it being the "definite version". ;-)

Janis

[johne...@gmail.com]

Very good thread, Mike. Thanks Mike and all those who contributed!

[Mike Sanders]

johne...@gmail.com <johne...@gmail.com> wrote:

> Very good thread, Mike. Thanks Mike and all those who contributed!

You're very welcome (hope it'll be useful for you) & I'm thankful
for the group's help too. I enjoy learning new things.

--
Mike Sanders
www.peanut-software.com