def sepdeg(a1,b1,a2,b2):
#
# find the separation of two points on the CS 
# inputs and out in degrees
# derived from a STARLINK SLA routine
#
    const = math.pi/180.
    ra1 = a1*const
    rb1 = b1*const
    ra2 = a2*const
    rb2 = b2*const
 
    v1_1 = math.cos(ra1)*math.cos(rb1)
    v1_2 = math.sin(ra1)*math.cos(rb1)
    v1_3 = math.sin(rb1)
 
    v2_1 = math.cos(ra2)*math.cos(rb2)
    v2_2 = math.sin(ra2)*math.cos(rb2)
    v2_3 = math.sin(rb2)
 
    w = ( (v1_1-v2_1)**2 + (v1_2-v2_2)**2 + (v1_3-v2_3)**2 )/4.0
 
    x = math.sqrt(w)
    y = math.sqrt(max(0.0,1.0-w))
    angle = 2.0*math.atan2(x,y)/const
    return angle