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

[Gil]

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.

Here is the code
h—>s.h
\<< "Before, run

P1P2\->D
 or 
P1\->P2
.

2 methods \-> 2 Results
as a complex#" DROP "Rough 
approximation!
" UNROT 'h2' STO 'h1' STO STD RAD \pi 2 / lat1 D\->RAD - SIN \pi 2 / lat2 D\->RAD - SIN * lon1 D\->RAD lon2 D\->RAD - COS * \pi 2 / lat1 D\->RAD - COS \pi 2 / lat2 D\->RAD - COS * + \->NUM ACOS \-> \Gm\Gh
\<< s \Gm\Gh / lat1 lat\->R DROP lat2 lat\->R DROP + 2 / 2 \->LIST \-> \GmR
\<< \GmR 1 GET h1 h2 MIN + \Gm\Gh * SQ h1 h2 - SQ + \v/ \GmR 2 GET DUP h1 h2 MIN + SWAP / s * SQ h1 h2 - SQ + \v/ R\->C DUP 's.h' STO "s(" h1 + "/" + h2 + ")" + \->TAG s "s(0,0)" \->TAG SWAP s.h s s R\->C - DUP '\GDs.h' STO "\GDs.h[m]" \->TAG
\>>
\>>
\>>

As usual, as there are some slight changes in the directory, I have included the full Directory to be copied, possibly with a shorter name, for example DIST.

Regards,
Gil Campart