|
|
|
#!/usr/bin/env rexx
/* ---------------------------------------------------------------- */
/* 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 */
/* ---------------------------------------------------------------- */
-- Uncomment next 2 lines for classic Rexx and comment ::rquires
--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
-- Comment next line and uncomment old RxFuncAdd stuff above for classic Rexx
::requires 'rxmath' LIBRARY
| If you feel inclined to make corrections, suggestions etc.,
please mail me any.
| |
|
All content © Ruurd Idenburg, 2007–,
except where marked otherwise. All rights reserved. This page is primarily for non-commercial use only.
The Idenburg website records no personal information and sets no ‘cookies’.
This site is hosted on a VPS(Virtual Private System) rented from Transip.nl, a Dutch company, falling under Dutch (privacy) laws (I think).
This page updated on by Ruurd Idenburg.
|
|