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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
#{ SKYBOT -- Find the minor planets in an image, display it and mark the
# asteroids.
procedure skybot (image)
string image { prompt = "Input image" }
bool display = yes { prompt = "Display result?" }
bool ned = no { prompt = "Overlay NED sources?" }
bool grid = no { prompt = "Overlay coordinate grid?" }
real nhours = 6.0 { prompt = "Number of hours of movement" }
bool verbose = yes { prompt = "Verbose output?" }
int status = 0 { prompt = "Service status code" }
begin
string img, date, res, coords, tvcoords, cmd
real ra, dec, size, mjd, nh
real radius, x1, x2, y1, y2, vmag, dra, ddec
bool verb, disp, do_ned, do_grid
img = image
verb = verbose
disp = display
nh = nhours
do_ned = ned
do_grid = grid
set clobber = yes
set imclobber = yes
if (imaccess (img)) {
# Compute the MJD from the DATE-OBS keyword
hselect (img, "DATE-OBS", yes) | scan (date)
if (date == "") {
errror (0, "No DATE-OBS keyword in the image.")
} else {
print ("print(julday(\""//date//"\"))") | \
astcalc (commands="STDIN") | scan (mjd)
}
iferr { wcsinfo (img) } then {
error (0, "Cannot determine image coords for `"//img//"'")
} else {
ra = wcsinfo.midx
dec = wcsinfo.midy
size = max (wcsinfo.width, wcsinfo.height) * 60.0 / 2.0
}
} else
error (0, "Image not found.")
# Create temp files for the output
res = mktemp ("tmp$imsky")
# Call the SKYBOT task to do the search.
sbquery(ra, dec, size, epoc=mjd, fields="ra,dec,vmag,dra,ddec,name", >& res)
# if (verbose || display == no)
# type (res)
if (display) {
cmd = mktemp ("tmp$imsky")
coords = mktemp ("tmp$imsky")
tvcoords = mktemp ("tmp$imsky")
display (img, 1, fill+, bpmask="", select+, >& "dev$null")
fields (res, "1-", >& "foo")
list = "foo"
while (fscan (list, ra, dec, vmag, dra, ddec, line) != EOF) {
printf ("%h %h %.1f\n", ra, dec, vmag) |& \
wcsctran ("STDIN", "STDOUT", img, "world", "logical",
units="h n", verbose-, >& tvcoords )
# Compute a radius representing 2-hours of movement.
x = dra # in "/hr
y = ddec # in "/hr
rad = max (10., sqrt((x*x)+(y*y)) * wcsinfo.scale) # in pixels
type (tvcoords) | scan (x, y)
printf ("%-11h %-12h %7.1f %7.1f %4.1f %6.4f %6.4f %4.1f %s\n",
ra, dec, x, y, vmag, dra, ddec, rad, line)
tvmark (1, tvcoords, mark="circle", radii="10", col=205,
nxoffset=12, nyoffset=2, lab+)
x1 = ra - (nh * (dra / 3600.) / 15.)
y1 = dec - (nh * (ddec / 3600.))
x2 = ra + (nh * (dra / 3600.) / 15.)
y2 = dec + (nh * (ddec / 3600.))
printf ("%h %h 101 s\n%h %h 101 s\n",
x1, y1, (ra-(dra/3600./15./2.)), (dec-(ddec/3600./2.))) |&
wcsctran ("STDIN", "STDOUT", img, "world", "logical",
units="h n", verbose-, >& tvcoords )
tvmark (1, coords="", commands=tvcoords, interact-, col=206, lab-)
printf ("%h %h 101 s\n%h %h 101 s\n",
(ra+(dra/3600./15./2.)), (dec+(ddec/3600./2.)), x2, y2) |&
wcsctran ("STDIN", "STDOUT", img, "world", "logical",
units="h n", verbose-, >& tvcoords )
tvmark (1, coords="", commands=tvcoords, interact-, col=204, lab-)
}
# Overlay the NED objects?
if (do_ned)
nedoverlay (img)
# Draw the coordinate overlay grid?
if (do_grid)
wcslab (img, 1, use-, fill+, overplot+, append+,
labout-, dev="imdy")
# Clean up
delete (cmd, verify-, >& "dev$null")
delete (coords, verify-, >& "dev$null")
delete (tvcoords, verify-, >& "dev$null")
}
# Clean up
delete (res, verify-, >& "dev$null")
end
|
