module cmdfit_subs

  use colteff_subs

  implicit none

  real, private :: extra_mag, high_tau=-1.0

  contains

  real function high_calc_tau()

    high_calc_tau=high_tau/2.0

  end function high_calc_tau

  subroutine natural_norm(mag_max, mag_min, data, axdata)
    
    real, intent(in) :: mag_max, mag_min
    real, intent(inout), dimension(:,:) :: data, axdata

    integer :: npix, imag
        
    ! Fint the number of pixels we must sum over (avoids the number
    ! of pixels changing by +/-1 as you move over pixel boundaries).
    npix = abs(nint(real(size(data,2)) * &
    (mag_max-mag_min)/(axdata(size(data,2),2)-axdata(1,2))))
    ! Find the nearest position in mag, avoiding problems if
    ! magnitude axes run in the wrong direction).
    imag=min( minval(minloc(abs(mag_min-axdata(1:size(data,2),2)))), &
    minval(minloc(abs(mag_max-axdata(1:size(data,2),2)))))
    if (imag < 1) then
      npix=imag+npix-1
      imag=1
    end if
    if (imag+npix-1 > size(data,2)) npix=size(data,2)+1-imag
    if (imag<1 .or. imag+npix-1>size(data,2)) then
      print*, 'Problem ', imag, imag+npix-1
      stop
    end if
    data=data/sum(data(:,imag:imag+npix-1))

  end subroutine natural_norm

    subroutine make_star(age, colflag, modflag, low_mass, up_mass, &
    pwr_mass, multiplicity, mag, col, p_mass, p_teff, flag_out)

    ! The flags indicating which age, colour and model to use.
    real, intent(in) :: age
    integer, intent(in) :: colflag, modflag
    ! For this model the range of mass and power-law distribution
    ! in mass to be simulated.
    real, intent(in) :: up_mass, low_mass, pwr_mass 
    ! dN/dlog(M) = M**(-1 + pwr_mass), or dN/dM = M**(-2 + pwr_mass), or
    ! N = M**(-1 + pwr_mass).
    ! Thus for a Salpeter mass function (-1+pwr_mass)=-1.35; pwr_mass=-0.35
    ! For a 40Myr DAM isochrone, the following numbers give roughly
    ! evenly spaced stars over an Mv range of 2 to 11.
    ! real, parameter :: up_mass=1.9, low_mass=0.36, pwr_mass=-1.5
    ! Flag indicating multiplicity (1=single, 2=binary, 3=triple etc).
    integer :: multiplicity

    ! The colour and magnitude of the stars created. 
    ! The mass and effective temperature of the primary star.
    real, intent(out) :: mag, col
    real, intent(out) :: p_mass, p_teff
    ! A character flag, 'OK' if things are.
    character(len=*), intent(out) :: flag_out

    ! Locals.
    logical, save :: first=.true.
    character(len=50) :: isofil
    real, save :: col1, col2, mag1, mag2
    real :: harvest, grad
    real, dimension(multiplicity) :: mass, uflux, bflux, vflux, iflux
    real :: umb, bmv, vmi, v
    integer :: i, iminus
    logical :: v_vmi_like, v_bmv_like, umb_bmv_like
    character(len=20) :: flag, work_flag

    if (modflag == 0) then

      ! This is test code, which creates a linear isochrone.

      if (multiplicity /= 1) then
        print*, 'Cannot create a binary sequence for a linear isochrone.'
        stop
      end if

      if (first) then
        print*, '> Give the name of the isochrone (.iso) file.'
        read(*,*) isofil
        open(unit=10, file=isofil, action='read')
        ! Rob's 6 lines of header.
        read(10,'(/,/,/,/,/)')
        read(10,*) col1, mag1
        read(10,*) col2, mag2
        close(10)
        ! Now, chop a magnitude off the start and end points to allow
        ! for a range of possible distance modulii when fitting.
        grad=(col1-col2)/(mag1-mag2)
        mag1=mag1+1.0
        col1=col1+grad
        mag2=mag2-1.0
        col2=col2-grad
        first=.false.
      end if

      call random_number(harvest)
      ! Calculate the start and range of the isochrone in colour space
      ! space.
      col=col1+(col2-col1)*harvest
      ! And now in magnitude (remembering to reverse the gradient).
      mag=mag1+(mag2-mag1)*harvest

    else

      ! Create a real isochrone.       

      ! First find out if the colours are of the form V vs V-I or V vs B-V.
      umb_bmv_like=.false.
      v_vmi_like=.false.
      v_bmv_like=.false.
      ! First is this colour/colour?
      iminus=index(mag_name(colflag), '-') 
      if (iminus > 0) then
        ! O.K., is it in the correct order?
        if (trim(mag_name(colflag)(iminus+1:len(mag_name(colflag)))) /= &
          col_name(colflag)(1:index(col_name(colflag), '-')-1)) then
          print*, 'Colour not of form U-B vs B-V.'
          stop
        else
          umb_bmv_like=.true.
        end if
      else
        iminus=index(col_name(colflag), '-') 
        if (iminus == 0) then
          print*, 'Cannot find the minus in ', col_name(colflag)
          stop
        else
          if (trim(mag_name(colflag)) == col_name(colflag)(1:iminus-1)) then
            v_vmi_like=.true.
          else if (trim(mag_name(colflag)) == &
            trim(col_name(colflag)(iminus+1:))) then
            v_bmv_like=.true.
          else
            print*, 'Cannot find ', trim(mag_name(colflag)), ' in ', &
            trim(col_name(colflag))
            stop
          end if
        end if
      end if


