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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
|
include <error.h>
# MAKEWCS -- Task to write the WCS information into a header of an image
# subset from the scanned plate catalog.
procedure t_makewcs()
int imtopenp(), imtgetim()
int list
char image[SZ_FNAME]
bool verbose, clgetb()
begin
list = imtopenp ("infiles")
verbose = clgetb("verbose")
while (imtgetim(list, image, SZ_FNAME) != EOF) {
call makecd(image)
if (verbose) {
call printf (" WCS for %s\n")
call pargstr(image)
call flush(STDOUT)
}
}
end
include <imhdr.h>
include <math.h>
define SZ_KEYWORD 8
# MAKECD -- Procedure to calculate the WCS values for a subset of a
# scanned plate that contains a general plate solution.
#
# Calculates CRPIX, CRVAL and CD matrix from the information in the GSSS
# header wich consist on the positions of the lower left corner of the
# subimage (CNPIX1, CNPIX2) w/r to the lower left corner of the scanned plate.
# From here create a grid of 100x100 points and do a local solution to
# the CD values and translate the reference position to the center of the
# subimage.
#
procedure makecd(image)
char image[SZ_FNAME]
double amdx[20], amdy[20] # Plate Soln.coeff.
double plate_cen_x, plate_cen_y # center position (microns)
double plate_cen_ra, plate_cen_dec # in radians
double x_pixel_size, y_pixel_size # step size (microns)
double plate_scale # arcsec/mm
int ra_h, ra_m # RA hours, minutes
double ra_s # RA seconds
char dec_sign[1] #
int dec_d, dec_m # DEC degrees, minutes
double dec_s # DEC seconds
double object_mag, object_col # Magnitud, color for object
double im_x_center_pix, im_y_center_pix # Pixel position of image center
# w/r to the original scanned plate
int xcorner, ycorner # Lower left position of image w/r to
# original scanned plate
int xsize, ysize # Naxis1, naxis2
real crpix1, crpix2
double crval1, crval2
double cdmatx[4]
double imgetd()
real imgetr()
int imgeti()
char parname[8]
int i
pointer im
int immap(), imaccf()
errchk imgetd, imgeti
begin
im = immap (image, READ_WRITE, 0)
# See if image header has the General Plate solution.
if (imaccf(im,"PPO3 ") == NO)
call error(0,"No Plate solution available in header")
# If we have an old DSS image (i.e. the one with the CRPIX), rename
# this keyword to CNPIX and proceed.
if (imaccf(im,"CRPIX1") == YES ) {
if (imaccf(im,"CRVAL1") == YES) {
if (imaccf(im,"CD1_1") == NO ){
# This is the case when we have CRPIX1, CRVAL1 and no CD1_1
# so, proceed to calculate the WCS again.
crpix1 = imgetr (im, "CRPIX1")
call imdelf (im, "CRPIX1")
call imaddr (im, "CNPIX1", crpix1)
crpix2 = imgetr (im, "CRPIX2")
call imdelf (im, "CRPIX2")
call imaddr (im, "CNPIX2", crpix2)
}
} else {
crpix1 = imgetr (im, "CRPIX1")
call imdelf (im, "CRPIX1")
call imaddr (im, "CNPIX1", crpix1)
crpix2 = imgetr (im, "CRPIX2")
call imdelf (im, "CRPIX2")
call imaddr (im, "CNPIX2", crpix2)
}
}
if (imaccf(im,"CNPIX1") == NO )
call error(0,"Lower left position CNPIX(1,2) not found")
plate_cen_x = imgetd(im, "PPO3 ")
plate_cen_y = imgetd(im, "PPO6 ")
x_pixel_size = imgetd(im, "XPIXELSZ")
y_pixel_size = imgetd(im, "YPIXELSZ")
plate_scale = imgetd(im, "PLTSCALE")
ra_h = imgeti (im, "PLTRAH ")
ra_m = imgeti (im, "PLTRAM ")
ra_s = imgetd (im, "PLTRAS ")
call imgstr(im, "PLTDECSN", dec_sign, 1)
dec_d = imgeti (im, "PLTDECD ")
dec_m = imgeti (im, "PLTDECM ")
dec_s = imgetd (im, "PLTDECS ")
plate_cen_ra = DEGTORAD ((ra_h + ra_m/60.0d0 + ra_s/3600.0)*15.0d0)
plate_cen_dec = DEGTORAD (dec_d + dec_m/60.0d0 + dec_s/3600.0d0)
if (dec_sign[1] == '-')
plate_cen_dec = -plate_cen_dec
# Get general plate solution coefficients
do i = 1, 20 {
call sprintf(parname, SZ_LINE, "AMDX%d")
call pargi(i)
amdx[i] = imgetd (im, parname)
}
do i = 1, 20 {
call sprintf(parname, SZ_LINE, "AMDY%d")
call pargi(i)
amdy[i] = imgetd (im, parname)
}
if (imaccf(im, "CNPIX1") == YES) {
xcorner = imgetr (im, "CNPIX1")
ycorner = imgetr (im, "CNPIX2")
} else
call error(0,"CNPIXs keyword found")
xsize = IM_LEN(im,1)
ysize = IM_LEN(im,2)
crpix1 = xsize/2.
crpix2 = ysize/2.
# Center of image w/r to original lower left corner of
# scanned plate
im_x_center_pix = xcorner + (xsize/2.) - 0.5
im_y_center_pix = ycorner + (ysize/2.) - 0.5
# Calculate equatorial coordinates for the center of subset giving
# the complete plate solution w/r to the original lower left corner
call ccgseq (plate_cen_ra,
plate_cen_dec,
plate_cen_x,
plate_cen_y,
x_pixel_size,
y_pixel_size,
plate_scale,
amdx,
amdy,
im_x_center_pix,
im_y_center_pix,
object_mag,
object_col,
crval1,
crval2)
# calculate CD matrix values for the input subset from the original
# plate solution.
call calcds (plate_cen_ra,
plate_cen_dec,
plate_cen_x,
plate_cen_y,
xcorner,
ycorner,
x_pixel_size,
y_pixel_size,
plate_scale,
xsize,
ysize,
crval1,
crval2,
amdx,
amdy,
cdmatx)
# Now update the image header
iferr (call addkf (im, "CTYPE1", 8, TY_CHAR,
"r.a. in tangent plane proyection"))
;
call impstr (im, "CTYPE1", "RA---TAN")
iferr (call addkf (im, "CTYPE2", 8, TY_CHAR,
"dec. in tangent plane proyection"))
;
call impstr (im, "CTYPE2", "DEC--TAN")
iferr (call addkf (im, "CRPIX1", 1, TY_REAL,
"reference pixel of first axis"))
;
call imputr (im, "CRPIX1", crpix1)
iferr (call addkf (im, "CRPIX2", 1, TY_REAL,
"reference pixel of second axis"))
;
call imputr (im, "CRPIX2", crpix2)
iferr (call addkf (im, "CRVAL1", 1, TY_REAL,
"r.a. at pixel crpix1 in decimal degrees"))
;
call imputd (im, "CRVAL1", RADTODEG(crval1))
iferr (call addkf (im, "CRVAL2", 1, TY_REAL,
"dec. at pixel crpix2 in decimal degrees"))
;
call imputd (im, "CRVAL2", RADTODEG(crval2))
# now add the new CD values
iferr (call addkf (im, "CD1_1", 1, TY_REAL,
"Change in r.a. per pixel along the first axis"))
;
else
call addkf (im, "", 1, TY_REAL,
"evaluated at the reference pixel")
call imputd (im, "CD1_1", cdmatx[1])
iferr (call addkf (im, "CD1_2", 1, TY_REAL,
"ditto along the second axis"))
;
call imputd (im, "CD1_2", cdmatx[2])
iferr (call addkf (im, "CD2_1", 1, TY_REAL,
"Change in dec. per pixel along the first axis"))
;
else
call addkf (im, "", 1, TY_REAL,
"evaluated at the reference pixel")
call imputd (im, "CD2_1", cdmatx[3])
iferr (call addkf (im, "CD2_2", 1, TY_REAL,
"ditto along the second axis"))
;
call imputd (im, "CD2_2", cdmatx[4])
call imunmap (im)
end
include <syserr.h>
include <imhdr.h>
include <imio.h>
include <mach.h>
define SZ_VALSTR 80
# ADDKF -- Add a user field to the image header. It is an error if the named
# field already exists.
procedure addkf (im, key, lenv, datatype, comment)
pointer im # image descriptor
char key[SZ_KEYWORD] # name of the new parameter
int lenv # keyword value length
int datatype # datatype of parameter
char comment[SZ_VALSTR] # comment describing new parameter
int max_lenuserarea, diff, slen
pointer sp, keyname, rp, ua, op
int idb_kwlookup(), idb_findrecord(), strlen()
errchk syserrs
begin
call smark (sp)
call salloc (keyname, SZ_FNAME, TY_CHAR)
call strcpy (key, Memc[keyname], SZ_FNAME)
# Check for a redefinition.
if ((idb_kwlookup (key) > 0) || (idb_findrecord (im, key, rp) > 0))
call syserrs (SYS_IDBREDEF, key)
# Open the user area string for appending. If the user area is not
# empty the last character must be the newline record delimiter,
# else the new record we add will be invalid.
max_lenuserarea = (LEN_IMDES + IM_LENHDRMEM(im) - IMU + 1) * SZ_STRUCT
ua = IM_USERAREA(im)
for (rp=ua; Memc[rp] != EOS; rp=rp+1)
;
if (rp - ua + SZ_VALSTR + 1 >= max_lenuserarea)
call syserrs (SYS_IDBOVFL, key)
if (rp > ua && Memc[rp-1] != '\n') {
Memc[rp] = '\n'
rp = rp + 1
}
# Append the new record with an uninitialized value field. Keyword
# value pairs are encoded in FITS format.
do op = rp, rp + SZ_VALSTR # blank fill card
Memc[op] = ' '
# Add the "= 'value' / comment".
slen = strlen(Memc[keyname])
call amovc (Memc[keyname], Memc[rp], slen)
if (slen != 0)
Memc[rp+9-1] = '='
if (datatype == TY_CHAR) {
Memc[rp+11-1] = '\''
Memc[rp+20-1] = '\''
}
# Add the comment field.
diff = lenv - 20
if (diff > 0) {
Memc[rp+34-1+diff] = '/'
call amovc (comment, Memc[rp+36-1+diff],
min (SZ_VALSTR-36+1-diff, strlen(comment)))
} else {
Memc[rp+32-1] = '/'
call amovc (comment, Memc[rp+34-1],
min (SZ_VALSTR-34+1, strlen(comment)))
}
# Terminate the card.
Memc[rp+SZ_VALSTR] = '\n'
Memc[rp+SZ_VALSTR+1] = EOS
call sfree (sp)
end
|
