1 files changed, 48 insertions, 0 deletions
diff --git a/src/slalib/dsep.f b/src/slalib/dsep.f
new file mode 100644
index 0000000..b1bd4be
--- /dev/null
+++ b/src/slalib/dsep.f
@@ -0,0 +1,48 @@
+      DOUBLE PRECISION FUNCTION sla_DSEP (A1, B1, A2, B2)
+*+
+*     - - - - -
+*      D S E P
+*     - - - - -
+*
+*  Angle between two points on a sphere (double precision)
+*
+*  Given:
+*     A1,B1    dp    spherical coordinates of one point
+*     A2,B2    dp    spherical coordinates of the other point
+*
+*  (The spherical coordinates are RA,Dec, Long,Lat etc, in radians.)
+*
+*  The result is the angle, in radians, between the two points.  It
+*  is always positive.
+*
+*  Called:  sla_DCS2C
+*
+*  P.T.Wallace   Starlink   April 1985
+*
+*  Copyright (C) 1995 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      DOUBLE PRECISION A1,B1,A2,B2
+
+      INTEGER I
+      DOUBLE PRECISION V1(3),V2(3),W
+
+
+
+*  Convert coordinates from spherical to Cartesian
+      CALL sla_DCS2C(A1,B1,V1)
+      CALL sla_DCS2C(A2,B2,V2)
+
+*  Modulus squared of half the difference vector
+      W=0D0
+      DO I=1,3
+         W=W+(V1(I)-V2(I))**2
+      END DO
+      W=W/4D0
+
+*  Angle between the vectors
+      sla_DSEP=2D0*ATAN2(SQRT(W),SQRT(MAX(0D0,1D0-W)))
+
+      END