### Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook

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

#### VTwin

• Hero Member
• Posts: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### jamie

• Hero Member
• Posts: 2896
##### 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.
Number 1 at blue screen app creations!

#### VTwin

• Hero Member
• Posts: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### jamie

• Hero Member
• Posts: 2896
##### 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 »
Number 1 at blue screen app creations!

#### VTwin

• Hero Member
• Posts: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### VTwin

• Hero Member
• Posts: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### sstvmaster

• Full Member
• Posts: 155
##### 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.

« Last Edit: September 15, 2019, 04:34:24 pm by sstvmaster »
OS: Windows 7 (32 bit) / Windows 10 (64 bit)
Lazarus: 2.0.6 x32 / 2.1.0 Trunk x32

#### VTwin

• Hero Member
• Posts: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### sstvmaster

• Full Member
• Posts: 155
##### Re: Vincenty's geodesy formulae
« Reply #8 on: September 15, 2019, 07:35:49 pm »
Thank you very much!!!
OS: Windows 7 (32 bit) / Windows 10 (64 bit)
Lazarus: 2.0.6 x32 / 2.1.0 Trunk x32

#### sstvmaster

• Full Member
• Posts: 155
##### 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
OS: Windows 7 (32 bit) / Windows 10 (64 bit)
Lazarus: 2.0.6 x32 / 2.1.0 Trunk x32

#### winni

• Hero Member
• Posts: 1341
##### 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: 882
• 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

macOS 10.13.6: Lazarus 2.0.7 fixes svn r62689 (64 bit Cocoa)
Ubuntu 18.04.3: Lazarus 2.0.6 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 2.0.6 (64 bit on VBox)
fpc 3.0.4

#### winni

• Hero Member
• Posts: 1341
##### 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

• Full Member
• Posts: 155
##### 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.
OS: Windows 7 (32 bit) / Windows 10 (64 bit)
Lazarus: 2.0.6 x32 / 2.1.0 Trunk x32

#### VTwin

• Hero Member
• Posts: 882
• 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.