1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
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
|
