Recent

Author Topic: Vincenty's geodesy formulae  (Read 3910 times)

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Vincenty's geodesy formulae
« on: September 14, 2019, 01:52:11 am »
I need a function to accurately calculate the distance, initial azimuth, and final azimuth between two points defined by latitude and longitude on a reference ellipsoid.

The correct solution is given by the inverse geodetic problem of Vincenty's formulae. I have found implementations in java and C#:

https://github.com/mgavaghan/geodesy

a Fortran version from NOAA:

https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html

and in C++, java, js, matlab, python, (legacy) Fortran, and (legacy) C:

https://geographiclib.sourceforge.io/

The first is likely the easiest to translate, the second is spaghetti code (apologies to the coder, but looks like a lot of work to translate), the third seems the most bombproof in terms of unusual cases, but also looks like a lot of work to translate.

Is anyone aware of solid Delphi/Pascal code for this problem? I don't want to dive in if it is already done.
« Last Edit: September 14, 2019, 02:17:02 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)

jamie

  • Hero Member
  • *****
  • Posts: 6091
Re: Vincenty's geodesy formulae
« Reply #1 on: September 14, 2019, 02:14:13 am »
The standard won't do for you?

 Sqrt(Sqr(number)+Sqr(otherNumber));

 Maybe I am thinking of something else.
The only true wisdom is knowing you know nothing

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Vincenty's geodesy formulae
« Reply #2 on: September 14, 2019, 02:20:32 am »
The standard won't do for you?

 Sqrt(Sqr(number)+Sqr(otherNumber));

 Maybe I am thinking of something else.

Thanks, but no. Pythagoras got it on a plane. Haversine figured it out for a sphere. Vincenty worked it out for an ellipsoid.
“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)

jamie

  • Hero Member
  • *****
  • Posts: 6091
Re: Vincenty's geodesy formulae
« Reply #3 on: September 14, 2019, 02:21:49 am »
Yes I was thinking of the curve .

 I write a Sat tracking app years ago In TO using the bga graphics etc. it was in my opinion a work of art!

These source files may help you, they are source code to Sat tracking and it has a lot o the math you need.

https://celestrak.com/software/tskelso-sw.php
« Last Edit: September 14, 2019, 02:32:07 am by jamie »
The only true wisdom is knowing you know nothing

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Vincenty's geodesy formulae
« Reply #4 on: September 14, 2019, 03:00:16 am »
Thanks, this looks very interesting, but I need an implementation of Vincenty's inverse geodetic problem.
« Last Edit: September 14, 2019, 03:44:13 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: Vincenty's geodesy formulae
« Reply #5 on: September 15, 2019, 12:37:50 am »
Thanks. I got the function I needed with a java translation.
“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)

sstvmaster

  • Sr. Member
  • ****
  • Posts: 299
Re: Vincenty's geodesy formulae
« Reply #6 on: September 15, 2019, 04:32:44 pm »
Hi VTWin,

here i have a simple programm, without reversing.

Can you show me your code, please?

« Last Edit: September 15, 2019, 04:34:24 pm by sstvmaster »
greetings Maik

Windows 10,
- Lazarus 2.2.6 (stable) + fpc 3.2.2 (stable)
- Lazarus 2.2.7 (fixes) + fpc 3.3.1 (main/trunk)

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Vincenty's geodesy formulae
« Reply #7 on: September 15, 2019, 07:25:55 pm »
Here it is. I had to remove some dependencies.

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)

sstvmaster

  • Sr. Member
  • ****
  • Posts: 299
Re: Vincenty's geodesy formulae
« Reply #8 on: September 15, 2019, 07:35:49 pm »
Thank you very much!!! :o
greetings Maik

Windows 10,
- Lazarus 2.2.6 (stable) + fpc 3.2.2 (stable)
- Lazarus 2.2.7 (fixes) + fpc 3.3.1 (main/trunk)

sstvmaster

  • Sr. Member
  • ****
  • Posts: 299
Re: Vincenty's geodesy formulae
« Reply #9 on: September 15, 2019, 10:22:57 pm »
Hi, i found a little bug, if lat or lon are negative floats without N,S,E,W?

