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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <syserr.h>
include <imhdr.h>
include <imio.h>
include <mach.h>
include "fxf.h"
# FXF_RHEADER -- Read a FITS header into the image descriptor and the
# internal FITS descriptor.
procedure fxf_rheader (im, group, acmode)
pointer im #I image descriptor
int group #I group number to read
int acmode #I access mode
long pixoff, mtime
pointer sp, fit, lbuf, poff
int compress, devblksz, i, impixtype
bool bfloat, lscale, lzero
bool fxf_fpl_equald()
int strncmp()
errchk fxf_rfitshdr, realloc, syserr, syserrs
begin
call smark (sp)
call salloc (lbuf, SZ_LINE, TY_CHAR)
fit = IM_KDES(im)
FIT_MAX(fit) = 0.0
FIT_MIN(fit) = 0.0
FIT_MTIME(fit) = 0.0
FIT_IM(fit) = im
FIT_OBJECT(fit) = EOS
IM_CLSIZE(im) = 0
# Read the header unit number 'group', setting the values of the
# reserved fields in the FIT descriptor saving it in the FITS cache.
call fxf_rfitshdr (im, group, poff)
IM_MIN(im) = FIT_MIN(fit)
IM_MAX(im) = FIT_MAX(fit)
IM_MTIME(im) = FIT_MTIME(fit)
call strcpy (FIT_OBJECT(fit), IM_TITLE(im), LEN_CARD)
# If there is no group specification in the filename, group is -1;
# new group number is in FIT_GROUP.
group = FIT_GROUP(fit)
IM_CLINDEX(im) = group
# Process the reserved keywords (set in the FIT descriptor) into the
# corresponding fields of the IMIO descriptor.
if (acmode != NEW_COPY) {
IM_NDIM(im) = FIT_NAXIS(fit) # IM_NDIM
do i = 1, IM_NDIM(im) { # IM_LEN
IM_LEN(im,i) = FIT_LENAXIS(fit,i)
if (IM_LEN(im,i) == 0) {
IM_NDIM(im) = 0
break
}
}
}
lscale = fxf_fpl_equald (1.0d0, FIT_BSCALE(fit), 1)
lzero = fxf_fpl_equald (0.0d0, FIT_BZERO(fit), 1)
# Determine if scaling is necessary.
bfloat = (!lscale || !lzero)
FIT_PIXTYPE(fit) = NULL
FIT_ZCNV(fit) = NO
switch (FIT_BITPIX(fit)) {
case 8:
FIT_PIXTYPE(fit) = TY_UBYTE
if (bfloat)
impixtype = TY_REAL
else
impixtype = TY_SHORT # convert from byte to short
FIT_ZCNV(fit) = YES
case 16:
FIT_PIXTYPE(fit) = TY_SHORT
if (bfloat) {
impixtype = TY_REAL
FIT_ZCNV(fit) = YES
} else
impixtype = TY_SHORT
if ((strncmp ("USHORT", FIT_DATATYPE(fit), 6) == 0) ||
(lscale && fxf_fpl_equald (32768.0d0, FIT_BZERO(fit),4))) {
impixtype = TY_USHORT
FIT_ZCNV(fit) = NO
}
case 32:
FIT_PIXTYPE(fit) = TY_INT
if (bfloat)
impixtype = TY_REAL
else
impixtype = TY_INT
case -32:
FIT_PIXTYPE(fit) = TY_REAL
impixtype = TY_REAL
case -64:
FIT_PIXTYPE(fit) = TY_DOUBLE
impixtype = TY_DOUBLE
default:
impixtype = ERR
}
IM_PIXTYPE(im) = impixtype
IM_NBPIX(im) = 0 # no. bad pixels
mtime = IM_MTIME(im)
if (IM_MAX(im) > IM_MIN(im))
IM_LIMTIME(im) = mtime + 1 # time max/min last updated
else
IM_LIMTIME(im) = mtime - 1 # Invalidate DATA(MIN,MAX)
IM_HISTORY(im) = EOS
# Call up IMIO to set up the remaining image header fields used to
# define the physical offsets of the pixels in the pixfile.
compress = YES # do not align image lines on blocks
devblksz = 1 # disable all alignment
pixoff = Memi[poff+group]
FIT_PIXOFF(fit) = pixoff
call imioff (im, pixoff, compress, devblksz)
call sfree (sp)
end
# FXF_FPL_EQUALD -- Compare 2 double precision quantities up to a precision
# given by a tolerance.
bool procedure fxf_fpl_equald (x, y, it)
double x, y #I Input quantities to be compare for equality
int it #I Tolerance factor of 10 to compare the values
int ex, ey
double x1, x2, normx, normy, tol
begin
# Check for the obvious first.
if (x == y)
return (true)
# We can't normalize zero, so handle the zero operand cases first.
# Note that the case 0 equals 0 is handled above.
if (x == 0.0D0 || y == 0.0D0)
return (false)
# Normalize operands and do an epsilon compare.
call fp_normd (x, normx, ex)
call fp_normd (y, normy, ey)
if (ex != ey)
return (false)
else {
tol = EPSILOND * 10.0D0 * it
x1 = 1.0D0 + abs (normx - normy)
x2 = 1.0D0 + tol
return (x1 <= x2)
}
end
|
