1 files changed, 48 insertions, 0 deletions
diff --git a/src/slalib/sep.f b/src/slalib/sep.f
new file mode 100644
index 0000000..c977985
--- /dev/null
+++ b/src/slalib/sep.f
@@ -0,0 +1,48 @@
+      REAL FUNCTION sla_SEP (A1, B1, A2, B2)
+*+
+*     - - - -
+*      S E P
+*     - - - -
+*
+*  Angle between two points on a sphere (single precision)
+*
+*  Given:
+*     A1,B1    real    spherical coordinates of one point
+*     A2,B2    real    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_CS2C
+*
+*  P.T.Wallace   Starlink   April 1985
+*
+*  Copyright (C) 1995 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      REAL A1,B1,A2,B2
+
+      INTEGER I
+      REAL V1(3),V2(3),W
+
+
+
+*  Convert coordinates from spherical to Cartesian
+      CALL sla_CS2C(A1,B1,V1)
+      CALL sla_CS2C(A2,B2,V2)
+
+*  Modulus squared of half the difference vector
+      W=0.0
+      DO I=1,3
+         W=W+(V1(I)-V2(I))**2
+      END DO
+      W=W/4.0
+
+*  Angle between the vectors
+      sla_SEP=2.0*ATAN2(SQRT(W),SQRT(MAX(0.0,1.0-W)))
+
+      END