Ok I'm stumped as to why my Spherical Law of Cosines calcs don't seem to work

150 views
Skip to first unread message

Paul LeChef

unread,
Jan 5, 2015, 2:40:13 PM1/5/15
to mitappinv...@googlegroups.com
HI All,

I'm trying to get the distance between GPS coordinates using the Spherical Law of Cosines since the points will rarely be more than a couple of kilometers from each other. 

Since I couldn't get it to work properly in my app, I decided to simplify things dramatically. I made another app that simply calculates the distance between 
  1. my house and my brother-in-law's house
  2. my house and my parents' house
According to http://www.movable-type.co.uk/scripts/latlong.html , distance 1 should be 2312 meters and distance 2 should be 18840 meters (ya, need it in meters although I would certainly settle if someone could tell me where I goofed in KM)

Results are 2831m for #1 and 26509m for #2

Any thoughts or suggestions would be welcomed
Cheers

Paul
distBlocks.pdf

Aubrey Colter

unread,
Jan 5, 2015, 3:19:33 PM1/5/15
to mitappinv...@googlegroups.com
Hi Paul,

I just took a quick look. I think you might have a simple math error. According to Wikipedia, the Spherical Law of Cosines is : 

\cos(c) = \cos(a) \cos(b) + \sin(a) \sin(b) \cos(C). \,


But in your code, you are doing:
cos(c) = sin(a)sin(b) + cos(a)cos(b)cos(c).

So try moving that last cos(c) block into your multiplication with the sin(a)sin(b) and see if that fixes it.

Good luck and happy apping!
Aubrey

Paul LeChef

unread,
Jan 5, 2015, 3:39:17 PM1/5/15
to mitappinv...@googlegroups.com
Hmmm, actually I don't think so, Aubrey.To work out distance using the law of cosines, as per http://www.movable-type.co.uk/scripts/latlong.htm

Spherical Law of Cosines

In fact, JavaScript (and most modern computers & languages) use ‘IEEE 754’ 64-bit floating-point numbers, which provide 15 significant figures of precision. By my estimate, with this precision, the simple spherical law of cosines formula (cos c = cos a cosb + sin a sin b cos C) gives well-conditioned results down to distances as small as a few metres on the earth’s surface. (Note that the geodetic form of the law of cosines is rearranged from the canonical one so that the latitude can be used directly, rather than thecolatitude).

This makes the simpler law of cosines a reasonable 1-line alternative to the haversine formula for many geodesy purposes (if not for astronomy). The choice may be driven by coding context, available trig functions (in different languages), etc – and, for very small distances an equirectangular approximation may be more suitable.

Law of cosines:d = acos( sin φ1 ⋅ sin φ2 + cos φ1 ⋅ cos φ2 ⋅ cos Δλ ) ⋅ R
JavaScript:
var φ1 = lat1.toRadians(), φ2 = lat2.toRadians(), Δλ = (lon2-lon1).toRadians(), R = 6371; // gives d in km
var d = Math.acos( Math.sin1)*Math.sin2) + Math.cos1)*Math.cos2) * Math.cos(Δλ) ) * R;
Excel:=ACOS( SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1) ) * 6371
(or with lat/lon in degrees):
=ACOS( SIN(lat1*PI()/180)*SIN(lat2*PI()/180) + COS(lat1*PI()/180)*COS(lat2*PI()/180)*COS(lon2*PI()/180-lon1*PI()/180) ) * 6371

I think I'll try out the excel formula to see if there is a difference 

Paul LeChef

unread,
Jan 5, 2015, 3:50:36 PM1/5/15
to mitappinv...@googlegroups.com
Well, I definitely have an error somewhere in the blocks. When I run the calculation in Excel as described in http://www.movable-type.co.uk/scripts/latlong.html, the numbers match up if rounded to an integer

Guess I'll keep hunting

Paul

Aubrey Colter

unread,
Jan 5, 2015, 4:05:03 PM1/5/15
to mitappinv...@googlegroups.com
Sorry about that! I was hoping it would be an easy fix. I'll poke around a little bit right now, too. 

Paul LeChef

unread,
Jan 5, 2015, 4:51:39 PM1/5/15
to mitappinv...@googlegroups.com
I went ahead and removed all passed in values then compared each piece with its equivalent excel result. I'm a little surprised by what I found. 

Every calc that refers to a lat is wrong. ( the sin's and the first two cos') I'm obviously not building these blocks properly.


distBlocks2.pdf

Aubrey Colter

unread,
Jan 5, 2015, 5:12:43 PM1/5/15
to mitappinv...@googlegroups.com
Try pulling your lat blocks out and putting the numbers in manually (using the number blocks), and see if it comes out right. Then we can see if it's the equation or the variables that are going wrong.

Paul LeChef

unread,
Jan 5, 2015, 5:25:02 PM1/5/15
to mitappinv...@googlegroups.com
Everything is fine right including the conversion to radians. Once it hits the trig function .... not so good at all.

Just had an epithany....(sp?) anyway... Is it possible that AI2 converts to radians in the background, so that my converting to radians is actually doing it twice?

Testing the epithany now


Paul LeChef

unread,
Jan 5, 2015, 5:29:16 PM1/5/15
to mitappinv...@googlegroups.com
OK, we can file this one under "undocumented features"

Removing the conversion to radians fixes the issue.

Thanks, for being a great sounding board, Aubry

Cheers

Aubrey Colter

unread,
Jan 5, 2015, 5:57:47 PM1/5/15
to mitappinv...@googlegroups.com
Thanks for patiently debugging. I'm glad you found the issue! I will figure out how to document it. What you're saying is that AI2 trig functions convert to radians in the background? That seems a little iffy. I'll take it from here!

M. Hossein Amerkashi

unread,
Jan 5, 2015, 6:07:53 PM1/5/15
to mitappinv...@googlegroups.com
@Paul, I checked the code and you are correct and under the hood, they are converted to radians.

-Hossein.

M. Hossein Amerkashi

unread,
Jan 5, 2015, 6:34:00 PM1/5/15
to mitappinv...@googlegroups.com
I might've spoke too soon and looked in the wrong place.
Use do-it on sin block and you see that its working as expected. Use scientific calculator to compare.

Paul LeChef

unread,
Jan 6, 2015, 1:11:06 PM1/6/15
to mitappinv...@googlegroups.com
@Hossein, the Do-It only works when I remove the convert-to-radians block.... Am I doing something wrong with the blocks?

At the SIN function:
Without Conversion: Do It Result: 0.71397
---
should be 0.713974972 (calculated in Excel)

With Conversion: Do It Result: 0.01388
---
should be 0.713974972

                        




Reply all
Reply to author
Forward
0 new messages