100   each_component: do i=1, multiplicity

        call random_number(harvest)
        if (i == 1) then
          mass(1) = (low_mass**pwr_mass - &
          harvest*(low_mass**pwr_mass-up_mass**pwr_mass))**(1.0/pwr_mass)
          ! Work out the primary's temperature.
          p_teff=10.0**isochr('mas', mass(i), 'tef', age, colflag, modflag, flag)
        else
          mass(i) = mass(i-1)*harvest
          ! Hack to cut off at 0.25Mo for simulation of effect of
          ! binary wedge for original tau^2 paper.
          !if (mass(i) < 0.25) then
          !  flux1(i)=0.0
          !  flux2(i)=0.0
          !  cycle each_component
          !end if
        end if
        if (v_vmi_like .or. v_bmv_like) then
          v=isochr('mas', mass(i), 'mag', age, colflag, modflag, flag)
        else if (umb_bmv_like) then
          ! In the colour-colour case get V from another isochrone.
          v=isochr('mas', mass(i), 'mag', age, colflag-1, modflag, flag)
        end if
        ! Now, if the flag is not OK, continue, but note the fact.
        work_flag=flag
        if (v_bmv_like .or. umb_bmv_like) then
          bmv=isochr('mas', mass(i), 'col', age, colflag, modflag, flag)
          vmi=0.0
        else if (v_vmi_like) then
          vmi=isochr('mas', mass(i), 'col', age, colflag, modflag, flag)
          bmv=0.0
        end if
        if (flag /= 'OK') work_flag=flag
        if (umb_bmv_like) then
          ! In the colour-colour case U-B is stored as a magnitude.
          umb=isochr('mas', mass(i), 'mag', age, colflag, modflag, flag)
          if (flag /= 'OK') work_flag=flag
        else
          umb=0.0
        end if
        if (i == 1) then
          flag_out=work_flag
        else
          ! This is a secondary star.  
          if (trim(work_flag) /= 'OK') then
            ! Only possible error is that its dropped 
            ! off the bottom of the isochrone.  
            ! Set fluxes already to zero.
            uflux(i)=0.0
            bflux(i)=0.0
            vflux(i)=0.0
            iflux(i)=0.0
            cycle each_component
          end if
        end if
        uflux(i)=10.0**(-0.4*(v+bmv+umb))
        bflux(i)=10.0**(-0.4*(v+bmv))
        vflux(i)=10.0**(-0.4*v)
        iflux(i)=10.0**(-0.4*(v-vmi))

      end do each_component

      ! Add up all the components.
      if (v_vmi_like) then
        mag=-2.5*log10(sum(vflux))
        col=-2.5*log10(sum(vflux)/sum(iflux))    
      else if (v_bmv_like) then
        mag=-2.5*log10(sum(vflux))
        col=-2.5*log10(sum(bflux)/sum(vflux))
      else if (umb_bmv_like) then
        mag=-2.5*log10(sum(uflux)/sum(bflux))
        col=-2.5*log10(sum(bflux)/sum(vflux))
        extra_mag=-2.5*log10(sum(vflux))
      end if
      p_mass=mass(1)

    end if

  end subroutine make_star

  real function v_for_ubv()

    v_for_ubv=extra_mag

  end function v_for_ubv

  subroutine iso_image_name(modflag, colflag, age, filnam)

    integer, intent(in) :: modflag, colflag
    real :: age
    character(len=*), intent(out) :: filnam

    character(len=8) :: age_name

    write(age_name, 30) int(age+tiny(age)), nint(100*(age-int(age+tiny(age))))
