      ! Mapccd 

      ! Purpose.
      ! Calculates the offset between 2 images by finding star positions in
      ! the second image.

      ! Note the Fortran 90 module structure means that the program follows
      ! the subroutines.

      ! History.
      ! (1) Originally by Robin Corbet, January 1984.
      ! (2) New pattern matching alogorithm added by Tim Naylor
      ! (3) Made to get more accurate offsets, and converted to Fortran 90
      !     by Tim Naylor, Christmas 1996.

      ! Version 1.1 
           ! Minor bug fixes, and made to write seeing.map
      ! Version 2.0
           ! The position angle of the Gaussian PSF allowed to run free.
      ! Version 2.1
           ! More comments to aid porting and module structure improved.
      ! Version 2.2
           ! Changed for new call to gfit.
      ! Version 2.3
           ! Changed to make sure the low and high are not silly.

      module mapccd_subs

        use flamsteed_subs
        use define_star
        use ark_file_io

        implicit none

      contains

        subroutine maps(data, low, high, approx_fwhm, &
        skycnt, skynos, fnsig, pix_flg, max_find, star, nstar, mstar)

        ! Finds stars - called by mapccd.

        ! Passed variables (inputs).
        integer, dimension(2), intent(in) :: low, high
        integer :: approx_fwhm
        real, dimension(:,:) :: data
        real :: skycnt, skynos, fnsig
        ! The pixel flag arrays.
        character, dimension(:,:), intent(in) :: pix_flg
        ! Stop searching when you've found how many stars?
        integer :: max_find

        ! Passed variables (outputs).
        integer :: nstar, mstar
        ! The peak flux (including sky), star position etc.
        type(a_star), dimension(mstar) :: star

        ! Other stuff.
        real :: up, good_max
      
        ! Counters.
        integer :: i, j, isig, istar, iaxis1, iaxis2

        ! For the significant events.
        integer, dimension(:), allocatable :: x_sig, y_sig
        integer :: n_sig, i_sig

        logical, dimension(:,:), allocatable :: bin_sig_map


        ! Just in case we're on the graphics screen, return to text.
        call pltcls

        isig=0
        nstar=0
        
        good_max=minval(data)
        do iaxis2=low(2), high(2)
          do iaxis1=low(1), high(1)
            if (pix_flg(iaxis1,iaxis2) == 'O') then
              good_max=max(good_max,data(iaxis1,iaxis2))
            end if
          end do
        end do
        print*, 'Highest good pixel has ', good_max, ' counts.'
        
        go_down: do

          ! Set the upper limit for an insignificant event.
          isig=isig+1
          up = skycnt + (good_max-skycnt)/real(2**isig)
          if (up-skycnt < fnsig*skynos) then
            if (isig == 1) then
              print*, 'Minimum sigma level is too low.'
              stop
            else
              exit go_down
            end if
          end if
          ! And set the counter for the number of events found.
          n_sig=0      
 
 
          allocate(bin_sig_map(high(1),high(2)))
          bin_sig_map=.false.
          ! Now hunt for significant events, within region specified.
          do j=low(2), high(2)
            do i=low(1), high(1)
              if (data(i,j) >= up) then
                if (pix_flg(i,j)=='O' .or. pix_flg(i,j)=='L' .or. &
                  pix_flg(i,j)=='S') then
                  bin_sig_map(i,j)=.true.
                end if
              end if
            end do  
          end do
          n_sig=count(bin_sig_map)
          allocate(x_sig(n_sig), y_sig(n_sig))
          i_sig=0
          do j=low(2), high(2)
            do i=low(1), high(1)
              if (bin_sig_map(i,j)) then
                i_sig=i_sig+1
                x_sig(i_sig)=i
                y_sig(i_sig)=j
                ! Move to the position of the significant event.
                call movxy(real(x_sig(i_sig)),real(y_sig(i_sig)))
                ! And draw a dot.
                call point(2)
              end if
            end do
          end do
          deallocate(bin_sig_map)
 
          ! Now change the list of significant events into one of stars.
          call filter(x_sig, y_sig, n_sig, high, data, pix_flg, &
          star, nstar)

          deallocate(x_sig, y_sig)
          
          print*, 'Pass ', isig, ' at ', (up-skycnt)/skynos, &
          ' sigma, gave ', n_sig, ' events and ', nstar, ' stars.'

          if (nstar >= max_find) then
            print*, nstar, ' found before completing all passes.'
            print*, 'This should be enough, so stopping search.'
            exit go_down
          end if

          if (nstar == size(star)) then
            print*, 'Fill the array star before completing all'
            print*, 'passes.'
            exit go_down
          end if
          
        end do go_down
        
        ! Make the fluxes sensible.
        do istar=1, nstar
          star(istar)%col(1)%data=0.0
          i=nint(star(istar)%x)
          j=nint(star(istar)%y)
          !print*, i,j
          do i=nint(star(istar)%x)-approx_fwhm, nint(star(istar)%x)+approx_fwhm
            if (i>=low(1) .and. i<=high(1)) then
              do j=nint(star(istar)%y)-approx_fwhm, nint(star(istar)%y)+approx_fwhm
                if (j>=low(2) .and. j<=high(2)) then
                  star(istar)%col(1)%data=star(istar)%col(1)%data+data(i,j)-skycnt
                  !print*, i,j, stsize(istar)
                end if
              end do
            end if
          end do
          ! write(44,*) stsize(istar), data(nint(xpos(istar)), &
          ! nint(ypos(istar)))-skycnt, xpos(istar), ypos(istar)
        end do

        end subroutine maps

      end module mapccd_subs


      program mapccd

        ! This module contains all the ARK I/O routines.
        use ark_file_io
        ! and the header routines
        use get_header
        ! This module contains the optimal extraction stuff.
        use opt_extr
        ! This is all the subroutines.
        use mapccd_subs
        use blnclc_mod
        ! The definition of the star structure.
        use define_star

        implicit none

        ! The maximum number of stars to be used in the pattern matching.
        ! This number largely depends on how much CPU you've got.
        integer, parameter :: max_match=100


        integer, dimension(2) :: naxis
        real, allocatable, dimension(:,:) :: data, axdata, noise_map
        ! The pixel flag arrays.
        character, dimension(:,:), allocatable :: pix_flg
        character :: ans
        character(len=50) :: ifname, ofname, answer
        integer :: iext
        integer :: iostat, iflag
        character :: cflag

        ! Significance above which something may be a star.
        real :: fnsig

        ! Counters.
        integer :: i, j

        ! The actual region to be used.
        integer, dimension(2) :: low, high
        logical :: first

        ! Arrays of fluxes in a star, an order of brilliance.
        integer :: mstar

        ! For the profile fitting.
        real, dimension(6) :: a_par
        real, dimension(2) :: e_pos
        real :: box

        ! For the sky fitting.
        real :: skycnt, skyerr, skynos
        real :: adu, approx_fwhm, fwhm, mean_noise

        ! All the stuff for blnclc and sort.
        integer :: ntoe, mmatch
        real :: difmag, diff, xoff, yoff, theoff, relmag
        type(a_star), allocatable, dimension(:) :: star2, star1
        integer :: nstar1, nstar2
        logical :: rotate

        ! Unit numbers for the offset files, MAPCCD.DAT, the magnitudes of
        ! the stars found, the unit number of the seeing file and the file 
        ! of frames not matched.
        integer :: ioff, nst, isee, ilost

        ! A general star counter.
        integer :: istar

        ! For finding the FWHM.
        real, dimension(3) :: shape_par
        integer :: ipsf, nfit
        real, allocatable, dimension(:) :: dpsf

        ! For the debugger.
        logical :: debug
        common /bug/ debug

        real :: peak

        ! Do we allow saturated stars?
        logical, parameter :: sat=.true.
        integer :: nrej

        call setbug()
        if (debug) call opt_debug_on()

        ! Create the three output files using the machine independent ARKOPN
        ! routine.  In each case ARKOPN returns the number of the unit
        ! that the file is opened on.
        call arkopn(ioff, 'offsets.dat', ' ', 'NEW', 'OVERWRITE', &
        'FORMATTED', ' ', 0)
        call arkopn(ilost, 'lost.dat', ' ', 'NEW', 'OVERWRITE', &
        'FORMATTED', ' ', 0)
        call arkopn(isee, 'seeing.map', ' ', 'NEW', 'OVERWRITE', &
        'FORMATTED', ' ', 0)        
 
        seeing: do 
          approx_fwhm=2.0
          print*, &
          '> Give approximate FWHM (in pixels) for stellar images <CR> = ', &
          approx_fwhm
          answer=' '
          read (*,'(a)') answer
          iostat=0
          if (answer .ne. ' ') read (answer,*,iostat=iostat) approx_fwhm
          if (iostat == 0) exit seeing
        end do seeing

        search: do 
          fnsig=3
          print*, '> What sigma events do you want to search for <CR> = ', &
          fnsig
          answer=' '
          read (*,'(a)') answer
          iostat=0
          if (answer .ne. ' ') read (answer,*,iostat=iostat) fnsig
          if (iostat /= 0) cycle search
          if (fnsig .le. 0.0) then
            print*,'* You cannot have a negative value.'
            cycle search
          end if               
          exit search
        end do search

        rotation: do
          print*, '> Might there be a rotation between the frames (<cr>=no)?' 
          ans=' '
          read(*,'(a)') ans
          if (ans.eq.'n' .or. ans.eq.'N' .or. ans.eq.' ') then
            rotate=.false.
            exit rotation
          else if (ans.eq.'y' .or. ans.eq.'Y') then
            rotate=.true.
            exit rotation
          else
            print*, '* Please answer yes or no.'
          end if
        end do rotation

        first=.true.

        frame: do

          print*,'**************************************************'

          ! Read the image.  The image array is data.
          call extark('fit')
          iflag=inpark(naxis, data, axdata)
          if (iflag == -1) then
            if (first) then
              exit frame
            else
              print*, 'Now for the next reference frame.'
              first=.true.
              deallocate(star1)
              cycle frame
            end if
          else if (iflag < -1) then
            print*, 'Failed to read file, INPARK returned = ', iflag
            cycle frame
          end if
          if (naxis(1)*naxis(2) == 0) cycle frame
          ! lstark returns the name of the last image file read.
          call lstark_in(ifname, iext)

          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 ', trim(ifname), ' extension ', iext
            print*, 'using the program pix_info.  Going to next file.'
            deallocate(data, axdata)
            cycle frame
          end if
          if (get_header_r('ARK-GAIN', adu) < 0) then
            print*, 'The gain needs setting for ', trim(ifname), ' extension ', iext
            print*, 'using the program pix_info.  Going to next file.'
            deallocate(data, axdata)
            cycle frame
          end if
          iflag=inpflg(naxis, pix_flg)

          !allocate(noise_map(naxis(1),naxis(2)))
          !call back(data, low, high, 200, noise_map, mean_noise, &
          !pix_flg)
          !deallocate(noise_map)

          ! Estimate the sky and the sigma on it.
          call skyfit(data, real(high(1)+low(1))/2.0, &
          real(high(2)+low(2))/2.0, approx_fwhm, &
          min(high(1)-low(1),high(2)-low(2)), low, high, &
          pix_flg, skycnt, skyerr, skynos, cflag)
          print*, 'Sky is at ', skycnt
          print*, 'Sigma is ', skynos
     
          ! Check there were some counts in the frame.
          if (skynos.eq.0.0 .or. skycnt.eq.0.0) then
            print*,'***WARNING***'
            print*,'NO COUNTS IN THAT FRAME, SKIPPING TO NEXT'
            write(ilost,*) trim(ifname), iext
            write(ilost,*) 'No counts in frame'
            deallocate(data, axdata, pix_flg)
            cycle frame
          end if

          ! Now allocate all the stuff for finding stars.  
          mstar=max(400, (high(1)-low(1))*(high(2)-low(2))/200)
          allocate(star2(mstar))
          call zero_star(star2)
          ! Give all the stars an extension number.
          star2%ccd=iext
          if (fnsig .gt. 0.0) then
            ! Start up the plotting package.
            call pltopn(2)
            ! Set the screen to data co-ordinate conversion.  First four
            ! arguments are x,y pairs in percentage of avialable screen, then
            ! the same points in data co-ordinates.
            call limit(5.0, 5.0, 50.0, 95.0, 0.0, 0.0, &
            real(naxis(1)+1), real(naxis(2)+1))
            ! Put the file name on the graphics screen above the image.
            call ptxtm('t', ifname)           
            ! Draw a box around the whole chip.
            call boxm(0.5,0.5,real(naxis(1))+0.5,real(naxis(2))+0.5)
            ! And now around the active region.
            call boxm(real(low(1))-0.5,real(low(2))-0.5, &
            real(high(1))+0.5,real(high(2))+0.5)
            call maps(data, low, high, nint(approx_fwhm), &
            skycnt, skynos, fnsig, pix_flg, 5*max_match, star2, nstar2, mstar)
          end if
          if (nstar2 > 0) then
            ! Now throw out stars that are saturated.
            if (.not. sat) then
              i=1
              junk_sat: do
                if (pix_flg(nint(star2(i)%x),nint(star2(i)%y)) == 'S') then
                  do j=i, nstar2-1
                    star2(j)=star2(j+1)
                  end do
                  nstar2=nstar2-1
                else
                  i=i+1
                end if
                if (nstar2 < i) exit junk_sat
              end do junk_sat
              print*, 'After removing saturated stars, number left =', nstar2
            end if
          end if
          if (nstar2 > 0) then
            i=1
            junk_neg: do
              if (star2(i)%col(1)%data <= 0.0) then
                do j=i, nstar2-1
                  star2(j)=star2(j+1)
                end do
                nstar2=nstar2-1
              else
                i=i+1
              end if
              if (nstar2 < i) exit junk_neg
            end do junk_neg
            print*, 'After removing stars with negative counts, number left =', nstar2
          end if
          if (nstar2 == 0) then
            print*,'***WARNING***'
            print*,'NO STARS FOUND IN THAT FRAME, SKIPPING TO NEXT'
            write(ilost,*) trim(ifname), iext
            write(ilost,*) 'No stars found.'
            deallocate(data, axdata, pix_flg)
            deallocate(star2)
            cycle frame
          end if

          ! Place the stars in order.
          do istar=1, nstar2
            call flux_to_mag(star2(istar)%col(1))
          end do
          call sort_star(star2, nstar2, 'increasing_col(1)%data')

          ! Find the seeing.
          allocate(dpsf(nstar2))
          dpsf=2.0*approx_fwhm
          call psf_calc(data, pix_flg, nstar2, star2%x, star2%y, &
          dpsf, adu, high, low, approx_fwhm, shape_par, ipsf, nfit)
          if (debug) print*, '@ Back from PSF calc.'
          deallocate(dpsf)
          if (ipsf < 0) then
            print*, '* WARNING, failed to measure fwhm, will use the guess.'
            fwhm=approx_fwhm
          else
            fwhm=1.665*sqrt(shape_par(1)*shape_par(2))
            print*, 'Seeing for the star at marked with triangle is ', &
            1.665*shape_par(1), ' by ',  1.665*shape_par(2)
            print*, 'pixels,  with a position angle of ', shape_par(3)*57.30, &
            ' degrees, at pixel ', star2(ipsf)%x, star2(ipsf)%y
            ! Do the next graphics in yellow.
            call colour('y')
            ! Move to the required position.
            call movxy(star2(ipsf)%x, star2(ipsf)%y)
            ! And draw a filled traingle.
            call point(12)
            ! Return the colour to white.
            call colour('w')
            ! And return to the text screen.
            call pltcls()
            ! And write out the seeing to file.
            write(isee,*) trim(ifname), iext
            write(isee,*) fwhm, 1.665*shape_par(1), 1.665*shape_par(2)
          end if

          ! Now fit each star properly.  
          i=1
          istar=1
          nrej=0
          fit_all: do i=1, nstar2
            a_par(1:3)=shape_par(1:3)
            a_par(4)=data(nint(star2(i)%x),nint(star2(i)%y))-skycnt
            a_par(5)=star2(i)%x
            a_par(6)=star2(i)%y
            if (debug) print*, '@ Start parameters ', a_par
            call gfit(data, .true., 3.0, .false., skycnt, skynos, adu, &
            low, high, pix_flg, a_par, e_pos, cflag)
            if (nint(a_par(5)) < 1) cycle fit_all
            if (nint(a_par(6)) < 1) cycle fit_all
            if (nint(a_par(5)) > naxis(1)) cycle fit_all
            if (nint(a_par(6)) > naxis(2)) cycle fit_all
            if ((cflag=='O' .or. cflag=='L' .or. (cflag=='S'.and.sat)) .and. a_par(4)>0.0) then
              star2(istar)=star2(i)
              star2(istar)%col(1)%data=sum_aper(data, fwhm, skycnt, skyerr, &
              skynos, a_par(5), a_par(6), adu, low, high, pix_flg, &
              star2(istar)%col(1)%err, cflag)
              peak=data(nint(a_par(5)), nint(a_par(6)))-skycnt
              !write(23,*) star2(istar)%col(1)%data/peak, star2(istar)%col(1)%data
              star2(istar)%x=a_par(5)
              star2(istar)%y=a_par(6)
              ! Now, lets apply a stellarity test.  This is based on the fact that
              ! the flux in a Gaussian is peak*sigma1*sigma2*pi.  We drop the
              ! factor pi to make it inclusive.
              !if (a_par(1)*a_par(2)*data(nint(a_par(5)), nint(a_par(6)))-skycnt < &
              !  star2(istar)%col(1)%data) then
                istar=istar+1
                ! Have we reached the maximum number we are going to pattern
                ! match.  If so, why fit any more?
                if (istar-1 == max_match) exit fit_all
              !else
              !  nrej=nrej+1
              !end if
            end if
          end do fit_all
          nstar2=istar-1
          !print*, 'Number rejected by Stuart''s filter was ', nrej
          print*, 'Number of stars sucessfully fitted was ', nstar2

          call doubles(fwhm, star2, nstar2)
          print*, 'Final number, after removing doubles is ', nstar2

          do istar=1, nstar2
            ! Move to the center of the star.
            call movxy(star2(istar)%x, star2(istar)%y)
            ! Change the graphics colour to red.
            call colour('r')
            ! And draw a small circle.
            call point(0)
          end do
          call pltcls()
          print*, 'Marked good stars with circles.'

          do istar=1, nstar2
            call flux_to_mag(star2(istar)%col(1))
          end do
          call sort_star(star2, nstar2, 'increasing_col(1)%data')
          ! Create an output file name.
          ofname=ifname
          i=index(ifname, '.')
          ofname(i:i)='_'
          write(ofname(i+1:i+2), '(i2.2)') iext
          ofname(i+3:50)='_map.cat'
          write(*,*) 'Sending postions and magnitudes to file ', trim(ofname)
          ! Create the file using the machine independent ARKOPN.
          call arkopn(nst, ofname, ' ', 'NEW', 'OVERWRITE', &
          'FORMATTED', ' ', 0)
          write(nst,*) '1 colour'
          write(nst,*) '?'
          write(nst,*)
          do i=1, nstar2
            call write_star(nst, star2(i), 1)
          end do
          close(nst)
          if (first) then
            nstar1=nstar2
            allocate(star1(nstar1))
            star1(1:nstar1) = star2(1:nstar1)
            write(ioff,*) trim(ifname), iext
            write(ioff,*) 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, skycnt, 0.0, 0, fwhm
            first=.false.
          else
            mmatch=min(nstar1,nstar2)
            ! The number should be set to the lowest number you need to be
            ! throughly convinced of a match.
            ntoe=min(40,max(nstar1,nstar2))
            difmag=3.0
            ! The difference allowable in pixels.  Differential atmospheric
            ! refraction across the CCD is the main problem here.  At 60
            ! degrees zenith distance its about 0.001 percent at V.  We allow
            ! double this, or two pixels, whichever is the greater.
            diff=max(2.0, 0.001*maxval(high-low))
            !call blnclc(ntoe, mmatch, difmag, diff, xpos1, ypos1, &
            !stsize1, nstar1, star2%x, star2%y, star2%col(1)%data, nstar2, &
            !rotate, xoff, yoff, theoff, relmag) 
            call blnclc(ntoe, mmatch, difmag, diff, star1, nstar1, &
            star2, nstar2, rotate, xoff, yoff, theoff, relmag) 
            if (ntoe .gt. 0) then
              write(*,*) 'Offset is ', xoff, yoff, ' pixels and '
              write(*,*) theoff*57.2957795, ' degrees from ', ntoe, &
              ' stars.'
              write(ioff,*) trim(ifname), iext
              write(ioff,*) cos(theoff), -sin(theoff), xoff, &
                            cos(theoff),  sin(theoff), yoff, &
                            skycnt, theoff, ntoe, fwhm
            else
              print*, 'Failed to match any stars.'
              write(ilost,*) trim(ifname), iext
              write(ilost,*) 'Failed to match any stars.'
            end if
          end if

          deallocate(data,axdata, pix_flg)
          deallocate(star2)

        end do frame
 
        print*,'Normal finish.'
        print*,'Offsets between frames written to "offsets.dat".'
        ! Finish with the graphics system.
        call pltend()

      end program mapccd
