49 50 Ver6.09.hp Geodesic distance & Earth Euclidean distance calculator, bearing

[Gil]

For HP49-HP50, programs to calculate exact distances (error less than 0.5 mm) on the earth surface.

Last version; 6.08, January 7th 2024.

Copy the file-directory ending by .doc or by .hp.

Inside are two programs:
1) P1P2—> (indirect method/problem)
2) P1—>P2 (direct method/problem)

Both programs give very accurate answers, for example the calculated distance is exact to 0.5 mm.

Let's assume you have two points P1 and P2 on the earth geoid.

1) To get the shortest distance between them, enter the following known 4 "objects" in the stack: lat1 and lon1 in the stack, for the 1st point P1, then for the 2nd point P2 lat2 and lon2 again in the stack, all of them in the format DD.mmsssss (be careful, not decimal degrees or radians).

Then press P1P2—>

After a few seconds you get the distance s, the forward azimuth alpha12 (initial bearing) and the backward azimuth alpha21.

2) The direct method/problem: you enter the coordinates of P1 in then stack, i. e. lat1 and lon2, then enter the chosen distance s in meters and the forward azimuth alpha12 (initial bearing) in DD.mmsssss.

Then press P1—>P2.

After a few seconds you get the coordinates of the new point P2, i.e. lat2 and lon2, as well as the backward azimuth alpha 21.

Iterations are made according to Vincenty's formulae.
Be careful: in this version of the indirect method/problem, the program P1P2—> might not find a solution for (near) antipodal points P1/P2; but those are special, rare cases.

You can change the radius a and b according to your preferred geoid; just modify the values in both programs P1P2—> and P1—>P2 and choose always meters (and not km or other units).

All "degrees" results are given in DD.mmssss, so that you can check the results of one program directly entering those found values for the other program.

Remarks are welcome.

Gil

[SlideRule]

A little more info, please: what format/font is the DOC file?

SlideRule

[Gil]

You can't read it without a HP49-50G calculator.

You just have to copy the doc/file/repertory (ending by doc) included in my first post into your computer, let us say in path "... document".

Then you use the EMU48 on your phone and execute "... load object".

And retrieve the "object" from your computer located in "path".

And there you are, with the programs written in RPN and algebraic mode.

Regards,
Gil

[Gil]

Here are the full codes of both programs included in the file/directory "... doc".

P1P2—>
\<< "lat1 lon1 lat2 lon2
in D.mmss
S < 0
W < 0

Vincenty fails for
nearly antipode pts
" DROP 'lon2' STO 'lat2' STO 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG lat2 "lat2 D.mmss" \->TAG lon2 "lon2 D.mmss" \->TAG RAD 6378137 6356752.3142 'b' STO 'a' STO lat1 D\->RAD 'lat1' STO lat2 D\->RAD 'lat2' STO lon1 D\->RAD 'lon1' STO lon2 D\->RAD 'lon2' STO a b - a / 'f' STO lon2 lon1 - DUP '\Gl' STO 'l' STO 2 \pi * \->NUM '\Gl\180' STO 'ATAN((1-f)*TAN(lat1))' \->NUM 'u1' STO 'ATAN((1-f)*TAN(lat2))' \->NUM 'u2' STO
WHILE \Gl \Gl\180 - ABS .000000000001 >
REPEAT '\v/((COS(u2)*SIN(\Gl))^2+(COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl))^2)' \->NUM 'SIN.\Gs' STO 'SIN(u1)*SIN(u2)+COS(u1)*COS(u2)*COS(\Gl)' \->NUM 'COS.\Gs' STO SIN.\Gs COS.\Gs ATAN2 '\Gs' STO 'COS(u1)*COS(u2)*SIN(\Gl)/SIN(\Gs)' EVAL 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO
IF COS\178.\Ga 0 \=/
THEN 'COS.\Gs-2*SIN(u1)*SIN(u2)/COS\178.\Ga' EVAL
ELSE 0
END 'COS.2\Gsm' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO \Gl '\Gl\180' STO 'l+(1-C)*f*SIN.\Ga*(\Gs+C*SIN.\Gs*(COS.2\Gsm+C*COS.\Gs*(-1+2*SQ(COS.2\Gsm))))' \->NUM '\Gl' STO
END 'COS\178.\Ga*(a^2-b^2)/b^2' EVAL 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' \->NUM 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 'B*SIN.\Gs*(COS.2\Gsm+B/4*(COS.\Gs*(-1+2*COS.2\Gsm^2)-B/6*COS.2\Gsm*(-3+4*SIN.\Gs^2)*(-3+4*COS.2\Gsm^2)))' \->NUM '\GD\Gs' STO 'b*A*(\Gs-\GD\Gs)' \->NUM 's' STO s 1000 / "km" \->TAG 'COS(u2)*SIN(\Gl)' \->NUM 'COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl)' \->NUM ATAN2 \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga1' STO "\Ga12 D.mmss\166 \|^\-> +90" \->TAG 'COS(u1)*SIN(\Gl)' \->NUM '-SIN(u1)*COS(u2)+COS(u1)*SIN(u2)*COS(\Gl)' \->NUM ATAN2 \pi + \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga2' STO "\Ga21 D.mmss\166 \|^\-> +90" \->TAG { \GD\Gs B A u\178 C COS.2\Gsm COS\178.\Ga SIN.\Ga \Gs COS.\Gs SIN.\Gs u2 u1 \Gl\180 l \Gl } PURGE lat1 RAD\->D 'lat1' STO lon1 RAD\->D 'lon1' STO lat2 RAD\->D 'lat2' STO lon2 RAD\->D 'lon2' STO
\>>

