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

[Gil]

Just a few words about the history structure change of the program.

Originally the program was intended to calculate only geodesic distances, i.e. distances on the theorical ellipsoid representing the Earth.

Soon it appeared that a large number of short distances to be checked were straight lines segments, and not geodesic "lines" (curves), so that the necessity to include Euclidean distances; initially - and wrongly - all on heights above the ellipsoid were supposed to be equal to zero metre, so that the entered points were only complex numbers (lat_i, long_i).

But seldom are the heights at 0 level on the ellipsoid (or on Mean Sea Level), so that the following, new notation for the points was adopted :
{lat_i long_i height_i}, besides already the first one with complex number, without the height (lat_i long_i height_i).

Now both ways (complex numbers () and lists {}) for calculating the distances are valid, even the list with only 2 objets (the lat and long, without including the height) like {lat_i long_i height_i}; in that latter case, and every time when the height does not appear clearly, the height above the ellipsoid will be set to 0.

Besides, everything is calculated automatically; no need anymore to enter in a second step (after P1P2—>) the heights and launch the program h12->s.h12, as that latter program is now integrated inside the main program P1P2->.

When calculating Euclidean distances between two registered points in the form A (30deg, 60deg), without the mention of its height, which is missing, and B {40deg, 60deg, 100m}, the height of A is, for practical use and calculation, set equal to 0m, as already alluded. Should it nevertheless not be the case for A, say it is 800 m above the ellipsoid, then you have to change the value of A:
{30 60 800} 'A' STO.

When adding new points, press first the DATA.DIST A-key that gives access to the file directory. Then enter your new points and save them with a specif name for each of them. Then order the following variables:
{ P1¨PNSs RESULT l0 —> —>—> —>—>—>
ls st ls0 s0t lsh12 sh12t }
LS MEM
DIR THE NEXT
ORDER.

For the meaning of the above variables
ls st ls0 s0t lsh12 sh12t, see further below (quickly: l=list, t=total, s=ell.distance, s.0=Euclidean distance at height 0, sh12=at height h1&h2).

Sometimes, it is easier to enter directly the values without () or {}, for example :
30 60 800 40 60 100 P1P2->
(always in the order lat1 long1 h1 lat2 long2 h3).

One problem was the case of (nearly) antipodal pair of points.
As soon the iterations number reaches 10, antipodal points are now detected and the program breaks down (FLAG 3 cleared) smoothly, deleting the intermediary, unnecessary variables, with a warning message, depending on the case.

