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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
include <math.h>
include <pkg/cq.h>
include <pkg/skywcs.h>
include "../../lib/astrom.h"
# AT_TQUERY -- Extract catalog objects from a text file stored in a results
# query structure.
pointer procedure at_tquery (at, cq, cres, hdrtext, nlines, fieldno)
pointer at #I the astrometry package descriptor
pointer cq #I the astrometric catalog descriptor
pointer cres #I the input catalog results descriptor
char hdrtext[ARB] #I the catalog header test
int nlines #I the number of lines in the header text
int fieldno #I the region number
double rac, decc, ra, dec, rawidth, decwidth, ra1, ra2, dec1, dec2, dist
double tra, trac
pointer sp, csystem, raname, decname, funits, tmpname, line
pointer res, ccoo, mw
int i, fd, strfd, stat, units
pointer cq_fquery()
int cq_rstati(), at_rcregion(), open(), strlen(), cq_hinfo(), sk_decwcs()
int stropen(), getline(), strdic(), cq_gvald(), cq_grecord(), sk_stati()
bool streq()
begin
# Return if the input catalog is undefined or contains no records.
if (cres == NULL)
return (NULL)
if (cq_rstati (cres, CQRNRECS) <= 0)
return (NULL)
# Return if the header is undefined.
if (nlines <= 0 || hdrtext[1] == EOS)
return (NULL)
# Get the region to be extracted.
if (at_rcregion (at, cres, fieldno, rac, decc, rawidth,
decwidth) == ERR)
return (NULL)
# Compute the ra and dec limits.
call at_rclimits (rac, decc, rawidth, decwidth, ra1, ra2, dec1, dec2)
# Get some working space.
call smark (sp)
call salloc (csystem, CQ_SZ_FNAME, TY_CHAR)
call salloc (tmpname, SZ_FNAME, TY_CHAR)
call salloc (raname, CQ_SZ_FNAME, TY_CHAR)
call salloc (decname, CQ_SZ_FNAME, TY_CHAR)
call salloc (funits, CQ_SZ_FUNITS, TY_CHAR)
call salloc (line, SZ_LINE, TY_CHAR)
# Open the catalog coordinate system.
if (cq_hinfo (cres, "csystem", Memc[csystem], SZ_FNAME) <= 0)
Memc[csystem] = EOS
if (Memc[csystem] == EOS || streq (Memc[csystem], "INDEF"))
call strcpy ("DEF_CATSYSTEM", Memc[csystem], SZ_FNAME)
# Open the query coordinate system data structure.
stat = sk_decwcs (Memc[csystem], mw, ccoo, NULL)
if (stat == ERR || mw != NULL) {
if (mw != NULL)
call mw_close (mw)
call sk_close (ccoo)
call sfree (sp)
return (NULL)
}
# Open the temporary results file.
call mktemp ("res", Memc[tmpname], SZ_FNAME)
fd = open (Memc[tmpname], NEW_FILE, TEXT_FILE)
# Write the file header to the temporary results file.
strfd = stropen (hdrtext, strlen(hdrtext), READ_ONLY)
call fprintf (fd, "# BEGIN CATALOG HEADER\n")
while (getline (strfd, Memc[line]) != EOF) {
call fprintf (fd, "# %s")
call pargstr (Memc[line])
}
call fprintf (fd, "# END CATALOG HEADER\n#\n")
call strclose (strfd)
# Determine the names of the ra and dec columns.
iferr (call at_stats (at, FIRA, Memc[raname], CQ_SZ_FNAME))
call strcpy ("ra", Memc[raname], CQ_SZ_FNAME)
iferr (call at_stats (at, FIDEC, Memc[decname], CQ_SZ_FNAME))
call strcpy ("dec", Memc[decname], CQ_SZ_FNAME)
# Determine the units of the ra and dec keywords.
call cq_funits (cres, Memc[raname], Memc[funits], CQ_SZ_QPUNITS)
units = strdic (Memc[funits], Memc[funits], CQ_SZ_FUNITS,
SKY_LNG_UNITLIST)
if (units > 0)
call sk_seti (ccoo, S_NLNGUNITS, units)
call cq_funits (cres, Memc[decname], Memc[funits], CQ_SZ_QPUNITS)
units = strdic (Memc[funits], Memc[funits], CQ_SZ_FUNITS,
SKY_LAT_UNITLIST)
if (units > 0)
call sk_seti (ccoo, S_NLATUNITS, units)
# Loop over the catalog records selecting those that match
# the region description.
do i = 1, cq_rstati (cres, CQRNRECS) {
# Decode the coordinates.
if (cq_gvald (cres, i, Memc[raname], ra) <= 0)
next
if (cq_gvald (cres, i, Memc[decname], dec) <= 0)
next
# Determine the coordinate units.
switch (sk_stati(ccoo, S_NLNGUNITS)) {
case SKY_HOURS:
ra = 15.0d0 * ra
case SKY_DEGREES:
;
case SKY_RADIANS:
ra = DRADTODEG(ra)
default:
;
}
switch (sk_stati(ccoo, S_NLATUNITS)) {
case SKY_HOURS:
dec = 15.0d0 * dec
case SKY_DEGREES:
;
case SKY_RADIANS:
dec = DRADTODEG(dec)
default:
;
}
# Test the limits
if (dec < dec1 || dec > dec2)
next
if (ra1 < ra2) {
if (ra < ra1 || ra > ra2)
next
} else {
if (ra > ra2 && ra < ra1)
next
}
# Check the longitude coordinate distance to remove pathologies
# in longitude or latitude strips involving the pole. This is
# an extra test of my own.
if (ra1 < ra2) {
dist = abs (ra - rac)
} else {
if (ra > ra1)
tra = ra - 360.0d0
else
tra = ra
if (rac > ra1)
trac = rac - 360.0d0
else
trac = rac
dist = abs (tra - trac)
}
if (abs (2.0d0 * dist *cos(DEGTORAD(dec))) > rawidth)
next
# Record has been selected.
if (cq_grecord (cres, Memc[line], SZ_LINE, i) <= 0)
next
call putline (fd, Memc[line])
}
# Close the tmeporary file.
call close (fd)
# Query the temporary file and then delete it.
res = cq_fquery (cq, Memc[tmpname], hdrtext)
call delete (Memc[tmpname])
# Clean up.
call sfree (sp)
return (res)
end
|