P1—>P2
\<< "lat1 lon1 s \Ga1
lat1 lon1 \Ga1 \166\|^\->90
in D.mmss
S < 0
W < 0
dist s in m
" DROP RAD '\Ga1' STO 's' STO 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG s "s [m]" \->TAG \Ga1 "\Ga12 D.mmss\166 \|^\-> +90" \->TAG 6378137 'a' STO 6356752.3142 'b' STO lat1 D\->RAD 'lat1' STO lon1 D\->RAD 'lon1' STO \Ga1 D\->RAD '\Ga1' STO a b - a / 'f' STO '(1-f)*TAN(lat1)' \->NUM 'TAN.u1' STO '1/\v/(1+TAN.u1^2)' \->NUM 'COS.u1' STO 'TAN.u1*COS.u1' \->NUM 'SIN.u1' STO TAN.u1 \Ga1 COS ATAN2 '\Gs1' STO 'COS.u1*SIN(\Ga1)' 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO 'COS\178.\Ga*(a^2-b^2)/b^2' \->NUM 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 's/(b*A)' \->NUM '\Gs' STO \Gs .01 - '\Gs\180' STO
WHILE \Gs \Gs\180 - ABS .000000000001 >
REPEAT 'COS(2*\Gs1+\Gs)' EVAL 'COS2.\Gsm' STO 'B*SIN(\Gs)*(COS2.\Gsm+B/4*(COS(\Gs)*(-1+2*COS2.\Gsm^2.)-B/6*COS2.\Gsm*(-3+4*SIN(\Gs)^2)*(-3+4*COS2.\Gsm^2)))' EVAL '\GD\Gs' STO \Gs '\Gs\180' STO 's/(b*A)+\GD\Gs' EVAL '\Gs' STO
END 'SIN.u1*COS(\Gs)+COS.u1*SIN(\Gs)*COS(\Ga1)' \->NUM '(1-f)*\v/(SIN.\Ga^2+(SIN.u1*SIN(\Gs)-COS.u1*COS(\Gs)*COS(\Ga1))^2)' \->NUM ATAN2 RAD\->D 'lat2' STO lat2 "lat2 D.mmss" \->TAG 'SIN(\Gs)*SIN(\Ga1)' \->NUM 'COS.u1*COS(\Gs)-SIN.u1*SIN(\Gs)*COS(\Ga1)' \->NUM ATAN2 '\Gl' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO '\Gl-(1-C)*f*SIN.\Ga*(\Gs+C*SIN(\Gs)*(COS2.\Gsm+C*COS(\Gs)*(-1+2*COS2.\Gsm^2)))' \->NUM 'l' STO lon1 l + RAD\->D 'lon2' STO lon2 "lon2 D.mmss" \->TAG SIN.\Ga '-SIN.u1*SIN(\Gs)+COS.u1*COS(\Gs)*COS(\Ga1)' \->NUM ATAN2 \pi + \->NUM RAD\->D HMS\-> 360 MOD \->HMS '\Ga2' STO \Ga2 "\Ga21 D.mmss\166 \|^\-> +90" \->TAG { \Gs\180 COS2.\Gsm \GD\Gs B A u\178 COS\178.\Ga SIN.\Ga \Gs1 SIN.u1 COS.u1 TAN.u1 l C \Gl \Gs } PURGE lat1 RAD\->D 'lat1' STO lon1 RAD\->D 'lon1' STO \Ga1 RAD\->D '\Ga1' STO
\>>

And the two subprograms :
RAD—>D
\<< 180 * \pi / \->NUM \->HMS
\>>

D—>RAD
\<< HMS\-> \pi * 180 / \->NUM
\>>

Regards,
Gil

[SlideRule]

Thanks all, most helpful!

SlideRule

[Gil]

Two more points

A) Iin the previous transcription (given by the program inout):

Gs stands for "Greek s", i. e. "sigma";
Ga stands for "Greek a", i. e. "alpha" ;
Gl stands for "Greek l", i. e. "lambda".

B) To work, when entering the four values in the "stack", the calculator should be in RPN mode.

Regards,
Gil

[Gil]

I forgot to include

the subprogram ATAN2, whose location is the utmost parent directory HOME :

ATAN2

\<< \-> y x
\<<
CASE x 0 >
THEN '2*ATAN(y/(\v/(x^2+y^2)+x))'
END x 0 \<= y 0 \=/ AND
THEN '2*ATAN((\v/(x^2+y^2)-x)/y)'
END x 0 < y 0 == AND
THEN '4*ATAN(1)'
END x 0 == y 0 == AND
THEN "Undef"
END
END \->NUM
\>>
\>>

