TheBackShed.com

Forum Index : Microcontroller and PC projects : GPS calculation / Electride / Fortran

Page 1 of 2

Author

Message

isochronic
Guru

Joined: 21/01/2012
Location: Australia
Posts: 689

Posted: 02:55am 31 Dec 2014

Copy link to clipboard

Print this post

What distance do you get, calculating it between Melbourne and Sydney ?

I get 712.8965 km but it seems a bit short, can anyone with a navsat etc confirm ?

Here is the data etc (rounded a little)

Distance = acos( sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) *

cos(long2 - long1) )

It has been constituted over several lines.

PROGRAM inpgpsdist;
C;
DOUBLE xo, yo, xd, yd, r ;
DOUBLE qf ;
WRITE "Start" ;
r = 57.29577951 ;
WRITE "What is the Origin latitude ?" ;
READ xo ;
WRITE ;
WRITE "What is the Origin longitude ?" ;
READ yo ;
WRITE ;
WRITE "Destination latitude ?" ;
READ xd ;
WRITE ;
WRITE "Destination longitude ?" ;
READ yd ;
WRITE ;
WRITE ;
WRITE "Origin latitude is ", xo ;
WRITE "Origin longitude is ", yo ;
WRITE "Destination latitude is ", xd ;
WRITE "Destination longitude is ", yd ;
xo = xo / r ;
yo = yo / r ;
xd = xd / r ;
yd = yd / r ;
qf = dist(xo,yo,xd,yd) ;
FORMAT (A,F6.1,A,/) ;
WRITE (6,0) "Distance is ", qf, " km" ;
END ;
C;
DOUBLE FUNCTION dist ( a, b, c, e ) ;
DOUBLE l, m, z, t, r ;
r = 57.29577951 ;
m = SIN a * SIN c ;
l = ( b - e );
l = COS l ;
t = COS a * COS c * l ;
z = m + t ;
z = ACOS z ;
z = z * 1.852 * 60.0 * r ;
dist = z ;
RETURN ;
END ;
\

Start
What is the Origin latitude ?
-37.814107

What is the Origin longitude ?
144.96328

Destination latitude ?
-33.8674869

Destination longitude ?
151.2069902

Origin latitude is -37.81410700
Origin longitude is 144.96328000
Destination latitude is -33.86748690
Destination longitude is 151.20699020
Distance is 712.9 km

Edited by chronic 2015-01-01

BobD

Guru

Joined: 07/12/2011
Location: Australia
Posts: 935

Posted: 04:45am 31 Dec 2014

Copy link to clipboard

Print this post

Google Earth gives it as 713.79 but I wasn't exactly on your co-ordinates. Close enough to verify your calcs.

matherp
Guru

Joined: 11/12/2012
Location: United Kingdom
Posts: 10310

Posted: 04:47am 31 Dec 2014

Copy link to clipboard

Print this post

This formula is supposed to be mathematically equivalent but less subject to rounding errors:

d=2*asin(sqrt((sin((lat1-lat2)/2))^2 +
cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))

cwilt
Senior Member

Joined: 20/03/2012
Location: United States
Posts: 147

Posted: 11:09am 31 Dec 2014

Copy link to clipboard

Print this post

What kind of precision are you looking for? There are several different formulas floating around.

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6283

Posted: 11:35am 31 Dec 2014

Copy link to clipboard

Print this post

These are the functions I use for distance and bearing:

function distance(lat1,long1,lat2,long2)
lon1=long1/57.29577951
lat1=lat1/57.29577951
lon2=long2/57.29577951
lat2=lat2/57.29577951
if (lat1=lat2) and (long1=long2) then
distance = 0
else
distance = acs(sin(lat1)*sin(lat2)+cos(lat1)* cos(lat2)*cos(lon1-lon2))* 6371
endif
end function

