Recent

Author Topic: Here is one I maybe can't find an answer to, Not easily.  (Read 2889 times)

JLWest

  • Hero Member
  • *****
  • Posts: 1293
Here is one I maybe can't find an answer to, Not easily.
« on: December 27, 2018, 04:44:18 am »
I have two points on the globe. Could be any two points. LAX and Hong Kong International.

I know the longitude and latitude of both places. I need to find the distance in miles or natural miles between the two.

I have been searching the web and all I get is calculators where you enter the data and they give you the distance. But I need it in code. And of course my data is, I think decimal format (47.449888889    -122.311777778) SeaTac Seattle.

Anyone know about this. Not holding out much hope and I realize it's kinda of topic.
FPC 3.2.0, Lazarus IDE v2.0.4
 Windows 10 Pro 32-GB
 Intel i7 770K CPU 4.2GHz 32702MB Ram
GeForce GTX 1080 Graphics - 8 Gig
4.1 TB

JLWest

  • Hero Member
  • *****
  • Posts: 1293
Re: It's called the Haversine formula
« Reply #1 on: December 27, 2018, 04:58:01 am »
Been around back to 1805. BC (Before Computers) they hand did the calculations for hundreds of points and published books for ships. Now its built into you Android.

I'll have to give it a go but I don't remember what a cosine is

Let the central angle Θ between any two points on a sphere be:

    Θ = d r {\displaystyle \Theta ={\frac {d}{r}}} {\displaystyle \Theta ={\frac {d}{r}}}

where:

    d is the distance between the two points (along a great circle of the sphere; see spherical distance),
    r is the radius of the sphere.

The haversine formula hav of Θ is given by:

    hav ⁡ ( Θ ) = hav ⁡ ( φ 2 − φ 1 ) + cos ⁡ ( φ 1 ) cos ⁡ ( φ 2 ) hav ⁡ ( λ 2 − λ 1 ) {\displaystyle \operatorname {hav} \left(\Theta \right)=\operatorname {hav} (\varphi _{2}-\varphi _{1})+\cos(\varphi _{1})\cos(\varphi _{2})\operatorname {hav} (\lambda _{2}-\lambda _{1})} {\displaystyle \operatorname {hav} \left(\Theta \right)=\operatorname {hav} (\varphi _{2}-\varphi _{1})+\cos(\varphi _{1})\cos(\varphi _{2})\operatorname {hav} (\lambda _{2}-\lambda _{1})}

where

    φ1, φ2: latitude of point 1 and latitude of point 2,
    λ1, λ2: longitude of point 1 and longitude of point 2.

Finally, the haversine function (half a versine) of an angle θ (applied above to the differences in latitude and longitude) is:

    hav ⁡ ( θ ) = sin 2 ⁡ ( θ 2 ) = 1 − cos ⁡ ( θ ) 2 {\displaystyle \operatorname {hav} (\theta )=\sin ^{2}\left({\frac {\theta }{2}}\right)={\frac {1-\cos(\theta )}{2}}} \operatorname {hav} (\theta )=\sin ^{2}\left({\frac {\theta }{2}}\right)={\frac {1-\cos(\theta )}{2}}

To solve for the distance d, apply the inverse haversine hav-1 to the central angle Θ or use the arcsine (inverse sine) function:

    d = r hav − 1 ⁡ ( h ) = 2 r arcsin ⁡ ( h ) {\displaystyle d=r\operatorname {hav} ^{-1}(h)=2r\arcsin \left({\sqrt {h}}\right)} d=r\operatorname {hav} ^{-1}(h)=2r\arcsin \left({\sqrt {h}}\right)

where h=hav(Θ), or more explicitly:

    d = 2 r arcsin ⁡ ( hav ⁡ ( φ 2 − φ 1 ) + cos ⁡ ( φ 1 ) cos ⁡ ( φ 2 ) hav ⁡ ( λ 2 − λ 1 ) ) {\displaystyle d=2r\arcsin \left({\sqrt {\operatorname {hav} (\varphi _{2}-\varphi _{1})+\cos(\varphi _{1})\cos(\varphi _{2})\operatorname {hav} (\lambda _{2}-\lambda _{1})}}\right)} {\displaystyle d=2r\arcsin \left({\sqrt {\operatorname {hav} (\varphi _{2}-\varphi _{1})+\cos(\varphi _{1})\cos(\varphi _{2})\operatorname {hav} (\lambda _{2}-\lambda _{1})}}\right)}
    = 2 r arcsin ⁡ ( sin 2 ⁡ ( φ 2 − φ 1 2 ) + cos ⁡ ( φ 1 ) cos ⁡ ( φ 2 ) sin 2 ⁡ ( λ 2 − λ 1 2 ) ) {\displaystyle =2r\arcsin \left({\sqrt {\sin ^{2}\left({\frac {\varphi _{2}-\varphi _{1}}{2}}\right)+\cos(\varphi _{1})\cos(\varphi _{2})\sin ^{2}\left({\frac {\lambda _{2}-\lambda _{1}}{2}}\right)}}\right)} {\displaystyle =2r\arcsin \left({\sqrt {\sin ^{2}\left({\frac {\varphi _{2}-\varphi _{1}}{2}}\right)+\cos(\varphi _{1})\cos(\varphi _{2})\sin ^{2}\left({\frac {\lambda _{2}-\lambda _{1}}{2}}\right)}}\right)}