You might choose the alternative formula and program:

ATAN2

\<< \-> y x
\<<
CASE x 0 >
THEN y x / ATAN
END x 0 < y 0 \>= AND
THEN y x / ATAN -17 FS?
IF
THEN \pi
ELSE 180
END +
END x 0 < y 0 < AND
THEN y x / ATAN -17 FS?
IF
THEN \pi
ELSE 180
END -
END x 0 == y 0 > AND
THEN -17 FS?
IF
THEN \pi 2 /
ELSE 90
END
END x 0 == y 0 < AND
THEN -17 FS?
IF
THEN \pi NEG 2 /
ELSE -90
END
END "Undef"
END \->NUM
\>>
\>>

Note
\v/ stands for square root.

Regards,
Gil

[Gil]

Shorter version for
ATAN2

\<< \-> y x
\<<
CASE x 0 >
THEN '2*ATAN(y/(\v/(x^2+y^2)+x))'
END x 0 \<= y 0 \=/ AND
THEN '2*ATAN((\v/(x^2+y^2)-x)/y)'
END x 0 < y 0 == AND
THEN '4*ATAN(1)'
END "Undef"
END \->NUM
\>>
\>>


And for the alternative formula, taking into account the three possible
"formats of angles", i. e. DEG/RAD/GRAD:
ATAN2


\<< \-> y x
\<<
CASE x 0 >
THEN y x / ATAN
END x 0 < y 0 \>= AND
THEN y x / ATAN -17 FS?
IF
THEN \pi
ELSE
IF -18 FC?
THEN 180
ELSE 200
END
END +
END x 0 < y 0 < AND
THEN y x / ATAN -17 FS?
IF
THEN \pi
ELSE
IF -18 FC?
THEN 180
ELSE 200
END
END -
END x 0 == y 0 > AND
THEN -17 FS?
IF
THEN \pi 2 /
ELSE
IF -18 FC?
THEN 90
ELSE 100
END
END
END x 0 == y 0 < AND
THEN -17 FS?
IF
THEN \pi NEG 2 /
ELSE
IF -18 FC?
THEN -90
ELSE -100
END
END
END "Undef"
END \->NUM
\>>
\>>

Regards,
Gil

[SammysHP]

atan2 is available as the built-in command ARG.

[Gil]

I did not know... and was surprised not to find that useful "function".

Thanks for your observation.

I think that I will adapt the program accordingly.

Regards,
Gil

[Gil]

RE: HP49-50G Geodesic distance calculator
Version 2

Added in this new version 2 the possibility to introduce the data relative to two points under the form of two complex numbers (lat1, lon1) (lat2, lon2), instead of only (in the first version) lat1 lon1 lat2 lon2.

Furthermore, it is now possible to calculate in one single step a journey between several points that you must have created and saved previously.
Suppose you want to calculate the distance between New York, Chicago and Berlin : 1) go to Directory DATA.DIST inside the current distance directory; 2) create there the three points in question (latitude first and longitude after latitude, both in DD.mmssss, minus sign for South-Latitude or for West-Longitude, R—>C to transform th lat_i lon_i into one point/complex number) 40.4651 -73.5838 R—>C 'N.York' STO 41.5255 -87.3723 R—>C 'Chicago' STO 52.3112 13.2418 R—>C 'Berlin' STO; 3) then select the three (less or more) points with {} and the way points names inside {N.York Chicago Berlin}; 4) then run/execute the program on the left screen P1..PN—>Σs.

As explained above, the DATA.DIST directory contains your way points A1, A2,... An, a way point Ai having the form (lat_i, lon_i), lat_i and lon_i in DD.mmssss, and the program P1..PN—>Σs to calculate the separate tracks/distances.

P1.. PN—>Sigma_s
\<< DUPDUP SIZE { } \-> l0 siz l.s
\<< l0 EVAL siz \->LIST 'l0' STO { } 'l.s' STO 1 siz 1 -
FOR i l0 i GET l0 i 1 + GET UPDIR P1P2\->D DATA.DIST 7 DROPN l.s s + 'l.s' STO
NEXT l.s DUP
IF siz 2 >
THEN \GSLIST
END "Total [m]" \->TAG 0 FIX
\>>
\>>

Used, as suggested by SammysHP, instead of my cumbersome and slow y x ATAN2 instructions (ATAN2 being a program to be created), the instructions (x y) ARC, as ARC as already built-in function in the calculator being about 30 times faster than the ATAN2 to be created in inserted.

The new codes are the following:

P1P2—>D
\<< STD "lat1 lon1 lat2 lon2
in D.mmss
S < 0
W < 0