function bearing(lat1,long1,lat2,long2)
if (lat1=lat2) and (long1=long2) then
bearing = 0
else
lon1=long1/57.29577951
lat1=lat1/57.29577951
lon2=long2/57.29577951
lat2=lat2/57.29577951
fa=SIN(lon2-lon1)*COS(lat2)
fb=COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1)
bearing=(atn(fa/fb))* 57.29577951
if fb<0 then bearing = 180+ bearing else if fa<0 then bearing = 360 + bearing
endif
end function

Lat and long are in decimal degrees
Degrees south are negative
Degrees West are negative

I have not tested them on recent releases of MMBasic but they should be OK still.
The routines are part of my Satellite location program from a few years ago.

Jim

VK7JH
MMedit

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6283

Posted: 12:30pm 31 Dec 2014

Copy link to clipboard

Print this post

I just realized that I am missing a function for ASC()

This breaks the distance function. All I have to do now is find where I had that function....

With all this type of calculations, the limits of 32 bit floating point maths will come into play.

I have added the missing function and changed the name of distance to dist so it doesn't conflict with a keyword in the micromite's

input "Latitude ";lat1
input "Longitude ";lon1
print
input "Latitude ";lat2
input "Longitude ";lon2

print "Distance = ";dist(lat1,lon1,lat2,lon2)

print "Bearing = ";bearing(lat1,lon1,lat2,lon2)

end

function dist(lat1,long1,lat2,long2)
lon1=long1/57.29577951
lat1=lat1/57.29577951
lon2=long2/57.29577951
lat2=lat2/57.29577951
if (lat1=lat2) and (long1=long2) then
dist = 0
else
dist = acs(sin(lat1)*sin(lat2)+cos(lat1)* cos(lat2)*cos(lon1-lon2))*6371
endif
end function

function bearing(lat1,long1,lat2,long2)
if (lat1=lat2) and (long1=long2) then
bearing = 0
else
lon1=long1/57.29577951
lat1=lat1/57.29577951
lon2=long2/57.29577951
lat2=lat2/57.29577951
fa=SIN(lon2-lon1)*COS(lat2)
fb=COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1)
bearing=(atn(fa/fb))* 57.29577951
if fb<0 then bearing = 180+ bearing else if fa<0 then bearing = 360 + bearing
endif
end function

function acs(c)
acs = 1.5708 - 2*atn(c/(1+sqr(1-c*c )))
end function

DOS MMBasic Version 4.5
Copyright 2011-2014 Geoff Graham

Latitude ? -37.814107
Longitude ? 144.96328

Latitude ? -33.8674869
Longitude ? 151.2069902
Distance = 713.402
Bearing = 57.7023
>

Jim

Edited by TassyJim 2015-01-01

VK7JH
MMedit

isochronic
Guru

Joined: 21/01/2012
Location: Australia
Posts: 689

Posted: 03:11am 01 Jan 2015

Copy link to clipboard

Print this post

Thanks all !
As an alternative I also used The formula (thanks matherp) which is called Haversine's I think, it is said to ok with single precision

DOUBLE FUNCTION dist ( a, b, c, e ) ;
DOUBLE l, m, z, t, r ;
r = 57.29577951 ;
m = SIN ( ( a - c ) / 2.0 ) ;
l = COS ( a ) * COS ( c ) ;
t = SIN ( ( e - b ) / 2.0 ) ;
z = ASIN ( SQRT ( m * m + l * t * t ) ) ;
dist = 2.0 * z * 1.852 * 60.0 * r ;
RETURN ;
END ;

and it gives the same result eg 712.8965 km
curiouser and curioser !! I havent tried tassys formula yet but I have no doubt the result would be the 713 result a-ok.. . The 6371 figure triggers a memory of being the mean earth radius, although it is actually oblate at the equator so the theoretical trig formulas are probably out a bit..i will try to find an "offical" figure , but i am relieved my arithmetic is ok !!
thanks again

