/* ---------------------------------------------------------------- */
/* 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 */