Vincenty fails for
nearly antipode pts
" DROP DUP TYPE 1 ==
IF
THEN OBJ\->
END 'lon2' STO 'lat2' STO DUP TYPE 1 ==
IF
THEN OBJ\->
END 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG lat2 "lat2 D.mmss" \->TAG lon2 "lon2 D.mmss" \->TAG RAD 6378137 6356752.3142 'b' STO 'a' STO lat1 D\->RAD 'lat1' STO lat2 D\->RAD 'lat2' STO lon1 D\->RAD 'lon1' STO lon2 D\->RAD 'lon2' STO a b - a / 'f' STO lon2 lon1 - DUP '\Gl' STO 'l' STO 2 \pi * \->NUM '\Gl\180' STO 'ATAN((1-f)*TAN(lat1))' \->NUM 'u1' STO 'ATAN((1-f)*TAN(lat2))' \->NUM 'u2' STO
WHILE \Gl \Gl\180 - ABS .000000000001 >
REPEAT '\v/((COS(u2)*SIN(\Gl))^2+(COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl))^2)' \->NUM 'SIN.\Gs' STO 'SIN(u1)*SIN(u2)+COS(u1)*COS(u2)*COS(\Gl)' \->NUM 'COS.\Gs' STO COS.\Gs SIN.\Gs R\->C ARG '\Gs' STO 'COS(u1)*COS(u2)*SIN(\Gl)/SIN(\Gs)' EVAL 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO
IF COS\178.\Ga 0 \=/
THEN 'COS.\Gs-2*SIN(u1)*SIN(u2)/COS\178.\Ga' EVAL
ELSE 0
END 'COS.2\Gsm' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO \Gl '\Gl\180' STO 'l+(1-C)*f*SIN.\Ga*(\Gs+C*SIN.\Gs*(COS.2\Gsm+C*COS.\Gs*(-1+2*SQ(COS.2\Gsm))))' \->NUM '\Gl' STO
END 'COS\178.\Ga*(a^2-b^2)/b^2' EVAL 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' \->NUM 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 'B*SIN.\Gs*(COS.2\Gsm+B/4*(COS.\Gs*(-1+2*COS.2\Gsm^2)-B/6*COS.2\Gsm*(-3+4*SIN.\Gs^2)*(-3+4*COS.2\Gsm^2)))' \->NUM '\GD\Gs' STO 'b*A*(\Gs-\GD\Gs)' \->NUM 's' STO s 1000 / "km" \->TAG 'COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl)' \->NUM 'COS(u2)*SIN(\Gl)' \->NUM R\->C ARG \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga1' STO "\Ga12 D.mmss\166 \|^\-> +90" \->TAG '-SIN(u1)*COS(u2)+COS(u1)*SIN(u2)*COS(\Gl)' \->NUM 'COS(u1)*SIN(\Gl)' \->NUM R\->C ARG \pi + \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga2' STO "\Ga21 D.mmss\166 \|^\-> +90" \->TAG { \GD\Gs B A u\178 C COS.2\Gsm COS\178.\Ga SIN.\Ga \Gs COS.\Gs SIN.\Gs u2 u1 \Gl\180 l \Gl } PURGE lat1 RAD\->D 'lat1' STO lon1 RAD\->D 'lon1' STO lat2 RAD\->D 'lat2' STO lon2 RAD\->D 'lon2' STO
\>>

P1—>P2
\<< STD "lat1 lon1 s \Ga1
lat1 lon1 \Ga1 \166\|^\->90
in D.mmss
S < 0
W < 0
dist s in m
" DROP RAD '\Ga1' STO 's' STO 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG s "s [m]" \->TAG \Ga1 "\Ga12 D.mmss\166 \|^\-> +90" \->TAG 6378137 'a' STO 6356752.3142 'b' STO lat1 D\->RAD 'lat1' STO lon1 D\->RAD 'lon1' STO \Ga1 D\->RAD '\Ga1' STO a b - a / 'f' STO '(1-f)*TAN(lat1)' \->NUM 'TAN.u1' STO '1/\v/(1+TAN.u1^2)' \->NUM 'COS.u1' STO 'TAN.u1*COS.u1' \->NUM 'SIN.u1' STO \Ga1 COS TAN.u1 R\->C ARG '\Gs1' STO 'COS.u1*SIN(\Ga1)' \->NUM 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO 'COS\178.\Ga*(a^2-b^2)/b^2' \->NUM 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 's/(b*A)' \->NUM '\Gs' STO \Gs .01 - '\Gs\180' STO
WHILE \Gs \Gs\180 - ABS .000000000001 >
REPEAT 'COS(2*\Gs1+\Gs)' EVAL 'COS2.\Gsm' STO 'B*SIN(\Gs)*(COS2.\Gsm+B/4*(COS(\Gs)*(-1+2*COS2.\Gsm^2.)-B/6*COS2.\Gsm*(-3+4*SIN(\Gs)^2)*(-3+4*COS2.\Gsm^2)))' EVAL '\GD\Gs' STO \Gs '\Gs\180' STO 's/(b*A)+\GD\Gs' EVAL '\Gs' STO
END '(1-f)*\v/(SIN.\Ga^2+(SIN.u1*SIN(\Gs)-COS.u1*COS(\Gs)*COS(\Ga1))^2)' \->NUM 'SIN.u1*COS(\Gs)+COS.u1*SIN(\Gs)*COS(\Ga1)' \->NUM R\->C ARG RAD\->D 'lat2' STO lat2 "lat2 D.mmss" \->TAG 'COS.u1*COS(\Gs)-SIN.u1*SIN(\Gs)*COS(\Ga1)' \->NUM 'SIN(\Gs)*SIN(\Ga1)' \->NUM R\->C ARG '\Gl' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO '\Gl-(1-C)*f*SIN.\Ga*(\Gs+C*SIN(\Gs)*(COS2.\Gsm+C*COS(\Gs)*(-1+2*COS2.\Gsm^2)))' \->NUM 'l' STO lon1 l + RAD\->D 'lon2' STO lon2 "lon2 D.mmss" \->TAG '-SIN.u1*SIN(\Gs)+COS.u1*COS(\Gs)*COS(\Ga1)' \->NUM SIN.\Ga R\->C ARG \pi + \->NUM RAD\->D HMS\-> 360 MOD \->HMS '\Ga2' STO \Ga2 "\Ga21 D.mmss\166 \|^\-> +90" \->TAG { \Gs\180 COS2.\Gsm \GD\Gs B A u\178 COS\178.\Ga SIN.\Ga \Gs1 SIN.u1 COS.u1 TAN.u1 l C \Gl \Gs } PURGE lat1 RAD\->D 'lat1' STO lon1 RAD\->D 'lon1' STO \Ga1 RAD\->D '\Ga1' STO
\>>

