program speed
!
!*******************************************************************************
!
!! SPEED measures the speed of writing and reading FITS files with FITSIO.
!
!
!  Local Parameters:
!
!    ascrow =  number of rows in the ASCII table.
!    binrow =  number of rows in the binary table.
!    dimsiz =  number of image pixels to write or read at one time.
!    xsize  =  number of pixels in each row of the image.
!    ysize  =  number of rows in the I*2 image.
!
  implicit none
!
  integer, parameter :: dimsiz = 20000
!
  integer array(dimsiz)
  integer, parameter :: ascrow = 100000
  integer, parameter :: binrow = 1250000
  character ( len = 80 ) :: fname = 'speedff.fit'
  integer i
  integer status
  integer unit
  integer, parameter :: xsize = 3000
  integer, parameter :: ysize = 3000
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SPEED'
  write ( *, '(a)' ) '  Measure the speed of reading and writing FITS files.'

  status = 0
!
!  Initialize the values in the array.
!
  do i = 1, dimsiz
    array(i) = i
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Create a large FITS file to test FITSIO'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '                                                    ', &
    '  SIZE / ELAPSE= RATE'
  write ( *, '(a)' ) '                                                    ', &
    '  (Mb)   (sec)  (Mb/s)'
!
!  Delete the file if it already exists, so we can then recreate it.
!
  call dfile ( fname, status )
!
!  Get an unused logical unit number to use to open the FITS file.
!
  call ftgiou ( unit, status )
!
!  Create the new empty file to test raw write and read speed.
!
  call ftinit ( unit, fname, 1, status )
!
!  Write than read a file using raw Fortran writes and reads.
!
  call raw_write ( unit, array, xsize, ysize, status )
  call raw_read ( unit, array, xsize, ysize, status )
!
!  Delete the raw file.
!
  call ftdelt ( unit, status )
!
!  Create the new empty file FITS file.
!
  call ftinit ( unit, fname, 1, status )

  call image_write ( unit, array, xsize, ysize, dimsiz, status )
  call bin_write ( unit, array, binrow, dimsiz, status )
  call ascii_write ( unit, array, ascrow, dimsiz, status )

  call image_read ( unit, array, xsize, ysize, dimsiz, status )
  call bin_read ( unit, array, binrow, dimsiz, status )
  call ascii_read ( unit, array, ascrow, dimsiz, status )
!
!  Close the file.
!
  call ftclos ( unit, status )
!
!  Delete the file.
!
  call dfile ( fname, status )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SPEED'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine raw_write ( unit, array, xsize, ysize, status )
!
!*******************************************************************************
!
!! RAW_WRITE tests the speed of raw Fortran writes.
!
  implicit none
!
  integer array(720)
  integer bitpix
  real elapse
  integer ios
  integer j
  integer naxes(2)
  integer naxis
  integer nloop
  real rate
  real size
  integer status
  integer unit
  integer xsize
  integer ysize
!
!  Write a dummy header, to prevent errors later when deleting the file.
!
  bitpix = 16
  naxis = 0
  call ftphpr ( unit, .true., bitpix, naxis, naxes, 0, 1, .true., status )
!
!  Flush the header to disk.
!
  call ftflus ( unit, status )
!
!  Number of loops required to write the same amount of data as
!  in the FITS image, writing 2880 bytes at a time.
!
  nloop = ( xsize * ysize ) / 720
!
!  Mark the starting time.
!
  call mktime
!
!  Write the array to the FITS file.
!
  do j = 2, nloop+1

    write ( unit, rec = j, iostat = ios ) array

    if ( ios /= 0 ) then
      status = 106
      exit
    end if

  end do
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )
!
!  Size of file, in megabytes.
!
  size = real ( 2880 * nloop ) / 1000000.0
  rate = size / elapse

  write ( *, 1000 ) size,elapse,rate