ed-
the calculation is the same but...
using 1.85200 (km/ nautical mile)
* 60.0 (minutes/degree)
* 57.29577951 degrees/radian....gives slightly less than 6371 as the factor

Using the 6371 as the factor instead gives 713.38...so that explains it. So i guess the technical spherical-arc distance is the 712, and the more practical oblate-average is the 713.Edited by chronic 2015-01-02

cwilt
Senior Member

Joined: 20/03/2012
Location: United States
Posts: 147

Posted: 07:06am 01 Jan 2015

Copy link to clipboard

Print this post

GIS calculations can be a very deep rabbit hole indeed. Thats why I asked what level of precision you were looking for. My current profession is as a pipeline surveyor and work with Trimble R10 GPS units. Our pipelines typically run over 500 miles and our target accuracy is 1/10 of a foot in that distance. We would often get some weather related down time which gave me time to try different calculation methods and compare to professional measurements.

Haversine is a fast approximation of distance between two points. In my tests it could have errors in excess of 100 feet. Somewhere on one my computers I have a mmbasic library of GPS calculation methods. I will look for it and post it.

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6283

Posted: 12:54pm 01 Jan 2015

Copy link to clipboard

Print this post

This is one of the sources I used when developing my routines:
http://www.movable-type.co.uk/scripts/latlong.html

The same formula used in Liberty Basic which has unlimited precision gives the same distance as MMBasic for DOS but different bearings.
MMBasic is out by a few degrees.

I will find out where the errors occur later - the lawn mower is calling me.

Jim

Edit: I have just been informed that we are taking the dogs out for a coffee, the lawns can wait.
Edited by TassyJim 2015-01-02

VK7JH
MMedit

mpep
Newbie

Joined: 09/11/2014
Location: New Zealand
Posts: 29

Posted: 01:07pm 01 Jan 2015

Copy link to clipboard

Print this post

As there appears to be only 1km difference between the 2, over 700km, you are looking at a variation of 0.14%.

Another approach is to use Vincenty's calculations. I am looking into this now for a project I'm working on.
Supposedly even closer than Haversine.

Your maths examples certainly give me hope that everything I want to do is indeed entirely possible. THANKS.

MPepEdited by mpep 2015-01-02

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6283

Posted: 03:53pm 01 Jan 2015

Copy link to clipboard

Print this post

I think better after a coffee...

This version has fixed the problem with the wrong bearing.
I was changing global variables within the functions.

I have also used the RAD() and DEG() functions where appropriate.

'input "Latitude ";latitude1
'input "Longitude ";longitude1
'print
'input "Latitude ";latitude2
'input "Longitude ";longitude2
latitude1 = -37.814107
longitude1 = 144.96328
latitude2 = -33.8674869
longitude2 = 151.2069902
PRINT "Distance = ";dist(latitude1,longitude1,latitude2,longitude2)
PRINT "Bearing = ";bearing(latitude1,longitude1,latitude2,longitude2)

END

FUNCTION dist(la1,lo1,la2,lo2)
LOCAL lat1,long1,lat2,long2
lat1 =RAD(la1)
long1=RAD(lo1)
lat2 =RAD(la2)
long2=RAD(lo2)
IF (lat1=lat2) AND (long1=long2) THEN
dist = 0
ELSE
dist = acs(SIN(lat1)*SIN(lat2)+COS(lat1)* COS(lat2)*COS(long1-long2))*6371
ENDIF
END FUNCTION

FUNCTION bearing(la1,lo1,la2,lo2)
LOCAL lat1,long1,lat2,long2
lat1 =RAD(la1)
long1=RAD(lo1)
lat2 =RAD(la2)
long2=RAD(lo2)
IF (lat1=lat2) AND (long1=long2) THEN
bearing = 0
ELSE
fa=SIN(long2-long1)*COS(lat2)
fb=COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(long2-long1)
bearing=DEG(ATN(fa/fb))
IF fb<0 THEN bearing = 180+ bearing ELSE IF fa<0 THEN bearing = 360 + bearing
ENDIF
END FUNCTION

