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