diff options
Diffstat (limited to 'noao/nproto/skysep.cl')
| -rw-r--r-- | noao/nproto/skysep.cl | 41 |
1 files changed, 41 insertions, 0 deletions
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 |