wrong calculation:
Code: Pascal  [Select][+][-]
  1.   lat1 = '38.88922'; // Washington DC
  2.   lon1 = '-77.04978';
  3.   lat2 = '48.85889'; // Eiffel Turm
  4.   lon2 = '2.29583';

right calculation:
Code: Pascal  [Select][+][-]
  1.   lat1 = '38.88922 N'; // Washington DC
  2.   lon1 = '-77.04978 W';
  3.   lat2 = '48.85889 N'; // Eiffel Turm
  4.   lon2 = '2.29583 E';

i have make changes to VLatLongUtils.pas.

for both functions LatitudeVal + LongitudeVal
Code: Pascal  [Select][+][-]
  1. function LatitudeVal(const s: string; var rad: double): boolean;
  2. ...
  3. t2: double;
  4. ...
  5.   if not StrToDMS(t, deg, min, sec) then
  6.     exit; // failed
  7.  
  8.   // new check
  9.   if FloatVal(t,t2) then
  10.     if t2 < 0 then
  11.       sign := -1;

is now it think its right
greetings Maik

Windows 10,
- Lazarus 2.2.6 (stable) + fpc 3.2.2 (stable)
- Lazarus 2.2.7 (fixes) + fpc 3.3.1 (main/trunk)

winni

  • Hero Member
  • *****
  • Posts: 3197
Re: Vincenty's geodesy formulae
« Reply #10 on: September 15, 2019, 10:39:55 pm »
From the technical point of view the N/S/E/W suffixes should  not be becessary, because:

Longitude: East is positiv, West is negative
Latitude: North is positive, South is negative



VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Vincenty's geodesy formulae
« Reply #11 on: September 15, 2019, 11:03:14 pm »
Thanks, dumb mistake on my part. Fix:

Code: Pascal  [Select][+][-]
  1. function IsLegal(c: string): boolean;
  2. begin
  3.   result := (c = '+') or (c = '-') or (c = '.') or // excludes comma
  4.     (c = '0') or (c = '1') or (c = '2') or (c = '3') or
  5.     (c = '4') or (c = '5') or (c = '6') or (c = '7') or
  6.     (c = '8') or (c = '9');
  7. end;  

Yes, W, E, N, and S are not necessary. Neither are minutes and seconds. All are annoying, but not every source or person uses signed decimal degrees.
« Last Edit: September 15, 2019, 11:10:30 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)

winni

  • Hero Member
  • *****
  • Posts: 3197
Re: Vincenty's geodesy formulae
« Reply #12 on: September 15, 2019, 11:25:04 pm »
Easier:

Code: Pascal  [Select][+][-]
  1.     function IsLegal(c: string): boolean;
  2.     begin
  3.       result := (Length(c) = 1 ) and ( c[1] in ['+', '-', '.', '0'..'9'] );
  4.     end;
  5.      
  6.  

Winni

sstvmaster

  • Sr. Member
  • ****
  • Posts: 299
Re: Vincenty's geodesy formulae
« Reply #13 on: September 16, 2019, 12:19:01 am »
No, thank you for your work.

...but not every source or person uses signed decimal degrees.

At the end you must use floating point to calculate. For me, i always use floating point.
greetings Maik

Windows 10,
- Lazarus 2.2.6 (stable) + fpc 3.2.2 (stable)
- Lazarus 2.2.7 (fixes) + fpc 3.3.1 (main/trunk)

VTwin

  • Hero Member
  • *****
  • Posts: 1215
  • Former Turbo Pascal 3 user
Re: Vincenty's geodesy formulae
« Reply #14 on: September 16, 2019, 01:32:39 am »
No, thank you for your work.

...but not every source or person uses signed decimal degrees.

At the end you must use floating point to calculate. For me, i always use floating point.

Glad it was useful.

I store all angles as double floating point radians, these routines are for reading and writing files and GUI controls only. If I have an unusual routine that takes degrees, I try to label it clearly with 'Deg' in the name and parameters.
“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)

 

TinyPortal © 2005-2018