     program opphot                         

        ! Carries out optimal extraction of photometry using the method
        ! described by Naylor (1998).  You need to run the code twice.
        ! The first pass will produce accurate positions of the stars
        ! (in stars.pos_new), which should then be fed into the second
        ! run.

        ! By Tim Naylor.

        ! Version 1.1
        !   Now fails elegantly when psf.pos or stars.pos is incorrectly
        !   formatted.
        ! Version 2.0
        !   The position angle of the Gaussian PSF allowed to run free.
        ! Version 2.1
        !   1) Small bug cured which referred to a non-existent element of
        !   dpsf when writing a comment to an output file, and to a 
        !   non-existant element of dpos when writing to the screen.
        !   2) comp reset to .false. for each new star.
        ! Version 2.2
        !   Changed format of .log file.
        ! Version 2.3
        !   More comments to aid porting, and default values of aperture and
        !   mask radius set using fwhm from first image.
        ! Version 2.4
        !   Added errors in position into .log file.
        ! Version 2.5
        !   Increased to 99,999 output file names for Monica
        ! Version 2.6
        !   Log file changed to imagename.opt file, and imagename.dat added.
        !   These are files used for colour-magnitude diagrams, and hence
        !   if there are more than 100 stars, the star????? files are not
        !   written.
        ! Version 2.7
        !   Minor changes so that it compiles with the NAG Fortran 95
        !   compiler, as well as with the DEC Fortran 90 compiler.
        ! Version 2.8
        !   1) Allows user to skip aperture photometry.
        !   2) Adaptively sets the clipping radius to be a mutliple of the 
        !      FWHM of the PSF star in each frame.
        !   3) Allows the user to skip the optimal photometry.
        !   4) Checks for a high skewness in the sky histogram.
        !   5) Signals problem measurements only through star%col(1)%flg, never by
        !      setting the error negative.
        ! Version 3.0
        !   General tidy-up of code, including changing many subroutine arguments.
        !   Introduction of pixel flags.
        !   Flags changed to characters rather than integers.

        ! This module contains all the ARK I/O routines.
        use ark_file_io
        ! This module contains the optimal extraction stuff.
        use opt_extr
        ! This contains the definition of the star structure.
        use define_star

        implicit none

        ! For the list of PSF stars.
        character(len=50) :: psfpos
        integer :: npsf, ipsf, nfit, start_ext, end_ext
        type(a_star), allocatable, dimension(:) :: psf
        real, allocatable, dimension(:) :: dpsf
        logical :: cluster_cat
                          
        ! For the list of target stars.
        character(len=50) :: strpos
        integer :: nstar, non
        real, allocatable, dimension(:) :: dpos
        type(a_star), allocatable, dimension(:) :: star
        integer :: nonstellar
        ! The output file for these stars.
	character(len=50), allocatable, dimension(:) :: ofopt, ofap
        real, allocatable, dimension(:) :: xpos_new, ypos_new
        integer, allocatable, dimension(:) :: npos_new
        character(len=80) :: data_str

        ! The positions of stars shifted from the current reference frame
        ! into the new one.
        real, allocatable, dimension(:) :: xpos0, ypos0

        integer :: i, j, istar
        ! integer :: jstar
        ! An answer string
        character(len=50) :: answer           
        logical :: there
        integer :: iunit1, iostat, iunit2, iunit4, iunit5
        real :: peak

        ! For the image.
        integer :: naxis(2)
        real, allocatable, dimension(:,:) :: data, axdata
        ! The active area of the detector.
        integer :: low(2), high(2)
        ! The pixel flag array.
        character, dimension(:,:), allocatable :: pix_flg

        ! Is there a companion, and what is its position?
        logical :: comp
        real :: xcomp, ycomp
                                        
        ! For skyfit.
        real :: skycnt, skyerr, skynos
        integer :: ibox
        logical :: sky_stuffed
        character :: cflag

        ! For the fitting routine.
        ! The first 3 parameters, which define the shape of the function.
        real, dimension(3) :: shape_par

        ! The gain.
        real :: adu
        ! The number of the star whose signal-to-noise is to be optimised.
        integer :: iopt
        ! And its peak counts.
        real :: optnrm

        ! For offsets.dat
        character(len=50) :: offset, ifname
        character(len=3) :: cext
        integer :: iext
        real :: sky, theoff, fwhm
        real, dimension(6) :: coff
        integer :: ntoe
        real :: work
                              
        ! For times.dat
        character(len=50) :: timfil, ctime
        real, allocatable, dimension(:) :: time, dtime
        character(len=50), allocatable, dimension(:) :: timnam
        integer ntim
        ! The times for the current frame.
        real :: timec, dtimec

        ! The shifted position of psf and target star
        real :: xpos, ypos
        ! And its error.
        real :: xerr, yerr
        ! The ouput optimal and aperture photometry stars.
        type(a_star) :: opt, ap

        ! For clipping the edge of the weight region.
        real :: cliprad
        ! How many times the FWHM cliprad is set.
        real :: clip_fwhm

        ! The aperture for aperture photometry.
        real :: aprad

        ! The distance to the nearest bad pixel.
        real :: badpix
                                   
        ! The log files for the optimal and aperture photometry.
        character(len=50) :: ifroot
        integer :: cluster_opt, cluster_ap, isky
        character(len=80), dimension(3) :: header

        real :: bad_sky_chi, bad_sky_skw

        integer :: ncol
        character(len=10), dimension(mcol) :: colstr
                            
        ! Variables for debugging                     
	logical :: debug
        common /bug/ debug

        ! ***************************************************************
        ! Actual program starts here.
        ! ***************************************************************
                     
        call setbug()

        if (debug) call opt_debug_on()


        ! ***************************************************************
        ! Read the co-ordinates for each psf star.
        ! ***************************************************************

        print*, ' '
        print*, '* First this program requires a list of possible PSF'
        print*, '* stars.  This should be a file where each line has'
        print*, '* a star number, followed by its X and Y co-ordinates'
        print*, '* in the reference data frame, and the radius to be'
        print*, '* searched for the center of the star.  If this last'
        print*, '* number is zero the star position is taken as fixed.'
        print*, '* (The first 3 lines of the file are ignored.)'
        print*, '* (Type end if you don''t want optimal photometry.)'

        psffil: do 
          print*, '> Give the name of this file (<cr>=psf.cat).'
          psfpos=' '
          read (*, '(a)') psfpos
          if (psfpos=='end' .or. psfpos=='END') then
            clip_fwhm = -1.0
            exit psffil
          else
            clip_fwhm=1.0
          end if
          if (psfpos .eq. ' ') psfpos='psf.cat'
          inquire(file=psfpos, exist=there)
          if (.not. there) then
            print*, '* File not found.'
          else
           exit psffil
          end if
        end do psffil
        cluster_cat=.false.
        if (index(psfpos, '.cat') > 0) then
          print*, '* Assuming this is a cluster format catalogue.'
          cluster_cat = .true.
        end if

        npsf=0
        if (clip_fwhm > 0.0) then

          if (.not. cluster_cat) then
            ! Open the file with the position(s) of the psf star(s) using a 
            ! machine independent ARK routine.  It returns the number of the
            ! unit the file is opened on.
            call arkopn(iunit1, psfpos, ' ', 'OLD', 'READONLY', &
            'FORMATTED', ' ', 0)
            read(iunit1,*,iostat=iostat)
            read(iunit1,*,iostat=iostat)
            read(iunit1,*,iostat=iostat)
            ctpsf: do
              read(iunit1, '(a)', iostat=iostat) psfpos
              if (iostat < 0) exit ctpsf
              npsf=npsf+1
            end do ctpsf
            allocate(psf(npsf))
            call zero_star(psf)
            allocate(dpsf(npsf))
            rewind (iunit1)
            read(iunit1,*,iostat=iostat)
            read(iunit1,*,iostat=iostat)
            read(iunit1,*,iostat=iostat)
            psfin: do i=1, npsf
              read(iunit1,*,iostat=iostat) psf(i)%id, psf(i)%x, psf(i)%y, dpsf(i)
              psf(i)%col(1)%flg='OO'
              if (iostat /= 0) then
                print*, 'Error ', iostat, ' in psf star position file at line ', i+3
                stop
              end if
            end do psfin
            close(iunit1)
          else
            call nxtcls_in(psfpos)
            iostat=read_cluster_file(npsf, ncol, colstr, psf) 
          end if
                                                  
          print*, '* ', npsf, ' psf star positions read.'

          if (cluster_cat) then 
            allocate(dpsf(npsf))
            print*, '> Give radius to be searched for each star (<0 to fix position).'
            read(*,*) work
            dpsf=work
          end if

        end if

        ! *******************************************************************
        ! Read the co-ordinates for each target star, and make a file for it.
        ! *******************************************************************

        print*, ' '
        print*, '* Next the program requires a list of the target'
        print*, '* stars.  This should be a file where each line has'
        print*, '* a star number, followed by its X and Y co-ordinates'
        print*, '* in the reference data frame, and the radius to be'
        print*, '* searched for the center of the star.  If this last'
        print*, '* number is zero the star position is taken as fixed.'
        print*, '* (The first 3 lines of the file are ignored.)'

        strfil: do
          print*, '> Give the name of this file (<cr>=stars.cat).'
          strpos=' '
          read (*, '(a)') strpos
          if (strpos == ' ') strpos='stars.cat'
          inquire(file=strpos, exist=there)
          if (.not. there) then
            print*, '* File not found.'
          else
           exit strfil
          end if
        end do strfil
        cluster_cat=.false.
        if (index(strpos, '.cat') > 0) then
          print*, '* Assuming this is a cluster format catalogue.'
          cluster_cat = .true.
        end if

        ! Open the file with the position(s) of the star(s) using a 
        ! machine independent ARK routine.  It returns the number of the
        ! unit the file is opened on.
        if (.not. cluster_cat) then
          call arkopn(iunit1, strpos, ' ', 'OLD', 'READONLY', &
          'FORMATTED', ' ', 0)
          ! Read three lines of header.
          read(iunit1, '(/,/)')
          ! And count the number of stars.
          nstar=0
          ctstr: do   
            read(iunit1, '(a)', iostat=iostat) strpos
            if (iostat < 0) exit ctstr
            nstar=nstar+1
          end do ctstr
          allocate(star(nstar))
          call zero_star(star)
          allocate(dpos(nstar))
          ! Now check to see if there's a non-stellar flag.
          rewind (iunit1)
          read(iunit1, '(/,/)')
          read(iunit1, '(a80)') data_str
          read(data_str,*, iostat=non) star(1)%id, star(1)%x, star(1)%y, &
          dpos(1), nonstellar
          rewind (iunit1)
          read(iunit1, '(/,/)')
	  strin: do i=1, nstar
            if (non < 0) then
	      read(iunit1,*,iostat=iostat) star(i)%id, star(i)%x, star(i)%y, &
              dpos(i)
              nonstellar=0
            else
	      read(iunit1,*,iostat=iostat) star(i)%id, star(i)%x, star(i)%y, &
              dpos(i), nonstellar
            end if
            if (nonstellar == 1) then
              star(i)%col(1)%flg='ON'
            else
              star(i)%col(1)%flg='OO'
            end if
            if (iostat /= 0) then
              print*, 'Error in star position file at line ', i+3
              stop
            end if
          end do strin
          close(iunit1)
        else
          call nxtcls_in(strpos)
          iostat=read_cluster_file(nstar, ncol, colstr, star)
          where(star%col(1)%flg /= 'ON') star%col(1)%flg='OO'
        end if

        print*, '* ', nstar, ' star positions read.'
        if (cluster_cat) then 
          allocate(dpos(nstar))
          print*, '> Give radius to be searched for each star (<0 to fix position).'
          read(*,*) work
          dpos=work
        else
          if (non ==0) print*, '* including non-stellar flags.'
        end if
        if (nstar >= 100) print*, '* Files will not be made for each star.'

        ! Allocate lots more things that need to know the number of stars.
        if (nstar < 100) then
          allocate(ofopt(nstar))
          allocate(ofap(nstar))
        end if
        allocate(xpos_new(nstar))
        allocate(ypos_new(nstar))
        allocate(npos_new(nstar))
        allocate(xpos0(npsf))
        allocate(ypos0(npsf))

        ! ***************************************************************
        !                   Open offsets.dat
        ! ***************************************************************

        print*, ' '
        print*, '* This program also requires a file which has an'
        print*, '* image file name on one line, followed by the' 
        print*, '* offsets of that image on the next.'
        offsets: do
          print*, '> Give this file name (<cr>=offsets.dat).'
          offset=' '
          read (*, '(a)') offset
          if (offset .eq. ' ') offset='offsets.dat'
          inquire(file=offset, exist=there)
          if (.not. there) then
            print*, '* File not found.'
          else
            exit offsets
          end if
        end do offsets

        ! Open the file of offsets with between the data frames using a
        ! machine independent ARK routine.  It returns the number of the
        ! unit the file is opened on.
        call arkopn(iunit1, offset, ' ', 'OLD', 'READONLY', &
        'FORMATTED', ' ', 0)
        ! read the first couple of lines so we get a clue as to the
        ! seeing.
	read(iunit1,'(a)',iostat=iostat) ifname
        if (iostat < 0) then
          print*, 'Error reading offsets.dat file'
          stop
        end if
        read(iunit1,*,iostat=iostat) coff, sky, theoff, ntoe, fwhm
        if (iostat < 0) then
          print*, 'Error reading offsets.dat file'
          stop
        end if
        ! But now rewind the file.
        rewind (iunit1)


        if (clip_fwhm > 0.0) then

        ! ***************************************************************
        !                   Open seeing.dat
        ! ***************************************************************

          ! Create the file which records the seeing in each frame using a
          ! machine independent ARK routine.  It returns the number of the
          ! unit the file is opened on.
          call arkopn(iunit5, 'seeing.dat', ' ', 'NEW', 'OVERWRITE', &
          'FORMATTED', ' ', 0)
          rewind (iunit5)
          write(iunit5,*) '! Column 1 is times from times.dat file.'
          write(iunit5,*) '! Column 2 is geometric mean of X and Y FWHM'
          write(iunit5,*) '! of PSF star, in pixels.'

        end if

        ! ***************************************************************
        !        Read the file of observation and exposure times.
        ! ***************************************************************

        print*, ' '
        print*, '* Next we require a file which has the image file'
        print*, '* name on one line, and the observation time and'
        print*, '* and half the exposure length on the next (in any'
        print*, '* units you like, they are just copied to the output).'
        times: do
          print*, &
          '> Give the name of this file ("end"=no file, <cr>=times.dat).'
          timfil=' '
          read (*, '(a)', iostat=iostat) timfil
          ntim=0
          if (iostat < 0) exit times
          if (timfil=='END' .or. timfil=='end') exit times
          if (timfil == ' ') timfil='times.dat'
          inquire(file=timfil, exist=there)
          if (there) then
            ! Open the file with the times that each image was taken.
            call arkopn(iunit2, timfil, ' ', 'OLD', 'READONLY', &
            'FORMATTED', ' ', 0)
            ntim=0
            cttimes:  do
              read(iunit2, '(a)', iostat=iostat) timfil
              if (iostat < 0) exit cttimes
              read(iunit2, '(a)', iostat=iostat) timfil
              if (iostat < 0) exit cttimes
              ntim=ntim+1
            end do cttimes
            allocate(time(ntim))
            allocate(dtime(ntim))
            allocate(timnam(ntim))
            rewind (iunit2)
            do i=1, ntim
              read(iunit2, '(a)') timnam(i)
              read(iunit2, *) time(i), dtime(i)
              ! Remove leading blanks.
              deblank1: do j=0, 50
                if (timnam(i)(1:1) .eq. ' ') then
                  timnam(i)=timnam(i)(2:50)
                else
                  exit deblank1
                end if
              end do deblank1
            end do
            close(iunit2)
            print*, '* ', ntim, ' times read.'
            exit times
          else
            print*, '* File not found.'
          end if
        end do times

        ! ***************************************************************
        !                   Set up all the options
        ! ***************************************************************

        bad_sky: do
          print*, '* The sky will be flagged as ill-determined if the'
          print*, '* skew or chi-sqaured exceeds certain values.  If you'
          print*, '* don''t want flagging in either value, type -1.'
          print*, '> Give the values of chi-sqaured and skew (<cr>=-1,-1).'
          answer=' '
          read (*,'(a)') answer
          iostat=0
          bad_sky_chi=-1.0
          bad_sky_skw=-1.0
          if (answer .ne. ' ') read (answer,*,iostat=iostat) bad_sky_chi, bad_sky_skw
          if (iostat == 0) then
            call set_bad_sky(bad_sky_chi, bad_sky_skw)
            exit bad_sky
          end if
        end do bad_sky

        if (clip_fwhm > 0.0) then
          clip: do 
            ! Most of the signal-to-noise is obtained by setting the clipping
            ! radius to be 2*fwhm.  But there is little gain in speed by
            ! making it less than 3 pixels.
            clip_fwhm=max(3.0/fwhm, 2.0)
            print*, &
            '> Give the radius for clipping the weight map in FWHM, <CR> = ', &
            clip_fwhm
            answer=' '
            read (*,'(a)') answer
            iostat=0
            if (answer .ne. ' ') read (answer,*,iostat=iostat) clip_fwhm
            if (iostat == 0) exit clip
          end do clip
        end if

        aper: do 
          ! Set the default aperture radius to the optimum (as shown in paper),
          ! but don't let it be less than 2 pixels.
          aprad=max(2.0, 2.0*fwhm/3.0)
          print*, &
          '> Give the radius for aperture photometry (<0 to skip), <CR> = ', &
          aprad
          answer=' '
          read (*,'(a)') answer
          iostat=0
          if (answer .ne. ' ') read (answer,*,iostat=iostat) aprad
          if (iostat == 0) exit aper
        end do aper     

        if (clip_fwhm > 0.0) then
          optstar: do 
            iopt=star(1)%id
            print*, '* The extraction can be optimised for a particular star'
            print*, '* or for sky limited stars.'
            print*, '> Give the star number (or -1 for sky limit), <CR> = ', &
            iopt
            answer=' '
            read (*,'(a)') answer
            iostat=0
            if (answer .ne. ' ') read (answer,*,iostat=iostat) iopt
            if (iostat == 0) then
              if (iopt == -1) exit optstar
              do i=1, nstar
                if (star(i)%id == iopt) then
                  iopt=i
                  exit optstar
                end if
              end do
              print*, '* That star was not in the file, try again.'
            end if
          end do optstar
        end if

        ! Begin some output files.
        if (nstar < 100) then
          do i=1, nstar
            ofopt(i)='star00.opt'
 	    write(ofopt(i)(5:6),'(i2.2)') star(i)%id
            if (clip_fwhm > 0.0) then
              ! Create the output file for the optimal photometry of star nstar.
              call arkopn(iunit2, ofopt(i), ' ', 'NEW', 'OVERWRITE', &
              'FORMATTED', ' ', 0) 
              write(iunit2,*) '! Weight mapped clipped at ', clip_fwhm,  &
              ' times the FWHM.'
              if (dpos(i) > 0) then
                write(iunit2,*) '! Radius searched for star ', dpos(i), ' pixels.'
              else
                write(iunit2,*) '! Star position fixed.'
              end if
              write(iunit2,*) '! '
              close(unit=iunit2)
            end if
            if (aprad > 0.0) then
              ofap(i)=ofopt(i)
	      ofap(i)(8:10)='dat'
              ! Create the output file for the aperture photometry of star nstar.
              call arkopn(iunit4, ofap(i), ' ', 'NEW', 'OVERWRITE', &
              'FORMATTED', ' ', 0)
	      write(iunit4,*) '! Aperture of ', aprad, ' pixels radius.'
              if (dpos(i) > 0) then
                write(iunit4,*) '! Radius searched for star ', dpos(i), ' pixels.'
              else
                write(iunit4,*) '! Star position fixed.'
              end if
              write(iunit4,*) '! '
	      close(unit=iunit4)
            end if
          end do
        end if

        xpos_new=0.0
        ypos_new=0.0
        npos_new=0

        ! ***************************************************************
        !                   Now begin processing data.
        ! ***************************************************************

        frame: do

          ! *************************************************************
          !                    Go around this loop for each file.
          ! *************************************************************

          !  Reads in filename and offsets.
	  read(iunit1,*,iostat=iostat) ifname, iext
          if (iostat < 0) exit frame
	  read(iunit1,*,iostat=iostat) coff, sky, theoff, ntoe, fwhm
          if (iostat /= 0) then
            print*, 'Error reading offsets.dat file.'
            stop
          end if
          theoff=theoff/57.2957795

          ! Reads in image.
	  print*, '********************************************'
          ! Tell sizark what the name of the next file to be read is.
          ifname=adjustl(ifname)
          call nxtark_in(ifname, iext)
          ! Read in the data.  The image is contained in the array data.
          i=inpark(naxis, data, axdata)
          if (i < 0) then
            print*, '* Going to next file.'
            cycle frame
          end if
          ! And any flags.
          i=inpflg(naxis, pix_flg)
          ! Save the root of the file name.
          ifroot=ifname
          i=index(ifroot, '.', .true.)
          if (i > 0) ifroot(i:len(ifroot))=' '
                      
          ! Try and match the file name to a known time.
          if (ntim /= 0) then
            do i=1, ntim
              if (ifname .eq. timnam(i)) then
                timec=time(i)
                dtimec=dtime(i)
                goto 560
              end if
            end do
            ! When this fails use the frame number.
            print*, 'Failed to find a time to match this frame.'
            print*, 'Will use the frame number.'
          end if
          ctime=ifname
          do i=1, 50
            if (ichar(ctime(i:i)).lt.48 .or. &
            ichar(ctime(i:i)).gt.57) ctime(i:i)=' '
          end do 
          timec=0.0
          read (ctime, *, iostat=iostat) timec
          dtimec=0.5
