|
/* REXX EXEC to calculate a great circle segment */
arg latd1 latm1 lond1 lonm1 latd2 latm2 lond2 lonm2
call RxFuncAdd 'MathLoadFuncs', 'RxMath', 'MathLoadFuncs'
call MathLoadFuncs
pi = 3.141592653589
if latd1 <> "" then do
lata = radians(dms(latd1 latm1 0))
lona = radians(dms(lond1 lonm1 0))
latb = radians(dms(latd2 latm2 0))
lonb = radians(dms(lond2 lonm2 0))
end
else do
say 'Usage is: gcvx.rex <latdeg1> <latmin1> <londeg1> <lonmin1> <latdeg2> <latmin2> <londeg2> <lonmin2>'
say 'Using Frankfurt (EDDF) - Singapore (WSSS) as an example now.'
/*
* Example: Frankfurt - Singapore
* EDDF - WSSS
*/
lona = radians(8.5589)
lonb = radians(103.9892)
lata = radians(50.0251)
latb = radians(1.3596)
end
/*
* Unity Vector to Position 'A'
*/
v1.x = cos(lata) * cos(lona)
v1.y = cos(lata) * sin(lona)
v1.z = sin(lata)
say "V1:" v1.x v1.y v1.z "(abs="||format(sqrt(v1.x*v1.x+v1.y*v1.y+v1.z*v1.z),,6)||")"
/*
* Unity Vector to Position 'B'
*/
v2.x = cos(latb) * cos(lonb)
v2.y = cos(latb) * sin(lonb)
v2.z = sin(latb)
say "V2:" v2.x v2.y v2.z "(abs="||format(sqrt(v2.x*v2.x+v2.y*v2.y+v2.z*v2.z),,6)||")"
/*
* Scalar Product to Obtain Angle between 'A' and 'B' on Disc
*/
ScalProdAB = v1.x * v2.x + v1.y * v2.y + v1.z * v2.z
DistRadians = acos(ScalProdAB)
DistNM = DistRadians * 60 / pi * 180
Distkm = DistNM * 1.852
/*
* Normal Vector on Disc
*/
n.x = v1.y * v2.z - v2.y * v1.z
n.y = -v1.x * v2.z + v2.x * v1.z
n.z = v1.x * v2.y - v2.x * v1.y
nabs = sqrt(n.x*n.x+n.y*n.y+n.z*n.z)
n0.x = n.x / nabs
n0.y = n.y / nabs
n0.z = n.z / nabs
say "N0:" n0.x n0.y n0.z "(abs="||format(sqrt(n0.x*n0.x+n0.y*n0.y+n0.z*n0.z),,6)||")"
/*
* Direction Vector
*/
n1.x = -n0.x
n1.y = -n0.y
n1.z = -n0.z
vabs = sqrt(v1.x*v1.x+v1.y*v1.y+v1.z*v1.z)
v0.x = v1.x / vabs
v0.y = v1.y / vabs
v0.z = v1.z / vabs
d.x = v0.y * n1.z - n1.y * v0.z
d.y = -v0.x * n1.z + n1.x * v0.z
d.z = v0.x * n1.y - n1.x * v0.y
say "D :" d.x d.y d.z "(abs="||format(sqrt(d.x*d.x+d.y*d.y+d.z*d.z),,6)||")"
/*
* Vector pointing to the North Pole along a Meridian
*/
npv.x = -cos(lona) * sin(lata)
npv.y = -sin(lona) * sin(lata)
npv.z = cos(lata)
npvabs = sqrt(npv.x*npv.x+npv.y*npv.y+npv.z*npv.z)
npv.x = npv.x / npvabs
npv.y = npv.y / npvabs
npv.z = npv.z / npvabs
say "NP:" npv.x npv.y npv.z "(abs="||format(sqrt(npv.x*npv.x+npv.y*npv.y+npv.z*npv.z),,6)||")"
/*
* Scalar Product of the North Pole Vector and the Direction Vector
* yields the cos() of the initial True Course
*/
scalprod = npv.x * d.x + npv.y * d.y + npv.z * d.z
say "ScalarProduct=" scalprod
CourseDeg = degrees(acos(scalprod))
if lonb < lona then CourseDeg = 360 - CourseDeg
RevCourseDeg = (CourseDeg + 180) // 360
say "Distance="||DistNM "NM or" Distkm "km"
say "TrueDepartureCourse="||CourseDeg "TrueArrivalCourse="||RevCourseDeg
/*
* Just a check, the Scalar Product of the Vector to 'A'
* with the North Pole Vector must be zero (perpendicular vectors)
*/
sp2 = npv.x * v1.x + npv.y * v1.y + npv.z * v1.z
say format(sp2,,6) "<-- this must be zero"
call MathDropFuncs
exit 0
degrees: procedure expose pi
arg rads
return (rads)/pi*180.0
radians: procedure expose pi
arg degs
return (degs)/180.0*pi
dms: procedure
arg deg min sec
deg = (sec + min * 60 + deg * 3600) / 3600
return deg
/*
* The following are needed when using RxMath with ooREXX.
* With OS/2 REXX RxMath had implementations with these names directly.
*/
sin: procedure
arg a
return RxCalcSin(a,16,'R')
cos: procedure
arg a
return RxCalcCos(a,16,'R')
acos: procedure
arg x
return RxCalcArcCos(x,16,'R')
sqrt: procedure
arg x
return RxCalcSqrt(x,16)
|