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

[Gil]

RE: HP49-50G Geodesic distance calculator
Version 2

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

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

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

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

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

The new codes are the following:

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

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

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

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

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

Nothing changed in

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

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

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

Remarks welcome.

Regards,
Gil