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

[Gil]

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 
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 -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 
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 -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