Compute distances between 2 points on earth

/* Calculate surface distance between centroids of two zip or postal codes, */
/* where centroids translated to latitude and longitude */
/* Inputs: 2 zipcodes, 2 latitudes, 2 longitudes */
/* Outputs: Distance between locations in miles or kilometers */
/* Simon Freidin 6/6/2003 */

* Following dummy data made up by Ray.
data list list /zip1 zip2 lat1 lat2 long1 long2.
begin data
1 2 125 126 98 100
1 3 125 126 99 99
end data.

/* only calculate distance if zip codes not identical and latitudes and */
/*  longitudes not identical */
do if ~(zip1 = zip2) and ~((lat1=lat2) and (long1=long2)).
/* calculate great circle distance */
/* http://codewalkers.com/getcode.php?id=247 */
compute theta = long2 - long1.
compute d1 = (sin((artan(1)/45)*(lat2)) * sin((artan(1)/45)*(lat1))) +
              (cos((artan(1)/45)*(lat2)) * cos((artan(1)/45)*(lat1)) * 
		cos((artan(1)/45)*(theta))).
if range(d1,-1,1) d2 = 2*artan(1) - arsin(d1) .
compute d3 = (45/artan(1))*(d2).
compute distmile = d3 * 60 * 1.1515 .
compute distkm = distmile * 1.609344.
else.
/* hasn't moved */
compute distmile=0.
compute distkm=0.
end if.

compute distkm=rnd(distkm).
compute distmile=rnd(distmile).

variable labels
  distmile 'DV: Distance (miles)'/
  distkm   'DV: Distance (kilometers)'/.

*Research Database Manager and Analyst
*Melbourne Institute of Applied Economic and Social Research
*The University of Melbourne
*MELBOURNE VIC 3010
*Tel: (03) 8344 3601 Fax: (03) 8344 5630.