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

[Gil]

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

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

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

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

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

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

Regards,
Gil