When using these formulae, one must ensure that h does not exceed 1 due to a floating point error (d is only real for h from 0 to 1). h only approaches 1 for antipodal points (on opposite sides of the sphere)—in this region, relatively large numerical errors tend to arise in the formula when finite precision is used. Because d is then large (approaching πR, half the circumference) a small error is often not a major concern in this unusual case (although there are other great-circle distance formulas that avoid this problem). (The formula above is sometimes written in terms of the arctangent function, but this suffers from similar numerical problems near h = 1.)

As described below, a similar formula can be written using cosines (sometimes called the spherical law of cosines, not to be confused with the law of cosines for plane geometry) instead of haversines, but if the two points are close together (e.g. a kilometer apart, on the Earth) you might end up with cos(d/R) = 0.99999999, leading to an inaccurate answer. Since the haversine formula uses sines, it avoids that problem.

I don't think I have much chance here.
WOW

FPC 3.2.0, Lazarus IDE v2.0.4
 Windows 10 Pro 32-GB
 Intel i7 770K CPU 4.2GHz 32702MB Ram
GeForce GTX 1080 Graphics - 8 Gig
4.1 TB

bytebites

  • Hero Member
  • *****
  • Posts: 640
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #2 on: December 27, 2018, 05:12:11 am »

engkin

  • Hero Member
  • *****
  • Posts: 3112
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #3 on: December 27, 2018, 05:15:02 am »
Not tested!

Derived from:
Cos(Shortest)=Sin(Lat1)Sin(Lat2)+Cos(Lat1)Cos(Lat2)Cos(Lon1-Lon2)

The params and result are in radians:
Code: Pascal  [Select][+][-]
  1. function ShortestPath(Lat1, Lon1, Lat2, Lon2: single): single;
  2. begin
  3.   Result:=Sin(Lat1)*Sin(Lat2)+Cos(Lat1)*Cos(Lat2)*Cos(Lon1-Lon2);
  4.   Result :=ArcCos(Result);
  5. end;

Multiply the result with earth radius to get the distance.

You can use RadToDeg/DegToRad from unit Math to convert between Radians and Degrees.

JLWest

  • Hero Member
  • *****
  • Posts: 1293
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #4 on: December 27, 2018, 05:43:34 am »
Wow That can't possibly work.

Found a site that wrote the haversine in several languages and one was Pascal.  And you are proposing two line of code.

Oh, another thing. There is a raging debate on the circumferences of the earth. Two camps with a difference of 1 mile. You would think they would just pace it off and get me an answer. I think it varies depending on the route. You go up Everest to go around it's gonna be longer.

Thanks

Code: Pascal  [Select][+][-]
  1. Program HaversineDemo(output);
  2.  
  3. uses
  4.   Math;
  5.  
  6. function haversineDist(th1, ph1, th2, ph2: double): double;
  7.   const
  8.    diameter = 2 * 6372.8;
  9.   var
  10.     dx, dy, dz: double;
  11.   begin
  12.     ph1 := degtorad(ph1 - ph2);
  13.     th1 := degtorad(th1);
  14.     th2 := degtorad(th2);
  15.  
  16.     dz := sin(th1) - sin(th2);
  17.     dx := cos(ph1) * cos(th1) - cos(th2);
  18.     dy := sin(ph1) * cos(th1);
  19.     haversineDist := arcsin(sqrt(dx**2 + dy**2 + dz**2) / 2) * diameter;
  20.   end;
  21.  
  22. begin
  23.   writeln ('Haversine distance: ', haversineDist(36.12, -86.67, 33.94, -118.4):7:2, ' km.');
  24. end.
FPC 3.2.0, Lazarus IDE v2.0.4
 Windows 10 Pro 32-GB
 Intel i7 770K CPU 4.2GHz 32702MB Ram
GeForce GTX 1080 Graphics - 8 Gig
4.1 TB

