import numpy as np
def spherdist(lon1, lat1, lon2, lat2):
   """
   Computes the distance between geo. pos.
   lon1,lat1 and lon2,lat2.
   input is in degrees.

   output is real*4 for better global consistancy,
   by truncating double precision roundoff errors.
   real*4 is not in f90, but is widely supported.

   Based on m_spherdist.F90 from Geir Evanson.
   """
   invradian=0.017453292
   rearth=6371001.0      # Radius of earth
   deg360 = 360.0 ; deg0 = 0.0
   deg90 = 90.0 ; one = 1.0
#
#     ensure that spherdist(ax,ay,bx,by) == spherdist(bx,by,ax,ay)
#
   dlon1 = lon1
   dlon1 = dlon1 % deg360
   if (dlon1 < deg0): dlon1 = dlon1 + deg360
   lat1 = lat1

   dlon2 = lon2
   dlon2 = dlon2 % deg360
   if (dlon2 < deg0) : dlon2 = dlon2 + deg360
   lat2 = lat2


   if (lat1 < lat2) :
     rlon1=dlon1*invradian            #lon1 in rad
     rlat1=(deg90-lat1)*invradian     #90-lat1 in rad
     rlon2=dlon2*invradian            #lon2 in rad
     rlat2=(deg90-lat2)*invradian     #90-lat2 in rad
   elif (lat1 == lat2) and (dlon1 <= dlon2) :
     rlon1=dlon1*invradian            #lon1 in rad
     rlat1=(deg90-lat1)*invradian     #90-lat1 in rad
     rlon2=dlon2*invradian            #lon2 in rad
     rlat2=(deg90-lat2)*invradian     #90-lat2 in rad
   else :
     rlon2=dlon1*invradian            #lon1 in rad
     rlat2=(deg90-lat1)*invradian     #90-lat1 in rad
     rlon1=dlon2*invradian            #lon2 in rad
     rlat1=(deg90-lat2)*invradian     #90-lat2 in rad




##
   x1= np.sin(rlat1)*np.cos(rlon1)        #x,y,z of pos 1.
   y1= np.sin(rlat1)*np.sin(rlon1)
   z1= np.cos(rlat1)
   ##
   x2= np.sin(rlat2)*np.cos(rlon2)        #x,y,z of pos 2.
   y2= np.sin(rlat2)*np.sin(rlon2)
   z2= np.cos(rlat2)
##
   dr=np.arccos(np.min([one,x1*x2+y1*y2+z1*z2]))  # Arc length
   spher=dr*rearth
      #spher = float(spher)
##
   return spher