program monte

  use ark_file_io
  use cmdfit_subs
  use colteff_subs

  implicit none

  integer :: iflag, colflag, modflag
  integer :: iage, nage
  real :: age, age_low, age_high, age_step
  integer, dimension(2) :: naxis
  real, allocatable, dimension(:,:) :: data, axdata
  character(len=50) :: ofname
  character(len=8) :: age_name
  integer :: i1, i2
  real :: col_start, col_end, mag_start, mag_end, low_mass, high_mass
  ! For a long time this was hard coded in each call as -0.1.
  real, parameter :: pwr_mass=-0.35
  integer :: useful_top, useful_bottom

  integer :: i, imag, icol, istar, multiple, icount
  real :: p_mass, p_teff
  ! The number of stars in the simulation.
  integer :: nstar
  ! The number of stars in the simple isochrone.
  integer, parameter :: niso=10000
  ! The size of the pixels in magnitude.
  real, parameter :: pix_size=0.0025
  ! The binary fraction.
  real :: binary_fraction
  ! More edges.
  real :: mag_min, mag_max, col_min, col_max

  real, dimension(:), allocatable :: star_mag, star_col
  character(len=50) :: isoflag
  logical, dimension(:), allocatable :: isgood
  ! Allow extraoplation in colour-Teff relationship?  You may want to do
  ! this for luminous evolved stars as the gravity grids run out.
  ! Its only allowed to extraoplate in gravity, and this currently only works
  ! for the Bessell atmospheres.
  logical, parameter :: extrap_col_teff=.false.
  !call setbug()

  ! 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 binary fraction do you want?'
  read(*,*) binary_fraction

  print*, 'What range of log10(ages) do you want?'
  read(*,*) age_low, age_high
  if (abs(age_low-age_high) > 2.0*tiny(age_low)) then    
    print*, 'What step in age do you want?'
    read(*,*) age_step
    nage=1+nint((age_high-age_low)/age_step)
    if (abs(age_high-age_low-real(nage-1)*age_step)/age_step > 0.0001) then
      print*, 'Age step does not make sense.  You require ', &
      nint(age_high-age_low)/age_step, ' age steps.' 
      stop
    end if
  else
    nage=1
  end if

  each_age: do iage=1, nage

    if (nage == 1) then
      age = age_low
    else
      age = age_low + real(iage-1)*(age_high-age_low)/real(nage-1)
    end if

    print*, 'Doing log10(age)', age


    ! Now calculate the number of stars required.
    nstar=1000000.0
    print*, 'With ', nstar, ' simulated objects.'
    allocate(star_mag(nstar), star_col(nstar), isgood(nstar))

    ! 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,*)

    call random_seed()

    create_each_star: do istar=1, nstar

      if (istar - 100000*nint(real(istar)/100000) == 0) &
      print*, 'Created ', 100.0*real(istar)/real(nstar), ' percent of stars.'

      multiple=1
      if (istar > nint(real(nstar)*(1.0-binary_fraction))) multiple=2

      call make_star(age, colflag, modflag, low_mass, high_mass, &
      pwr_mass, multiple, star_mag(istar), star_col(istar), p_mass, &
      p_teff, isoflag)
      !call make_star_ubv(age, modflag, low_mass, high_mass, &
      !pwr_mass, multiple, mag, col)
      isgood(istar)=.false.
      if (trim(isoflag)=='OK') isgood(istar)=.true.
      if (trim(isoflag)=='Colour-Teff (grav)' .and. extrap_col_teff) &
      isgood(istar)=.true.
      if (istar <= niso) write(12,*) star_col(istar), star_mag(istar), &
      p_mass, p_teff, trim(isoflag), isgood(istar)
      !print*, istar, star_mag(istar), star_col(istar), isgood(istar)
    end do create_each_star
    close(12)

    ! Pack the arrays.
    icount=0
    do istar=1, nstar
      if (isgood(istar)) then
        icount=icount+1
        star_mag(icount)=star_mag(istar)
        star_col(icount)=star_col(istar)
      end if
    end do
    nstar=icount
    print*, 'There were ', nstar, ' usable stars.'

    ! Find the size of image needed.    
    mag_start=maxval(star_mag(1:nstar))
    mag_end  =minval(star_mag(1:nstar))
    ! Its often best to be a little wider in colour.
    col_start=minval(star_col(1:nstar))-0.2
    col_end  =maxval(star_col(1:nstar))+0.2
    ! Round these to "nice" values.  This makes the grids look nicer, since
    ! the aliasing effects are minimised.
    mag_start=real(int(mag_start/pix_size)  )*pix_size
    mag_end  =real(int(mag_end  /pix_size)+1)*pix_size
    col_start=real(int(col_start/pix_size)+1)*pix_size
    col_end  =real(int(col_end  /pix_size)  )*pix_size
    naxis(1)=nint((col_end-col_start)/pix_size)+1
    naxis(2)=nint((mag_start-mag_end)/pix_size)+1
    print*, 'Grid will be ', naxis
    allocate(data(naxis(1),naxis(2)), axdata(maxval(naxis),2))
    print*, 'Covering mags ', mag_start, ' to ', mag_end
    print*, 'and colours ', col_start, ' to ', col_end
    do i=1, naxis(1)
      axdata(i,1)=col_start+(col_end-col_start)*real(i-1)/real(naxis(1)-1)
    end do
    do i=1, naxis(2)
      axdata(i,2)=mag_start+(mag_end-mag_start)*real(i-1)/real(naxis(2)-1)
    end do

    do istar=1, nstar

      if (istar - 100000*nint(real(istar)/100000) == 0) &
      print*, 'Placed ',  100.0*real(istar)/real(nstar), &
      ' percent of stars in grid.'

      ! Find the pixel number.
      imag=minloc(abs(axdata(1:naxis(2),2)-star_mag(istar)),1)
      icol=minloc(abs(axdata(1:naxis(1),1)-star_col(istar)),1)  

      ! Put the star into the image.
      data(icol,imag)=data(icol,imag)+1.0

    end do

    ! Make the sum over all area (including that not on the grid) one.
    data=data/((1.0+binary_fraction)*real(nstar))

    ! First write out the image file.
    call put_header_i('MODFLAG', modflag, 'Whose isochrones' , 1)
    call put_header_i('COLFLAG', colflag, 'Which colour' , 1)
    call put_header_r('BINFRAC', binary_fraction, 'Binary star fraction' , 1)
    call nxtark_out(ofname)
    iflag=makark(naxis, data, axdata)

    deallocate(data, axdata, star_mag, star_col, isgood)

  end do each_age

end program monte
