m_lldist_p <- function(long,lat,long0,lat0){
# M_LLDIST Spherical earth distance between points in long/lat coordinates. 
#   RANGE<-M_LLDIST_P(LONG,LAT,LONG0,LAT0) gives the distance in kilometers between
#   points in the vectors LONG and LAT and the point with coordinates LONG0 and LAT0,computed
#   using the Haversine formula on a spherical earth of radius
#   6378.137km. Distances are probably good to better than 1# of the
#   "true" distance on the ellipsoidal earth
#


# Rich Pawlowicz (rich@ocgy.ubc.ca) 6/Nov/00
# This software is provided "as is" without warranty of any kind. But
# it's mine, so you can't sell it.
#
# 30/Dec/2005 - added n-point geodesic computations, based on an algorithm
#               coded by Jeff Barton at Johns Hopkins APL in an m-file
#               I looked at at mathworks.com.


pi180 <- pi/180;
earth_radius<-6378.137;

m<-length(long);

long1<-long*pi180;
long2<-matrix(long0,m,1)*pi180;
lat1<- lat*pi180;
lat2<- matrix(lat0,m,1)*pi180;

dlon <- long2 - long1; 
dlat <- lat2 - lat1; 
a <- (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2;
angles <- 2 * atan2( sqrt(a), sqrt(1-a) );
dist <- earth_radius * angles;
dist
}