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 335 336 337 338 339 340 341 342 343 344 345 346 347
|
program f77iterate_c
C
C This example program illustrates how to use the CFITSIO iterator function.
C
C This program creates a 2D histogram of the X and Y columns of an event
C list. The 'main' routine just creates the empty new image, then executes
C the 'writehisto' work function by calling the CFITSIO iterator function.
C
C 'writehisto' opens the FITS event list that contains the X and Y columns.
C It then calls a second work function, calchisto, (by recursively calling
C the CFITSIO iterator function) which actually computes the 2D histogram.
C external work function to be passed to the iterator
external writehisto
integer ncols
parameter (ncols=1)
integer units(ncols), colnum(ncols), datatype(ncols)
integer iotype(ncols), offset, n_per_loop, status
character*70 colname(ncols)
integer naxes(2), ounit, blocksize
character*80 fname
logical exists
C include f77.inc -------------------------------------
C Codes for FITS extension types
integer IMAGE_HDU, ASCII_TBL, BINARY_TBL
parameter (
& IMAGE_HDU = 0,
& ASCII_TBL = 1,
& BINARY_TBL = 2 )
C Codes for FITS table data types
integer TBIT,TBYTE,TLOGICAL,TSTRING,TSHORT,TINT
integer TFLOAT,TDOUBLE,TCOMPLEX,TDBLCOMPLEX
parameter (
& TBIT = 1,
& TBYTE = 11,
& TLOGICAL = 14,
& TSTRING = 16,
& TSHORT = 21,
& TINT = 31,
& TFLOAT = 42,
& TDOUBLE = 82,
& TCOMPLEX = 83,
& TDBLCOMPLEX = 163 )
C Codes for iterator column types
integer InputCol, InputOutputCol, OutputCol
parameter (
& InputCol = 0,
& InputOutputCol = 1,
& OutputCol = 2 )
C End of f77.inc -------------------------------------
C**********************************************************************
C Need to make these variables available to the 2 work functions
integer xsize,ysize,xbinsize,ybinsize
common /histcomm/ xsize,ysize,xbinsize,ybinsize
C**********************************************************************
status = 0
xsize = 480
ysize = 480
xbinsize = 32
ybinsize = 32
fname = 'histoimg.fit'
ounit = 15
C delete previous version of the file if it exists
inquire(file=fname,exist=exists)
if( exists ) then
open(ounit,file=fname,status='old')
close(ounit,status='delete')
endif
99 blocksize = 2880
C create new output image
call ftinit(ounit,fname,blocksize,status)
naxes(1) = xsize
naxes(2) = ysize
C create primary HDU
call ftiimg(ounit,32,2,naxes,status)
units(1) = ounit
C Define column as TINT and Output
datatype(1) = TINT
iotype(1) = OutputCol
C force whole array to be passed at one time
n_per_loop = -1
offset = 0
C execute the function to create and write the 2D histogram
print *,'Calling writehisto iterator work function... ',status
call ftiter( ncols, units, colnum, colname, datatype, iotype,
& offset, n_per_loop, writehisto, 0, status )
call ftclos(ounit, status)
C print out error messages if problem
if (status.ne.0) then
call ftrprt('STDERR', status)
else
print *,'Program completed successfully.'
endif
stop
end
C--------------------------------------------------------------------------
C
C Sample iterator function.
C
C Iterator work function that writes out the 2D histogram.
C The histogram values are calculated by another work function, calchisto.
C
C--------------------------------------------------------------------------
subroutine writehisto(totaln, offset, firstn, nvalues, narrays,
& units_out, colnum_out, datatype_out, iotype_out, repeat,
& status, userData, histogram )
integer totaln,offset,firstn,nvalues,narrays,status
integer units_out(narrays),colnum_out(narrays)
integer datatype_out(narrays),iotype_out(narrays)
integer repeat(narrays)
integer histogram(*), userData
external calchisto
integer ncols
parameter (ncols=2)
integer units(ncols), colnum(ncols), datatype(ncols)
integer iotype(ncols), rowoffset, rows_per_loop
character*70 colname(ncols)
integer iunit, blocksize
character*80 fname
C include f77.inc -------------------------------------
C Codes for FITS extension types
integer IMAGE_HDU, ASCII_TBL, BINARY_TBL
parameter (
& IMAGE_HDU = 0,
& ASCII_TBL = 1,
& BINARY_TBL = 2 )
C Codes for FITS table data types
integer TBIT,TBYTE,TLOGICAL,TSTRING,TSHORT,TINT
integer TFLOAT,TDOUBLE,TCOMPLEX,TDBLCOMPLEX
parameter (
& TBIT = 1,
& TBYTE = 11,
& TLOGICAL = 14,
& TSTRING = 16,
& TSHORT = 21,
& TINT = 31,
& TFLOAT = 42,
& TDOUBLE = 82,
& TCOMPLEX = 83,
& TDBLCOMPLEX = 163 )
C Codes for iterator column types
integer InputCol, InputOutputCol, OutputCol
parameter (
& InputCol = 0,
& InputOutputCol = 1,
& OutputCol = 2 )
C End of f77.inc -------------------------------------
C**********************************************************************
C Need to make these variables available to the 2 work functions
integer xsize,ysize,xbinsize,ybinsize
common /histcomm/ xsize,ysize,xbinsize,ybinsize
C**********************************************************************
if (status .ne. 0) return
C name of FITS table
fname = 'iter_c.fit'
iunit = 16
C do sanity checking of input values
if (totaln .ne. nvalues) then
C whole image must be passed at one time
status = -1
return
endif
if (narrays .ne. 1) then
C number of images is incorrect
status = -2
return
endif
if (datatype_out(1) .ne. TINT) then
C input array has wrong data type
status = -3
return
endif
C open the file and move to the table containing the X and Y columns
call ftopen(iunit,fname,0,blocksize,status)
call ftmnhd(iunit, BINARY_TBL, 'EVENTS', 0, status)
if (status) return
C both the columns are in the same FITS file
units(1) = iunit
units(2) = iunit
C desired datatype for each column: TINT
datatype(1) = TINT
datatype(2) = TINT
C names of the columns
colname(1) = 'X'
colname(2) = 'Y'
C leave column numbers undefined
colnum(1) = 0
colnum(2) = 0
C define whether columns are input, input/output, or output only
C Both input
iotype(1) = InputCol
iotype(1) = InputCol
C take default number of rows per iteration
rows_per_loop = 0
rowoffset = 0
C calculate the histogram
print *,'Calling calchisto iterator work function... ', status
call ftiter( ncols, units, colnum, colname, datatype, iotype,
& rowoffset, rows_per_loop, calchisto, histogram, status )
call ftclos(iunit,status)
return
end
C--------------------------------------------------------------------------
C
C Iterator work function that calculates values for the 2D histogram.
C
C--------------------------------------------------------------------------
subroutine calchisto(totalrows, offset, firstrow, nrows, ncols,
& units, colnum, datatype, iotype, repeat, status,
& histogram, xcol, ycol )
integer totalrows,offset,firstrow,nrows,ncols,status
integer units(ncols),colnum(ncols),datatype(ncols)
integer iotype(ncols),repeat(ncols)
integer histogram(*),xcol(*),ycol(*)
C include f77.inc -------------------------------------
C Codes for FITS extension types
integer IMAGE_HDU, ASCII_TBL, BINARY_TBL
parameter (
& IMAGE_HDU = 0,
& ASCII_TBL = 1,
& BINARY_TBL = 2 )
C Codes for FITS table data types
integer TBIT,TBYTE,TLOGICAL,TSTRING,TSHORT,TINT
integer TFLOAT,TDOUBLE,TCOMPLEX,TDBLCOMPLEX
parameter (
& TBIT = 1,
& TBYTE = 11,
& TLOGICAL = 14,
& TSTRING = 16,
& TSHORT = 21,
& TINT = 31,
& TFLOAT = 42,
& TDOUBLE = 82,
& TCOMPLEX = 83,
& TDBLCOMPLEX = 163 )
C Codes for iterator column types
integer InputCol, InputOutputCol, OutputCol
parameter (
& InputCol = 0,
& InputOutputCol = 1,
& OutputCol = 2 )
C End of f77.inc -------------------------------------
integer ii, ihisto, xbin, ybin
C**********************************************************************
C Need to make these variables available to the 2 work functions
integer xsize,ysize,xbinsize,ybinsize
common /histcomm/ xsize,ysize,xbinsize,ybinsize
C**********************************************************************
if (status .ne. 0) return
C --------------------------------------------------------
C Initialization procedures: execute on the first call
C --------------------------------------------------------
if (firstrow .eq. 1) then
C do sanity checking of input values
if (ncols .ne. 2) then
C number of arrays is incorrect
status = -4
return
endif
if (datatype(1).ne.TINT .or. datatype(2).ne.TINT) then
C wrong datatypes
status = -5
return
endif
C initialize the histogram image pixels = 0, including null value
do 10 ii = 1, xsize * ysize + 1
histogram(ii) = 0
10 continue
endif
C ------------------------------------------------------------------
C Main loop: increment the 2D histogram at position of each event
C ------------------------------------------------------------------
do 20 ii=2,nrows+1
xbin = xcol(ii) / xbinsize
ybin = ycol(ii) / ybinsize
ihisto = ( ybin * xsize ) + xbin + 2
histogram(ihisto) = histogram(ihisto) + 1
20 continue
return
end
|