91 lines
4.9 KiB
Rexx
Executable File
91 lines
4.9 KiB
Rexx
Executable File
|
|
/* ---------------------------------------------------------------- */
|
|
/* Calculates the distance between (possibly a series) of */
|
|
/* 2 latitude/longitude pairs. May be used in a pipe. */
|
|
/* ---------------------------------------------------------------- */
|
|
/* */
|
|
/* Originally by Ruurd J. Idenburg */
|
|
/* */
|
|
/* */
|
|
/* No copyright, no licence, no guarantees or warrantees, be it */
|
|
/* explicit, implicit or whatever. Usage is totally and completely */
|
|
/* at the users own risk, the author shall not be liable for any */
|
|
/* damages whatsoever, for any reason whatsoever. */
|
|
/* */
|
|
/* Please keep this comment block intact when modifying this code */
|
|
/* and add a note with date and a description. */
|
|
/* */
|
|
/* ---------------------------------------------------------------- */
|
|
/* Parameter(s): */
|
|
/* */
|
|
/* None */
|
|
/* */
|
|
/* Result: */
|
|
/* */
|
|
/* For each pair of latitude/longitude the distance in thousands */
|
|
/* of a kilometer to the previous pair is going to STDOUT */
|
|
/* (usually the console). */
|
|
/* */
|
|
/* ---------------------------------------------------------------- */
|
|
/* 2007/12/31 - Initial version approximately. */
|
|
/* 2011/12/31 - Isolate calculus to allow use outside a pipe. */
|
|
/* 2021/05/20 - Updated shebang and ::requires iso RxMathLoadFunc */
|
|
/* ---------------------------------------------------------------- */
|
|
Call RxFuncAdd "MathLoadFuncs","rxmath","MathLoadFuncs"
|
|
Call MathLoadFuncs
|
|
|
|
/* ---------------------------------------------------------------- */
|
|
/* The following formula's can be used to calculate the distance */
|
|
/* between two points on earth: */
|
|
/* */
|
|
/* d=R*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2)) */
|
|
/* */
|
|
/* or */
|
|
/* */
|
|
/* d=R*2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)* */
|
|
/* (sin((lon1-lon2)/2)^2)) */
|
|
/* */
|
|
/* where */
|
|
/* */
|
|
/* R -- the radius of the earth, commonly taken to be: 6378,137 km */
|
|
/* lat1/lat2 -- latitude in degrees (NOT degrees minutes seconds) */
|
|
/* lon1/lon2 -- longitude in degrees (NOT degrees minutes seconds */
|
|
/* */
|
|
/* The first formula is more subject to possible rounding errors */
|
|
/* for short distances and since the route description produced by */
|
|
/* Google Earth/Maps often has points within a couple of meters */
|
|
/* this code uses the second formula. */
|
|
/* */
|
|
/* ---------------------------------------------------------------- */
|
|
|
|
/* -oldradius = 6366710 */
|
|
googleradius = 6378137 /* -- Equatorial earth radius used by Google */
|
|
|
|
/* Doesn't seem necessary for Windows, might be necessary for Linux */
|
|
Do Until Lines()>0
|
|
Call SysSleep 0.1
|
|
End
|
|
|
|
Parse Value Linein() With flat flon tlat tlon
|
|
If (tlat<>'' & tlon<>'') Then Say Distance(flat,flon,tlat,tlon)
|
|
Do While Lines()>0
|
|
Parse Value Linein() With tlat tlon .
|
|
Say Distance(flat,flon,tlat,tlon)
|
|
flat = tlat
|
|
flon = tlon
|
|
End
|
|
Exit
|
|
|
|
Distance: Procedure Expose googleradius
|
|
Parse Arg flat, flon, tlat, tlon
|
|
sinlat = RxCalcSin((flat-tlat)/2)
|
|
cosplat = RxCalcCos(flat)
|
|
coslati = RxCalcCos(tlat)
|
|
sinlon = RxCalcSin((flon-tlon)/2)
|
|
sinlatp = (sinlat)**2
|
|
sinlonp = (sinlon)**2
|
|
radians = 2*RxCalcArcSin(RxCalcSqrt(sinlatp + cosplat*coslati*sinlonp),,'R')
|
|
distance = (radians*googleradius)~format(,0)/1000 'km'
|
|
Return distance /* -- returns distance in thousands of a Km */
|
|
/* ::requires 'rxmath' LIBRARY */
|