|
HP49-50G Geodesic distance calculator +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: HP49-50G Geodesic distance calculator (/thread-15243.html) |
49 50 Ver7.hp Geodesic distance & Earth Euclidean distance calculator, bearing - Gil - 06-22-2020 03:10 PM For HP49-HP50, programs to calculate exact distances (error less than 0.5 mm) on the earth surface. Last version; 6.08, January 7th 2024. Copy the file-directory ending by .doc or by .hp. Inside are two programs: 1) P1P2—> (indirect method/problem) 2) P1—>P2 (direct method/problem) Both programs give very accurate answers, for example the calculated distance is exact to 0.5 mm. Let's assume you have two points P1 and P2 on the earth geoid. 1) To get the shortest distance between them, enter the following known 4 "objects" in the stack: lat1 and lon1 in the stack, for the 1st point P1, then for the 2nd point P2 lat2 and lon2 again in the stack, all of them in the format DD.mmsssss (be careful, not decimal degrees or radians). Then press P1P2—> After a few seconds you get the distance s, the forward azimuth alpha12 (initial bearing) and the backward azimuth alpha21. 2) The direct method/problem: you enter the coordinates of P1 in then stack, i. e. lat1 and lon2, then enter the chosen distance s in meters and the forward azimuth alpha12 (initial bearing) in DD.mmsssss. Then press P1—>P2. After a few seconds you get the coordinates of the new point P2, i.e. lat2 and lon2, as well as the backward azimuth alpha 21. Iterations are made according to Vincenty's formulae. Be careful: in this version of the indirect method/problem, the program P1P2—> might not find a solution for (near) antipodal points P1/P2; but those are special, rare cases. You can change the radius a and b according to your preferred geoid; just modify the values in both programs P1P2—> and P1—>P2 and choose always meters (and not km or other units). All "degrees" results are given in DD.mmssss, so that you can check the results of one program directly entering those found values for the other program. Remarks are welcome. Gil RE: HP49-50G Geodesic distance calculator - Gil - 06-23-2020 04:02 PM 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 RE: HP49-50G Geodesic distance calculator - Gil - 06-23-2020 05:02 PM Two more points A) Iin the previous transcription (given by the program inout): Gs stands for "Greek s", i. e. "sigma"; Ga stands for "Greek a", i. e. "alpha" ; Gl stands for "Greek l", i. e. "lambda". B) To work, when entering the four values in the "stack", the calculator should be in RPN mode. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 06-23-2020 05:39 PM I forgot to include the subprogram ATAN2, whose location is the utmost parent directory HOME : ATAN2 \<< \-> y x \<< CASE x 0 > THEN '2*ATAN(y/(\v/(x^2+y^2)+x))' END x 0 \<= y 0 \=/ AND THEN '2*ATAN((\v/(x^2+y^2)-x)/y)' END x 0 < y 0 == AND THEN '4*ATAN(1)' END x 0 == y 0 == AND THEN "Undef" END END \->NUM \>> \>> You might choose the alternative formula and program: ATAN2 \<< \-> y x \<< CASE x 0 > THEN y x / ATAN END x 0 < y 0 \>= AND THEN y x / ATAN -17 FS? IF THEN \pi ELSE 180 END + END x 0 < y 0 < AND THEN y x / ATAN -17 FS? IF THEN \pi ELSE 180 END - END x 0 == y 0 > AND THEN -17 FS? IF THEN \pi 2 / ELSE 90 END END x 0 == y 0 < AND THEN -17 FS? IF THEN \pi NEG 2 / ELSE -90 END END "Undef" END \->NUM \>> \>> Note \v/ stands for square root. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 06-25-2020 03:40 AM 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 RE: HP49-50G Geodesic distance calculator - Gil - 07-07-2020 10:23 AM New version 2.3 As the circumference calculation (lat—>R) on a constant latitude was done with the wrong formula. Now it matches when checking with Visenty's formula. Regards RE: HP49-50G Geodesic distance calculator - Gil - 02-24-2021 02:29 PM Version 2.5 For that earth distance calculator — ccurateness within 1 mm —, I just added here the possibility to modulate the radius equatorial radius a, the pole radius b and flatness f (or INV.f), knowing that (a-b)/a = INV.f (Note that you cannot enter the value of the flatness f, but instead its inverse value INV. f (=1/f). If you know f — instead of 1/f —, then put f in the stack and press the key 1/x (juste above the 9-key).) For that, if changes really required, just run before everything the small "solver" program pressing 1) C.abc ENTER 2) Then give two values: - for a and b, then press Left-Shift INV.f (to get INV.f = 1/f) ; - or for a and INV.f, then press LEFT-SHIFT b (to get the pole radius b) ; - or for b and INV.f, then press Left-Shift a (to get the equatorial radius a). 3) Then run normally (as many times you want) - P1P2—> - or P1—>P2 - or lat—>R without having to repeat anymore the steps 1 and 2 above. Code for C.abc \<< "WGS84 a:6378137 INV.f:298.257223563 " DROP '(a-b)/a-1/INV.f' STEQ 30 MENU \>> Enjoy and regards, Gil Campart RE: HP49-50G Geodesic distance calculator - Gil - 02-24-2021 08:01 PM Small correction in prog lat—>R \<< "1 single Argument: lat in D.mmss To change a, b, INV.f run before all C.abf " DROP 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 Centre-Earth(" lat + ")" + \->TAG 'b/a*TAN(la)' \->NUM ATAN COS a * 2 * \pi * \->NUM "Circumf(" lat + ")" + \->TAG \>> \>> You should have only one DROP instruction before DEG. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 02-25-2021 03:19 AM Slight 'appearance" change of the order of the variables. Here is the full Directory. Copy it in your calculator under the file/Directory name DISTANCE à. Then with your calculator enter your new directory DISTANCE pressing DISTANCE ENTER. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 02-26-2021 04:29 PM New version 2.8b When entering, in D.mmsss Lat1 Long1 Lat2 Long2 and pressing P1P2—> Then the above program used to convert the given lat/long into decimal degrees, then into radian. And at the end it converted back the radian lat/long into the initial D.mmsss format, losing, sometimes, some of the original lat/long values. Now this has been corrected. The same applied to the P1—>P2 program with the given input lat1, lon1 and alpha degree-values. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 02-27-2021 03:33 PM Version 2.9 Changed the subprogram lat—> to calculate the radius from Earth centre to the given surface latitude. It was already correct, but, for special latitude = +/-90 (at the poles), it did not give the expected exact radius b =6356752.31425, but 6356752.3142, due to roundings in the formulae. The same occurred with latitude = 0 (at equator line), reckoned 6378137.0001 in the program, instead of the true value of the radius a = 6378137. Regards, Gil Campart RE: HP49-50G Geodesic distance calculator - Gil - 02-28-2021 04:16 PM New version 3 Geodesic Distance Calculator Including now an attempt to take into account the height/altitude above sea level of two points. For multiple point calculation of journeys, the height option is not available. The new program is called h—>s.h It's short for "from h (actually 2 heights: height A in stack level 2; height B in stack level) you get the distance s.h (distance above sea level, h standing for height)". To use it, you have before to have run once the program P1P2—> or P1—>P2. Observe that running P1P2—> or P1—>P2 will set h1 and h2 to 0 ; however it keeps the values of s.h and delta.s.h The calculation given here is a really rough approximation, giving wrong results. No maths or geodesic theory behind it, contrary to the main program P1P2—> or P1—>P2, based on solid Vincenty's formulae. The purpose is to remember that very accurate formulae like Vincenty's on the Earth surface might diverge from the reality, as most of the time the effective points are not on sea level. If you have better ideas please feel to divide them with the community. I adopted the following points: 0) Point A(lat 30, long 0) and B(lat 60, long 120). First Approach 1) Find effective distance s on Earth Surface without taking into account the height. s = 8635843.51351 2) Approximation 1 (suppose a sphere) cos AB= sin(pi/2-lat1)*sin(pi/2-lat2)*cos(lon1-lon2) +cos(pi/2-lat1)*cos(pi/2-lat2). theta = arcos(AB) =1.35256181244 R. mean × theta_in_radian = s 2b)R. mean × theta_in_radian = s R. mean =s / theta_in_radian =6384805.06701 3) Use 2), but, instead of R.mean, calculating new distance with "R.mean + height_min". —> (R.mean + height_min) × theta_in_radian = s_height_min = 8636384.53824 where height_min = min(height point A, height point B). 4) Use Pythagore requested_final_s_in_altitude = (s_height_min² + (height point A — height point B)²)^.5. Final distance in altitude =8636384.54751 It seems to me that it is being a "fair" (reasonable) approximation for tackling the problem. Second "solution" 1) R from centre to point A or B (((a^2*COS(la))^2+(b^2*SIN(la))^2) /((a*COS(la))^2+(b*SIN(la))^2)) Calculate then radius for A and B, R.A and R.B. R.A 6372824.4203 R.B 6362132.22441 2) Calculate the mean radius R.mean (R.A + R.B) / 2 = 6367478.32235 3) R.New = R.Mean + MIN (alt.A, alt B) =6367878.32235 4) As the angle theta has not changed S.New.horizontal = s × R.New/R.mean. 8635843.51351 × 6367878.32235/6367478.32235 =8636386.01046 5) Final distance s_altitude = ((S.New.horizontal² +(alt.A—alt B)²)^.5 = 8636386.01972 Both results are very similar and seem to make sense. Thanks for your commentaries. Regards, Gil Campart RE: HP49-50G Geodesic distance calculator - Gil - 03-01-2021 09:47 AM For HP49-HP50, programs basically to calculate exact distances (error less than 0.5 mm) on the Earth surface at theorical sea level Present Version 3.2 "corrects" my previous reasoning for approximated real distance above sea level (altitudes h1 & h2). Explanation From the two points P1 & P2, I calculate a first "possible" angle (mean.thets) between them, then a mean radius (=s/mean.theta), considering the ellipsoid to be a sphere. I calculate also radii R.1 & R.2 from Earth centre to point P1 and to P2 (with program lat—>R), then I build a second mean radius.2= (R1+R2)/2. From those two mean radii (meanR12), I have to consider a new, mean, greater radius at altitude (h1+h2)/2. In other words, s.h (distance between points P1 & P2 at altitude h1 & h2) =s (original distance at sea level) × [meanR12+ (h1+h2)/2]/meanR12. Inside are the five user programs : A) Four very accurate programs 1) C.abh (to calibrate first your ellipsoid) 2) P1P2—> (indirect method/problem, to find the distance) 3) P1—>P2 (direct method/problem, to find point P2) ) 4) lat—> (gives - distance from Earth centre to latitude at sea level - & Earth circumference along that latitude) B) Essay to give an "idea" of the "real", effective distance above sea level. No real maths or geodesic, scentifical consideration. Just a possible approximation to give a rough idea of the difference when calculating very "accurately", but only at Earth sea level. 5) h—>s.h. Observations or commentaries welcome. Regards, Gil RE: HP49-50G Geodesic distance calculator - Gil - 03-11-2021 09:53 AM 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. 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 RE: HP49-50G Geodesic distance calculator - Gil - 03-14-2021 08:46 AM RE: HP49-50G Geodesic distance calculator Version 4.2 New Very slight, "cosmetical" changes. 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. Regards, Gil Campart RE: HP49-50G Geodesic distance calculator - Gil - 03-15-2021 09:25 AM 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. Regards, Gil Campart RE: HP49-50G Geodesic distance calculator - Gil - 03-19-2021 05: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. Regards, Gil RE: (49) (50) Geodesic distance calculator - Gil - 03-20-2021 06: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 Regards, Gil RE: (49) (50) Geodesic distance calculator - Gil - 03-25-2021 09:17 AM 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. Regards, Gil Campart RE: (49) (50) Geodesic distance calculator - Gil - 03-27-2021 06: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 RE: (49) (50) Geodesic distance calculator - Gil - 03-27-2021 06: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 RE: (49) (50) Geodesic distance calculator - Gil - 03-30-2021 05: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] Regstds RE: (49) (50) Geodesic distance calculator - Gil - 04-02-2021 05: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 RE: (49) (50) Geodesic distance calculator - Gil - 04-04-2021 05: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 RE: (49) (50) Geodesic distance calculator - Gil - 04-07-2021 01: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... RE: (49) (50) Geodesic distance calculator - Gil - 05-14-2021 05: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). RE: (49) (50) Geodesic distance calculator - Gil - 05-14-2021 06: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 RE: 49 50 Ver06b.hp Geodesic distance & Earth Euclidean distance calculator, beari... - Gil - 03-29-2023 05:54 PM Correction in program P1 —>P2. Regards, Gil RE: 49 50 Ver06b.hp Geodesic distance & Earth Euclidean distance calculator, beari... - Gil - 01-05-2024 06:42 PM 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. RE: 49 50G Ver7 Geodesic distance & Earth Euclidean distance calculator, bearing - Gil - 12-27-2024 08:35 AM New version 7 Includes now Andoyer 2nd order formula for distance calculation. RE: 49 50 Ver8.hp Geodesic distance & Earth Euclidean distance calculator, bearing - Gil - 2025-02-01 07:03 AM Version 8 As before. Just used, for increased accurateness, the built-in D->R and R—>D functions of the HP50G (instead of the manual conversion by 180/pi). |