module ken_hart

    ! Given the colour, this interpolates
    ! the Kenyon and Hartmann relationship to get the bolometric correction
    ! and colour.

    ! Error flag.
    ! 00 = success.
    ! 12 = out of range for colour.

  use quad

  implicit none

  integer, parameter :: nteff=55
  real, private, save, dimension(nteff) :: teff, bc_v, u_b, b_v, &
      v_rc, v_rj, v_i, v_ij, v_j, v_h, v_k, v_l 
  logical, private, save :: first=.true.

  contains

  subroutine readfiles_kh()

    integer :: iteff
    character(len=5) :: sp 

    character(len=30) :: data_dir

    first=.false.
    ! Open up table 5.
    call get_environment_variable('CMDDATA', data_dir)
    open(unit=1, file=trim(data_dir)//'/KenyonHartmann/table5.dat', &
    action='read')
    read(1,'(/,/)')
    each_teff: do iteff=1, nteff
      read(1,*) sp, teff(iteff), bc_v(iteff), u_b(iteff), b_v(iteff), &
      v_rc(iteff), v_rj(iteff), v_i(iteff), v_ij(iteff), v_j(iteff), &
      v_h(iteff), v_k(iteff), v_l(iteff) 
    end do each_teff
    teff=log10(teff)
    close(1)

  end subroutine readfiles_kh

  real function ken_hart_col_vi(teff_in, cflag)

    real, intent(in) :: teff_in
    character(len=*), intent(out) :: cflag

    if (first) call readfiles_kh()
    ken_hart_col_vi=quadint(teff, v_i, teff_in, cflag)

  end function ken_hart_col_vi


  real function ken_hart_bc_v(teff_in, iflag)

    real, intent(in) :: teff_in
    integer, intent(out) :: iflag

    character(len=2) :: cflag

    if (first) call readfiles_kh()
    ken_hart_bc_v=quadint(teff, bc_v, teff_in, cflag)
    iflag=0
    if (cflag /= 'OK') iflag=1

  end function ken_hart_bc_v

  real function ken_hart_col_jh(teff_in, cflag)

    real, intent(in) :: teff_in
    character(len=*), intent(out) :: cflag

    if (first) call readfiles_kh()
    ken_hart_col_jh=quadint(teff, v_h-v_j, teff_in, cflag)

  end function ken_hart_col_jh


  real function ken_hart_bc_j(teff_in, iflag)

    real, intent(in) :: teff_in
    integer, intent(out) :: iflag

    character(len=2) :: cflag

    if (first) call readfiles_kh()
    ! Think about the sign here.  As the star becomes redder, it becomes
    ! brighter in Mj, since the BC appears as negative in the expression
    ! for Mj.
    ken_hart_bc_j=quadint(teff, bc_v+v_j, teff_in, cflag)
    iflag=0
    if (cflag /= 'OK') iflag=1

  end function ken_hart_bc_j

end module ken_hart