Besides, in that file/directory dist a new program lat—>R has been created to determinate the Earth Radius and the Earth circumference in function of the latitude:

lat—>R
\<< DEG DUP 'lat' STO HMS\-> \-> la
\<< '\v/(((a^2*COS(la))^2+(b^2*SIN(la))^2)/((a*COS(la))^2+(b*SIN(la))^2))' \->NUM "R Earth(" lat + ")" + \->TAG DUP la COS * 2 * \pi * \->NUM "Circumf(" lat + ")" + \->TAG
\>>
\>>

Nothing changed in

RAD—>D
\<< 180 * \pi / \->NUM \->HMS
\>>

D—>RAD
\<< HMS\-> \pi * 180 / \->NUM
\>>

Notes
All the latitudes, longitudes and angles have, as defore, always to be in DD.mmssss.
And programs are to be used/run in RPN mode.

Remarks welcome.

Regards,
Gil

[Gil]

Just added explanations for the arguments to be entered.

[Gil]

New version 2.3

As the circumference calculation
(lat—>R)
on a constant latitude
was done with the wrong formula.

Now it matches when checking with Visenty's formula.

Regards

[Gil]

Version 2.5

For that earth distance calculator — ccurateness within 1 mm —,
I just added here the possibility to modulate the radius equatorial radius a, the pole radius b and flatness f (or INV.f),
knowing that (a-b)/a = INV.f

(Note that you cannot enter the value of the flatness f,
but instead its inverse value INV. f (=1/f).
If you know f — instead of 1/f —, then put f in the stack and press the key 1/x (juste above the 9-key).)

For that, if changes really required, just run before everything
the small "solver" program pressing
1) C.abc ENTER

2) Then give two values:
- for a and b, then press Left-Shift INV.f
(to get INV.f = 1/f) ;
- or for a and INV.f, then press LEFT-SHIFT b
(to get the pole radius b) ;
- or for b and INV.f, then press Left-Shift a
(to get the equatorial radius a).

3) Then run normally (as many times you want)
- P1P2—>
- or P1—>P2
- or lat—>R

without having to repeat anymore the steps 1 and 2 above.

Code for C.abc
\<< "WGS84
a:6378137
INV.f:298.257223563
" DROP '(a-b)/a-1/INV.f' STEQ 30 MENU
\>>

Enjoy and regards,

Gil Campart

[Gil]

Small correction
in prog lat—>R

\<< "1 single Argument:

lat
 in D.mmss

To 
change a, b, INV.f


run
 before all 
C.abf

" DROP DEG DUP 'lat' STO HMS\-> \-> la
\<< '\v/(((a^2*COS(la))^2+(b^2*SIN(la))^2)/((a*COS(la))^2+(b*SIN(la))^2))' \->NUM "R Centre-Earth(" lat + ")" + \->TAG 'b/a*TAN(la)' \->NUM ATAN COS a * 2 * \pi * \->NUM "Circumf(" lat + ")" + \->TAG
\>>
\>>
You should have only one DROP instruction before DEG.

Regards,
Gil

[Gil]

Slight 'appearance" change of the order of the variables.

Here is the full Directory.

Copy it in your calculator
under the file/Directory name DISTANCE à.

Then with your calculator enter your new directory DISTANCE pressing DISTANCE ENTER.

Regards,
Gil

[Gil]

New version 2.8b

When entering, in D.mmsss
Lat1 Long1 Lat2 Long2

and pressing P1P2—>

Then the above program
used to convert the given lat/long
into decimal degrees, then into radian.

And at the end it converted back the radian lat/long into the initial D.mmsss format, losing, sometimes, some of the original lat/long values.

Now this has been corrected.

The same applied to the P1—>P2 program with the given input lat1, lon1 and alpha degree-values.

Here are the new codes:

P1P2—>
\<< "4 Arguments

lat1 lon1 lat2 lon2


or 2 Complex Arguments

(lat1,lon1)(lat2,lon2)


\[] all in D.mmss
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

\[]To 
change a, b, INV.f