30  format(i5.5, '.', i2.2)

    filnam=trim(mod_name(modflag))//'_'// &
           trim(mag_name(colflag)) //'_'// &
           trim(col_name(colflag)) //'_'// &
                age_name          //'.fit'

  end subroutine iso_image_name

  real function likelihood(data, axdata, grad, icol, imag, sig_col, &
    sig_mag, correlated, flag, sum_gauss)
    
    ! The unsmoothed image and axes.
    real, dimension(:,:), intent(in) :: data, axdata
    real, dimension(:), allocatable :: grad
    ! The co-ordinates of the position being considered.
    integer, intent(in) :: icol, imag
    real, intent(in) :: sig_col, sig_mag
    ! If correlated is zero the colour and magnitude uncertainties are taken
    ! to be uncorrelated.  If its 1, they are correlated as would
    ! be expected in I vs R-I; if its -1 as for V vs V-I (i.e. its whether 
    ! you add or subtract the colour to make the other magnitude).  
    ! If the uncertainties are correlated, make sure sig_col > sig_mag.
    integer, intent(in) :: correlated
    character, intent(out) :: flag

    real, optional, intent(out) :: sum_gauss

    ! We cut off the Gaussian smoothing the 2D distribution at some sigma.
    real, parameter :: sig_cut=3.0
    real, parameter :: pi=3.141592654

    ! Gradient in mag-mag space, and uncertainty in second magnitude.
    real :: mag_grad, mag_sig
    integer :: iicol, iimag, i1, i2
    real :: col, mag
    real :: work
    real, save :: area, pix1, pix2
    real :: d_col, d_mag, gauss, rho, gauss_norm
    logical, save :: first=.true.
    
    if (first) then
      ! Find the area of a pixel.
      pix1=(axdata(size(data,1),1)-axdata(1,1))/real(size(data,1)-1)
      pix2=(axdata(size(data,2),2)-axdata(1,2))/real(size(data,2)-1)
      area=abs(pix1*pix2)
      first=.false.
    end if

    ! Find the nearest positions in X and Y.
    col=axdata(icol,1)
    mag=axdata(imag,2)

    ! How far out do we need to go? 
    iicol=nint(sig_cut*sig_col/abs(pix1))
    iimag=nint(sig_cut*sig_mag/abs(pix2))

    ! For correlated uncertainties calculate mag_sig.
    if (correlated == 0) then
      mag_sig=sig_col
    else
      if (sig_mag > sig_col) then
        mag_sig=sig_mag/2.0
      else
        mag_sig=max(sqrt(sig_col**2.0 - sig_mag**2.0), sig_mag/2.0)
      end if
    end if

    ! Now we normalise the model so it is one within the magnitude
    ! range.


    if (present(sum_gauss)) then
      sum_gauss=0.0
      
      fudge: do i1=icol-iicol, icol+iicol

        ! This sums up the volume under the uncertainty ellipse which
        ! lies on the area of the model.  If it falls outside the colour
        ! range of the image we still add up the volume.  This is because
        ! monte trims in colour to save space.  In magnitude its a
        ! different story.  The model may well be non-zero outside the
        ! magnitude range, so we need to calculate what fraction of this 
        ! error ellipse falls within the monte grid.

        ! This works for all values of i1, even if they lie outside the grid.
        work=axdata(1,1) + (pix1*real(i1-1)) - col

        do i2=max(imag-iimag, 1), min(imag+iimag, size(data,2))
          
          d_mag=(axdata(i2,2)-mag)/sig_mag
          if (correlated == 0) then
            d_col=work/sig_col
          else if (correlated == 1) then
            d_col=(axdata(i2,2)-mag)+work
            d_col=d_col/mag_sig
          else if (correlated == -1) then
            d_col=(axdata(i2,2)-mag)-work
            d_col=d_col/mag_sig
          else 
            print*, 'Unacceptable value of ', correlated
            stop
          end if
          gauss=exp(-0.5*(d_mag**2.0 + d_col**2.0))
          sum_gauss=sum_gauss+(area*gauss)
          
        end do
      end do fudge

      sum_gauss=sum_gauss/(2.0*pi*mag_sig*sig_mag)

    end if


    ! Removed this if.  It used to check you didn't skip off the edge of the
    ! image.  Now the assmumption is that if you do, the values of the image
    ! pixels should be zero.
    !if (icol-iicol>0       .and. imag-iimag>0 .and. &
    !  icol+iicol<=size(data,1) .and. imag+iimag<=size(data,2)) then

      work=0.0

    
      outer: do i1=max(icol-iicol, 1), min(icol+iicol, size(data,1))

        if (allocated(grad)) then
          if (grad(i1) < tiny(grad(i1))) cycle outer
          ! I.e. if there are any non-zero data points at this colour.
          ! Calculate the normalisation.
          if (correlated == 0) then
            rho=(1.0/sig_mag)**2.0 + 1.0/(sig_col*grad(i1))**2.0
          else 
            mag_grad=1.0/(1.0-(1.0/grad(i1)))
            ! Avoid a numberical error if the two uncertainties are equal.
            mag_sig=max(sqrt(sig_col**2.0 - sig_mag**2.0), sig_mag/2.0)
            ! write(22,*) 'Got ', grad(i1), mag_grad, i1, mag_sig
            rho=(1.0/sig_mag)**2.0 + 1.0/(mag_sig*mag_grad)**2.0
            ! write(22,*) 'Which gave ', rho
          end if
          ! rho=sqrt(rho/pi)
          ! Correct the factor two error.
          rho=sqrt(0.5*rho/pi)
          gauss_norm=1.0
        else
          rho=1.0
          gauss_norm=1.0/(2.0*pi*mag_sig*sig_mag)
        end if



        do i2=max(imag-iimag, 1), min(imag+iimag, size(data,2))
          if (data(i1,i2) >= tiny(data(i1,i2))) then
            ! Find the colour difference from this point to the
            ! data point, and put it in terms of sigma.
            if (correlated == 0) then
              d_col=(axdata(i1,1)-col)/sig_col
            else if (correlated == 1) then
              ! But now comes the swindle to account for correlated
              ! uncertainties.  The colour distance we calculate, is 
              ! actually the distance in the other magnitude.
              d_col=(axdata(i2,2)-mag)+(axdata(i1,1)-col)
              ! Thus we have to divide by the uncertainty in the other 
              ! magnitude, reconstructed from the given uncertainties.
              d_col=d_col/mag_sig
            else if (correlated == -1) then
              d_col=(axdata(i2,2)-mag)-(axdata(i1,1)-col)
              d_col=d_col/mag_sig
            else 
              print*, 'Unacceptable value of ', correlated
              stop
            end if
            ! The same for the magnitude.  
            d_mag=(axdata(i2,2)-mag)/sig_mag
            ! gauss=exp(-1.0*(d_mag**2.0 + d_col**2.0))
            ! Correct the factor two error.
            gauss=exp(-0.5*(d_mag**2.0 + d_col**2.0))
            ! Now multiply the Gaussian at this point by the value of the
            ! model, the area of the pixel, and apply the normalisation.
            ! write(22,*) work, rho, area, gauss, data(i1,i2), i1, i2
            work=work+rho*area*gauss_norm*gauss*data(i1,i2)
          end if
        end do
      end do outer

      likelihood=work
      flag = 'O'

    !else

    !  likelihood=0.0
    !  flag = 'A'

    !end if

  end function likelihood


  real function func2d(a_par, dat, colflag, modflag, correlated, &
    distrib, distrib_flag, rnpts)

    use ark_file_io
    use define_star

    real, intent(in), dimension(:) :: a_par
    type(a_star), dimension(:), intent(in) :: dat
    integer, intent(inout) :: colflag
    integer, intent(in) :: modflag
    ! The correlation flag for s/r likelihood.
    integer, intent(in) :: correlated
    real, intent(out), dimension(:) :: distrib
    character, intent(out), dimension(:) :: distrib_flag
    real, optional, intent(out) :: rnpts

    integer :: i
    real, save :: last_age=-1
    real, dimension(:,:), allocatable, save :: data, axdata
    integer, dimension(2), save :: naxis
    real, dimension(:), allocatable, save :: grad, axgrad
    integer, save :: n_grad
    integer :: iflag

    integer :: ipos, icol, imag
    real :: mag, col
    character(len=50) :: ifname

    real :: sum_gauss

    ! The nomrmalisation for the model.
    real :: rho_norm
    integer :: npix
    real :: mag_min, mag_max
    logical, save :: natural

    ! In this example case the function is simply a straight line with 
    ! a width given by the uncertainty in the second co-ordinate.

    if (abs(last_age-a_par(1)) > tiny(last_age)) then
      call iso_image_name(modflag, colflag, a_par(1), ifname)
      call nxtark_in(ifname)
      iflag=inpark(naxis, data, axdata)
      if (iflag < 0) then
        print*, 'Error reading 2D isochrone file ', ifname
        stop
      end if
      ! At the moment colflag could be for the wrong flavour of the filters 
      ! (e.g. Bessell not Jeffries/Naylor).  Lets make sure its right.
      if (get_header_i('COLFLAG', colflag) < 0) then
        print*, 'Cannot find colour number in header.'
        stop
      end if

      last_age=a_par(1)

      ! Now read the gradient file.
      natural=.false.
      ipos=index(ifname, '.fit')
      ifname(ipos:ipos+3)='.grd'
      call nxtark_in(ifname)
      iflag=inpark(n_grad, grad, axgrad)
      if (iflag == -2) then
        print*, 'No gradient file found, using natural normalisation.'
        natural=.true.
      else if (iflag < 0) then
        print*, 'Error reading gradient file file ', ifname
        stop
      end if

    end if

    if (natural) then
      ! Find the range of magnitudes in the data.
      find_range: do i=1, size(dat,1)
        ! Change to reddening-free absolute magnitude.
        mag=dat(i)%col(1)%data
        col=dat(i)%col(2)%data
        call reddening(colflag, -a_par(2), -a_par(3), col, mag)
        if (i ==1) then
          mag_min=mag
          mag_max=mag
        else
          mag_min=min(mag_min, mag)
          mag_max=max(mag_max, mag)
        end if
      end do find_range
      call natural_norm(mag_max, mag_min, data, axdata)
    else
      rho_norm=1.0
      data=data/rho_norm
    end if

    func2d = 0.0
    if (present(rnpts)) rnpts=0.0


    each_point: do i=1, size(dat,1)

      ! Change to reddening-free absolute magnitude.
      mag=dat(i)%col(1)%data
      col=dat(i)%col(2)%data
      call reddening(colflag, -a_par(2), -a_par(3), col, mag)

      ! Find the nearest positions in X and Y.  (S/R likelihood will spot
      ! if we are close to the array edge).
      icol=minval(minloc(abs(col-axdata(1:size(data,1),1))))
      imag=minval(minloc(abs(mag-axdata(1:size(data,2),2))))

      if (present(rnpts)) then
        distrib(i)=likelihood(data, axdata, grad, icol, imag, &
        dat(i)%col(2)%err, dat(i)%col(1)%err, correlated, distrib_flag(i), &
        sum_gauss)
        rnpts=rnpts+sum_gauss
      else
        distrib(i)=likelihood(data, axdata, grad, icol, imag, &
        dat(i)%col(2)%err, dat(i)%col(1)%err, correlated, distrib_flag(i))
      end if

      if (distrib_flag(i) /= 'O') then

        ! Probably outside grid range, so set value to that of the soft
        ! clipping.
        distrib(i) = 0.0

      end if

      func2d=0.0

    end do each_point

    data=data*rho_norm

  end function func2d


  real function isochr(axis_in, val_in, axis_out, age, colflag, modflag, flag)

    ! Given one parameter from an isochrone, this function interpolates 
    ! another parameter from the table.

    ! The axis the input parameter is from, should be 'mag, 'col', or 'mass'.
    character(len=3), intent(in) :: axis_in
    ! The value of the input parameter.
    real, intent(in) :: val_in
    ! The axis the output parameter is from, should be 'mag, 'col', or 'mass'.
    character(len=3), intent(in) :: axis_out
    ! The age, colour and model to be used.
    real :: age
    integer :: colflag, modflag
    ! And a flag.  
    character(len=*), intent(out) :: flag

    real :: val_out, dx
    integer :: i, iostat, j

    logical, save :: first=.true.

    integer :: npt
    real, allocatable, dimension(:) :: newcol, temp

    type chrone 
      ! Now the specification of the iscohrone.  Its age, colours and model.
      real :: age
      integer :: colflag, modflag, npts
      real, dimension(:), allocatable :: col, mag, mass, teff
      character(len=50), dimension(:), allocatable :: &
      colflg, magflg, massflg, teffflg
    end type chrone

    type(chrone), save :: a, b, c

    if (colflag > 20 .and. colflag<=30) then
      isochr=isochr_p(axis_in, val_in, axis_out, age, colflag, modflag, flag)
    else

      if (first) then
        b%age=0.0
        c%age=0.0
        b%colflag=-1
        c%colflag=-1
        first=.false.
      end if

      if (abs(a%age-age)<=2.0*tiny(age) .and. a%colflag==colflag) then
        ! Still using same isochrone.
      else if (abs(b%age-age)<=2.0*tiny(age) .and. b%colflag==colflag) then
        a=b
      else if (abs(c%age-age)<=2.0*tiny(age) .and. c%colflag==colflag) then
        a=c
      else
        a%age=age
        a%colflag=colflag
        a%modflag=modflag
        if (colflag < 10) then
          ! Read the colour-Teff relationship
          open(unit=11, file=&
          trim(data_dir())//'/'//trim(mod_name(modflag))//'_'//&
          trim(mag_name(colflag))//'_'//trim(col_name(colflag))//'.dat', &
          status='old', iostat=iostat)      
          if (iostat /= 0) then
            print*, 'Error opening file ', &
            trim(data_dir())//'/'//trim(mod_name(modflag))//'_'//&
            trim(mag_name(colflag))//'_'//trim(col_name(colflag))//'.dat'
            stop
          end if
          read(11,'(//)')
          npt=0
          count_teff: do
            read(11,*,iostat=iostat)
            if (iostat < 0) exit count_teff
            npt=npt+1
          end do count_teff
          rewind(11)
          allocate(newcol(npt), temp(npt))
          read(11,'(//)')
          do i=1, npt
            read(11,*) newcol(i), temp(i)
          end do
          temp=log10(temp)
          close(11)
        end if
 
        ! Calculate the isochrone
        !print*, 'Calling iso_calc.'
        call iso_calc(a%age, a%colflag, a%modflag, &
        npt, newcol, temp, a%npts, a%col, a%colflg, a%mag, a%magflg, &
        a%mass, a%massflg, a%teff, a%teffflg)

        !do i=1, a%npts
        !  print*, a%col(i), a%colflg(i), a%mag(i), a%magflg(i)
        !end do

        if (b%colflag == -1) then
          b=a
        else if (c%colflag == -1) then
          c=a
        end if

      end if

      if (axis_in=='mag' .and. axis_out=='col') then
        val_out=quadintflg(a%mag, a%magflg, a%col, a%colflg, val_in, flag)
        if (trim(flag) == '1') flag='outside magnitude range'
      else if (axis_in=='col' .and. axis_out == 'mag') then
        val_out=quadintflg(a%col, a%colflg, a%mag, a%magflg, val_in, flag)
        if (trim(flag) == '1') flag='outside colour range'
      else if (axis_in=='mas' .and. axis_out == 'mag') then
        val_out=quadintflg(a%mass, a%massflg, a%mag, a%magflg, val_in, flag)
        if (trim(flag) == '1') flag='outside mass range'
      else if (axis_in=='mas' .and. axis_out == 'col') then
        val_out=quadintflg(a%mass, a%massflg, a%col, a%colflg, val_in, flag)
        if (trim(flag) == '1') flag='outside mass range'
      else if (axis_in=='mas' .and. axis_out == 'tef') then
        val_out=quadintflg(a%mass, a%massflg, a%teff, a%teffflg, val_in, flag)
        if (trim(flag) == '1') flag='outside mass range'
      else if (axis_in=='mag' .and. axis_out == 'mas') then
        val_out=quadintflg(a%mag, a%magflg, a%mass, a%massflg, val_in, flag)
        if (trim(flag) == '1') flag='outside magnitude range'
      else
        print*, 'Illegal combination of axis_in ', axis_in
        print*, 'and axis_out ', axis_out, ' in function isochr.'
        stop
      end if

      isochr=val_out

    end if

  end function isochr

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

      use define_star

      ! Performs a 2d grid search.

      ! The model parameters.
      real, dimension(:) :: a_par
      ! The data points.
      type(a_star), dimension(:), intent(in) :: dat
      ! The model and colours being used.
      integer, intent(inout) :: colflag
      integer, intent(in) :: modflag
      ! The correlation flag for s/r likelihood.
      integer, intent(in) :: correlated
      ! The numbers of the parameters to be searched.
      integer, dimension(2), intent(in) :: var_par
      ! The ranges of parameters to be searched.
      real, dimension(2), intent(in) :: start, end
      real, intent(in) :: soft_clip
      ! The output grid.
      real, dimension(:,:), intent(out) :: data, axdata
      ! And an optional fraction of the stars on the grid.
      real, dimension(:,:), intent(out), optional :: rnpts

      real, intent(out) :: best_value
      logical :: first_value
      real, intent(out), dimension(:) :: best_par
      real, intent(out), dimension(:) :: best_distrib

      ! Locals.
      character, allocatable, dimension(:) :: distrib_flag
      real, allocatable, dimension(:) :: distrib
      integer :: i1, i2, i3
      real, dimension(size(data,1), size(data,2), size(dat,1)) :: distrib_full
      real :: dummy, work, low_tau
      integer, dimension(2) :: loc

      allocate(distrib_flag(size(dat,1)), distrib(size(dat,1)))

      ! Set up the values of best_par to be those of a_par.  We'll
      ! the correct the ones we fit.
      best_par=a_par

      first_value=.true.
      do i1=1, size(data,1)
        if (size(data,1) == 1) then
          a_par(var_par(1))=start(1)
        else
          a_par(var_par(1))=start(1)+(end(1)-start(1))*real(i1-1) &
          /real(size(data,1)-1)
        end if
        axdata(i1,1)=a_par(var_par(1))
        do i2=1, size(data,2)
          if (size(data,2) == 1) then 
            a_par(var_par(2))=start(2)
          else
            a_par(var_par(2))=start(2)+(end(2)-start(2))*real(i2-1) &
            /real(size(data,2)-1)
          end if
          axdata(i2,2)=a_par(var_par(2))
          if (present(rnpts)) then
            dummy=func2d(a_par, dat, colflag, modflag, correlated, &
                 distrib_full(i1,i2,:), distrib_flag, rnpts(i1, i2))
          else 
            dummy=func2d(a_par, dat, colflag, modflag, correlated, &
                 distrib_full(i1,i2,:), distrib_flag)
          end if
        end do
      end do

      ! Find the highest probability in the whole fit.
      low_tau=-2.0*log(maxval(distrib_full))
      print*, 'Lowest tau-squared for any point is ', low_tau
      if (soft_clip > 0.0) then
        work=exp(-0.5*(soft_clip+low_tau))
      else
        work=0.0
      end if
      ! Soft clip the probabilities.
      distrib_full=distrib_full+work

      data=0.0
      high_tau=-2.0*log(2.0*tiny(1.0))*size(dat,1)
      do i1=1, size(data,1)
        do i2=1, size(data,2)
          ! Construct the tau**2 for this grid point by adding
          ! the tau^2 of each point.
          each_point: do i3=1, size(dat,1)
            ! Protect against infinities.
            if (distrib_full(i1,i2,i3) > 2.0*tiny(1.0)) then
              data(i1,i2)=data(i1,i2)-2.0*log(distrib_full(i1,i2,i3))
            else
              ! At least one data point is returning a probablilty so low
              ! we can't cope, so set the entire grid point to a very
              ! high tau^2.
              data(i1,i2)=high_tau
              exit each_point
            end if
          end do each_point
        end do
      end do

      ! Find the best fit.
      loc=minloc(data)
      best_distrib=-2.0*log(distrib_full(loc(1),loc(2),:))
      best_par(var_par(1))=axdata(loc(1),1)
      best_par(var_par(2))=axdata(loc(2),2)
      best_value=sum(best_distrib)

      deallocate(distrib_flag, distrib)

    end subroutine grid2d

    subroutine tau_sort(tau_data)

      real, intent(inout), dimension(:)  :: tau_data

      integer :: k, l, m
      real :: safe

      outer: do k=2, size(tau_data)
        safe=tau_data(k)
        do l=1,k-1
          if (safe < tau_data(l)) then
            do m=k, l+1, -1
              tau_data(m)=tau_data(m-1)
            end do
            tau_data(l)=safe
            cycle outer
          endif
        end do
      end do outer

    end subroutine tau_sort

    real function isochr_p(axis_in, val_in, axis_out, age, colflag, &
         modflag, flag)

      implicit none

      ! The axis the input parameter is from, should be 'mag, 'col', or 'mass'.
      character(len=3), intent(in) :: axis_in
      ! The value of the input parameter.
      real, intent(in) :: val_in
      ! The axis the output parameter is from, should be 'mag, 'col', or 'mass'.
      character(len=3), intent(in) :: axis_out
      ! Now the specification of the iscohrone.  Its age, colours and model.
      real, intent(in) :: age
      integer, intent(in) :: colflag, modflag
      ! And a flag.
      character(len=*), intent(out) :: flag

      real :: dmass, dage, val_out
      integer :: iage, jpt
      logical, save :: first=.true.

      type(isochron), dimension(3), save :: isoc, last_but_one_isoc, safe
      real, dimension(3) :: work_out
      real, save :: last_age
      integer, save :: last_colflag, last_but_one_colflag
      real, save :: min_mass, max_mass

      if (.not. first) then
        if (abs(age-last_age) > 2.0*tiny(age)) then
          first=.true.
          last_colflag=colflag
          last_but_one_colflag=0
        else
          if (colflag /= last_colflag) then
            if (colflag == last_but_one_colflag) then
              ! This is the last-but-one colour used, so swap isoc and
              ! but-one isoc.
              safe=isoc
              isoc=last_but_one_isoc
              last_but_one_isoc=safe
              last_but_one_colflag=last_colflag
              last_colflag=colflag
            else
              ! Never used before.  So put the current isoc into the
              ! last-but-one, and read a new one.
              first=.true.
              last_but_one_colflag=last_colflag
              last_but_one_isoc=isoc
              last_colflag=colflag
            end if
          else
            ! Using the current isoc.
          end if          
        end if
      end if
      last_age=age
      last_colflag=colflag

      if (first) then
        if (modflag == 8) then
          call read_padova_col(age, colflag, isoc)
        else if (modflag == 9) then
          call read_geneva_col(age, colflag, isoc)
        else
          print*, 'Only the Padova and Geneva isochrones are available'
          print*, 'with their own supplied colours.'
          stop
        end if
        first=.false.
      end if

      min_mass=maxval(isoc%mass(1))
      max_mass=min(isoc(1)%mass(isoc(1)%npts), &
      min(isoc(2)%mass(isoc(2)%npts), &
      isoc(3)%mass(isoc(3)%npts)))

      isochr_p=0.0
      if (val_in < min_mass .and. axis_in=='mas') then
        flag='Below minimum mass'
      else if (val_in > max_mass .and. axis_in=='mas') then
        flag='Above maximum mass'
      else
        ! Interpolate the three isochrones to the appropriate mass.
        do iage=1, 3
          ! Find the nearest point to the mass required.
          if (axis_in == 'mas') then
            call locate(isoc(iage)%mass, isoc(iage)%npts, val_in, jpt)
          else if (axis_in == 'mag') then
            call locate(isoc(iage)%magnitude, isoc(iage)%npts, val_in, jpt)
          end if
          ! Avoid running off the end.
          jpt = min(max(jpt-1,1), isoc(iage)%npts-2)
          ! Now interpolate (parabolic so n=3), returns colour_out and 
          ! dmass (error)   
          if (axis_in == 'mas') then
            if (axis_out == 'col') then
              call polint(isoc(iage)%mass(jpt:jpt+2), &
              isoc(iage)%colour(jpt:jpt+2), 3, val_in, work_out(iage), dmass)
            else if (axis_out == 'mag') then
              call polint(isoc(iage)%mass(jpt:jpt+2), &
              isoc(iage)%magnitude(jpt:jpt+2), 3, val_in, work_out(iage), dmass)
            else if (axis_out == 'tef') then
              work_out(iage)=0.0
            else
              print*, 'Error.  Requested output axis is ', axis_out
              print*, 'For input axis ', axis_in
              stop
            end if
          else if (axis_in == 'mag') then
            if (axis_out == 'mas') then
              call polint(isoc(iage)%magnitude(jpt:jpt+2), &
              isoc(iage)%mass(jpt:jpt+2), 3, val_in, work_out(iage), dmass) 
            else
              print*, 'Error.  Requested output axis is ', axis_out
              print*, 'For input axis ', axis_in
              stop
            end if
          else
            print*, 'axis_in must be mag or mas'
            stop
          end if
        end do

        ! Interpolate to the appropriate age.
        call polint(isoc%age, work_out, 3, age, val_out, dage)

        isochr_p=val_out
        flag='OK'
      end if

    end function isochr_p

    end module cmdfit_subs
