Initial commit

1 files changed, 347 insertions, 0 deletions

diff --git a/pkg/tbtables/cfitsio/iter_c.f b/pkg/tbtables/cfitsio/iter_c.f
new file mode 100644
index 00000000..f9abeaa8
--- /dev/null
+++ b/pkg/tbtables/cfitsio/iter_c.f
@@ -0,0 +1,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
+