run
 before all 
C.abf

" DROP STD DUP TYPE 1 ==
IF
THEN OBJ\->
END 'lon2' STO 'lat2' STO DUP TYPE 1 ==
IF
THEN OBJ\->
END 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG lat2 "lat2 D.mmss" \->TAG lon2 "lon2 D.mmss" \->TAG RAD INV.f INV 'f' STO lat1 D\->RAD lon1 D\->RAD lat2 D\->RAD lon2 D\->RAD \-> lat1 lon1 lat2 lon2
\<< lon2 lon1 - DUP '\Gl' STO 'l' STO 2 \pi * \->NUM '\Gl\180' STO 'ATAN((1-f)*TAN(lat1))' \->NUM 'u1' STO 'ATAN((1-f)*TAN(lat2))' \->NUM 'u2' STO
WHILE \Gl \Gl\180 - ABS .000000000001 >
REPEAT '\v/((COS(u2)*SIN(\Gl))^2+(COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl))^2)' \->NUM 'SIN.\Gs' STO 'SIN(u1)*SIN(u2)+COS(u1)*COS(u2)*COS(\Gl)' \->NUM 'COS.\Gs' STO COS.\Gs SIN.\Gs R\->C ARG '\Gs' STO 'COS(u1)*COS(u2)*SIN(\Gl)/SIN(\Gs)' EVAL 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO
IF COS\178.\Ga 0 \=/
THEN 'COS.\Gs-2*SIN(u1)*SIN(u2)/COS\178.\Ga' EVAL
ELSE 0
END 'COS.2\Gsm' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO \Gl '\Gl\180' STO 'l+(1-C)*f*SIN.\Ga*(\Gs+C*SIN.\Gs*(COS.2\Gsm+C*COS.\Gs*(-1+2*SQ(COS.2\Gsm))))' \->NUM '\Gl' STO
END 'COS\178.\Ga*(a^2-b^2)/b^2' EVAL 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' \->NUM 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 'B*SIN.\Gs*(COS.2\Gsm+B/4*(COS.\Gs*(-1+2*COS.2\Gsm^2)-B/6*COS.2\Gsm*(-3+4*SIN.\Gs^2)*(-3+4*COS.2\Gsm^2)))' \->NUM '\GD\Gs' STO 'b*A*(\Gs-\GD\Gs)' \->NUM 's' STO s 1000 / "km" \->TAG 'COS(u1)*SIN(u2)-SIN(u1)*COS(u2)*COS(\Gl)' \->NUM 'COS(u2)*SIN(\Gl)' \->NUM R\->C ARG \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga1' STO "\Ga1 D.mmss\166 \|^\-> +90" \->TAG '-SIN(u1)*COS(u2)+COS(u1)*SIN(u2)*COS(\Gl)' \->NUM 'COS(u1)*SIN(\Gl)' \->NUM R\->C ARG \pi + \->NUM RAD\->D HMS\-> 360 MOD 360 MOD \->HMS DUP '\Ga2' STO "\Ga2 D.mmss\166 \|^\-> +90" \->TAG { \GD\Gs B A u\178 C COS.2\Gsm COS\178.\Ga SIN.\Ga \Gs COS.\Gs SIN.\Gs u2 u1 \Gl\180 l \Gl } PURGE
\>>
\>>


P1—>P2
\<< "4 Arguments:

lat1 lon1 s \Ga1

\[] lat1 lon1 \Ga1 \166\|^\->90
in D.mmss
\[] S < 0
W < 0
\[] dist s in m

\[] 
Change a, b, INV.f
:

run
 before all 
C.abf

