program sim

  ! Creates a simulated cluster CMD.

  use define_star
  use quad
  use cmdfit_subs

  implicit none

  type(a_star) :: star

  ! First set up the uncertainties to be used in each colour.  icorr is 
  ! between the errors.  Allowed values are:
  ! icorr =  1, for a CMD such as V, V-I;
  ! icorr = -1, for a CMD such as B, U-B;
  ! icorr =  0, gives you uncorrelated errors.
  ! The algebra is set up so that the uncertainties you supply are the total 
  ! uncertainty ignoring correlation (i.e. what most people quote).
  ! This means that for any icorr other than zero, the uncertainty in 
  ! magnitude must be less than that in colour.
  ! There are two ways to set up the uncertainties.  Choose which by
  ! setting this logical.
  logical, parameter :: fixed_uncer = .false.
  ! If this is set true, the uncertainties are equal for all data points, 
  ! and set by these two parameters.
  real, parameter :: fix_col_err=0.043, fix_mag_err=0.03
  ! Otherwise the magnitude uncertainty is taken to be of the form
  ! sqrt(mag_ind**2.0 (norm_10*10.0**(mag_norm*(true_mag-mag_shift)))**2.0),
  ! where mag_ind is a magnitude independent uncertainty which dominates at
  ! bright magnitudes, norm_10 and mag_norm are normalisation terms for
  ! what amounts to s sky-limit term, and mag_shift indicates at what 
  ! magnitude the sky noise begins to dominate.
  ! real, parameter :: mag_ind=0.005, norm_10=0.03, mag_norm=0.3, mag_shift=11.5
  real, parameter :: mag_ind=0.01, norm_10=0.03, mag_norm=0.4, mag_shift=9.5
  ! Finally, the uncertainty in colour is taken to be a constant times this.
  real, parameter :: col_ratio=1.3
  real :: col_err, mag_err
  integer :: icorr=0
  integer :: iband, iminus
  real :: age, dm, low_mass, up_mass

  integer :: istar, i
  real :: col, mag
  real :: true_col, true_mag
  real :: harvest
  real, external :: rnd_gauss
  character(len=50) :: ofname
  real :: v_err, i_err
  integer :: multiplicity
  integer :: imodel, icol
  real :: binary_fraction
  integer :: nstars
  real :: p_mass, p_teff
  character(len=20) :: flag
  real :: delta_b, delta_v
  real, parameter :: ubv_uncer=0.02

  print*, 'The available models are as follows.'
  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(*,*) imodel

  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(*,*) icol

  print*, '> Give binary fraction.'
  read(*,*) binary_fraction

  print*, '> Give output catalogue file name.'
  read(*,*) ofname
  open(unit=3, file=ofname)
  if (icol == 15) then
    write(3,*) '3 colours were created'
    write(3,*) 'V', ' ', trim(mag_name(icol)), ' ', trim(col_name(icol))
    write(3,*)
  else
    write(3,*) '2 colours were created'
    write(3,*) trim(mag_name(icol)), ' ', trim(col_name(icol))
    write(3,*)
  end if

  print*, '> Give the required age (in Myr).'
  read(*,*) age

  print*, '> Give the required distance modulus.'
  read(*,*) dm

  ! The range used in the Naylor and Jeffries paper was 0.36 1.9Mo
  print*, '> Give range of masses (low, high).'
  read(*,*) low_mass, up_mass

  ! The number of stars used in Naylor and Jeffries was 100.
  print*, 'How many stars do you want?'
  read(*,*) nstars

  ! Now lets sort out whether you add or subtract the colour from the
  ! magnitude to get the other magnitude.
  iminus=index(col_name(icol), '-')
  if (iminus == 0) then
    print*, 'Cannot find the minus sign in colour ', col_name(icol)
    stop
  end if
  iband=index(col_name(icol), trim(mag_name(icol)))
  if (iband == 0) then
    print*, 'Cannot find the magnitude ', trim(col_name(icol)), &
    ' in the colour ', trim(col_name(icol))
    print*, 'So assuming the colours are uncorrelated.'
    icorr=0
  else if (iband > iminus) then
    print*, 'Will add ', trim(mag_name(icol)), ' to ', trim(col_name(icol)), &
    ' to create ', col_name(icol)(1:iminus-1)
    icorr = -1
  else
    print*, 'Will subtract ', trim(col_name(icol)), ' from ', &
    trim(mag_name(icol)), &
    ' to create ', trim(col_name(icol)(iminus+1:len(col_name(icol))))
    icorr = 1
  end if

  open(unit=21, file='sim.iso')
  write(21,*) 'Columns are color, mag, primary mass, flag and multiplicity.'
  write(21,*)
  write(21,*)

  do istar=1, nstars

     call random_seed()
     call random_number(harvest)
     if (harvest < binary_fraction) then
       multiplicity=2
     else
       multiplicity=1
     end if

10   call make_star(age, icol, imodel, low_mass, up_mass, -0.35, &
     multiplicity, true_mag, true_col, p_mass, p_teff, flag)

     write(21,*) true_col, true_mag, p_mass, trim(flag), multiplicity
     if (trim(flag) /= 'OK') then
       print*, 'Star ', istar, ' of mass ', p_mass, ' has flag ', trim(flag)
       print*, 'Trying again.'
       goto 10
     end if

     if (icol == 15) then

       delta_b=rnd_gauss(0.0, ubv_uncer)
       delta_v=rnd_gauss(0.0, ubv_uncer)

       star%col(1)%data= v_for_ubv()+delta_v
       star%col(1)%err = ubv_uncer
       star%col(1)%flg = 'OO'

       star%col(2)%data= true_mag+delta_b-delta_v
       star%col(2)%err = ubv_uncer
       star%col(2)%flg = 'OO'

       star%col(3)%data= true_col+rnd_gauss(0.0, ubv_uncer)-delta_b
       star%col(3)%err = ubv_uncer
       star%col(3)%flg = 'OO'

       call write_star(3, star, 3)

     else

       ! Now, assume the CMD is of the form V vs V-I.  Then we must 
       ! correlate the uncertainties in the following way.
       ! First get the error in each magnitude.
       if (fixed_uncer) then
         mag_err=fix_mag_err
         col_err=fix_col_err
       else
         mag_err=sqrt(mag_ind**2.0 &
         + (norm_10*10.0**(mag_norm*(true_mag-mag_shift)))**2.0)
         col_err=col_ratio*mag_err
       end if
       v_err=rnd_gauss(0.0, mag_err)
       i_err=rnd_gauss(0.0, sqrt(col_err**2.0-(mag_err*real(icorr))**2.0))
       ! Now for the magnitude this is easy.
       mag=true_mag+v_err
       ! And for the colour.
       if (icorr == 0) then
         col=true_col+i_err
       else
         col=true_col+real(icorr)*(v_err-i_err)
       end if
       
       call zero_star(star)
 

       star%col(1)%data=mag+dm
       star%col(1)%err =mag_err
       star%col(1)%flg ='OO'
     
       star%col(2)%data=col
       star%col(2)%err =col_err
       star%col(2)%flg='OO'

     call write_star(3, star, 2)

   end if


   end do
   close(21)
   close(3)

end program sim
     