engkin

  • Hero Member
  • *****
  • Posts: 3112
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #5 on: December 27, 2018, 06:17:50 am »
I think it varies depending on the route. You go up Everest to go around it's gonna be longer.
I guess using the average of the altitude of the trip would give correct distance. Of course: AverageAltitude+EarthRadius.

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #6 on: December 27, 2018, 07:02:15 am »
Maybe already answered, but here is the code I use:

Code: Pascal  [Select][+][-]
  1. { LatLongToDistance
  2.   Uses the Haversine formula to calculate the great-circle distance between
  3.   two points. Bearing is the initial bearing. Uses mean spherical earth radius.
  4.   Angles are in radians.
  5.   Ref: http://www.movable-type.co.uk/scripts/latlong.html }
  6. procedure LatLongToDistance(const lat0, long0, lat1, long1: double;
  7.   var dist, bearing: double);
  8. var
  9.   R: double = 6371000; // earth radius in meters
  10.   dlat, dlong, slat, slong, a, c, x, y: double;
  11. begin
  12.   // dist
  13.   dlat := lat1 - lat0;
  14.   dlong := long1 - long0;
  15.   slat := Sin(0.5 * dlat);
  16.   slong := Sin(0.5 * dlong);
  17.   a := slat * slat + Cos(lat0) * Cos(lat1) * slong * slong;
  18.   c := 2.0 * ArcTan2(Sqrt(a), Sqrt(1.0-a));
  19.   dist := R * c;
  20.   // bearing
  21.   y := Sin(long1 - long0) * Cos(lat1);
  22.   x := Cos(lat0) * Sin(lat1) - Sin(lat0) * Cos(lat1) * Cos(long1 - long0);
  23.   bearing := ArcTan2(y, x);
  24. end;

Distance in SI units, meters, angles in radians. Note that initial bearing is different than final, reverse start and end points and add pi for final bearing.
« Last Edit: December 27, 2018, 07:54:25 am by VTwin »
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 12.1: Lazarus 2.2.6 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 2.2.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.2.6 (64 bit on VBox)

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #7 on: December 27, 2018, 07:47:21 am »
An alternate, analytical geometry, solution is to convert the spherical coordinates to direction cosines, take the dot product, and use the spherical earth radius to get the arc length. That might be faster than the trigonometric solution I gave.
« Last Edit: December 27, 2018, 07:52:08 am by VTwin »
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 12.1: Lazarus 2.2.6 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 2.2.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.2.6 (64 bit on VBox)

dsyrios

  • Jr. Member
  • **
  • Posts: 57
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #8 on: December 27, 2018, 10:44:08 am »
Take a look at a mine old  simple sample.
It's just a demo with no many explanations.
Maybe take some idea
Laz 2.0.2/FPC 3.0.4/ win10/64
Laz 2.0.0/FPC 3.0.4/ mint-mate 19.1

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #9 on: December 27, 2018, 05:36:14 pm »
Code: Pascal  [Select][+][-]
  1. const
  2.   LatLAX  = 33.9425222;
  3.   LongLAX = -118.4071611;
  4.   LatHKG  = 22.3088888889;
  5.   LongHKG = 113.9144444444;
  6.  
  7. function Calc: double;
  8. var
  9.   d, b: double;
  10.   lat0, long0, lat1, long1: double;
  11. begin
  12.   lat0 := DegToRad(latLAX);
  13.   long0 := DegToRad(longLAX);
  14.   lat1 := DegToRad(latHKG);
  15.   long1 := DegToRad(longHKG);
  16.   LatLongToDistance(lat0, long0, lat1, long1, d, b);
  17.   result := 0.001 * d  // 11664.569 km for spherical Earth radius = 6371.000 km
  18. end;
  19.  

http://edwilliams.org/gccalc.htm

has calculations for various Earth models. The model used here is FAI (Fédération Aéronautique Internationale) Sphere.

This model gives 11664.569 km from LAX to HKG.
« Last Edit: December 27, 2018, 06:35:17 pm by VTwin »
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 12.1: Lazarus 2.2.6 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 2.2.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.2.6 (64 bit on VBox)

JLWest

  • Hero Member
  • *****
  • Posts: 1293
Re: Here is one I maybe can't find an answer to, Not easily.
« Reply #10 on: December 27, 2018, 11:27:28 pm »
Well I guess I got an answer to that one. You know back in the day BC (Before Computers) they should have invented computers and saved all that time do the horrendous calc.
FPC 3.2.0, Lazarus IDE v2.0.4
 Windows 10 Pro 32-GB
 Intel i7 770K CPU 4.2GHz 32702MB Ram
GeForce GTX 1080 Graphics - 8 Gig
4.1 TB

 

TinyPortal © 2005-2018