" DROP STD RAD '\Ga1' STO 's' STO 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG s "s [m]" \->TAG \Ga1 "\Ga1 D.mmss\166 \|^\-> +90" \->TAG lat1 D\->RAD lon1 D\->RAD \Ga1 D\->RAD \-> lat1 lon1 \Ga1
\<< INV.f INV 'f' STO '(1-f)*TAN(lat1)' \->NUM 'TAN.u1' STO '1/\v/(1+TAN.u1^2)' \->NUM 'COS.u1' STO 'TAN.u1*COS.u1' \->NUM 'SIN.u1' STO \Ga1 COS TAN.u1 R\->C ARG '\Gs1' STO 'COS.u1*SIN(\Ga1)' \->NUM 'SIN.\Ga' STO 1 SIN.\Ga SQ - 'COS\178.\Ga' STO 'COS\178.\Ga*(a^2-b^2)/b^2' \->NUM 'u\178' STO '1+u\178/16384*(4096+u\178*(-768+u\178*(320-175*u\178)))' 'A' STO 'u\178/1024*(256+u\178*(-128+u\178*(74-47*u\178)))' \->NUM 'B' STO 's/(b*A)' \->NUM '\Gs' STO \Gs .01 - '\Gs\180' STO
WHILE \Gs \Gs\180 - ABS .000000000001 >
REPEAT 'COS(2*\Gs1+\Gs)' EVAL 'COS2.\Gsm' STO 'B*SIN(\Gs)*(COS2.\Gsm+B/4*(COS(\Gs)*(-1+2*COS2.\Gsm^2.)-B/6*COS2.\Gsm*(-3+4*SIN(\Gs)^2)*(-3+4*COS2.\Gsm^2)))' EVAL '\GD\Gs' STO \Gs '\Gs\180' STO 's/(b*A)+\GD\Gs' EVAL '\Gs' STO
END '(1-f)*\v/(SIN.\Ga^2+(SIN.u1*SIN(\Gs)-COS.u1*COS(\Gs)*COS(\Ga1))^2)' \->NUM 'SIN.u1*COS(\Gs)+COS.u1*SIN(\Gs)*COS(\Ga1)' \->NUM R\->C ARG RAD\->D 'lat2' STO lat2 "lat2 D.mmss" \->TAG 'COS.u1*COS(\Gs)-SIN.u1*SIN(\Gs)*COS(\Ga1)' \->NUM 'SIN(\Gs)*SIN(\Ga1)' \->NUM R\->C ARG '\Gl' STO 'f/16*COS\178.\Ga*(4+f*(4-3*COS\178.\Ga))' \->NUM 'C' STO '\Gl-(1-C)*f*SIN.\Ga*(\Gs+C*SIN(\Gs)*(COS2.\Gsm+C*COS(\Gs)*(-1+2*COS2.\Gsm^2)))' \->NUM 'l' STO lon1 l + RAD\->D 'lon2' STO lon2 "lon2 D.mmss" \->TAG '-SIN.u1*SIN(\Gs)+COS.u1*COS(\Gs)*COS(\Ga1)' \->NUM SIN.\Ga R\->C ARG \pi + \->NUM RAD\->D HMS\-> 360 MOD \->HMS '\Ga2' STO \Ga2 "\Ga2 D.mmss\166 \|^\-> +90" \->TAG { \Gs\180 COS2.\Gsm \GD\Gs B A u\178 COS\178.\Ga SIN.\Ga \Gs1 SIN.u1 COS.u1 TAN.u1 l C \Gl \Gs } PURGE
\>>
\>>

And variable
Version
"2.8b - 2026.02.26
campart@hotmail.com

check also:
https://geographiclib.sourceforge.io/cgi-bin/GeodSolve
by charles@karney.com"

In any case here is included the new, full directory.

Regards,
Gil

[Gil]

Version 2.9

Changed the subprogram lat—> to calculate the radius from Earth centre to the given surface latitude.

It was already correct, but, for special latitude = +/-90 (at the poles), it did not give the expected exact radius b =6356752.31425,
but 6356752.3142, due to roundings in the formulae. The same occurred with latitude = 0 (at equator line), reckoned 6378137.0001 in the program, instead of the true value of the radius a = 6378137.

Here are the codes
for lat—>:
\<< "1 single Argument:

lat
 in D.mmss

To 
change a, b, INV.f


run
 before all 
C.abf

" DROP DEG DUP 'lat' STO HMS\-> \-> la
\<<
IF la 0 \=/
THEN
IF la ABS 90 \=/
THEN '\v/(((a^2*COS(la))^2+(b^2*SIN(la))^2)/((a*COS(la))^2+(b*SIN(la))^2))' \->NUM
ELSE b
END
ELSE a
END "R Centre-Earth(" lat + ")" + \->TAG 'b/a*TAN(la)' \->NUM ATAN COS a * 2 * \pi * \->NUM "Circumf(" lat + ")" + \->TAG
\>>
\>>

As usual here is attached the full Directory with all the programs and variables.

Regards,
Gil Campart

[Gil]

New version 3 Geodesic Distance Calculator
Including now an attempt to take into account the height/altitude above sea level of two points.

For multiple point calculation of journeys, the height option is not available.

The new program is called
h—>s.h
It's short for "from h (actually 2 heights: height A in stack level 2; height B in stack level) you get the distance s.h (distance above sea level, h standing for height)".

To use it, you have before to have run once the program P1P2—> or P1—>P2.

Observe that running P1P2—> or P1—>P2
will set h1 and h2 to 0 ; however it keeps the values of s.h and delta.s.h

The calculation given here is a really rough approximation, giving wrong results.
No maths or geodesic theory behind it, contrary to the main program P1P2—> or P1—>P2, based on solid Vincenty's formulae.

The purpose is to remember that very accurate formulae like Vincenty's on the Earth surface might diverge from the reality, as most of the time the effective points are not on sea level.

If you have better ideas please feel to divide them with the community.

I adopted the following points:

0) Point A(lat 30, long 0) and B(lat 60, long 120).

First Approach
1) Find effective distance s on Earth Surface without taking into account the height.
s = 8635843.51351

2) Approximation 1 (suppose a sphere)

cos AB=
sin(pi/2-lat1)*sin(pi/2-lat2)*cos(lon1-lon2)
+cos(pi/2-lat1)*cos(pi/2-lat2).
theta = arcos(AB) =1.35256181244
R. mean × theta_in_radian = s

2b)R. mean × theta_in_radian = s
R. mean =s / theta_in_radian
=6384805.06701

3) Use 2), but, instead of R.mean, calculating new distance with "R.mean + height_min".
—> (R.mean + height_min) × theta_in_radian = s_height_min = 8636384.53824
where height_min = min(height point A, height point B).

