From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 noao/nproto/skysep.cl | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100644 noao/nproto/skysep.cl

(limited to 'noao/nproto/skysep.cl')

diff --git a/noao/nproto/skysep.cl b/noao/nproto/skysep.cl
new file mode 100644
index 00000000..d70742d5
--- /dev/null
+++ b/noao/nproto/skysep.cl
@@ -0,0 +1,41 @@
+# SEP -- Separation between two celestial coordinates.
+
+procedure sep (ra1, dec1, ra2, dec2)
+
+real	ra1			{ prompt="RA (hr|deg)"}
+real	dec1			{ prompt="DEC (deg)"}
+real	ra2			{ prompt="RA (hr|deg)"}
+real	dec2			{ prompt="DEC (deg)"}
+string	raunit = "hr"		{ prompt="RA unit (hr|deg)", enum="hr|deg" }
+bool	verbose = no		{ prompt="Verbose?"}
+real	sep			{ prompt="Separation (arcsec)"}
+
+begin
+	real	r1, d1, r2, d2
+	real	c1, c2, x, y, z
+
+	if (raunit == "hr") {
+	    r1 = ra1 * 15.
+	    d1 = dec1
+	    r2 = ra2 * 15.
+	    d2 = dec2
+	} else {
+	    r1 = ra1
+	    d1 = dec1
+	    r2 = ra2
+	    d2 = dec2
+	}
+
+	c1 = dcos(d1)
+	c2 = dcos(d2)
+	x = dcos(r1) * c1 - dcos(r2) * c2
+	y = dsin(r1) * c1 - dsin(r2) * c2
+	z = dsin(d1) - dsin(d2)
+	c1 = (x*x + y*y + z*z) / 4.
+	c2 = max (0., 1.-c1)
+	sep = 2 * datan2(sqrt(c1),sqrt(c2)) * 3600
+
+	if (verbose)
+	    printf ("%.2f arcsec = (%.2H, %.1h) - (%.2H, %.1h)\n",
+		sep, r1, d1, r2, d2)
+end
-- 
cgit