program grid

  use define_star
  use cmdfit_subs
  use ark_file_io
  use colteff_subs

  implicit none

  integer, parameter :: n_par=3
  type(a_star), dimension(:), allocatable :: star
  real, dimension(n_par) :: a_par, best_par
  integer :: iflag, n_term, colflag, icol, modflag
  integer :: iostat, i
  integer :: points
  integer, dimension(2) :: naxis, jaxis
  integer, dimension(2) :: var_par
  real, dimension(2) :: start, end
  real, dimension(:,:), allocatable :: data, axdata, rnpts
  character(len=80), dimension(2) :: ctype
  real :: best_value, age_step
  type(a_star) :: test
  integer :: icol1, icol2
  integer :: ncol
  character(len=10), dimension(:), allocatable :: bands
  character(len=50) :: ofname, ifname
  real, dimension(:), allocatable :: best_distrib

  type(a_star), dimension(:), allocatable :: dat

  integer :: nleft
  real, allocatable, dimension(:) :: sorted_distrib
  real :: clipped_tau
  real :: soft_clip

  integer :: iminus, iband, correlated

  real :: mag_shift, colour_shift, global_mag_shift, pix_size, &
  global_colour_shift
  integer :: i_mag_shift
  real :: sys_mag, sys_col
  real :: low_tau

  ! Mode is 1 for searching in age and/or distance modulus.
  ! Anything else for searching in distance modulus and/or extinction.
  real, parameter :: mode=1

  print*, '> Give input catalogue name.'
  read(*,'(a)') ifname
  open(unit=1, file=ifname, action='read')
  read(1,*) ncol
  allocate(bands(ncol))
  read(1,*) bands
  read(1,*)

  print*, '* Will fit magnitude ', bands(1)
  if (ncol == 2) then
    print*, '* Against colour ', bands(2)
    icol1=1
    icol2=2
  else
    print*, '* Available magnitudes are '
    do icol=1, ncol
      print*, '* ', icol, bands(icol) 
    end do
    print*, '> Give the number of the one you want.'
    read(*,*) icol1
    print*, '* Available colours are '
    do icol=1, ncol
      print*, '* ', icol, bands(icol) 
    end do
    print*, '> Give the number of the one you want.'
    read(*,*) icol2
  end if
     

  ! O.K., lets sort out the colours being fitted.
  colflag=-1
  find_col: do icol=1, n_names
    if (trim(col_name(icol)) == trim(bands(icol2))) then
      if (trim(mag_name(icol)) == trim(bands(icol1))) then
        print*, 'Fitting magnitude ', trim(mag_name(icol)), &
        ' against colour ', trim(col_name(icol)), &
        '  ('//trim(pht_comm(icol))//')'
        colflag=icol
        exit find_col
      end if
    end if
  end do find_col
  if (colflag == -1) then
    print*, 'Cannot find colours ', bands(icol1), ' ', bands(icol2)
    stop
  end if
  ! print*, 'Fitting colour flag ', colflag

  ! Now lets sort out whether you add or subtract the colour from the
  ! magnitude to get the other magnitude.
  if (index(bands(icol1), '-') == 0) then
    ! This is truly a magnitude.
    iminus=index(bands(icol2), '-')
    if (iminus == 0) then
      print*, 'Cannot find the minus sign in colour ', bands(icol2)
      stop
    end if
    iband=index(bands(icol2), trim(bands(icol1)))
    if (iband == 0) then
      print*, 'Cannot find the magnitude ', trim(bands(icol2)), &
      ' in the colour ', trim(bands(icol2))
    else if (iband > iminus) then
      print*, 'Will add ', trim(bands(icol1)), ' to ', trim(bands(icol2)), &
      ' to create ', bands(icol2)(1:iminus-1)
      correlated = 1
    else
      print*, 'Will subtract ', trim(bands(icol2)), ' from ', &
      trim(bands(icol)), &
      ' to create ', trim(bands(icol2)(iminus+1:len(bands(icol2))))
      correlated = -1
    end if
  else
    ! Its colour-colour data.
    correlated=0
  end if

  print*, '> Give the extra uncertainty to be added to each magnitude.'
  read(*,*) sys_mag
  print*, '> Give the extra uncertainty to be added to each colour.'
  read(*,*) sys_col

  ! And the model.
  print*, 'The available models are '
  do i=1, n_models
    if (len_trim(mod_name(i)) > 0) print*, i, trim(mod_name(i))
  end do
  print*, '> Give the number for the model.'
  read(*,*) modflag

  points=0
  do
    iostat=read_star(1, test, ncol)
    if (iostat < 0) exit
    if (iostat > 0) then
      print*, 'Error reading file ', trim(ifname)
      stop
    end if
    if (test%col(icol1)%flg /= 'OO') cycle
    if (test%col(icol2)%flg /= 'OO') cycle
    points=points+1
  end do
  allocate(star(points), best_distrib(points))

  rewind(1)
  read(1,'(/,/)')
  i=1
  do 
    iostat=read_star(1, test, ncol)
    if (iostat < 0) exit
    if (test%col(icol1)%flg /= 'OO') cycle
    if (test%col(icol2)%flg /= 'OO') cycle
    star(i)=test
    i=i+1
  end do
  print*, 'Read ', points, ' of ', trim(bands(icol1)), ' vs ', &
  trim(bands(icol2))

  if (mode == 1) then
    ! Age.
    var_par(1)=1
    print*, '> Give range of log10(ages) to be searched.'
    read(*,*) start(1), end(1)
    if (abs(start(1)-end(1)) > 2.0*tiny(end(1))) then
      print*, '> And the log10(age) step to be used.'
      read(*,*) age_step
      naxis(1)=1+nint((end(1)-start(1))/age_step)
      if (abs(end(1)-start(1)-real(naxis(1)-1)*age_step)/age_step > 0.0001) then
        print*, 'Age step does not make sense.  You require ', &
        nint(end(1)-start(1))/age_step, ' age steps.' 
        print*, end(1)-start(1), real(naxis(1)-1)*age_step
        stop
      end if
    else
      naxis(1)=1
    end if
  else 
    ! The extinction in the colour in question.
    var_par(1)=3
    print*, '> Give the range of reddening to search in ', &
    trim(col_name(colflag))
    read(*,*) start(1), end(1)
    print*, '> And the number of reddening points in the grid.'
    read(*,*) naxis(1)
  end if
  ! The distance modulus.
  print*, '> Give the range in distance modulus to be searched.'
  read(*,*) start(2), end(2)
  print*, '> And the number of distance modulus points in the grid.'
  read(*,*) naxis(2)
  var_par(2)=2

  if (mode == 1) then
    ! The extinction in the colour in question.
    print*, '> Give reddening in ', trim(col_name(colflag))
    read(*,*) a_par(3)
  else
    ! Age
    print*, '> Give the log10(age).'
    read(*,*) a_par(1)
  end if

  ! The tau^2 at which data points get clipped out.
  print*, '> Give maximum allowed range of tau^2 (i.e. value of the soft clipping).'
  read(*,*) soft_clip

  allocate(data(naxis(1),naxis(2)), axdata(maxval(naxis),2), &
  rnpts(naxis(1),naxis(2)))

  n_term=2

  allocate(dat(points))

  do i=1, points
    ! Increase the uncertainties by the systematic.
    star(i)%col(icol1)%err  = sqrt(star(i)%col(icol1)%err**2.0 + sys_mag**2.0)
    ! Bug corrected here August 11 2006. 
    ! Used to be...
    ! star(i)%col(2)%err  = sqrt(star(i)%col(icol2)%err**2.0 + 2.0*(sys**2.0))
    ! The changed in the final stages of the NGC2169 paper (1st September 2006)
    ! to allow different uncertainties in mag and colour.
    star(i)%col(icol2)%err  = sqrt(star(i)%col(icol2)%err**2.0 + sys_col**2.0)
    ! Copy the star structure accross, which deals with colour 1.
    dat(i)=star(i)
    ! And now make sure colours are right.
    dat(i)%col(1)%data = star(i)%col(icol1)%data
    dat(i)%col(1)%err  = star(i)%col(icol1)%err
    dat(i)%col(2)%data = star(i)%col(icol2)%data
    dat(i)%col(2)%err  = star(i)%col(icol2)%err
  end do

  if (sum(dat%col(1)%err) >= 0.99*sum(dat%col(2)%err)) then
    print*, 'The uncertanties in magnitude are larger than in colour so'
    print*, 'I will take them as uncorrelated.'
    correlated = 0
  end if

  !  write(20,*) a_par
  !  write(20,*) colflag
  !  write(20,*) modflag
  !  write(20,*) correlated
  !  write(20,*) start
  !  write(20,*) end
  !  write(20,*) soft_clip

  call grid2d(a_par, dat, colflag, modflag, correlated, &
  var_par, start, end, &
  soft_clip, data, axdata, best_value, best_par, best_distrib, rnpts)

  low_tau=minval(best_distrib)

  ofname='grid.fit'
  call typark(2)
  call nxtark_out(ofname)
  jaxis=naxis
  ctype='NOT SET'
  if (mode == 1) then
    if (naxis(1) > 1) ctype(1)='LOG10(AGE)'
  else
    if (naxis(1) > 1) ctype(1)='EXTINCTION'
  end if
  if (naxis(2) > 1) ctype(2)='DISTANCE MODULUS'
  if (trim(ctype(1)) /= 'NOT SET') call put_header_s('CTYPE1', ctype(1), ' ', 1)
  if (trim(ctype(2)) /= 'NOT SET') call put_header_s('CTYPE2', ctype(2), ' ', 1)

  call put_header_r('HI_TAU', high_calc_tau(), 'HIGHEST CALCULABLE TAU^2', 1)

  if (naxis(2) == 1) then
    open(unit=11, file=trim(ofname)//'.asc')
    write(11,'(/,/)')
    do i=1, naxis(1)
      write(11,*) axdata(i,1), data(i,1)
    end do
    close(11)
  else if (naxis(1) == 1) then
    open(unit=11, file=trim(ofname)//'.asc')
    write(11,'(/,/)')
    do i=1, naxis(2)
      write(11,*) axdata(i,2), data(1,i)
    end do
    close(11)
    jaxis(1)=0
  end if
  i=makark(jaxis, data, axdata)

  ofname='grid_npts.fit'
  call nxtark_out(ofname)
  i=makark(jaxis, rnpts, axdata)

  print*, 'Lowest value ', best_value
  do i=1, 2
    print*, ' '
    if (ctype(i)(1:5) == 'LOG10') then
      print*, ' For parameter ', trim(ctype(i)(6:len(ctype)))//'/10^6'
      print*, ' Range searched ', 10.0**(axdata(1,i)-6.0), &
      10.0**(axdata(naxis(i),i)-6.0)
      print*, ' Best value ', 10.0**(best_par(var_par(i))-6.0)
    else
      print*, ' For parameter ', trim(ctype(i))
      print*, ' Range searched ', axdata(1,i), axdata(naxis(i),i)
      print*, ' Best value ', best_par(var_par(i))
    end if
  end do

  call iso_image_name(modflag, colflag, best_par(1), ifname)
  call nxtark_in(ifname)
  iflag=inpark(naxis, data, axdata)

  ! Write a file of the tau^2s
  open (unit=66, file='tau_point.out')
  write(66,*) 'Columns are field, id, colour, magnitude and tau^2'
  write(66,'(/)')
  do i=1, points
    ! If the value of tau is more than half that of a fully clipped point,
    ! clip it.
    write(66,*) dat(i)%field, dat(i)%id, dat(i)%col(2)%data, &
    dat(i)%col(icol1)%data, best_distrib(i)
  end do
  close(66)

  ! Find the number of unclipped point.
  if (soft_clip > 0.0) then
    nleft=0
    clipped_tau=0.0
    do i=1, points
      ! If the value of tau is more than half that of a fully clipped point,
      ! clip it.
      if (best_distrib(i) < low_tau+soft_clip/2.0) then
        nleft=nleft+1
        clipped_tau=clipped_tau+best_distrib(i)
      end if
    end do

    print*, 'After removing clipped points, tau is ', clipped_tau
    print*, 'For ', nleft, ' data points.'
    print*, '(There were ', points, ' before clipping.)'

    ! Now allocate an array, and fill it with the values of tau.
    allocate(sorted_distrib(nleft))
    nleft=0
    do i=1, points
      if (best_distrib(i) < low_tau+soft_clip/2.0) then
        nleft=nleft+1
        sorted_distrib(nleft)=best_distrib(i)
      end if
    end do
  else
    nleft=points
    allocate(sorted_distrib(nleft))
    sorted_distrib=best_distrib
  end if

  ! Now sort them.
  call tau_sort(sorted_distrib)

  ! Write out the unclipped points in apparent magnitude.
  open(unit=32, file='unclipped.cat')
  write(32,*) ncol
  write(32,*) bands(1:ncol)
  write(32,*)
  do i=1, points
    if (soft_clip > 0.0) then
      if (best_distrib(i) < low_tau+soft_clip/2.0) &
      call write_star(32, star(i), ncol)
    else
      call write_star(32, star(i), ncol)
    end if
  end do
  close(32)


  ! Change to absolute magnitude and write out the unclipped points.
  bands(icol1)='('//trim(bands(icol1))//')o'
  bands(icol2)='('//trim(bands(icol2))//')o'
  open(unit=32, file='unclipped_abs.cat')
  write(32,*) ncol
  write(32,*) bands(1:ncol)
  write(32,*)
  do i=1, points
    call reddening(colflag, -best_par(2), -best_par(3), &
    star(i)%col(icol2)%data, star(i)%col(icol1)%data)
    if (soft_clip > 0.0) then
      if (best_distrib(i) < low_tau+soft_clip/2.0) &
      call write_star(32, star(i), ncol)
    else
      call write_star(32, star(i), ncol)
    end if
  end do
  close(32)


  ! Change to absolute magnitude and write out the unclipped points.
  ! This is complex as the magnitude shift is colour dependent.
  ! The following only does the job to the nearest half pixel.
  ! Find the magnitude shift at the middle pixel.
  global_colour_shift=(axdata(1,1)+axdata(naxis(1),1))/2.0
  global_mag_shift=0.0
  call reddening(colflag, best_par(2), best_par(3), &
  global_colour_shift, global_mag_shift)
  ! Apply it.
  axdata(1:naxis(2),2)=axdata(1:naxis(2),2)+global_mag_shift
  do i=1, naxis(1)
    ! For each colour, work out the shift from nominal.
    colour_shift=axdata(i,1)
    mag_shift=0.0
    call reddening(colflag, best_par(2), best_par(3), &
    colour_shift, mag_shift)
    mag_shift=mag_shift-global_mag_shift
    pix_size=(axdata(naxis(2),2)-axdata(1,2))/real(naxis(2)-1)
    i_mag_shift=nint(mag_shift/pix_size)
    if (i_mag_shift > 0) then
      data(i,1+i_mag_shift:naxis(2))=data(i,1:naxis(2)-i_mag_shift)
      data(i,1:i_mag_shift)=0.0
    else if (i_mag_shift < 0) then
      data(i,1:naxis(2)+i_mag_shift)=data(i,1-i_mag_shift:naxis(2))
      data(i,naxis(2)+i_mag_shift+1:naxis(2))=0.0
    end if
  end do
  ! Apply the colour shift..
  global_colour_shift=global_colour_shift-(axdata(1,1)+axdata(naxis(1),1))/2.0
  axdata(1:naxis(1),1)=axdata(1:naxis(1),1)+global_colour_shift

  call nxtark_out('best_model.fit')
  iflag=makark(naxis, data, axdata)

  open(unit=22, file='distrib.tau')
  write(22, '(/,/)')
  write(22,*) 0.0, 1.0
  do i=1, nleft
    write(22,*) sorted_distrib(i), real(nleft-i+1)/real(nleft)
    write(22,*) sorted_distrib(i), real(nleft-i)/real(nleft)
  end do
  close(22)

end program grid
