program iso

  use ark_file_io
  use cmdfit_subs
  use colteff_subs

  implicit none

  integer :: iflag, colflag, modflag, ipos
  integer :: iage, nage
  real :: age, age_low, age_high
  integer, dimension(2) :: naxis
  real, allocatable, dimension(:,:) :: data, axdata
  real, allocatable, dimension(:) :: grad
  character(len=50) :: ofname
  character(len=8) :: age_name
  integer :: i1, i2
  real :: mass, low_mass, high_mass
  integer :: useful_top, useful_bottom
  ! The image, and the associated values in axis1 (colour) and axis2 (mag).
  type a_grader
    real, dimension(2) :: low, high
    integer :: npts
  end type a_grader
  type(a_grader), allocatable, dimension(:) :: grader

  integer :: i, imag, icol, istar
  real :: col, mag, p_mass, p_teff
  character(len=50) :: isoflag
  ! The number of stars in the simulation.
  integer, parameter :: nstar=10000
  ! Are we to use the natural normalisation?  If so make sure the integral
  ! of the model is one and don't write out the gradient file.  The code
  ! will force the use of the natural normalisation for colour-colour work.
  ! If you set this false you get the old chi^2-like normalisation.
  logical :: natural=.true.
  ! More edges.
  real :: mag_min, mag_max, col_min, col_max


  ! This used to be hard coded 0.04 3.0
  print*, 'What range of masses do you want (low, high)?'
  read(*,*) low_mass, high_mass
  
  print*, 'The following models are available.'
  do i=1, n_models
    if (len_trim(mod_name(i)) > 0) print*, i, trim(mod_name(i))
  end do
  print*, 'What model number do you want?'
  read(*,*) modflag

  print*, 'The available CMDs are '
  do i=1, n_names
    if (len_trim(mag_name(i)) > 0) & 
    print*, i, trim(mag_name(i)), ' vs ', trim(col_name(i)), &
    '   ('//trim(pht_comm(i))//')'
  end do
  print*, 'What colour number do you want?'
  read(*,*) colflag

  print*, 'What age do you want?'
  read(*,*) age


    ! Write out a simple isochrone file.
    call iso_image_name(modflag, colflag, age, ofname)
    open (unit=12, file=trim(ofname)//'.iso')
    write(12,*) 'Colour, magnitude, mass, teff, flag.'
    write(12,*)
    write(12,*)

    ! We need to calculate the grader array, so we create an array of
    ! single stars whatever the binary fraction is.
    each_single_star: do istar=1, nstar

      mass=low_mass + (high_mass-low_mass)*real(istar)/real(nstar)

      call make_star(age, colflag, modflag, mass, mass, -0.1, &
      1, mag, col, p_mass, p_teff, isoflag)
      write(12,*) col, mag, p_mass, p_teff, trim(isoflag)

    end do each_single_star
    close(12)

end program iso