1000  format(' Raw Fortran Writes (2880 bytes/record)              ', &
   '  ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine raw_read ( unit, array, xsize, ysize, status )
!
!*******************************************************************************
!
!! RAW_READ tests the speed of raw Fortran reads
!
  implicit none
!
  integer array(720)
  real elapse
  integer j
  integer nloop
  real rate
  real size
  integer status
  integer unit
  integer xsize
  integer ysize
!
!  Number of loops required to read the same amount of data as
!  in the FITS image, reading 2880 bytes at a time.
!
  nloop = ( xsize * ysize ) / 720
!
!  Mark the starting time.
!
  call mktime
!
!  Read the array from the FITS file.
!
  do j = 2, nloop+1

    read ( unit, rec = j, iostat = status ) array

    if ( status /= 0 ) then
      status = 108
      exit
    end if

  end do
!
!  Get the elapsed time, in seconds
!
  call gttime ( elapse )
!
!  Size of file, in megabytes.
!
  size = real ( 2880 * nloop ) / 1000000.0

  rate = size / elapse

  write(*,1000) size,elapse,rate
1000  format(' Raw Fortran Reads                                   ', &
   '  ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine image_write ( unit, array, xsize, ysize, dimsiz, status )
!
!*******************************************************************************
!
!! IMAGE_WRITE writes a FITS primary array containing a 2-D image.
!
  implicit none
!
  integer xsize
!
  integer array(xsize)
  integer bitpix
  integer dimsiz
  real elapse
  integer j
  integer naxes(2)
  integer naxis
  integer nsize
  real rate
  real size
  integer status
  integer unit
  integer ysize
!
!  Initialize parameters about the FITS image.
!
  bitpix = 32
  naxis = 2
  naxes(1) = xsize
  naxes(2) = ysize
!
!  Write the required header keywords.
!
  call ftphpr ( unit, .true., bitpix, naxis, naxes, 0, 1, .true., status )

  nsize = xsize * ysize
!
!  Mark the starting time.
!
  call mktime
!
!  Write the array to the FITS file.
!
  do j = 1, nsize, dimsiz
    call ftpprj ( unit, 1, j, dimsiz, array, status )
  end do
!
!  Flush the file to disk.
!
  call ftflus ( unit, status )
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )
!
!  Size of file, in megabytes.
!
  size = real ( 4 * xsize * ysize ) / 1000000.0
  rate = size / elapse

  write(*,1000) xsize,ysize,dimsiz,size,elapse,rate
1000  format(' Write ',i5,' x ',i4,' row I*4 image (',i6, &
   ' pixels/loop) ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine bin_write ( unit, array, nrows, dimsiz, status )
!
!*******************************************************************************
!
!! BIN_WRITE writes a binary table.
!
  implicit none
!
  integer dimsiz
!
  integer array(dimsiz)
  real elapse
  character ( len = 16 ) :: extnam = 'speed'
  integer frow
  integer ngroup
  integer nrows
  integer ntodo
  real rate
  integer remain
  real size
  integer status
  character ( len = 16 ), dimension ( 2 ) :: tform = (/ '1J',   '1J'   /)
  character ( len = 16 ), dimension ( 2 ) :: ttype = (/ 'int1', 'int2' /)
  character ( len = 16 ), dimension ( 2 ) :: tunit = (/ ' ',    ' '    /)
  integer unit
  integer width
!
  width = 8
!
!  Write the required header parameters for the binary table.
!
  call ftibin ( unit, nrows, 2, ttype, tform, tunit, extnam, 0, status )

  call ftgrsz ( unit, ngroup, status )

  ngroup = min ( ngroup, dimsiz )
  remain = nrows
  frow = 1
!
!  Mark the starting time.
!
  call mktime

  do

    ntodo = min ( ngroup, remain )

    call ftpclj ( unit, 1, frow, 1, ntodo, array, status )
    call ftpclj ( unit, 2, frow, 1, ntodo, array, status)
    frow = frow + ntodo
    remain = remain - ntodo

    if ( remain <= 0 ) then
      exit
    end if

  end do
!
!  Flush the file to disk.
!
  call ftflus ( unit, status )
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )
!
!  Size of table in megabytes.
!
  size = real ( nrows * width ) / 1000000.0
  rate = size / elapse

  write(*,1000)nrows,ngroup,size,elapse,rate
1000  format(' Write ',i7,' row x 2 col BINTABLE  (',i5, &
   ' rows/loop) ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine ascii_write ( unit, array, nrows, dimsiz, status )
!
!*******************************************************************************
!
!! ASCII_WRITE writes an ASCII table.
!
  implicit none
!
  integer dimsiz
!
  integer array(dimsiz)
  real elapse
  character ( len = 16 ) :: extnam = 'aspeed'
  integer frow
  integer ngroup
  integer nrows
  integer ntodo
  integer remain
  real rate
  real size
  integer status
  integer, dimension(2) :: tbcol = (/ 1, 8 /)
  character ( len = 16 ), dimension(2) :: tform = (/ 'I6', 'I6' /)
  character ( len = 16 ), dimension(2) :: ttype = (/ 'int1', 'int2' /)
  character ( len = 16 ), dimension(2) :: tunit = (/ ' ', ' ' /)
  integer unit
  integer width
!
  width = 13
!
!  Write the required header parameters for the binary table.
!
  call ftitab ( unit, 13, nrows, 2, ttype, tbcol, tform, tunit, extnam, status )

  call ftgrsz ( unit, ngroup, status )

  ngroup = min ( ngroup, dimsiz )
  remain = nrows
  frow = 1
!
!  Mark the starting time.
!
  call mktime

  do

    ntodo = min ( ngroup, remain )

    call ftpclj ( unit, 1, frow, 1, ntodo, array, status )

    call ftpclj ( unit, 2, frow, 1, ntodo, array, status )

    frow = frow + ntodo
    remain = remain - ntodo

    if ( remain <= 0 ) then
      exit
    end if

  end do
!
!  Flush the file to disk.
!
  call ftflus ( unit, status )
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )
!
!  Size of table in megabytes.
!
  size = real ( nrows * width ) / 1000000.0
  rate = size / elapse

  write(*,1000)nrows,ngroup,size,elapse,rate
1000  format(' Write ',i7,' row x 2 col ASCII tab (',i5, &
   ' rows/loop) ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine image_read ( unit, array, xsize, ysize, dimsiz, status )
!
!*******************************************************************************
!
!! IMAGE_READ reads a FITS primary array containing a 2-D image.
!
  implicit none
!
  integer xsize
!
  logical anynul
  integer array(xsize)
  integer dimsiz
  real elapse
  integer hdutype
  integer j
  integer nsize
  integer nulval
  real rate
  real size
  integer status
  integer unit
  integer ysize
!
!  Move to the primary array.
!
  call ftmahd ( unit, 1, hdutype, status )

  nulval = 0
  nsize = xsize * ysize
!
!  Mark the starting time.
!
  call mktime
!
!  Read the array from the FITS file.
!
  do j = 1, nsize, dimsiz
    call ftgpvj ( unit, 1, j, dimsiz, nulval, array, anynul, status )
  end do
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )

  size = real ( 4 * xsize * ysize ) / 1000000.0
  rate = size / elapse

  write(*,1000)dimsiz,size,elapse,rate
1000  format(' Read back I*4 image  (',i6,' pixels/loop)      ', &
   '       ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine bin_read ( unit, array, nrows, dimsiz, status )
!
!*******************************************************************************
!
!! BIN_READ reads a binary table.
!
  implicit none
!
  integer dimsiz
!
  logical anynul
  integer array(dimsiz)
  real elapse
  integer frow
  integer hdutyp
  integer ngroup
  integer nrows
  integer ntodo
  integer nulval
  real rate
  integer remain
  real size
  integer status
  integer unit
  integer width
!
!  Move to the next extension.
!
  call ftmrhd ( unit, 1, hdutyp, status )

  call ftgrsz ( unit, ngroup, status )
  ngroup = min ( ngroup, dimsiz )
  remain = nrows
  frow = 1

  nulval = 0
  width = 8
!
!  Mark the starting time.
!
  call mktime

  do

    ntodo = min ( ngroup, remain )

    call ftgcvj ( unit, 1, frow, 1, ntodo, nulval, array, anynul, status )

    call ftgcvj ( unit, 2, frow, 1, ntodo, nulval, array, anynul, status )

    frow = frow + ntodo
    remain = remain - ntodo

    if ( remain <= 0 ) then
      exit
    end if

  end do
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )

  size = real ( nrows * width ) / 1000000.0
  rate = size / elapse

  write(*,1000)ngroup,size,elapse,rate
1000  format(' Read back Binary table  (',i5,' rows/loop)       ', &
   '      ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine ascii_read ( unit, array, nrows, dimsiz, status )
!
!*******************************************************************************
!
!! ASCII_READ reads an ASCII table.
!
  implicit none
!
  integer dimsiz
!
  logical anynul
  integer array(dimsiz)
  real elapse
  integer frow
  integer hdutyp
  integer ngroup
  integer nrows
  integer ntodo
  integer nulval
  real rate
  integer remain
  real size
  integer status
  integer unit
  integer width
!
!  Move to the next extension.
!
  call ftmrhd ( unit, 1, hdutyp, status )

  call ftgrsz ( unit, ngroup, status )
  ngroup = min ( ngroup, dimsiz )
  remain = nrows
  frow = 1

  nulval = 0
  width = 13
!
!  Mark the starting time.
!
  call mktime

  do

    ntodo = min ( ngroup, remain )
    call ftgcvj ( unit, 1, frow, 1, ntodo, nulval, array, anynul, status )

    call ftgcvj ( unit, 2, frow, 1, ntodo, nulval, array, anynul, status )

    frow = frow + ntodo
    remain = remain - ntodo

    if ( remain <= 0 ) then
      exit
    end if

  end do
!
!  Get the elapsed time, in seconds.
!
  call gttime ( elapse )

  size = real ( nrows * width ) / 1000000.0
  rate = size / elapse

  write(*,1000)ngroup,size,elapse,rate
1000  format(' Read back ASCII table   (',i5,' rows/loop)       ', &
   '      ',f5.2,1x,f6.2,1x,f6.2)
!
!  Check for any error.
!
  call perror ( status )

  return
end
subroutine perror ( status )
!
!*******************************************************************************
!
!! PERROR prints out the FITSIO error messages to the user.
!
  implicit none
!
  character ( len = 80 ) errmessage
  character ( len = 30 ) errtext
  integer status
!
!  Check if status is OK (no error); if so, simply return.
!
  if ( status <= 0 ) then
    return
  end if
!
!  Get the text string which describes the error.
!
  call ftgerr ( status, errtext )
  write ( *, '(a)' ) ' '
  write ( *, '(a,i9)' ) 'FITSIO Error Status =', status
  write ( *, '(a)' ) trim ( errtext )
!
!  Read and print out all the error messages on the FITSIO stack.
!
  call ftgmsg ( errmessage )

  do while ( len_trim ( errmessage ) /= 0 )
    write ( *, '(a)' ) trim ( errmessage )
    call ftgmsg ( errmessage )
  end do

  return
end
subroutine dfile ( filename, status )
!
!*******************************************************************************
!
!! DFILE deletes a FITS file.
!
  implicit none
!
  integer block
  character ( len = * )  filename
  integer status
  integer unit
!
!  Simply return if status is greater than zero.
!
  if ( status > 0 ) then
    return
  end if
!
!  Get an unused logical unit number to use to open the FITS file.
!
  call ftgiou ( unit, status )
!
!  Try to open the file, to see if it exists.
!
  call ftopen ( unit, filename, 1, block, status )

  if ( status == 0 ) then
!
!  File was opened, so now delete it.
!
    call ftdelt ( unit, status )

  else if ( status == 103 ) then
!
!  File doesn't exist, so just reset status to zero and clear errors.
!
    status = 0
    call ftcmsg

  else
!
!  There was some other error opening the file; delete the file anyway.
!
    status = 0
    call ftcmsg
    call ftdelt ( unit, status )

  end if
!
!  Free the unit number for later reuse.
!
  call ftfiou ( unit, status )

  return
end
subroutine mktime
!
!*******************************************************************************
!
!! MKTIME marks the start time.
!
  implicit none
!
  integer dd
  integer hh
  integer iarray(8)
  integer mm
  integer ss
  integer stemp
!
  common/times/hh,mm,ss,dd
!
!  Start recording the elapsed time on a seconds tick.
!  This provides more consistent measurements if
!  the times are only accurate to the nearest second.
!
  call date_and_time ( values = iarray )

  stemp = iarray(7)
!
!  Keep getting the time, until the seconds change.
!
  do

    call date_and_time ( values = iarray )

    if ( iarray(7) /= stemp ) then
      exit
    end if

  end do

  hh = iarray(5)
  mm = iarray(6)
  ss = iarray(7)
  dd = iarray(8)

  return
end
subroutine gttime ( elapse )
!
!*******************************************************************************
!
!! GTTIME gets the time that has elapsed since MKTIME was last called.
!
  implicit none
!
  integer dd
  real elapse
  real frac
  integer hh
  integer iarray(8)
  integer itemp
  integer jarray(8)
  integer mm
  integer nloop1
  integer nloop2
  integer ss
!
  common/times/hh,mm,ss,dd

  call date_and_time ( values = iarray )

  if ( dd == 0 .and. iarray(8) == 0 ) then
!
!  Clock is not reporting the milliseconds, so use an
!  alternate method to estimate the factional part of the seconds
!
    itemp = iarray(7)
    nloop1 = 0
!
!  Keep getting the time, until the seconds change.
!
    do

      call date_and_time ( values = jarray )
      nloop1 = nloop1 + 1

      if ( jarray(7) /= itemp ) then
        exit
      end if

    end do
!
!  Now count the number of loops in 1 second.
!
    itemp = jarray(7)
    nloop2 = 0
!
!  Keep getting the time, until the seconds change.
!
    do

      call date_and_time ( values = jarray )
      nloop2 = nloop2 + 1

      if ( jarray(7) /= itemp ) then
        exit
      end if

    end do

    frac = real ( nloop1 ) / real ( nloop2 )
    frac = min ( frac, 1.0 )

  else

    frac = real ( iarray(8) - dd ) / 1000.0

  end if

  elapse = real ( &
      ( iarray(5) - hh ) * 3600 &
    + ( iarray(6) - mm ) * 60 &
    + ( iarray(7) - ss ) ) &
    + frac

  return
end