FUNCTION acs(c)
acs = 1.5708 - 2*ATN(c/(1+SQR(1-c*c )))
END FUNCTION

Tested with MMBasic DOS, Maximite 4.5 and micromite 4.6

Jim

VK7JH
MMedit

cwilt
Senior Member

Joined: 20/03/2012
Location: United States
Posts: 147

Posted: 05:02pm 01 Jan 2015

Copy link to clipboard

Print this post

mpep said As there appears to be only 1km difference between the 2, over 700km, you are looking at a variation of 0.14%.

Another approach is to use Vincenty's calculations. I am looking into this now for a project I'm working on.
Supposedly even closer than Haversine.

Your maths examples certainly give me hope that everything I want to do is indeed entirely possible. THANKS.

MPep

You are correct, Vincenty formula will be more accuracte.

isochronic
Guru

Joined: 21/01/2012
Location: Australia
Posts: 689

Posted: 05:34pm 01 Jan 2015

Copy link to clipboard

Print this post

Man... the vincenty site quotes half-a-millimeter !!
I can't even see that detail level..Nice though.
I only tried originally as a test of some substantial arithmetic - then was puzzled as to why the result didn't seem right. ( I have driven from Melbourne to Sydney, it is a day trip, kind of boring, the road distance is quite a bit longer )
Incidentally the single-precision calculation gives a very similar result, it is out by a few meters only, unsurprising - six significant digits is six meters in 6000 km of course but it could add up if repeated I guess.
thanks all, and nice work tassy - that acos function is a new one for me
I'll give my eyes a rest now for a while !!!

isochronic
Guru

Joined: 21/01/2012
Location: Australia
Posts: 689

Posted: 07:43pm 03 Jan 2015

Copy link to clipboard

Print this post

Ok...a site with a vincenty calc using an elliptical earth
indeed gives the 713.78 google earth figure, so I'll take that !! (thanks bobd)

After that it gets wacky - the surface depends on which geoid ,
definition of a surface etc...and some very high level stuff I guess NASA uses,
beautiful but I will definitely leave it to the rocket scientists..eesh.

A last thing, in the function (thanks tassy)

FUNCTION acs(c)
acs = 1.5708 - 2*ATN(c/(1+SQR(1-c*c )))
END FUNCTION

the acos is derived from atan as normal, and I guess the subtraction from pi/2 is a chart-rotation -to - trig-rotation conversion...but where is the 2 from ?Edited by chronic 2015-01-05

TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 6283

Posted: 08:33pm 03 Jan 2015

Copy link to clipboard

Print this post

chronic said
A last thing, in the function (thanks tassy)

FUNCTION acs(c)
acs = 1.5708 - 2*ATN(c/(1+SQR(1-c*c )))
END FUNCTION

the acos is derived from atan as normal, and I guess the subtraction from pi/2 is a chart-rotation -to - trig-rotation conversion...but where is the 2 from ?

Somewhere in the dusty text books, I will have the proof but my head hurts when I try to do too much geometry!

I have a very useful book called "The Basic Handbook 2nd edition" by David A. Lien, published in 1981.
It covers all the popular dialects of Basic and, for most functions, has an entry called "If your computer doesn't have it"

I have been referring to it for the last 30 years and strongly recommend it if you can find a second hand copy.

It is not the same as "The IBM Basic handbook" by the same author.
The 3rd edition is available secondhand at good prices.

Jim

VK7JH
MMedit

BobD

Guru

Joined: 07/12/2011
Location: Australia
Posts: 935

Posted: 08:39pm 03 Jan 2015

Copy link to clipboard

Print this post

chronic
when I did the Google Earth measurement before I wasn't spot on your co-ordinates. I just fixed that (to 6 decimal places) and measured again and it came out to 713.81 at a bearing of 53.92. GE seems to be more accurate than I thought.
Bob

