  module define_star

    implicit none

    integer, parameter :: mcol=10

    type a_colour
      real :: data = 0.0
      real :: err = 0.0
      logical :: neg_flux = .false.
      character(len=2) :: flg = 'AA'
    end type a_colour
  
    type a_star
      integer :: field = 0
      integer :: id    = 0
      integer :: ccd   = 0
      integer :: ra_h  = 0
      integer :: ra_m  = 0
      real    :: ra_s  = 0.0
      integer :: dc_d  = 0 
      integer :: dc_m  = 0
      real    :: dc_s  = 0.0
      ! At the moment the rule is that if a declination is negative one of
      ! the above three should be set negative.  Eventually it will all be
      ! done via this variable.
      character(len=1) :: dc_sign = '+'
      real :: x = 0.0
      real :: y =0.0
      type(a_colour), dimension(mcol) :: col
    end type a_star

    logical, private :: pmt_nxtcls_out=.true., pmt_nxtcls_in=.true.
    character(len=100) :: nam_nxtcls_out, nam_nxtcls_in, nam_lstcls_in

    character(len=80), dimension(3), public, save :: cluster_header=' '

    interface zero_star
      
      module procedure zero_star_one
      module procedure zero_star_array

    end interface

    contains

    subroutine write_star(iunit, star, ncol)
      
      ! Writes out a star.  If any of the dc_d, dc_m or dc_s are negative 
      ! the declination is written out with a minus sign.
      ! Otherwise its given a plus sign.
      ! Eventually this code will use dc_sign, but at the moment not enough
      ! code sets it for its value to be reliable.

      integer, intent(in) :: iunit
      type(a_star), intent(in) :: star
      integer, optional :: ncol

      integer :: icol, jcol
      character(len=15) :: form
      type(a_star) :: st
      real :: flux, ferr
      
      jcol=mcol
      if (present(ncol)) jcol=ncol

      st=star

      if (st%x<-9999.99 .or. st%x>99999.99 &
     .or. st%y<-9999.99 .or. st%y>99999.99) then
        st%x=0.0
        st%y=0.0
      end if

      ! Deal with any lingering sign issues.
      if (st%dc_sign/='+' .and. st%dc_sign/='-') st%dc_sign='+'
      if (st%dc_d < 0) st%dc_sign='-'
      if (st%dc_m < 0) st%dc_sign='-'
      if (st%dc_s < 0.0) st%dc_sign='-'

      ! Now negative fluxes, noted in the code by neg_flux.
      do icol=1, jcol
        if (st%col(icol)%neg_flux .and. st%col(icol)%flg(2:2)=='O') then
          ! These are flagged in the output file by M.
          st%col(icol)%flg(2:2)='M'
          ! First get the flux from the magnitude.
          flux=10.0**(-0.4*st%col(icol)%data)
          ! And now uncertainty in flux.
          ferr=10.0**((st%col(icol)%err-st%col(icol)%data)/2.5) - flux
          ! Now the horrible way these are encoded in cluster files.
          ! The magnitude is placed in the uncertainty column with the sign 
          ! changed.
          st%col(icol)%err=-st%col(icol)%data
          ! And the magnitude column is the flux uncertainty converted to a 
          ! magnitude.
          st%col(icol)%data=-2.5*log10(ferr)
        end if
      end do
      
      if (st%id>999999 .or. st%id<-99999) then
        print*, 'Format statement cannot cope as star%id is ', st%id
      end if
      
      inquire(unit=iunit, form=form)
      if (form == 'UNFORMATTED') then
        write(iunit) st%field, st%ccd, st%id, st%ra_h, &
        st%ra_m, st%ra_s, st%dc_sign, abs(st%dc_d), abs(st%dc_m), &
        abs(st%dc_s), st%x, st%y, (st%col(icol)%data, st%col(icol)%err, &
        st%col(icol)%flg,icol=1,jcol)
      else
        write(iunit,10) st%field+real(st%ccd)/100.0, st%id, st%ra_h, &
        st%ra_m, st%ra_s, st%dc_sign, abs(st%dc_d), abs(st%dc_m), &
        abs(st%dc_s), st%x, st%y, (st%col(icol)%data, st%col(icol)%err, &
        st%col(icol)%flg, icol=1,jcol)
      end if
      
 10   format(1x,f6.2,2x,i6,2x,i2.2,1x,i2.2,1x,f6.3,1x, &
                           a1,i2.2,1x,i2.2,1x,f5.2,2x, &
      2(f9.3,2x),4(f9.3,2x,f9.3,2x,a2))
      
    end subroutine write_star

    integer function read_star(iunit, star, ncol)

      integer, intent(in) :: iunit
      type(a_star), intent(inout) :: star
      integer, optional :: ncol

      integer :: icol, jcol, iostat
      character(len=4) :: dc_d
      real :: field_ccd
      character(len=15) :: form
      real :: ferr, flux

      jcol=mcol
      if (present(ncol)) jcol=ncol

      inquire(unit=iunit, form=form)
      if (trim(form) == 'FORMATTED') then
        read(iunit,*,iostat=iostat) field_ccd, star%id, star%ra_h, &
        star%ra_m, star%ra_s, dc_d, star%dc_m, star%dc_s, star%x, &
        star%y, (star%col(icol)%data, star%col(icol)%err, &
        star%col(icol)%flg,icol=1,jcol)
      else if (trim(form) == 'UNFORMATTED') then
        read(iunit,iostat=iostat) star%field, star%ccd, star%id, &
        star%ra_h, star%ra_m, star%ra_s, star%dc_sign, star%dc_d, star%dc_m, &
        star%dc_s, star%x, star%y, (star%col(icol)%data, star%col(icol)%err, &
        star%col(icol)%flg,icol=1,jcol)
      else
        print*, 'Attempting to read file of form ', form
        stop
      end if
      
      if (iostat == 0) then
        
        if (jcol < mcol) then
          do icol=jcol+1, mcol
            star%col(icol)%data=0.0
            star%col(icol)%err=0.0
            star%col%neg_flux=.false.
            star%col(icol)%flg='AA'
          end do
        end if

        do icol=1, ncol
          if (star%col(icol)%flg(2:2) == 'M') then
            ! This is where we sort out the ghastly way negative fluxes are
            ! are stored in cluster files.
            ! The data column stores the uncertainty as a magnitude, get
            ! it as a flux.
            ferr=10.0**(-0.4*star%col(icol)%data)
            ! The uncertaintly stores the flux as a magnitude, but with the
            ! wrong sign.  So first get the flux from it.
            flux=10.0**(0.4*star%col(icol)%err)
            ! And now set the magnitude.
            star%col(icol)%data=-star%col(icol)%err
            ! And is a negative flux.
            star%col(icol)%neg_flux=.true.
            ! Which is unflagged.
            star%col(icol)%flg(2:2)='O'
            ! Finally create an uncertainty in magnitude space.
            star%col(icol)%err= 2.5*log10((ferr+flux)/flux)
          else
            star%col(icol)%neg_flux=.false.
          end if
        end do

        if (trim(form) == 'FORMATTED') then 
          ! Sort out the field and ccd numbers.
          star%field=int(field_ccd)
          if (100*star%field - nint(100.0*field_ccd) == 0) then
            star%ccd=0
          else
            star%ccd=nint(100.0*(field_ccd-real(star%field)))
          end if
      
          read(dc_d,*) star%dc_d
      
          ! Now find all the ways a negative sign declination could have been set.
          star%dc_sign='+'
          if (star%dc_d < 0) star%dc_sign='-'
          if (star%dc_m < 0) star%dc_sign='-'
          if (star%dc_s < 0.0) star%dc_sign='-'
          if (dc_d(1:1) == '-') star%dc_sign='-'
      
          ! When all cluster programs flag negative declination through 
          ! star%dc_sign, we won't need this bit.
          if (star%dc_sign == '-') then
            star%dc_d=-1*abs(star%dc_d)
            star%dc_m=-1*abs(star%dc_m)
            star%dc_s=-1.0*abs(star%dc_s)
          end if

          ! Convert from old-style flags.
          do icol=1, ncol
            if (star%col(icol)%flg(2:2) == ' ') then
              star%col(icol)%flg(2:2)=star%col(icol)%flg(1:1)
              star%col(icol)%flg(1:1)='O'
            end if
            call flagconv(star%col(icol)%flg)
          end do

          ! Get rid of aged values which may crash formats.
          star%col%data=min(star%col%data, 99.999)
          star%col%err =min(star%col%err,  99.999)
          star%col%data=max(star%col%data, -99.999)
          star%col%err =max(star%col%err,  -99.999)

        else

          if (star%dc_sign == '-') then
            star%dc_d=-star%dc_d
            star%dc_m=-star%dc_m
            star%dc_s=-star%dc_s
          end if

        end if

      end if

      read_star=iostat

    end function read_star

    subroutine zero_star_array(star)

      type(a_star), intent(inout), dimension(:) :: star

      integer :: icol

      star%field=0
      star%ccd=0
      star%id=0
      star%ra_h=0
      star%ra_m=0
      star%ra_s=0.0
      star%dc_d=0
      star%dc_m=0
      star%dc_s=0.0
      star%x=0.0
      star%y=0.0
      do icol=1, mcol
        star%col(icol)%data=0.0
        star%col(icol)%err=0.0
        star%col(icol)%neg_flux=.false.
        star%col(icol)%flg='AA'
      end do
    
    end subroutine zero_star_array


    subroutine zero_star_one(star)

      type(a_star), intent(inout) :: star

      integer :: icol

      star%field=0
      star%ccd=0
      star%id=0
      star%ra_h=0
      star%ra_m=0
      star%ra_s=0.0
      star%dc_d=0
      star%dc_m=0
      star%dc_s=0.0
      star%x=0.0
      star%y=0.0
      do icol=1, mcol
        star%col(icol)%data=0.0
        star%col(icol)%err=0.0
        star%col(icol)%neg_flux=.false.
        star%col(icol)%flg='AA'
      end do
    
    end subroutine zero_star_one

    subroutine flux_to_mag(col)

      type(a_colour), intent(inout) :: col
      
      real :: flux, ferr, tinyflux

      ! The smallest absolute value of the flux that we can represent in 
      ! in magnitude space.  (If its negative it uses the M flag.) We never 
      ! want a magnitude as faint as 100th (since it will overrun our format 
      ! statements), which implies a flux greater than 10**-100/2.5. Nor do 
      ! we want a flux so near zero that when we take the log of it we get a 
      ! floating point error.
      tinyflux=max(100.0*tiny(1.0), 1.0e-40)
    
      flux=col%data
      ferr=col%err

      if (flux < 0.0) then
        col%neg_flux = .true.
      else
        col%neg_flux = .false.
      end if

      flux=abs(flux)
      flux=max(tinyflux, flux)
      ferr=max(tinyflux, ferr)

      col%data=-2.5*log10(flux)
      col%err = 2.5*log10((ferr+flux)/flux)

    end subroutine flux_to_mag


    subroutine mag_to_flux(col)

      type(a_colour), intent(inout) :: col
      
      real :: mag, err
    
      mag=col%data
      err=col%err

      if (-mag/2.5 >= log10(huge(col%err)/10.0)) then 
        col%data=huge(col%err)/10.0
      else
        col%data=10.0**(-mag/2.5)
      end if
      ! Evaluating it this way avoids foloating point overflows
      ! if err and mag are both large.  Such a condition happens
      ! when the flux is small.
      if ((err-mag)/2.5 >= log10(huge(col%err)/10.0)) then
        col%err=huge(col%err)/10.0
      else
        col%err= 10.0**((err-mag)/2.5) - col%data
      end if

      if (col%flg =='OM' .or. col%neg_flux) col%data=-1.0*col%data
      ! And remove the flag.
      if (col%flg(2:2) == 'M') col%flg(2:2)='O'

    end subroutine mag_to_flux

    subroutine flagconv(aflag)

      ! Converts the old numerical flags into the new character ones,
      ! assuming the old flags have been read as characters.

      character(len=2), intent(inout) :: aflag

      character, dimension(0:9) :: convert=(/'O', 'N', 'E', 'B', 'S', &
      'I', 'V', 'A', 'F', 'M'/)
      integer :: i, ichr
      character :: test

      do ichr=1, 2
        do i=0, 9
          write(test,'(i1)') i
          if (aflag(ichr:ichr) == test) aflag(ichr:ichr)=convert(i)
        end do
      end do

    end subroutine flagconv


    subroutine sort_star(star, nstars, sort)

      type(a_star), dimension(:), intent(inout) :: star
      integer, intent(in) :: nstars
      character(len=*), intent(in) :: sort

      integer :: i, j, l, ir
      type(a_star) :: swap
      real, dimension(:), allocatable :: sort_key
      real :: sort_swap

      if (nstars < 1) then

        print*, ' Error in s/r sort.  Number of stars is ', nstars
        stop

      else if (nstars > 1) then

        ! If nstars==1, then the data are already sorted!
        
        allocate(sort_key(nstars))
        if (sort == 'increasing_col(1)%data') then
          sort_key=star%col(1)%data
        else if (sort == 'increasing_ra') then
          sort_key = real(star%ra_h) + real(star%ra_m)/60.0 + star%ra_s/3600.0
        else if (sort == 'increasing_dec') then
          sort_key = real(star%dc_d) + real(star%dc_m)/60.0 + star%dc_s/3600.0
          where(star%dc_sign == '-') sort_key=-1.0*abs(sort_key)
        else
          print*, 'Error in s/r sort.  Undefined sort key ', trim(sort)
          stop
        end if

        ir=nstars
        l=ir/2+1
        do
          if (l > 1)then
            l=l-1
            swap=star(l)
            sort_swap = sort_key(l)
          else
            swap=star(ir)
            sort_swap=sort_key(ir)
            star(ir)=star(1)
            sort_key(ir)=sort_key(ir)
            IR=IR-1
            if (ir == 1) then
              star(1)=swap
              sort_key=sort_swap
              exit
            end if
          end if
          i=l
          j=l+l
