Threaded Mode | Linear Mode

Gil · 03-11-2021, 02:53 PM

New version 4.1
New features and corrections for straight line distance calculation (Euclidean distance); the latter is no more an approximation, but the exact, Euclidean solution (easily found transforming the GeoDetic latitude, the GeoDetic longitude and the ellipsoidal altitude h into Cartesian coordinates).

Note
- All the entries are mandatory in D.mmsss (for example 3d2'7.59'' should be entered 3.020759, with a dot after the 3 to separate the degrees from the minutes; instead, 30d can be entered directly, simply as 30, witout a dot after 30, but 30. or 30.0 are also valid).
- Most calculations are done in radian, but the resulting values shown and saved are all again in D.mmsss.
- According to the context, the label (attributed by ->tag) « D.mmsss » is called « o.'s » (upper « o » stands for degrees).
- All the program have a name with an arrow (->) ; though not soon visible, the same applies for the prog P1..P2->Ss (program for multiple ellipsoidal distances calculation, to be found in DATA-DIST directory, to be accessed by the first A-key).
- Ellipsoidal distance s is calculated with Vincenty's formulae with a precision of 0,1 mm.
- Nevertheless, it fails to get a distance in its iterative process for antipodal (opposite) or nearly antipodal points, as mentionend by its author. For calculations of any pair of points, with nanometre precision, see https://geographiclib.sourceforge.io/cgi.bin/GeodSolve by Charles Karney.

Page 1
Main, upper directory.
The six A, B, C, D, E and F upper keys have the following labels :
DATA.DIST af->b P1P2-> h12->s.h12 P1->P2

DATA.DIST
If you enter that directory, you have a set of points on Earth surface already saved in D.mmss format as a set of complex numbers.
Put the chosen points inside {}, writing for instance {P1 P2 P3} ENTER
or 'P1', 'P2', 'P3' 3 ->LST (LS PRG TYPE 3 ->LIST)
Then pressing the program P1..P2->Ss gives the ellipsoidal distances between P1-P2, P2-P3 and P3-P4, together with the total distance-
All that last results & inputs are saved, in that directory, under a matrix name RESULT, the first column of which describing the travel journey, the second column giving the corresponding distance.

Example
Let P1 = (10, 0), P2 = (20, -20) and P3= (30,30).
{Then after P1 P2 P3} P1..P2->Ss
you should get
[['Total_m' 7'542'625.15014]
['|P1_P2' 2'415318.01827]
['|P2_P3' 5'127'307.13187]].

Note
- If instead ot the location number P1, P2 and P3 you gave complex numbers, all of them in a single list : {10, 0), (20, -20) (30,30)} ENTER.
After pressing P1..P2->Ss,
you should get something of the kind :
[['Total_m' 7'542'625.15014]
['|10.|0.->►20.|¬20.' '415318.01827 ]
['|20.|¬20.►30.|30.' 5'127'307.13187]].
Special signs have replaced the usual one, in order to be allowed in the matrix format.
- If you prefer the results in form of a list, P1..P2->Ss suppress the command ->M.RESULT at the end of program P1..P2->Ss.
Then the stack will show :
{P1 P2 P3} or {(10, 0) (20, -20)(30,30)} on level 3
Total [m¨: 7 '542'625.15014 on level 2
{415318.0182 5'127'307.13187} on level 1.

In that directory, you can create your own points, latitude first, then longitude.
Always in the format D.mmssss.
'For instance, for point (45.5d 120.75d) :
45.30 ENTER (for latitude)
120.45 ENTER (for longitude)
For more complicated transformations of decimal degrees, use :
45.5 ->HMS
120.75 ->HMS.
And transform this pair of values into a « complex number » by R->C (LS MATH COMPLEX).
You get then (45.30 120.45).
Give it a name, for instance 'Pnew' STO.

Note
Better, but not compulsory, before performing her the multiple ellipsoid distance calculation s, to put
the var-prog P1..P2->Ss
and the variable RESULT
both at top left position under label-key A and B, respectively.
Proceed as follows :
Put P1..P2->Ss and RESULT both inside {}
- 'P1..P2->Ss' ENTER 'RESULT' 2 -> LST (LS PRG TYPE 2 ->LIST)
- or {} (LS +), then press succesively P1..P2->Ss RESULT and ENTER-key.
Having { P1..P2->Ss } on stack-level 1,
press the command ORDER (LS PRG MEM DIR NXT ORDER).

You can write
(25 30) ENTER (point A)
(35 40) ENTER (point B)
(45 50) ENTER (point C)
3 -> LIST

Then you get
{(25,30) (35,40) (45,50)}

Press then
the prog P1..P2->Ss

If you want just a meter precision, write 0 FIX and press ENTER.

Note
- If you have saved those values under A, B and C
(25,30) 'A' STO (35,40) 'B' STO (45,50) 'C' STO,
you could have done {A, B, C} P1..P2->Ss and got the same results.
- When program P1..P2->Ss ends, it goes automatically to the main, upper directory.
Suppose you have cleared meanwhile your screen and want to get back the result of your last « journey ». Then go down to DATA.DIST and recover your journey data by pressing the variable RESULT (B-key).

Back to first page, upper, main directory :

Program af->B (previously named C.abf).
The letter a stands for the equatorial radius.
The letter f is for flattning, but the given value for ellispoids is f^-1, called here inv.f.
The letter b is the polar radius, normally calculated after a and b according to (a-b)/a=f.
Press the B-key corresponding the program af->B
Then appears the Solver-Menu.
Then enter value for the equatorial radius a and press a (A-key)
Press the value of f inv.f (=f^-1) and inv.f (B-key).
(You could have entered first inv.f, then a.)
Then solve – always at the end – for your searched variable ; here, for b : LS b (i.e LS C-key)
When you have finished with the solver, end the program by pressing LS CONT :
e^2 and e'^2 will be calculated and all the parameters, with the exception of e', will appear on the stack.

Note
If you had a and b given and wanted to calcute inv.f, proceed similarly :
a-value A-key
b-value C-key
LS B-key
Then, as mentioned before, finish the program with LS CONT.

Normally, you don't have to use that program or change its values.

Example
Calculate the ellipsoidal distance between A (0, 30d30') and B (10d, 31d0'30'').
Enter (0, 30.3) and (10,31.0030) and press P1P2->
You get the required distance
s 1107.28708093 km.

Bellow that value is also automatically given Euclidian distance
at (ellipoidal) height level 0, called here s.0
(straight line/segment, even if goes through the Earth), = 1105.87858377 km.
Then forward bearing alpha1 and backward bearing alpha2.

Note
s.0 should always be < s.
But because of roundings in small distances (and limitation of Vincenty's formalae to 0.1 mm accurateness), it could happen – wrongly – that calculated s.0 > s (when effective s.0 <0).
A concrete case of that strange, « wrong calculation » is for
lat1 = 46.3404832 ( 46f34'04.832'')
long1 = 6.4201738 (6d42'01.738''
lat2 = 46.3335241 (46d33'35.241'')
long2 = 6.4152569 (6d41'52.569'')
The results calculated are s=934.358 186 768
ans s0=934.358 190 371 -> error of 1/100 mm !
s.h12 is generally < s, but, according to the requested height, it is perfectly possible that s.h12>s.

Example
Suppose now that A (0, 30d30'), according to your GPS, is 700 m above the reference ellipsoid and B (10d, 31d0'30'') 800 m.

Mandatory : having calculated previously the distance s with the program P1P2->
Write then 700 ENTER
Write 800
Then press program h12->s.h12 (Euclidian distance at levels h1=700 m for A and h2=8000 m for B)
Then you will get 3 distances :
- the previous ellipsoidal distance s = 1107.28708093.
- the previous Euclidean distance at height 0
above ellipsoid s.0 = 1105.87858377.
- the requested Euclidean distance at height 700 for A and 800 for B s.h12 = 1106.00948808.

From point A, suppose that you want to travel 50 km with initial bearing = 10.5 degrees
Enter 0 (for lat of A) ENTER
30.30 (for long of B) ENTER
50000 (in meters)
10.30 (not 10.5)
Then press [b]P1->P2[/b]
{ :lat2 D.mmss: .264060543708 :lon2 D.mmss: 30.3454674799 }
You get lat of point B = 0d26'40.60543708" and a long of 30d34'54. 674799".
In that case, the Euclidean distance is never calculated ; to avoid wrong interpretation, the previous values of s.0 and s.h12 will be set to 0 ; accordingly, h1 and h2 will also be set to 0.

Note
The height is to be carefully handled.
Check that the given height is really the GPS, geodetic height (height above the ellipsoide = segment perpendicular to the ellipsoide ; it's prolongated line does not go through the Earth centre).
The GPS, ellipsoidal height to be entered here is not the same as the orthometric height ; they differ from +/- 100 metres.

Pages 2, 3 and 4 : a summary of the entered values.

Note
If you do not enter the altitudes of A and B, both points are supposed to be at ellipsoidal height = 0.
Then h1=0 and h2=0 and s.h12= s.0.

Page 5
lat->R
With the usual latitude – GeoDetic latitude, not GeoCentric latitude – the program calculates the radius from Earth Centre to the given latitude (GPS, GeoDetic latitude) on the ellipsoidal surface.

Note
Fort the usual coordinates, i.e. GeodDetic coodinates, the calculated radius here that goes right to the Earth cenre is not the normal line from the ellipsoidal surface at GeoDetic latitude.

Page 6
laGD-> gives, from GeoDetic latitude, the GeoCentric latitude called laGC (angle to the Earth centre).
laGC-> gives, from GeoCentric latitude, the GeoDetic latitude called laGD (a somewhat greater angle than the angle to the Earth centre).

Page 7
->XYZ gives the Cartesian coordinates when having GeoDetic latitude, GeoDetic longitude and GPS height h above the ellipsoide.
XYZ-> gives the reverse : from the Cartesian coordinates, you get the GeoDetic latitude, the GeoDetic longitude and the GPS height h above the ellipsoide.

Page 8
RAD->D : converts radian to D.mmsss
D->RAD : reverse process that converts D.mmsss to radian.

Checking
Point (40,30), 40 being the GPS, geodetic latitude.
GPS, GeoDetic altitude height h = 0.
Radius from that point to Earth Centre ?
40 -lat-> gives 6369344.86324 m (result 1)
Is it correct ?
Let's calculate the Cartesian coordinates.
40 30 0 ->XYZ gives :
X = 4'237'209.07494
Y= 2'446'353.80003
Z = 4'077'985.57221
Distance of that point to centre :
[ (4'237'209.07494 - 0)^2 + (2'446'353.80003 – 0)^2 + (4'077'985.57221- 0)^2 ]^0.5
=6369344.86324 (result 2).
We see that result 1 = result 2.

Below are the codes.

af—>b
\<< "\[]Run this SOLVER-prog
only if necessary &
before P1P2\->D, 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 '(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" \->TAG f "f" \->TAG 2 f * f SQ - DUP 'e\178' STO "e\178" \->TAG 'e\178/(1-e\178)' \->NUM 'e\180\178' STO 2.01 MENU
\>>

P1P2—>
\<< "2 Complex Arguments
(lat1,lon1)(lat2,lon2)

or 4 Arguments
lat1 lon1 lat2 lon2

\[] all GeoDetic [D.'s]
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

\[]To change a, b, INV.f
run before all af\->b
" 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 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 / "s 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
\>> 0 0 h12\->s.h12 DROP 6 ROLLD 3 DROPN
\>>

h12—>sh12

\<< "\[]2 arg: ellipsoid
heights h1 h2 [m]

\[]But before, run once
P1P2\->D or P1\->P2.

h: GPS ellipsoid alt
\=/MSL (h=H+N)
H: orthomet alt (from
Point P to MSL-geoid)" DROP -105 CF -3 SF STD DEG 2 \->ARRY [ 0 0 ] SWAP 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 s 1000 / "s [km] Ellipsoid" \->TAG s.0 1000 / "s.0 [km] Eucl(0\1660)" \->TAG s.h12 1000 / "s.h12 Euc(" h1 + "\166" + h2 + ")" + \->TAG lat2 lon2 R\->C 4 ROLLD lat1 lon1 R\->C 5 ROLLD
\>>
\>>

P1—>P2

\<< "\[]2 arg: ellipsoid
heights h1 h2 [m]

\[]But before, run once
P1P2\->D or P1\->P2.

h: GPS ellipsoid alt
\=/MSL (h=H+N)
H: orthomet alt (from
Point P to MSL-geoid)" DROP -105 CF -3 SF STD DEG 2 \->ARRY [ 0 0 ] SWAP 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 s 1000 / "s [km] Ellipsoid" \->TAG s.0 1000 / "s.0 [km] Eucl(0\1660)" \->TAG s.h12 1000 / "s.h12 Euc(" h1 + "\166" + h2 + ")" + \->TAG lat2 lon2 R\->C 4 ROLLD lat1 lon1 R\->C 5 ROLLD
\>>
\>>

lat—>R
\<< "\[]2 arg: ellipsoid
heights h1 h2 [m]

\[]But before, run once
P1P2\->D or P1\->P2.

h: GPS ellipsoid alt
\=/MSL (h=H+N)
H: orthomet alt (from
Point P to MSL-geoid)" DROP -105 CF -3 SF STD DEG 2 \->ARRY [ 0 0 ] SWAP 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 s 1000 / "s [km] Ellipsoid" \->TAG s.0 1000 / "s.0 [km] Eucl(0\1660)" \->TAG s.h12 1000 / "s.h12 Euc(" h1 + "\166" + h2 + ")" + \->TAG lat2 lon2 R\->C 4 ROLLD lat1 lon1 R\->C 5 ROLLD
\>>
\>>

laGD—>
\<< "2 arg
\[] latGeoDetic \^o.'s
\[] h [above ellipsoid]

" DROP DTAG "h[m]" \->TAG SWAP DTAG "latGeoDetic \^o.\180s" \->TAG SWAP DUP2 RAD 'h' STO DUP 'laGD' STO D\->RAD \-> laGD
\<< 'a/\v/(1-e\178*SIN(laGD))' \->NUM \-> R.N
\<< 'ATAN((1-e\178*R.N/(R.N+h))*TAN(laGD))' \->NUM RAD\->D
\>>
\>> "\->latGeoCentric \^o.'s" \->TAG DUP 'laGC' STO
\>>

laGC—>
\<< "2 arg
\[] latGeoCentric \^o.'s
\[] h [above ellipsoid]

" DROP DTAG "h[m]" \->TAG SWAP DTAG "latGeoCentric \^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 "\->latGeDetic \^o's" \->TAG DUP 'laGD' STO { R.N } PURGE
\>>

XYZ—>
\<< "3 Arg:
X Y Z
" 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)
" DROP DTAG "[m]" \->TAG ROT DTAG "laGD \^o's" \->TAG ROT DTAG "loGD \^o's" \->TAG ROT 3 DUPN RAD 'h' STO 'loGD' STO 'laGD' STO laGD D\->RAD loGD D\->RAD \-> laGD loGD
\<< '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
\>>
\>>
\>>

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

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

Version
\<< "4.1 - 2021.03.11
campart@hotmail.com

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

EXPL

\<< "\[]lat/lon/h=here always
GPS GeoDetic \-> Normal
radius does not go
through Earth centre
\[]GeoCentric \-> radius
through Earth centre
\[]h=H+N H:Geoid to Pt 1
h:Ellipsoid to Pt1 1
\[]Always D.mmss (\^o\180s)" 1 DISP 7 FREEZE
\>>

P1.. PN—>Ss
\<< "Enter min. 2 Points
Names w/ {}:
{Pt1 Pt2} P1\168PN\->\GSs
w/ (lat_i, lon_i)
'Pt_i' STO &
lat_i lon_i in D.mmsss
" DROP 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 8 DROPN l.s s + 'l.s' STO
NEXT l.s DUP
IF siz 2 >
THEN \GSLIST
ELSE OBJ\-> DROP
END "Total [m]" \->TAG SWAP \->M.RESULT UPDIR
\>>
\>>

—>M.RESULT
\<< DUP SIZE 1 + DUP 2 2 \->LIST 0 CON 4 ROLL \-> l.vr l.v s M t
\<< M { 1 1 } 'Total\175m' PUT { 1 2 } t PUT 'M' STO 1 s 1 -
FOR i M i 1 + 1 2 \->LIST l.vr i GET \->STR "\|>" + l.vr i 1 + GET \->STR + "'" "" SREPL DROP "," "\166" SREPL DROP "(" "" SREPL DROP ")" "" SREPL DROP "-" "\172" SREPL DROP "\166" SWAP + OBJ\-> PUT 'M' STO M i 1 + 2 2 \->LIST l.v i GET PUT 'M' STO
NEXT M DUP 'RESULT' STO
\>>
\>>

Should you find possible discrepancies with known data or errors, do not hesitate to communicate them to me — though, I confess, I am completely alien to the GeoDetic, survey field.

Commentaries are also welcome.

Regards,
Gil Campart

Gil · 03-14-2021, 01:46 PM

RE: HP49-50G
Geodesic distance calculator
Version 4.2
New

Very slight, "cosmetical" changes.

The main one is
the set of instructions
SWAP DUP TYPE 1 ==
IF
THEN OBJ—> ROT
ELSE SWAP
END

at the beginning of
—>XYZ, whose new, full code is:
\<< "\[] 3 Arg:
laGD loGD h(ellipsoid)

or 2 Arg:
(laGD loGD) h(ellips.)

\[]To change a, b, INV.f
\-> run af\->b before!
" DROP DTAG "[m]" \->TAG SWAP DUP TYPE 1 ==
IF
THEN OBJ\-> ROT
ELSE SWAP
END ROT DTAG "laGD \^o's" \->TAG ROT DTAG "loGD \^o's" \->TAG ROT 3 DUPN RAD 'h' STO 'loGD' STO 'laGD' STO laGD D\->RAD loGD D\->RAD \-> laGD loGD
\<< '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
\>>
\>>
\>>

Instead of entering 3 arguments:
laGD
loGD
h

you can now, alternatively, enter also the equivalent 2 arguments:
(laGD, loGD)
h

and execute the —>XYZ program.

Included here, as usual, the full Directory.
Copy and save it into your calculator under the simple name DIST4.2, for instance, and enjoy it.

Regards,
Gil Campart

Gil · 03-15-2021, 02:25 PM

RE: HP49-50G
Geodesic distance calculator
Version 4.3
New

Very slight, "cosmetical" changes.

You can use your previous versions as they are.

But now I added 2 columns in the Result-Matrix for the multiple points distance calculation.

I thought nice to have, when entering in one step several points, to get also the Euclidean, straight line distances s.0 between them, as I did previously for the ellipsoidal distances s.

Meaning of Matrix Result (in DATA.DIST directory).
The 1st column describes "from where to where."
The 2nd column gives the corresponding s ellipsoidal distances.
The 3rd, new column gives the Euclidean, straight lines s.0 distances.
The 4th, last and new column checks that calculated ellipsoidal s distances are really > than Euclidean, straight line s.0 distances. Note: the formulae are correct, but because of rounding errors, it may happen that calculated s.0 is just — unduly — somewhat > calculated s (normally by less than 0.1 mm) ; in that case, both results are almost correct (as always to 0.1 mm), or no one is absolutely correct (remember that Vincenty's gives accurateness within 0.1 mm).

Go to DATA.DIST directory and try for instance
{ (45 0) (45.01 0 }
and press program P1.. P2—>Ss (S stands for capital, Greek sigma).
Then wait and edit the Result-Matrix appearing on the stack, go to last column and see the commentaries. Examine then last line/2nd column (s value) and last line/3rd column (s.0 value) ; you will see that the calculated ellipsoidal distance s 1852.19895819 is — unduly — < the Euclidean distance s.0 1852.1990012. In that case, the difference is -.04301 mm : a tiny value... that of course should be positive, and is positive in reality.

Here are the 2 full codes relative to the multiple distance calculation:

P1.. P2—>Ss
\<< "Enter min. 2 Points

\[] w/ their Names in {}
Example: {Pt1 Pt2}

Note: 'Pt_i' RCL
(lat_i, lon_i)

\[] or directly Complex
# values in {}
{ (lat_1 lon1)
(lat_2 lon2) }

\[] lat_i lon_i in [\^o's]
" DROP DUP SIZE 'siz' STO 'l0' STO { } 'ls' STO { } 'ls0' STO l0 EVAL siz \->LIST 'lc' STO 1 siz 1 -
FOR i lc i GET lc i 1 + GET UPDIR P1P2\->s s.0 s >
IF
THEN DROP
END 8 DROPN DATA.DIST 'ls' s STO+ 'ls0' s.0 STO+
NEXT
IF siz 2 >
THEN ls \GSLIST 'lst' STO ls0 \GSLIST 'ls0t' STO
ELSE ls 1 GET 'lst' STO ls0 1 GET 'ls0t' STO
END \->RESULT { siz l0 lc ls lst ls0 ls0t } PURGE UPDIR
\>>

—>RESULT
\<< 1 CF siz 2 + 4 2 \->LIST 0 CON { 1 1 } 'm' PUT { 1 2 } 's' PUT { 1 3 } 's0' PUT { 2 1 } 'Total' PUT { 2 2 } lst PUT { 2 3 } ls0t PUT { 3 1 } '\175' PUT { 3 2 } '\175' PUT { 3 3 } '\175' PUT { 3 4 } '\175' PUT { 1 4 } '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 "-" "\172" SREPL DROP "\166" SWAP + OBJ\-> PUT i 3 + 2 2 \->LIST ls i GET PUT i 3 + 3 2 \->LIST ls0 i GET PUT i 3 + 4 2 \->LIST ls i GET ls0 i GET <
IF
THEN 'Tiny.Error' 1 SF
ELSE '\175'
END PUT { 2 4 } 1 FS?
IF[/b]
THEN 'Attention!'
ELSE 'All.dist\166OK'
END PUT
NEXT DUP 'RESULT' STO
\>>

Included here, as usual, the full Directory.
Copy and save it into your calculator under the simple name DIST4.3, for instance, and enjoy it.

Regards,
Gil Campart

Gil · 03-19-2021, 10:03 AM

RE: HP49-50G
Geodesic distance calculator
Version 4.4
New

Very slight, "cosmetical" changes.

The main one is a DOERR instruction for nearly antipodal points distance calculation (with program P1P2—>s).

I thought nice to introduce such a case, in order for the program to stop after 10 unsuccessful iterations and delete by itself the intermediary, useless variables.

The full, new code for this variable-program

P1P2—>s :
\<< "\[] 2 Complex Arg:
(lat1,lon1)(lat2,lon2)

or 4 Arg:
lat1 lon1 lat2 lon2

\[] all GeoDetic [\^o.'s]
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

\[]To change a, b, INV.f
\-> run af\->b before!
" 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 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 0 'SUM' STO
WHILE \Gl \Gl\180 - ABS .000000000001 >
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 { 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 '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 / "s 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 { 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
\>> 0 0 h12\->s.h12 s.0 s >
IF
THEN DROP
END DROP 6 ROLLD 3 DROPN s.0 s >
IF
THEN "s.0:" s.0 + "
>? s:" + s +
END
\>>

Attached here, as always, the full Directory with its ordered variables and programs.

Regards,
Gil

Gil · 03-20-2021, 11:16 PM

Gil - Yesterday 11:03 AM
RE: HP49-50G
Geodesic distance calculator
Version 5
New

Allows now also distance calculation between "identical" points such as:

a) point 1(lat 90, long x) and point 2 (lat 90, long y) or
point 1(lat -90, long x) and point 2 (lat -90, long y),
—> point 1 = point 2 —> distance = 0

b) point 1(lat c, long d) and point 2 (lat c, long d)
—> point 1 = point 2 —> distance = 0

Here is the new code for
P1P2—>s
\<< "\[] 2 Complex Arg:
(lat1,lon1)(lat2,lon2)

or 4 Arg:
lat1 lon1 lat2 lon2

\[] all GeoDetic [\^o.'s]
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

\[]To change a, b, INV.f
\-> run af\->b before!
" 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 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 "Any value" { \Ga1 \Ga2 } STO s "s" \->TAG \Ga1 "\Ga1" \->TAG \Ga2 "\Ga2" \->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 >
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 { 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 '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 / "s 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 { 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
\>> 0 0 h12\->s.h12 s.0 s >
IF
THEN DROP
END DROP 6 ROLLD 3 DROPN s.0 s >
IF
THEN "s.0:" s.0 + "
>? s:" + s +
END
\>>

As always, just copy / load the full file/Directory included here into your (phone) HP49-50G calculator and enjoy it.

Regards,
Gil

Gil · 03-25-2021, 02:17 PM

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 GeoDetic [\^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 GeoDetic [\^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:
GeoDetic 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:
GeoDetic 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:
latGeoCentric [\^o.'s]
h [above ellipsoid]

\[]To change a, INV.f, b
\-> run af\->b before!
" DROP DTAG "h[m]" \->TAG SWAP DTAG "latGeoCentric \^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 "\->latGeDetic \^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 GeoDetic [\^o.'s]
\->Normal radius doesn't
go through EarthCentre
\[]laGC: latitude GC
GC=GeoCentric \->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

Jim Horn · 03-27-2021, 01:29 AM

Wow, Gil, nice work! Back in the early 1980s I needed to compute the great-circle distances involved in the CAFE250 and CAFE400 air races in California. Fortunately, the RAND report (Calculator Programs for Staff Officers) included a Vincenti algorithm program for the then-new HP-67 which I adopted to the HP-41C. I added correction for each waypoint's altitude as well, assuming a constant rate of climb or descent by scaling the distance with a ratio of mid-segment earth radius at altitude to earth radius. A secondary addition of the increase due to the slope itself was made as well. Probably not as accurate as your algorithms but it was fun to do. And it allowed us to modify the race course over the years and get consistent results.

I'll have to bookmark this and both try it on my much-used HP 50g and perhaps port it to my more recent Prime. Should be fun!

Thank you for all the good work!

Gil · 03-27-2021, 11:18 PM

Thanks for your feedback, Jim.

Regarding the HP67, it was then my dream to possess such a calculator, that I never could accomplish as a teen.
Later, however, I got an HP28S, then the 48SX and GX, before using the 50g... and losing it.
Now I loaded the great emu48 Android Application on my phone; it is an absolute flawless must : stable and fast, very easy to transfer files from and into the phone calculator. Indeed a pure marvel.

Back to my geodesic distance calculator, the main goal was:
- clarity in input arguments;
- clarity in outputs of the results.

Here are three main changes, that we might feel sorry not to have included here.

- When executing the P1P2->s, with normal font it was not possible to see at once, in one single, quick glance, or without manipulating the calculator keys, the principal result s.
Now, changing the order of the alpha1 variable, the 3 distances s, s.0 and s.h12 appear together, normally on stack levels 3, 2 and 1, respectively (or 4, 3 and 2, when s is induly calculated < s.0).

- When going into subdirectory DATA.DIST* for multiple way points calculation and executing P1..PN->Ss, it appeared only the results in form of a matrix; it was and still is quite complete, but somewhat difficult to compare horizontal cells, in particular s, s.0 ans s.h12, if the columns were or are stretched;
now, I gave to that form (format) of result the name of RES.M (M for Matrix), that appears on level 1, and added the results in form (format) of a list named RES.L (L for List) that you can easily watch by a SWAP touch (right arrow, i.e. the key upper the back-space key), as it is given on stack-level 2;

- On stack level 3, you find the list of the way points you last gave for the P1..PN->Ss execution, argument list that you can also see again pressing its (new) name LAST.LIST;

Another addition:
For some results, the textbook mode or formal (flat -79 on OFF) might hide part of the results, so that now, when necessary, it was set on ON to avoid that issue.

* Note
When entering and saving new point in the directory DATA.DIST, for example;
(45.3059, 79) 'Pt.a' STO (meaning 45d 30'59''North, 79d long East),
{-46.0020, -79} 'Pt.b' STO (meaning 46d0'20'' lat South, 79d long West),
{46, -12, 100} 'Ptc' STO (meaning 46d lat North, 12d l'on West, 100m above the ellipsoid),

a smart reflex is to think of "cleaning" your first page of that directory DATA.LIST, leaving it exactly how it was originally, executing:

{ P1..PN->SS RES.M RES.L
LAST.LIST "DOUBLE-TRIANGLE LOOKING RIGHT" "TRIPLE-TRIANGLE LOOKING RIGHT" }
ORDER.
Then all your variables will appear starting from page 2, leaving the page 1 only for calculation and results.

Enjoy
and regards,
Gil

Gil · 03-27-2021, 11:21 PM

Thanks for your feedback, Jim.

Regarding the HP67, it was then my dream to possess such a calculator, that I never could accomplish as a teen.
Later, however, I got an HP28S, then the 48SX and GX, before using the 50g... and losing it.
Now I loaded the great emu48 Android Application on my phone; it is an absolute flawless must : stable and fast, very easy to transfer files from and into the phone calculator. Indeed a pure marvel.

Back to my geodesic distance calculator, the main goal was:
- clarity in input arguments;
- clarity in outputs of the results.

Here are three main changes, that we might feel sorry not to have included here.

- When executing the P1P2->s, with normal font it was not possible to see at once, in one single, quick glance, or without manipulating the calculator keys, the principal result s.
Now, changing the order of the alpha1 variable, the 3 distances s, s.0 and s.h12 appear together, normally on stack levels 3, 2 and 1, respectively (or 4, 3 and 2, when s is induly calculated < s.0).

- When going into subdirectory DATA.DIST* for multiple way points calculation and executing P1..PN->Ss, it appeared only the results in form of a matrix; it was and still is quite complete, but somewhat difficult to compare horizontal cells, in particular s, s.0 ans s.h12, if the columns were or are stretched;
now, I gave to that form (format) of result the name of RES.M (M for Matrix), that appears on level 1, and added the results in form (format) of a list named RES.L (L for List) that you can easily watch by a SWAP touch (right arrow, i.e. the key upper the back-space key), as it is given on stack-level 2;

- On stack level 3, you find the list of the way points you last gave for the P1..PN->Ss execution, argument list that you can also see again pressing its (new) name LAST.LIST;

Another addition:
For some results, the textbook mode or formal (flat -79 on OFF) might hide part of the results, so that now, when necessary, it was set on ON to avoid that issue.

* Note
When entering and saving new point in the directory DATA.DIST, for example;
(45.3059, 79) 'Pt.a' STO (meaning 45d 30'59''North, 79d long East),
{-46.0020, -79} 'Pt.b' STO (meaning 46d0'20'' lat South, 79d long West),
{46, -12, 100} 'Ptc' STO (meaning 46d lat North, 12d l'on West, 100m above the ellipsoid),

a smart reflex is to think of "cleaning" your first page of that directory DATA.LIST, leaving it exactly how it was originally, executing:

{ P1..PN->SS RES.M RES.L
LAST.LIST "DOUBLE-TRIANGLE LOOKING RIGHT" "TRIPLE-TRIANGLE LOOKING RIGHT" }
ORDER.
Then all your variables will appear starting from page 2, leaving the page 1 only for calculation and results.

Enjoy
and regards,
Gil

Gil · 03-30-2021, 10:06 PM

New version 6.3
Previous versions are OK in principle.

Here are nevertheless some minor changes in this new version regarding, for most of them, the case of North Pole to South Pole or reverse.

I was not convinced to adopt here the conventions of other programs that consider the values of the longitudes when computing the bearings.
For me, when you are on the Pole and want to go to the other Pole, you can chose any meridian line, so that I decided - for the time being - to set, in that special case, alpha 1 and alpha 2 equal to 'free'.
Indeed, being on the Pole, what is the longitude on that special position? Clearly, for me, any value should be the appropriate answer - though, perhaps, depending on the situation, not a genuine geodesic one.

Take, for instance, (90, 0) or (90,25) on the stack (level 3).
Give the shortest geodesic (ellipsoidal) distance to the South pole on the stack (level 2), i.e. 20'003'391.4586 (with the HP50G calculator and WGS84).
Then chose any value for the bearing alpha 1 (in level 1) and press P1..PN->Ss.
The answer should always give the South Pole - and will of course always do so in that (calculator) program, but - once again for the time being - with a free longitude, whatever the introduced longitude in level 1 for the North Pole.

Note that's my special view as "amateur" or understanding on the matter.
Netherless, I still did consider the longitude in the North (or South) Pole case when the other point is not the other Pole.
Not convinced however that this traditional or geodesic approach is the best one for the "common" man, as this case is somewhat similar to the other above first Pole to Pole case.

Here are again the code
[/b]

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 GeoDetic [\^o.'s]
\[] S < 0
W < 0

Vincenty fails for
nearly antipode pts

\[]To change a, INV.f, b
\-> run af\->b before!
" DROP -79 SF STD \-> p
\<<
CASE p TYPE 5 ==
THEN p OBJ\-> 2 == 0 IFT
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 == 0 IFT
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 DUP ABS 'la' STO D\->RAD lon1 D\->RAD lat2 D\->RAD lon2 D\->RAD \-> lat1 lon1 lat2 lon2
\<< lat1 lat2 == la 90 == AND lat1 lat2 == lon1 lon2 == AND OR
IF
THEN 'la' PURGE 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 4 ROLLD { 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 la } 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 lat1 lat2 + 0 == la 90 == AND
IF
THEN 'Any.value' { \Ga1 \Ga2 } STO \Ga1 "\Ga1\1662" \->TAG
ELSE '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 360 MOD '\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 360 MOD '\Ga2' STO "\Ga2 D.mmss\166 \|^\-> +90" \->TAG
END
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 la } PURGE
END
\>> s '?' SAME NOT
IF
THEN s 1000 / "s [km]" \->TAG 1 2
START 4 \Ga1 'Any.value' SAME 0 1 IFTE + ROLL
NEXT s s.0 <
IF
THEN 's<?s.0'
END
END 3 CF
\>>

P1..PN->Ss
\<< "
\[]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 GeoDetic [\^o.'s]
\[] S < 0
W < 0
\[] dist s in m

\[]To change a, INV.f, b
\-> run af\->b before!
" DROP -79 SF 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 DUP 'la1' 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 'la2' STO la1 ABS 90 == la2 ABS 89.5959999999 == la2 ABS 90 == OR AND
IF
THEN 'Free' { \Ga2 lon2 } STO la1 NEG s 0 == -1 1 IFTE * 'lat2' STO lat2 "lat2" \->TAG lon2 "lon2" \->TAG \Ga2 "\Ga2" \->TAG
ELSE la2 '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 360 MOD \->HMS
IF la1 ABS 90 ==
THEN DROP la1 SIGN -90 * 90 +
END '\Ga2' STO \Ga2 "\Ga2 D.mmss\166 \|^\-> +90" \->TAG
END { \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 la1 la2 } PURGE
\>>
\>>

Regstds

Gil · 04-02-2021, 10:10 AM

Version 6.4

The latitudes should always be in the interval [-90 +90].

Now, if is not the case (lat < -90 or lat > 90), a DOERR instruction was introduced.
Note also that, in that case (lat < -90 or lat > 90), when using multiple points distances calculation P1.. PN—>SS, the whole process will — unfortunately — be interrupted and the intermediary variables will consequently be deleted (in order to leave an original, "clean" DATA.DIST directory).

When a faulty error appear, according to the program, the introduced inputs are left in their initial order — or not, when less than 3 values appear in the stack.

In that latter case, the wrong latitude will appear on stack level 1. Depending on the program, you then have to do either:
- DROP, write new value in place of faulty latitude and run again the program;
- DROP, write new value in place of faulty latitude, SWAP and run again the program;
- DROP, write new value in place of faulty latitude, ROT, ROT and run again the program.

For the time being, I left the Poles treatment as already mentioned in the previous versions; possibly, however, I will have to adapt it back to the usual, geodesic approach.

Regards,
Gil

Gil · 04-04-2021, 10:01 PM

Suppose you are not on the Poles and want to go to North Pole, the forward azimuth should be exactly 0 degree, and not a very tiny value ≠ 0 ; and exactly 180 degrees to go to South Pole.

Suppose now your are on North Pole on go somewhere not to the South Pole, the backward azimuth should be exactly 0 degree, and not a very tiny value ≠ 0 ; and 180 degrees when starting from South Pole.

Small changes in the programs have been made in order to comply fully with the above.

Regards,
Gil

Gil · 04-07-2021, 06:39 PM

Nearly antipodal points

For distance and forward / backward calculation
P1P2—>s, or inverse problem, of nearly antipodal points, the basic idea would be to use the direct solution P1—>P2, in a modified version, of course.

We give then the initial point 1 (lat1, lon1) and final antipodal point 2 (lat 2, lon 2) and try to reach that last point, via a modified P1—>P2 program, just making vary at the same time the distance s and alpha 1.

That should give something of the kind:

'f1mod(P1—>P2(lat1, long1, s, alpha1))=lat2'
'f2mod(P1—>P2(lat1, long1, s, alpha1))=lon2'
2 —>ARRY
's'
'Alpha1'
2 —>ARRY
20000000
Pi/2 —>NUM
2 —>ARRY
MSLV

Still to be done...

Gil · 05-14-2021, 10:51 PM

New version 6.6

Like Version 6.4 (or 6.5?) , but corrected one
result appearance
(the unit d's was not clearly seen).

Gil · 05-14-2021, 11:00 PM

New version 6.6

Like Version 6.4 (or 6.5?) , but corrected one
result appearance
(the letter s of seconds in the unit d's).

Regards,
Gil

MS_VieuxGaminHP · 08-18-2022, 03:43 PM

(06-23-2020 08:41 PM)Gil Wrote: 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

I'm very interested in the work you've done (I already was "trying" to choose to do it my self on my DM42, 50G or 15C and at first trying to re-learn Spheric trig)
I would very much love to simply use it on my 50g, but i don't get why you do not provide either the code for us to enter with our own fingers or the file necessary to put on a SD card, as it should be.
Don't take it bad, i'm just starting to re-learn my HP's and certainly lakc the necessary knowledge to dive into it.

I will try what you gave as the code. At least (thank you) it will force my old brain to remember not only the maths, but how to program my 50g.

Have a wonderful whatever Wink

Pierre (from Belgium)

Gil · 08-20-2022, 09:05 PM

If I am not mistaken, I gave all the codes in the several posts regarding this program.

Please inform me if it is not so and I will forward the missing subprograms.

The easiest should be to load the whole program and read then the code directly from your calculator or the emu48 application on a phone.

As mentioned previously, you can give a full journey with several way points (intermediate positions) and all will be calculated in one single step, with of course the intermediate distances.

Regards

Gil · 01-12-2023, 01:18 PM

Version 6b
with the ending .hp
(instead of the confusing .DOC ending).

Gil · 03-29-2023, 10:54 PM

Correction in program P1 —>P2.

Regards,
Gil

Gil · (This post was last modified: 01-07-2024 11:57 AM by Gil.)

I was dealing with coordinates (lat_dd.mmss, lon_dd.mmss) introduced like complex numbers (a, b).

But I was somewhat disturbed, by pressing ENTER, to see the result appear with (r, angle), with r=sqrt(a²+b²).

All normal, of course, when flag -16 is set.

For clarity, I just added -16CF whenever needed.

Sorry, I adapted an old version.

The last version to be used is 6.09 (or 6d).
My apologises for the inconvenience.