isochronic
Guru

Joined: 21/01/2012
Location: Australia
Posts: 689

Posted: 08:44pm 05 Jan 2015

Copy link to clipboard

Print this post

I used some coordinates in degrees-minutes-seconds instead of decimals when I tried the vincenty method, it was probably a little different. And I typed ...78 when I meant ...79, so it all a-ok...maybe the 6371 radius approximation should be nudged a little for aus. The vincenty javascript is available but I'll leave it I think.
Incidentally the last time I looked in the few Sydney bookshops that have anything remotely technical, the cooking section was about three times larger than then all I.T.!..only 1 book even on C...the net has taken over I think.

cdeagle
Senior Member

Joined: 22/06/2014
Location: United States
Posts: 265

Posted: 05:02am 06 Feb 2015

Copy link to clipboard

Print this post

Attached is a zipped file containing an MM BASIC 4.6 implementation of yet another way to calculate the geodesic distance between two locations. It's based on Sodano's algorithm which also determines the azimuth to and from each location.

2015-02-06_145620_demo_geodesic.zip

The following is a typical user interaction with the program for the example given earlier.

program demo_geodesic

first location

please input the geographic latitude
(degrees [-90 to +90], minutes [0 - 60], seconds [0 - 60])
(north latitudes are positive, south latitudes are negative)
? -37,48,50

please input the geographic longitude
(degrees [0 - 360], minutes [0 - 60], seconds [0 - 60])
(east longitude is positive, west longitude is negative)
? 144,57,47

second location

please input the geographic latitude
(degrees [-90 to +90], minutes [0 - 60], seconds [0 - 60])
(north latitudes are positive, south latitudes are negative)
? -33,52,3

please input the geographic longitude
(degrees [0 - 360], minutes [0 - 60], seconds [0 - 60])
(east longitude is positive, west longitude is negative)
? 151,12,25

program demo_geodesic

azimuth 1 to 2 54.0432 degrees
azimuth 2 to 1 230.383 degrees

geodesic distance 713.802 kilometers

redrok

Senior Member

Joined: 15/09/2014
Location: United States
Posts: 209

Posted: 04:44am 13 Feb 2015

Copy link to clipboard

Print this post

Hi TassyJim;
Thanks for the code!!!
Just what I was looking for.
However, not for distance on the Earth.

I want to use it for controlling a Heliostat.

Changing Mean_Radius = 180/PI it gives distance in degrees.
In this case the spheroid calculations are perfect for my needs.

Would it be to much to ask you to add another function for
finding the midpoint between the end points.

Do you have a link for the source of these calculations?
Maybe they have the Bisector equations.

Yes, I am aware of the movable-type reference.
However, they use the "atan2(y,x)" function.

Maybe you can show us how to implement the atan2(y,x) function in MMBasic.

This would eliminate the Divide by 0 problems.

Thanks

redrok

TassyJim said I think better after a coffee...

This version has fixed the problem with the wrong bearing.
I was changing global variables within the functions.

I have also used the RAD() and DEG() functions where appropriate.

Tested with MMBasic DOS, Maximite 4.5 and micromite 4.6

Jim

cdeagle
Senior Member

Joined: 22/06/2014
Location: United States
Posts: 265

Posted: 05:14am 13 Feb 2015

Copy link to clipboard

Print this post

Here's an atan3 function.

function atan3(a, b)

' four quadrant inverse tangent function

' input

' a = sine of angle
' b = cosine of angle

' output

' c = angle ( 0 =< c <= 2 * pi : radians )

pidiv2 = 1.57079632

if (abs(a) < 0.00000001) then

atan3 = (1.0 - sgn(b)) * pidiv2

exit function

else

c = (2.0 - sgn(a)) * pidiv2

endif

if (abs(b) < 0.0000001) then

atan3 = c

exit function

else

atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2)

endif

end function

Page 1 of 2

Print this page