560       print*, 'Frame time is ', timec, '+/-', dtimec

          if ( & 
            get_header_i('ARK-LO1', low(1))  < 0 .or.   &
            get_header_i('ARK-LO2', low(2))  < 0 .or.   &
            get_header_i('ARK-HI1', high(1)) < 0 .or.  &
            get_header_i('ARK-HI2', high(2)) < 0) then
            print*, 'The active area needs setting for ', ifname
            print*, 'using the program ark_key.  Going to next file.'
            deallocate(data, axdata, pix_flg)
            cycle frame
          end if
          if (get_header_r('ARK-GAIN', adu) < 0) then
            print*, 'The gain needs setting for ', ifname
            print*, 'using the program ark_key.  Going to next file.'
            deallocate(data, axdata, pix_flg)
            cycle frame
          end if
          
          if (clip_fwhm > 0.0) then
            xpos0 = psf%x*coff(1) + psf%y*coff(2) + coff(3)
            ypos0 = psf%y*coff(4) + psf%x*coff(5) + coff(6)
            do i=1, npsf              
              if (debug) print*, 'Estimated position of PSF star ', psf(i)%id, &
              ' is ', xpos0(i), ypos0(i)
            end do

            ! Now fit the PSF star.

            ! First find the sub-part of the list.
            do i=1, npsf
              if (psf(i)%ccd == iext) then
                start_ext=i
                exit
              end if
            end do
            do i=start_ext, npsf
              if (psf(i)%ccd /= iext) exit
              end_ext=i
            end do

            call psf_calc(data, pix_flg, end_ext-start_ext+1, xpos0(start_ext:end_ext), &
            ypos0(start_ext:end_ext), dpsf(start_ext:end_ext), &
            adu, high, low, fwhm, shape_par, ipsf, nfit)
            if (ipsf == -1) then
              print*, 'No good PSF star could be found within frame.'
              deallocate(data, axdata, pix_flg)
              cycle frame
            end if
            ipsf=ipsf+start_ext-1
            print*, 'Fitted PSF star ', psf(ipsf)%id, ' from ccd ', psf(ipsf)%ccd, &
            ' which gave FWsHM ', 1.665*shape_par(1), 1.665*shape_par(2)
            print*, 'Rotated at an angle of ', 57.29*shape_par(3), &
            ' degrees from the vertical.'
            print*, 'Chosen using ', nfit, ' stars.'
            cliprad = clip_fwhm*1.665*sqrt(shape_par(1)*shape_par(2))
            print*, 'Will use a clipping radius of ', cliprad, ' pixels.'
            !write(logopt_unit,*) '! PSF star ', psf(ipsf)%id, ' had FWsHM ', &
            !1.665*shape_par(1), 1.665*shape_par(2)
            !write(logopt_unit,*) '! Rotated at an angle of ', 57.29*shape_par(3), &
            !' degrees from the vertical.'
            !write(logopt_unit,*) '! Clipping radius was ', cliprad, &
            !' pixels.  Star chosen from', nfit
            ! Write to seeing.dat.
            write(iunit5,*) timec, sqrt(1.665*shape_par(1)*1.665*shape_par(2))
            ! We now have a better estimate of fwhm, so use this instead.
            fwhm=1.665*sqrt(shape_par(1)*shape_par(2))

            ! Now get a handle on the flux in the star to be extracted
            ! optimally.
            if (iopt > 0) then

              if (star(iopt)%ccd /= iext) then
                print*, 'Star to be optimised is not on this extension.'
                stop
              end if

              xpos = star(iopt)%x*coff(1) + star(iopt)%y*coff(2) + coff(3)
              ypos = star(iopt)%y*coff(4) + star(iopt)%x*coff(5) + coff(6)
              ! Call the extraction, with the normalisation explicitly set to one, and
              ! return the peak flux in optnrm.
              call extr(data, pix_flg, xpos, ypos, dpos(iopt), &
              adu, high, low, fwhm, cliprad, shape_par, &
              0.0, .false., xcomp, ycomp, star(iopt)%col(1)%data, star(iopt)%col(1)%err, &
              star(iopt)%x, star(iopt)%y, xerr, yerr, optnrm, cflag, skynos)
              if (cflag /= 'O') then
                print*, 'Star to be optimised has flag ', cflag
                print*, 'Going to next frame.'
                deallocate(data, axdata, pix_flg)
                cycle frame
              end if
              print*, 'Extractions will be optimised for stars with ', &
              optnrm, ' peak counts.'
              optnrm=optnrm/(skynos*skynos)
            else                            
              optnrm=0.0
            end if

          else
            
            ! Don't do optimal photometry, but we will need a FWHM for fitting
            ! the star to find its position.
            shape_par(1)=fwhm/1.665
            shape_par(2)=fwhm/1.665
            shape_par(3)=0.0
            ! And a cliprad.
            cliprad=2.0*fwhm
            ! And (believe it or not) a star flux to be optimised.
            optnrm = 0.0
          end if

          ! This is a new file name, not merely an extension of the last
          ! one, so let's open some new files.
          if (iext <= 1) then
            header(2)=' '
            write(header(3),*) timec, dtimec, ' are the TDB and half exposure time.' 
            colstr(1)='?'
            if (clip_fwhm > 0.0) then
              write(header(1),*) 'Photometry parameters ', &
                   shape_par(1:3), cliprad
              call set_cluster_header(header)
              call nxtcls_out(trim(ifroot)//'_opt.cat')
              call start_cluster_file(1, colstr, cluster_opt)
            end if
            if (aprad > 0.0) then
              write(header(1),*) 'Aperture of ', aprad
              call set_cluster_header(header)
              call nxtcls_out(trim(ifroot)//'_ap.cat')
              call start_cluster_file(1, colstr, cluster_ap)
            end if
            call arkopn(isky, trim(ifroot), 'sky', 'NEW', &
                 'OVERWRITE', 'FORMATTED', ' ', 0)
            write(isky,*) 'Columns are CCD and star number, reduced chi-squared'
            write(isky,*) 'skew, sky counts, star counts, star uncertainty and flag.'
            write(isky,*)             
          else
            if (clip_fwhm > 0.0) then
              open(unit=cluster_opt, file=trim(ifroot)//'_opt.cat', &
              position='append')
            end if
            if (aprad > 0.0) then
              open(unit=cluster_ap, file=trim(ifroot)//'_ap.cat', &
              position='append')
            end if
            open(unit=isky, file=trim(ifroot)//'.sky', position='append')
          end if

          each_star: do istar=1, nstar                   
            ! ***********************************************************
            !                Go around this loop for each star.
            ! ***********************************************************

            ! Is it in this extension?
            if (iext /= star(istar)%ccd) cycle each_star 


            ! Calculate approximate position of star.
            xpos = star(istar)%x*coff(1) + star(istar)%y*coff(2) + coff(3)
            ypos = star(istar)%y*coff(4) + star(istar)%x*coff(5) + coff(6)

            comp=.false.

            opt=star(istar)
            !print*, 'Inheriting flag ', opt%col(1)%flg

            call extr(data, pix_flg, xpos, ypos, dpos(istar), &
            adu, high, low, fwhm, cliprad, shape_par, &
            optnrm, comp, xcomp, ycomp, opt%col(1)%data, opt%col(1)%err, &
            opt%x, opt%y, xerr, yerr, peak, cflag)
            
            if (opt%col(1)%flg=='OO') opt%col(1)%flg='O'//cflag

            ! If there's a bad pixel, how far away is it?
            ! First set nearest bad pixel to
            ! be closest edge of the frame.
            badpix=min(opt%x-real(low(1)), real(high(1))-opt%x, &
            opt%y-real(low(2)), real(high(2))-opt%y )            
            if (opt%col(1)%flg == 'OF') then
              do i=int(opt%x-cliprad), int(opt%x+cliprad)+1
                do j=int(opt%y-cliprad), int(opt%y+cliprad)+1
                  if (pix_flg(i,j) == 'F') then
                    badpix = min (sqrt( &
                    (opt%x-real(i))*(opt%x-real(i)) + &
                    (opt%y-real(j))*(opt%y-real(j)) ), badpix)  
                  end if
                end do
              end do
            end if

            if (clip_fwhm > 0.0) then
              !if (opt%col(1)%flg == 'OO') then
              !  print*, 'Target star ', opt%id, ' gave ', &
              !  opt%col(1)%data, ' counts +/- ', &
              !  100.0*opt%col(1)%err/abs(opt%col(1)%data), '%'
              !else
              !  print*, 'Target star ', opt%id, &
              !  ' was flagged with ', opt%col(1)%flg
              !end if
              if (100*(istar/100) == istar) then
                print*, 'Done ', istar, ' out of ', nstar, ' stars.'
              end if
              if (nstar < 100) then
                ! Write it to the appropriate star file.
                call arkopn(iunit2, ofopt(istar), ' ', 'OLD', &
                'APPEND', 'FORMATTED', ' ', 0)
                if (opt%col(1)%flg == 'OF') then
                  write(iunit2,*) timec, opt%col(1)%data, opt%col(1)%err, & 
                  dtimec, max(nint(badpix),1)
                else
                  write(iunit2,*) timec, opt%col(1)%data, opt%col(1)%err, & 
                  dtimec, opt%col(1)%flg
                end if
                close(iunit2)
              end if
              ! And the params.sky file.
              write(isky,*) opt%ccd, opt%id, sky_chi, sky_par(1),  sky_par(2), &
              opt%col(1)%data, opt%col(1)%err, opt%col(1)%flg
              ! And to a cluster file.
              call flux_to_mag(opt%col(1))
              call write_star(cluster_opt, opt, 1)
            end if

            if (dpos(istar) > 0.0) then
              print*, 'Fitted position of star ', opt%id, &
              sqrt((xpos-opt%x)**2.0 + (ypos-opt%y)**2.0), &
              'from start position.'
            end if

            if (opt%col(1)%flg == 'OO') then
              ! Keep a track on the mean position of the star in all frames.
              work =           ((opt%x-coff(3)) - &
              (coff(2)/coff(4))*(opt%y-coff(6)))/ &
              (coff(1) - (coff(5)*coff(2)/coff(4)))
              xpos_new(istar) = xpos_new(istar) + work
              ypos_new(istar) = ypos_new(istar) + &
              ((opt%y - coff(5)*work - coff(6))/coff(4))
              ! xpos_new(istar)=xpos_new(istar)+opt%x-xoff
              ! ypos_new(istar)=ypos_new(istar)+opt%y-yoff
              npos_new(istar)=npos_new(istar)+1
              if (debug) print*, '@ Star is at ', opt%x, opt%y
            end if

            ! And the aperture photometry. 
            ap=star(istar)
            ! Use the centroided positions from the optimal photometry.
            ap%x=opt%x
            ap%y=opt%y
            if (aprad > 0.0) then
              sky_stuffed = .false.
              ! The box size is determined in the way described in the paper.
              ibox=int(sqrt(628.4*aprad*aprad+ 4.0*fwhm*fwhm))
              call skyfit(data, ap%x, ap%y, fwhm, ibox, low, high, &
              pix_flg, skycnt, skyerr, skynos, cflag)
              if (cflag == 'I') sky_stuffed = .true.
              if (cflag == 'B') then
                ap%col(1)%data = 0.0
                ap%col(1)%err = 0.0
              else
                ap%col(1)%data=sum_aper(data, aprad, skycnt, skyerr, skynos, &
                ap%x, ap%y, adu, low, high, pix_flg, ap%col(1)%err, cflag)
                if (sky_stuffed) cflag = 'I'
                if (ap%col(1)%flg/='ON') ap%col(1)%flg='O'//cflag
                if (clip_fwhm < 0.0) then
                  if (ap%col(1)%flg == 'OO') then
                    print*, 'Target star ', ap%id, ' gave ', &
                    ap%col(1)%data, ' counts +/- ', &
                    100.0*ap%col(1)%err/abs(ap%col(1)%data), '%'
                  else
                    ! print*, 'Target star ', ap%id,  &
                    ! ' was flagged with ', ap%col(1)%flg
                  end if
                else  
                  ! if (ap%col(1)%flg /= 'OO') &
                  ! print*, 'Aperture photometry was flagged with ', ap%col(1)%flg
                end if
              end if
            end if

            if (aprad > 0.0) then
              if (nstar < 100) then
                call arkopn(iunit4, ofap(istar), ' ', 'OLD', &
                'APPEND', 'FORMATTED', ' ', 0)
                if (ap%col(1)%flg == 'OF') then 
                  write(iunit4,*) timec, ap%col(1)%data, ap%col(1)%err, &
                  dtimec, max(nint(badpix),1)
                else
                  write(iunit4,*) timec, ap%col(1)%data, ap%col(1)%err, &
                  dtimec, ap%col(1)%flg
                end if
                close(iunit4)
              end if
              call flux_to_mag(ap%col(1))
              call write_star(cluster_ap, ap, 1)
            end if

          end do each_star

          deallocate(data, axdata, pix_flg)
          if (clip_fwhm > 0.0) close(cluster_opt)
          if (aprad > 0.0) close(cluster_ap)
          close(isky)

        end do frame

        close(iunit1)
        if (clip_fwhm > 0.0) close(iunit5)

        ! Finally, write a file of updated star positions.
        call arkopn(iunit1, 'stars.pos_new', ' ', 'NEW', 'OVERWRITE', &
        'FORMATTED', ' ', 0)
        write(iunit1,*)
        write(iunit1,*)
        write(iunit1,*)
        do istar=1, nstar              
          if (npos_new(istar) > 0) then
            write(iunit1,*) star(istar)%id, &
            xpos_new(istar)/real(npos_new(istar)), &
            ypos_new(istar)/real(npos_new(istar)), 0.0
          end if
        end do        

      end program opphot