No warning message in two cases:
- when multiple distance calculating with
{ () () ()...} P1.. P2—>Ss
(because, for
{A B P1 P2 (40, 60) } P1..PN->Ss,
setting FLAG 3 at the beginning of each loop of the program P1..PN->Ss was/is a trick to permit to carry on the whole process — (A—B, B—P1, P1—P2, P2—(40,60) distances — for all the points in the given list {}, even when finding antipodal points, avoiding the warning message and the DOERR execution) ;
- when having set once "incidentally" the flag 3 and running for the first time the program P1P2—>; note that, after each execution of the program P1P2—> (or also for the other program P1..PN->SS), the flag 3 will be cleared, so that the subsequent runs of P1P2—> will duly produce a warning message when in presence of nearly antipodal pairs of points.

No more any restrictions for the points and for the following possible formats given below in examples:

(60,0) (60,0) P1P2-> or
{60,0, 0} {60,0,0} P1P2-> or
{60,0} {60,0,0} P1P2-> or
(60,0) {60,0,0} P1P2-> or
60 0 0 60 0 0 P1P2->
6 arguments when direct values in that latter case — and not only 4 or 5 when heights are unknown/missing (set then equal to zero to complete 6 arguments).

Though ridiculous, you should get the expected 0 result.

The same applies for

{ (60,0) (65,10) } P1..P2->Ss or
{ {60,0, 0} {65,10,0} }P1..P2->Ss or
{ {60,0} {65,10,0} } P1..2->Ss or
{ (60,0) {60,10,0} } P1..P2->Ss or
{60 0 0 60 10 0} P1..P2->Ss
In the latter case, a LIST {} with 6 argument inside — and not only 4 or 5 values inside the {} when heights are unknown (mandatory to write down here the missing heights h1 and h2, for instance equal to zero, in order to get 6 values/arguments in the single {} before launching P1P2—>).

In the case of multiple distances calculations, it's easier to enter points names, for instance:

{ A B C P1 }, then execute P1..P2—>Ss,
where A is:
- (latA, longA) or
- {latA longA hA} or
- {latA longA}.

Besides general rules, here below some signs and variables / calculation explanations popping up in the RESULT Matrix (after executing P1..P2—>Ss) or in the DATA.DIST directory.

l0= list zero, list original {}, as entered just before P1..P2->SS; remember that for multiple distance calculation P1..P2->SS the points are to be included in a list {} (in the above given example, its content/value is { A B C P1 } for l0)

ls= list of the different ellipsoidal distances s calculated with P1..P1->Ss;
st= total sum of the different ellipsoidal distances s calculated with P1..P2->Ss;

ls0= list of the different Euclidean distances s.0 calculated with P1..P1->Ss;
s0t= total sum of the different Euclidean distances s.0 calculated with P1..P1->Ss;

lsh12= list of the different Euclidean distances calculated with P1..P1->Ss;
sh12t= total sum of the different Euclidean distances calculated with P1..P1->Ss;

> Means that the given total sum of the several ellipsoid distances s was calculated setting all the antipodal point distances equal to zero ; consequently, the effective total real ellispoid distance s is greater (>) than the given value;

~ Has two meanings, the first one being to fill/represent the blanks in the point description of the Matrix result (first column) in the directory Data.DIST and the second signification being to remember that, due to roundings, the calculated ellipsoid distance s is, for the mentioned specific pair or points, unduly smaller than the Euclidean distance s.0 (the .0 in s.0 means at height zero, in opposition to s.h12, for the calculation at heights h1 and h2);

"Bold arrow pointing to the right": separates two points A & B when calculating the distance between them (a kind of A—B);

"Curved minus": the minus sign, as the real minus sign is not accepted inside the following format 'name';

"Two vertical lines, one under the other: replaces the comma when having complex numbers, as the real comma sign is not accepted inside the following format 'name';

"Tiny error" in second column of the Matrix RESULT means that the ellipsoid calculated distance s is "unduly" less than Euler distance s.0; the case may appear for points not too far away from each other (a not clear-cut notion), for instance s between (45, 0) and (45d01'00", 0) = 1852.19895819,
< s.0 = 1852.1990012.

Note also that, for simplification, generally because of lack of place, D.mmsss is sometimes written o.'s.

Only have to be used the latter format D.mmsss, and never decimal degrees.
For instance, 7.5d should be entered as 7.30 (meaning 7 degrees and 30').
Moreover, a value like 7.304589 means everywhere, even in the saved variables, 7d30'45.89".

Only use GeoDetic latitudes (laGD), which are the latitudes usually given for the place locations (Google, maps and books , GPS).
Should you need to convert them to Geocentric latitudes, use laGD—>; for the reverse problem, use instead laGC—>.

All the programs have an arrow (—>) in their name.

Now as usual here are the full codes of those programs.

af—>b
\<< "\[]Run this 
SOLVER
-prog

only if necessary
 &

before
 P1P2\->s, P1\->P2
or other prog.

\[] Finish it executing

CONT
 (LS 'ON-KEY').

For WGS84
a:6378137
INV.f:298.257223563

For IERS 2003
a:6378136.6
INV.f:298.25642
" DROP 3 CF '(a-INV.f/INV.f*b)/a-1/INV.f' STEQ 30 MENU HALT DEPTH 1 MIN DROPN INV.f INV 'f' STO a "a" \->TAG b "b" \->TAG INV.f "INV.f 
&
 f" \->TAG f 2 f * f SQ - DUP 'e\178' STO "e\178 
&
 e\180\178" \->TAG 'e\178/(1-e\178)' \->NUM DUP 'e\180\178' STO 2.01 MENU
\>>

P1P2—>s
\<< "\[] With ellipsoid
heights \=/ 0

2 Arg, each in 
{}
:

{lat1,lon1,h1}
{lat2,lon2,h2}


or 6 Arg:

lat1 lon1 h1
lat2 lon2 h2


\[] Or without heights
\-> heights set = 0

2 Arg, each in 
{}
:

{lat1,lon1}
{lat2,lon2}


or 2 Arg, each in 
()
:

(lat1,lon1)
(lat2,lon2)


\[] all 
G
eo
D
etic [\^o.'s]
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

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

\-> run 
af\->b before
!
" DROP STD \-> p
\<<
CASE p TYPE 5 ==
THEN p OBJ\-> 2 ==
IF
THEN 0
END
END p TYPE 1 ==
THEN p OBJ\-> 0
END p
END
\>> 'h2' STO 'lon2' STO 'lat2' STO \-> p
\<<
CASE p TYPE 5 ==
THEN p OBJ\-> 2 ==
IF
THEN 0
END
END p TYPE 1 ==
THEN p OBJ\-> 0
END p
END
\>> 'h1' STO 'lon1' STO 'lat1' STO lat1 "lat1 D.mmss" \->TAG lon1 "lon1 D.mmss" \->TAG h1 "h1" \->TAG lat2 "lat2 D.mmss" \->TAG lon2 "lon2 D.mmss" \->TAG h2 "h2" \->TAG h12\->s.h12 s.0 1000 / "s.0 [km] Eucl(0\1660)" \->TAG s.h12 1000 / "s.h12 Euc(" h1 + "\166" + h2 + ")" + \->TAG RAD lat1 D\->RAD lon1 D\->RAD lat2 D\->RAD lon2 D\->RAD \-> lat1 lon1 lat2 lon2
\<< lat1 lat2 == lat1 ABS 90 D\->RAD == AND lat1 lat2 == lon1 lon2 == AND OR
IF
THEN 0 's' STO h2 h1 - ABS 's.h12' STO 'Any.value' { \Ga1 \Ga2 } STO DROP s.h12 1000 / "s.h12 Euc(" h1 + "\166" + h2 + ")" + \->TAG \Ga1 "\Ga1\1662" \->TAG
ELSE 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 0 'SUM' STO
WHILE \Gl \Gl\180 - ABS .000000000001 > SUM 10 < AND
REPEAT 1 'SUM' STO+ '\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 SUM 10 ==
IF
THEN '?' { s \Ga1 \Ga2 } STO s "s" \->TAG UNROT '\Ga1?\166\Ga2?' 3 FC?
IF
THEN { SUM \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 "2 nearly antipodal points:
as expected,
Vincenty's algor.failed!" DOERR
END
END
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 SUM 10 ==
IF
THEN DROP
ELSE 's' STO '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
END { SUM \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
END
\>> s '?' SAME NOT
IF
THEN s 1000 / "s [km]" \->TAG 4 s 0 == 0 1 IFTE + ROLLD s s.0 <
IF
THEN "s.0:" s.0 + "
>? s:" + s +
END
END 3 CF
\>>

P1—>P2

\<< "
\[]4 Arg:
lat1 lon1 s \Ga1


or 3: 
(lat1 lon1) s \Ga1


or 
 {lat1 lon1} s \Ga1

or
 {lat1 lon1
 
h1} s \Ga1


\[]
Note
: heights are 
not

considered in program
\->
Pt1
 & 
2
: supposed on
ellips surface 
at 0 
m

\[]lat1 lon1 \Ga1 \166\|^\->90
all 
G
eo
D
etic [\^o.'s]
\[] S < 0
W < 0
\[] dist s in m

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

\-> run 
af\->b before
!
 

" DROP STD 0 { s.0 s.h12 } STO RAD '\Ga1' STO 's' STO \-> p1
\<<
CASE p1 TYPE 1 ==
THEN p1 OBJ\->
END p1 TYPE 5 ==
THEN p1 OBJ\-> 2 - DROPN
END p1
END
\>> 'lon1' STO 'lat1' STO "Pt1 & 2:at 0 m on" "ellipsoid surface" 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
\<< '(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
\>>
\>>

lat—>R
\<< "\[] 1 single Arg:

G
eo
D
etic
 lat
 [\^o.'s]

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

\-> run 
af\->b before
!
" DROP DEG "From laGD [\^o.'s]" \->TAG DUPDUP '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 "to Earth-centre, R" \->TAG DUP 'R.\GW\175latGD' STO
IF la ABS 90 ==
THEN 0
ELSE 'b/a*TAN(la)' \->NUM ATAN COS a * 2 * \pi * \->NUM
END "Circumf(" lat + ")" + \->TAG DUP 'Circ' STO
\>>
\>>

laGD—>
\<< "\[] 1 single Arg:

G
eo
D
etic
 lat
 [\^o.'s]

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

\-> run 
af\->b before
!
" DROP DEG "From laGD [\^o.'s]" \->TAG DUPDUP '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 "to Earth-centre, R" \->TAG DUP 'R.\GW\175latGD' STO
IF la ABS 90 ==
THEN 0
ELSE 'b/a*TAN(la)' \->NUM ATAN COS a * 2 * \pi * \->NUM
END "Circumf(" lat + ")" + \->TAG DUP 'Circ' STO
\>>
\>>

laGC—>
\<< "\[]2 Arg:

la
t
G
eo
C
entric [\^o.'s]

h
 [above ellipsoid]

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

\-> run 
af\->b before
!
" DROP DTAG "h[m]" \->TAG SWAP DTAG "
la
t
G
eo
C
entric \^o.'s" \->TAG SWAP DUP2 RAD 'h' STO D\->RAD 'laGC' STO
\<< 'a/\v/(1-e\178*SIN(laGD))' \->NUM 'R.N' STO 'TAN(laGD)=TAN(laGC)/(1-e\178*R.N/(R.N+h))' EVAL
\>> 'laGD' laGC ROOT RAD\->D "\->
la
t
G
e
D
etic \^o.'s" \->TAG DUP 'laGD' STO { R.N } PURGE
\>>

XYZ—>
\<< "3 Arg:Cartesian 
X Y Z


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

\-> run 
af\->b before
!
" DROP DTAG "Z" \->TAG ROT DTAG "X" \->TAG ROT DTAG "Y" \->TAG ROT 3 DUPN 'Z' STO 'Y' STO 'X' STO RAD 'TAN(LA)=(Z+a/\v/(1-e\178*SIN(LA)^2)*e\178*SIN(LA))/\v/(X^2+Y^2)' 'LA' 'ATAN(Z*(1+e\180\178)/\v/(X^2+Y^2))' \->NUM ROOT 'laGD' STO X Y R\->C ARG 'loGD' STO X SQ Y SQ + \v/ laGD COS / 'a/\v/(1-e\178*SIN(laGD)^2)' \->NUM - "\->h [m]" \->TAG laGD RAD\->D DUP 'laGD' STO "\->laGD \^o.'s" \->TAG loGD RAD\->D DUP 'loGD' STO "\->loGD \^o.'s" \->TAG ROT 'LA' PURGE
\>>

—>XYZ
\<< "\[] 3 Arg:

laGD loGD h
(ellipsoid)

or 2 Arg, in 
()
&
h
:

(laGD loGD) h
(ellips.)

or 1 Arg, in 
{}
:

{laGD loGD h
(ellips)
}


\[]
laGD 
&
 loGD
 in [\^o.'s]

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

\-> run 
af\->b before
!
" DROP \-> p
\<< p TYPE 5 ==
IF
THEN p OBJ\-> DROP
ELSE \-> q
\<< q TYPE 1 ==
IF
THEN q OBJ\-> p
ELSE q p
END
\>>
END
\>> 'h' STO 'loGD' STO 'laGD' STO laGD "laGD \^o's" \->TAG loGD "loGD \^o's" \->TAG h "h [m]" \->TAG RAD 3 DUPN DROP D\->RAD SWAP D\->RAD \-> loGD laGD
\<< 'a/\v/(1-e\178*SIN(laGD)^2)' \->NUM \-> N
\<< '(N+h)*COS(laGD)*COS(loGD)' \->NUM DUP 'X' STO "\->X" \->TAG '(N+h)*COS(laGD)*SIN(loGD)' \->NUM DUP 'Y' STO "\->Y" \->TAG '(N*(1-e\178)+h)*SIN(laGD)' \->NUM DUP 'Z' STO "\->Z" \->TAG
\>>
\>>
\>>

h12—>s.h12
\<< -105 CF -3 SF STD DEG [ 0 0 ] h1 h2 2 \->ARRY 2 \->LIST \-> l.h
\<< 1 2
FOR j l.h j GET OBJ\-> DROP 'h2' STO 'h1' STO 1 2
FOR i a 1 e\178 "lat" i + OBJ\-> HMS\-> SIN SQ * - \v/ / \-> N
\<< N "h" i + OBJ\-> + "lat" i + OBJ\-> HMS\-> COS * "lon" i + OBJ\-> HMS\-> COS * N "h" i + OBJ\-> + "lat" i + OBJ\-> HMS\-> COS * "lon" i + OBJ\-> HMS\-> SIN * 'b^2/a^2' \->NUM N * "h" i + OBJ\-> + "lat" i + OBJ\-> HMS\-> SIN * 3 \->ARRY
\>>
NEXT - DUP DOT \v/
NEXT
\>> 's.h12' STO 's.0' STO
\>>

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

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

P1.. PN—>Ss
\<< "Arg: list of
min. 2 Points in 
{}


\[] with their Names
Example: 
{
Pt1 Pt2
}


Format of Pt_i is:
. 
(
lat_i, lon_i
)
 or
. 
{
lat_i, lon_i
}
 or
. 
{
lat_i, lon_i, 
hi}


\[] or with complex #

{
 
(
lat_1 lon1
)


(
lat_2 lon2
)
 
}


\[] or { {}{} }:

{
 
{
lat_1 lon1
}


{
lat_2 lon2
}
 
}


\[] or { {}{} }:

{
 
{
lat_1 lon1 
h1}


{
lat_2 lon2 
h2}
 
}


\[] when 
no heights

\-> heights set = 
0 
m

\[]
lat_i lon_i 
in [\^o.'s]


" DROP 1 CF 2 CF DUPDUP SIZE 'siz' STO 'l0' STO { } 'ls' STO { } 'ls0' STO { } 'lsh12' STO 0 { st s0t sh12t } STO l0 EVAL siz \->LIST 'lc' STO 1 siz 1 -
FOR i lc i GET lc i 1 + GET UPDIR 3 SF P1P2\->s DUP TYPE 2 == { DROP } IFT 10 \Ga1 '?' SAME \Ga1 'Any.value' SAME OR 0 1 IFTE + DROPN DATA.DIST 'ls' s STO+ 'ls0' s.0 STO+ 'lsh12' s.h12 STO+ 'st' s '?' SAME
IF
THEN 0 1 SF
ELSE s
END STO+ 's0t' s.0 STO+ 'sh12t' s.h12 STO+ s '?' SAME NOT
IF
THEN s s.0 <
IF
THEN 2 SF
END
END
NEXT s0t sh12t SAME 2 FS? AND
IF
THEN "~" sh12t + OBJ\-> 'sh12t' STO
END 2 FS?
IF
THEN "~" s0t + OBJ\-> 's0t' STO
END 1 FS?
IF
THEN "'\166>"
ELSE 2 FS?
IF
THEN "~"
ELSE ""
END
END st + OBJ\-> 'st' STO \->RESULT { siz lc } PURGE UPDIR 3 CF
\>>

—>RESULT
\<< 2 CF siz 2 + 5 2 \->LIST 0 CON { 1 3 } 's' PUT { 1 4 } 's0' PUT { 1 5 } 's.h12' PUT { 2 1 } 'Total\175m' PUT { 2 3 } st PUT { 2 4 } s0t PUT { 2 5 } sh12t PUT 1 5
FOR i { 3 i } '\175' PUT
NEXT { 1 2 } 'TEST\166s>?s.0' PUT 1 siz 1 -
FOR i i 3 + 1 2 \->LIST l0 i GET \->STR "\|>" + l0 i 1 + GET \->STR + "'" "" SREPL DROP "," "\166" SREPL DROP "(" "" SREPL DROP ")" "" SREPL DROP "{ " "" SREPL DROP " }" "" SREPL DROP "-" "\172" SREPL DROP " " "~" SREPL DROP "\166" SWAP + OBJ\-> PUT i 3 + 3 2 \->LIST ls i GET PUT i 3 + 4 2 \->LIST ls0 i GET PUT i 3 + 5 2 \->LIST lsh12 i GET PUT i 3 + 2 2 \->LIST ls i GET \-> s
\<<
CASE s '?' SAME
THEN 1 SF 'Antip.Pt'
END s 0 ==
THEN 'Same.Pt'
END s ls0 i GET <
THEN 2 SF 'Tiny.Error'
END '\175'
END
\>> PUT
NEXT { 2 2 } 1 FS? 2 FS? OR
IF
THEN 'Attention!'
ELSE 'All.dist\166OK'
END DUP 4 ROLLD PUT { 1 1 } ROT PUT DUP 'RESULT' STO ls "ls" \->TAG st "st" \->TAG ls0 "ls0" \->TAG s0t "s0t" \->TAG lsh12 "lsh12" \->TAG sh12t "sh12t" \->TAG 7 ROLL
\>>


And the program-variables:

Note1
\<< "lat[i]:latit always 
GD


GD
=GPS 
G
eo
D
etic [\^o.'s]
\->Normal radius doesn't
go through EarthCentre
\[]la
GC
: latitude 
GC


GC=G
eo
C
entric \->radius
through Earth centre
\[]Only 
h
[m]=H+N:Ellip\->P
\=/H(orthom):Geoid \-> PtP" 1 DISP 7 FREEZE
\>>

Note2
\<< "\[]
s
: geodesic distance

ON
 ellipsoid \-> 
h=0

\[]
s.0
: Euclidean dist

NOT
 on ellipsoid, but

through
/inside Earth,
P1 & P2 at height 
h=0

\[]
s.h12
: Euclidean dist
at heights 
h1
 & 
h2


above
 the ellipsoid
" 1 DISP 7 FREEZE
\>>

Version
\<< "6 - 2021.03.25
Last published:
5 - 2021.03.20
campart@hotmail.com

check also:
https://geographiclib.
sourceforge.io/cgi-bin
/GeodSolve
by charles@karney.com
" 1 DISP 7 FREEZE
\>>

Just to convert all those text codes into HP49-50G compatible codes, use the following program:

INOUT to be saved in your calculator
\<< RCWS RCLF \-> ws f
\<< 3 TRANSIO DUP TYPE 2
IF ==
THEN \->STR f SIZE 3 > # 2F34Dh # 3016Bh IFTE SYSEVAL + STR\->
ELSE STD 64 STWS \->STR f SIZE 3 > # 2F34Fh # 2FEC9h IFTE SYSEVAL
END ws STWS f STOF
\>>
\>>

Regards,
Gil Campart