20        IF(J.LE.IR)THEN
            IF(J.LT.IR)THEN
              IF(sort_key(j) < sort_key(j+1)) j=j+1
            ENDIF
            IF(sort_swap < sort_key(j))THEN
              star(i)=star(j)
              sort_key(i)=sort_key(j)
              I=J
              J=J+J
            ELSE
              J=IR+1
            ENDIF
          GO TO 20
          ENDIF
          star(i)=swap
          sort_key(i)=sort_swap
        end do

        deallocate(sort_key)

      end if

    end subroutine sort_star

    subroutine nxtcls_in(ifname)

      character(len=*) :: ifname

      pmt_nxtcls_in=.false.
      nam_nxtcls_in=ifname

    end subroutine nxtcls_in


    integer function start_reading_cluster_file(iunit, nstars, ncol, colstr)

      ! Opens a cluster format file, and reads the header.
      ! The header ends up in cluster_header.

      ! The unit to be used (one day this should be made clever -- se below).
      integer, intent(out) :: iunit
      ! The number of stars (set to -1 if its a formatted file), and the
      ! number of colours.
      integer, intent(out) :: nstars, ncol
      ! The names of the colours.
      character(len=*), dimension(:), intent(out) :: colstr

      logical :: there
      integer :: iostat, i, itrim
      character(len=10), dimension(mcol) :: unformatted_colstr

      character(len=120) :: work_header

      ! Set iunit.
      iunit=21

      if (pmt_nxtcls_in) then
        get_name: do 
          print*, '> Give input cluster file name.'
          read(*,'(a)') nam_nxtcls_in
          inquire(file=nam_nxtcls_in, exist=there)
          if (there) exit get_name
          print*, 'Cannot find file ', trim(nam_nxtcls_in)
          print*, 'Try again.'
        end do get_name
      else
        pmt_nxtcls_in=.true.
      end if

      open(iunit,file=nam_nxtcls_in,status='old',iostat=iostat)
      start_reading_cluster_file=iostat
      if (iostat /= 0) return
      
      ! Three lines of header.
      read(iunit,*,iostat=iostat) ncol

      ! Now we can establish of this is an unformatted file or not.
      if (iostat/=0 .or. ncol<0 .or. ncol>mcol) then
        ! Try unformatted.
        print*, 'Trying unformatted read for file ', trim(nam_nxtcls_in)
        close(iunit)
        open(iunit,file=nam_nxtcls_in,status='old',form='unformatted')
        read(iunit) cluster_header(2)
        read(iunit) cluster_header(1)
        read(iunit) cluster_header(3)
        read(iunit) ncol
        read(iunit) (unformatted_colstr(i),i=1,ncol)
        colstr(1:ncol)=unformatted_colstr(1:ncol)
        read(iunit) nstars
      else
        rewind(iunit)
        read(iunit,'(a120)') work_header
        read(work_header,*,iostat=iostat) ncol
        if (iostat < 0) ncol=4
        ! Now clear out any comments about the number of colours.
        itrim=index(work_header, '#')
        if (itrim < 1) then
          itrim=index(work_header, 'colours were created.')+21
        end if
        if (itrim < 0) itrim=0
        cluster_header(1)=work_header(itrim+1:len(work_header))
        ! The next line should have the names of the colours on it, but
        ! it may not.
        read(iunit,'(a80)') cluster_header(2)
        read(cluster_header(2),*,iostat=iostat) (colstr(i),i=1,ncol)
        if (iostat < 0) colstr=' '
        read(iunit,'(a80)') cluster_header(3)
        nstars=-1
      end if

      nam_lstcls_in=nam_nxtcls_in

    end function start_reading_cluster_file

    subroutine lstcls_in(ifname)

      character(len=*), intent(out) :: ifname
      
      ifname=nam_lstcls_in

    end subroutine lstcls_in


    integer function read_cluster_file(nstars, ncol, colstr, star) 

      ! Reads in a cluster format file.

      integer, intent(out) :: nstars, ncol
      character(len=*), dimension(:), intent(out) :: colstr
      type(a_star), dimension(:), allocatable, intent(out) :: star
      type(a_star) :: test_star

      integer :: iostat, istar, iunit
      character(len=15) :: form

      read_cluster_file=0

      iostat=start_reading_cluster_file(iunit, nstars, ncol, colstr)

      if (nstars == -1) then
        ! Its a formatted file, so we have to count the hard way.
        rewind(iunit)
        read(iunit,'(/,/)')
        nstars=0
        count_them: do
          iostat=read_star(iunit, test_star, ncol)
          if (iostat < 0) exit count_them
          nstars=nstars+1
        end do count_them
        rewind(iunit)
        read(iunit, '(/,/)')
      end if

      if (allocated(star)) deallocate(star)
      allocate(star(nstars))
      do istar=1, nstars
        iostat=read_star(iunit, star(istar), ncol)
        if (iostat /= 0) then
          print*, 'Error reading that file, iostat is', iostat
          print*, 'For star record ', istar
        end if
      end do

      close(iunit)

    end function read_cluster_file

    subroutine nxtcls_out(ifname)

      character(len=*) :: ifname

      pmt_nxtcls_out=.false.
      nam_nxtcls_out=ifname

    end subroutine nxtcls_out

    subroutine write_cluster_file(ncol, colstr, star, nstars)

      integer, intent(in) :: ncol
      character(len=*), dimension(:), intent(in) :: colstr
      type(a_star), dimension(:), intent(in) :: star
      integer, intent(in) :: nstars

      integer :: iunit, istar

      call start_cluster_file(ncol, colstr, iunit, nstars) 

      do istar=1, nstars
        call write_star(iunit, star(istar), ncol)
      end do

      close(iunit)

    end subroutine write_cluster_file

    subroutine start_cluster_file(ncol, colstr, iunit, nstars) 

      ! Prompts for a cluster file name and writes out the header.

      ! The file will be created as unformatted, if the optional argument
      ! nstars is there, and larger than the operating system variable 
      ! CLUSTER_FORMATTED_MAX_STARS.

      ! Need to add the ability to write extra information in lines
      ! 1 and 3 (optional arguments).

      integer, intent(in) :: ncol
      character(len=*), dimension(:), intent(in) :: colstr
      integer, intent(out) :: iunit
      ! If nstars is not given, it will never write the file unformatted.
      integer, optional :: nstars

      ! Locals.
      logical :: open
      character(len=8) :: date
      character(len=10) :: time
      character(len=80) ::date_time
      integer :: i, iostat
      character(len=80) :: outstr
      logical :: unformatted
      character(len=10), dimension(mcol) :: unformatted_colstr
      integer :: cfms

      ! Get a free unit number.
      find: do i=12, 100
        if (i == 100) then
          print*,'s/r start_cluster_file cannot find a free unit number.'
          stop
        end if
        iunit=i
        inquire(unit=iunit, opened=open)
        if (.not. open) exit find
      end do find

      if (pmt_nxtcls_out) then
        print*, '> Give output file name for catalogue.'
        read(*,'(a)') nam_nxtcls_out
      else
        pmt_nxtcls_out=.true.
      end if
      call date_and_time(date, time)
      date_time=' # File written at '//date(1:4)//'/'//date(5:6)//'/'//date(7:8) &
      //' '//time(1:2)//':'//time(3:4)//':'//time(5:10)

      call getenv('CLUSTER_FORMATTED_MAX_STARS', outstr)
      read(outstr, *, iostat=iostat) cfms
      if (iostat /= 0) cfms=huge(cfms)
      unformatted=.false.
      if (present(nstars)) then
        if (nstars > cfms) unformatted=.true.
      end if
          
      if (unformatted) then
        open(unit=iunit, file=nam_nxtcls_out, status='unknown', &
        form='unformatted')
        outstr=' '
        write(iunit) date_time
        write(iunit) cluster_header(1)
        write(iunit) cluster_header(3)
        write(iunit) ncol
        unformatted_colstr(1:ncol)=colstr(1:ncol)
        write(iunit) unformatted_colstr(1:ncol)
        write(iunit) nstars
      else
        open(unit=iunit, file=nam_nxtcls_out, status='unknown')
        write(iunit,*) ncol, 'colours were created. # ', cluster_header(1)
        write(iunit,*) (trim(colstr(i))//' ',i=1,ncol), date_time
        write(iunit,*) cluster_header(3)
      end if

    end subroutine start_cluster_file

    subroutine set_cluster_header(cluster_header_in)

      character(len=*), dimension(3), intent(in) :: cluster_header_in

      cluster_header=cluster_header_in

    end subroutine set_cluster_header

  end module define_star