4) Use Pythagore
requested_final_s_in_altitude =
(s_height_min² +
(height point A — height point B)²)^.5.
Final distance in altitude =8636384.54751

It seems to me that it is being a "fair" (reasonable) approximation for tackling the problem.

Second "solution"
1) R from centre to point A or B
(((a^2*COS(la))^2+(b^2*SIN(la))^2)
/((a*COS(la))^2+(b*SIN(la))^2))

Calculate then radius for A and B, R.A and R.B.

R.A 6372824.4203
R.B 6362132.22441

2) Calculate the mean radius
R.mean (R.A + R.B) / 2 = 6367478.32235

3) R.New = R.Mean + MIN (alt.A, alt B)
=6367878.32235

4) As the angle theta has not changed
S.New.horizontal
= s × R.New/R.mean.
8635843.51351 ×
6367878.32235/6367478.32235
=8636386.01046

5) Final distance s_altitude
= ((S.New.horizontal² +(alt.A—alt B)²)^.5
= 8636386.01972

Both results are very similar and seem to make sense.

Thanks for your commentaries.

Here is the code
h—>s.h
\<< "Before, run

P1P2\->D
 or 
P1\->P2
.

2 methods \-> 2 Results
as a complex#" DROP "Rough 
approximation!
" UNROT 'h2' STO 'h1' STO STD RAD \pi 2 / lat1 D\->RAD - SIN \pi 2 / lat2 D\->RAD - SIN * lon1 D\->RAD lon2 D\->RAD - COS * \pi 2 / lat1 D\->RAD - COS \pi 2 / lat2 D\->RAD - COS * + \->NUM ACOS \-> \Gm\Gh
\<< s \Gm\Gh / lat1 lat\->R DROP lat2 lat\->R DROP + 2 / 2 \->LIST \-> \GmR
\<< \GmR 1 GET h1 h2 MIN + \Gm\Gh * SQ h1 h2 - SQ + \v/ \GmR 2 GET DUP h1 h2 MIN + SWAP / s * SQ h1 h2 - SQ + \v/ R\->C DUP 's.h' STO "s(" h1 + "/" + h2 + ")" + \->TAG s "s(0,0)" \->TAG SWAP s.h s s R\->C - DUP '\GDs.h' STO "\GDs.h[m]" \->TAG
\>>
\>>
\>>

As usual, as there are some slight changes in the directory, I have included the full Directory to be copied, possibly with a shorter name, for example DIST.

Regards,
Gil Campart

[Gil]

For HP49-HP50, programs basically to calculate exact distances (error less than 0.5 mm) on the Earth surface at theorical sea level

Present Version 3.2 "corrects" my previous reasoning for approximated real distance above sea level (altitudes h1 & h2).

Explanation
From the two points P1 & P2, I calculate a first "possible" angle (mean.thets) between them, then a mean radius (=s/mean.theta), considering the ellipsoid to be a sphere.

I calculate also radii R.1 & R.2 from Earth centre to point P1 and to P2 (with program lat—>R), then I build a second mean radius.2= (R1+R2)/2.

From those two mean radii (meanR12), I have to consider a new, mean, greater radius at altitude (h1+h2)/2.

In other words,
s.h (distance between points P1 & P2 at altitude h1 & h2) =s (original distance at sea level)
× [meanR12+ (h1+h2)/2]/meanR12.

The new code for that latter program:
h—>s.h
\<< "Before, run

P1P2\->D
 or 
P1\->P2
.

2 methods \-> 2 Results
as a complex#" DROP "Rough 
approximation!
" UNROT 'h2' STO 'h1' STO STD RAD \pi 2 / lat1 D\->RAD - SIN \pi 2 / lat2 D\->RAD - SIN * lon1 D\->RAD lon2 D\->RAD - COS * \pi 2 / lat1 D\->RAD - COS \pi 2 / lat2 D\->RAD - COS * + \->NUM ACOS s SWAP / lat1 lat\->R DROP lat2 lat\->R DROP + 2 / 2 \->LIST '\GmR12' STO 1 2
FOR i \GmR12 i GET DUP h1 h2 + 2 / + SWAP / s *
NEXT R\->C DUP 's.h' STO "s(" h1 + "/" + h2 + ")" + \->TAG s "s(0,0)" \->TAG SWAP s.h s s R\->C - DUP '\GDs.h' STO "\GDs.h[m]" \->TAG
\>>

Copy the file-directory ending by .doc.

Inside are the five user programs :
A) Four very accurate programs
1) C.abh (to calibrate first your ellipsoid)
2) P1P2—> (indirect method/problem, to find the distance)
3) P1—>P2 (direct method/problem, to find point P2) )
4) lat—>
(gives
- distance from Earth centre to latitude at sea level
- & Earth circumference along that latitude)

B) Essay to give an "idea" of the "real", effective distance above sea level.
No real maths or geodesic, scentifical consideration.
Just a possible approximation to give a rough idea of the difference when calculating very "accurately", but only at Earth sea level.
5) h—>s.h.

Observations or commentaries welcome.

Regards,
Gil