1: #!/usr/bin/env rexx
   2: /* ---------------------------------------------------------------- */
   3: /* Calculates the distance between (possibly a series) of           */
   4: /* 2 latitude/longitude pairs. May be used in a pipe.               */
   5: /* ---------------------------------------------------------------- */
   6: /*                                                                  */
   7: /* Originally by Ruurd J. Idenburg                                  */                                                             
   8: /*                                                                  */
   9: /*                                                                  */
  10: /* No copyright, no licence, no guarantees or warrantees, be it     */
  11: /* explicit, implicit or whatever. Usage is totally and completely  */
  12: /* at the users own risk, the author shall not be liable for any    */ 
  13: /* damages whatsoever, for any reason whatsoever.                   */
  14: /*                                                                  */
  15: /* Please keep this comment block intact when modifying this code   */
  16: /* and add a note with date and a description.                      */
  17: /*                                                                  */
  18: /* ---------------------------------------------------------------- */
  19: /*  Parameter(s):                                                   */ 
  20: /*                                                                  */
  21: /*    None                                                          */
  22: /*                                                                  */
  23: /*  Result:                                                         */
  24: /*                                                                  */
  25: /*    For each pair of latitude/longitude the distance in thousands */
  26: /*    of a kilometer to the previous pair is going to STDOUT        */
  27: /*    (usually the console).                                        */
  28: /*                                                                  */
  29: /* ---------------------------------------------------------------- */
  30: /* 2007/12/31 - Initial version approximately.                      */
  31: /* 2011/12/31 - Isolate calculus to allow use outside a pipe.       */
  32: /* 2021/05/20 - Updated shebang and ::requires iso RxMathLoadFunc   */
  33: /* ---------------------------------------------------------------- */
  34: -- Uncomment next 2 lines for classic Rexx and comment ::rquires
  35: --Call RxFuncAdd "MathLoadFuncs","rxmath","MathLoadFuncs"
  36: --Call MathLoadFuncs
  37: 
  38: /* ---------------------------------------------------------------- */
  39: /* The following formula's can be used to calculate the distance    */
  40: /* between two points on earth:                                     */ 
  41: /*                                                                  */
  42: /* d=R*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2)) */
  43: /*                                                                  */
  44: /* or                                                               */ 
  45: /*                                                                  */
  46: /* d=R*2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)*    */
  47: /*   (sin((lon1-lon2)/2)^2))                                        */
  48: /*                                                                  */
  49: /* where                                                            */
  50: /*                                                                  */
  51: /* R -- the radius of the earth, commonly taken to be: 6378,137 km  */
  52: /* lat1/lat2 -- latitude in degrees (NOT degrees minutes seconds)   */
  53: /* lon1/lon2 -- longitude in degrees (NOT degrees minutes seconds   */
  54: /*                                                                  */
  55: /* The first formula is more subject to possible rounding errors    */
  56: /* for short distances and since the route description produced by  */
  57: /* Google Earth/Maps often has points within a couple of meters     */
  58: /* this code uses the second formula.                               */
  59: /*                                                                  */
  60: /* ---------------------------------------------------------------- */
  61: 
  62: --oldradius = 6366710
  63: googleradius = 6378137 -- Equatorial earth radius used by Google
  64: 
  65: -- Doesn't seem necessary for Windows, might be necessary for Linux
  66: --Do Until Lines()>0
  67: --  Call SysSleep 0.1
  68: --End
  69: 
  70: Parse Value Linein() With flat flon tlat tlon
  71: If (tlat<>'' & tlon<>'') Then Say Distance(flat,flon,tlat,tlon)
  72: Do While Lines()>0
  73:   Parse Value Linein() With tlat tlon .
  74:   Say Distance(flat,flon,tlat,tlon)
  75:   flat = tlat
  76:   flon = tlon
  77: End
  78: Exit
  79: 
  80: Distance: Procedure Expose googleradius
  81:   Parse Arg flat, flon, tlat, tlon
  82:   sinlat = RxCalcSin((flat-tlat)/2)
  83:   cosplat = RxCalcCos(flat)
  84:   coslati = RxCalcCos(tlat)
  85:   sinlon = RxCalcSin((flon-tlon)/2)
  86:   sinlatp = (sinlat)**2
  87:   sinlonp = (sinlon)**2
  88:   radians = 2*RxCalcArcSin(RxCalcSqrt(sinlatp + cosplat*coslati*sinlonp),,'R')
  89:   distance = (radians*googleradius)~format(,0)/1000 'km'
  90: Return distance -- returns distance in thousands of a Km
  91: -- Comment next line and uncomment old RxFuncAdd stuff above for classic Rexx 
  92: ::requires 'rxmath' LIBRARY