module phoenix_bessell

  use quad

  implicit none

  integer, parameter, private :: m_grav=13, n_teff=66
  type a_teff
    real, dimension(m_grav) :: logg
  end type a_teff
  type a_corr
    type(a_teff), dimension(n_teff) :: teff
  end type a_corr

  ! Things read in from the data file.
  type(a_corr), save, private :: bc_u, bc_b, bc_v, bc_r, bc_i
  real, save, dimension(n_teff), private :: temp
  real, save, dimension(m_grav), private :: grav
  integer, save, dimension(n_teff), private :: n_grav
  logical, save, private :: first=.true.

contains

  subroutine readfiles()

    ! Allows table of bolometric corrections to be read in.
    integer, parameter :: m_grav=13
    integer, parameter :: n_teff=66
    type a_teff
      real, dimension(m_grav) :: logg
    end type a_teff
    type a_corr
      type(a_teff), dimension(n_teff) :: teff
    end type a_corr
    integer :: i, j, int_grav, int_temp, iostat, j_grav
    real :: u_in, b_in, v_in, r_in, i_in
    real :: grav_in, temp_in

    real, dimension(m_grav) :: col_grav, work_grav
    integer :: k_grav
    character(len=20) :: flag, iflg
    character(len=30) :: data_dir

    ! Define the gravity and temperature arrays.
    grav(1)=5.5
    do i=2, m_grav
      grav(i)=grav(i-1)-0.5
    end do

    temp(1)=2000.
    do i=2, n_teff
      temp(i)=temp(i-1)+100.
      if (i > 51.) then
        temp(i)=temp(i-1)+200.
      end if
    end do

    temp=log10(temp)

    call get_environment_variable('CMDDATA', data_dir)
    open(10, file=trim(data_dir)//'/phoenix_ubv_new1.dat', action='read')
    read(10,*)
    read(10,*)

    n_grav=0
    do 
      read(10,*,iostat=iostat)  temp_in, grav_in, u_in, b_in, v_in, r_in, i_in
      if (iostat < 0) exit
      temp_in=log10(temp_in)
      int_grav=minloc(abs(grav-grav_in),1)
      int_temp=minloc(abs(temp-temp_in),1)
      bc_u%teff(int_temp)%logg(int_grav)=u_in
      bc_b%teff(int_temp)%logg(int_grav)=b_in
      bc_v%teff(int_temp)%logg(int_grav)=v_in
      bc_r%teff(int_temp)%logg(int_grav)=r_in
      bc_i%teff(int_temp)%logg(int_grav)=i_in
      n_grav(int_temp)=n_grav(int_temp)+1
    end do

    close(10)

  end subroutine readfiles


  type(a_corr) function subtract_corr(bc1, bc2)

    type(a_corr), intent(in) :: bc1, bc2

    integer :: i_teff, i_grav

    subtract_corr = bc1
    do i_teff=1, n_teff
      do i_grav=1, m_grav
        subtract_corr%teff(i_teff)%logg(i_grav) = &
        bc1%teff(i_teff)%logg(i_grav) - bc2%teff(i_teff)%logg(i_grav)
      end do
    end do

  end function subtract_corr




  ! Functions to create colours from effective temperatures from 
  ! tabulated ones bolometric corections.
  real function phoenix_col_bv(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_col_bv=interp(grav_in, temp_in, subtract_corr(bc_v,bc_b), flg)

  end function phoenix_col_bv

  real function phoenix_col_vi(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_col_vi=interp(grav_in, temp_in, subtract_corr(bc_i,bc_v), flg)

  end function phoenix_col_vi
  


  ! Functions to create bolometric corrections from tabulated ones.
  real function phoenix_blc_u(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_blc_u=interp(grav_in, temp_in, bc_u, flg)

  end function phoenix_blc_u

  real function phoenix_blc_b(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_blc_b=interp(grav_in, temp_in, bc_b, flg)

  end function phoenix_blc_b

  real function phoenix_blc_v(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_blc_v=interp(grav_in, temp_in, bc_v, flg)

  end function phoenix_blc_v

  real function phoenix_blc_r(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_blc_r=interp(grav_in, temp_in, bc_r, flg)

  end function phoenix_blc_r

  real function phoenix_blc_i(grav_in, temp_in, flg)

    real, intent(in) :: grav_in, temp_in
    character(len=*), intent(out) :: flg

    if (first) call readfiles()
    phoenix_blc_i=interp(grav_in, temp_in, bc_i, flg)

  end function phoenix_blc_i

  real function interp(grav_in, temp_in, col, flg)

    real, intent(in) :: grav_in, temp_in
    type(a_corr), intent(in) :: col
    character(len=*), intent(out) :: flg

    real, dimension(m_grav) :: col_grav, work_grav
    real, dimension(n_teff) :: col_temp, work_temp
    character(len=20) :: flag

    integer :: j_grav, k_temp, j_temp

    flg='OK'

    ! At each gravity, interpolate the bolometric correction at the
    ! appropriate temperature, to make an array at a fixed temperature
    ! but a range of gravities.
    !do j_grav=1, m_grav
    !  col_grav(k_grav+1)=quadint(temp(:), col%teff(:)%logg(j_grav), &
    !       temp_in, flag)
    !  work_grav(k_grav+1)=grav(j_grav)
    !print*, grav(j_grav)
    !  k_grav=k_grav+1
    !end do
    !if (k_grav > 1) then
    !  interp=quadint(work_grav(1:k_grav), col_grav(1:k_grav), grav_in, flag)
    !end if

    ! At each temperature, interpolate the bolometric correction at
    ! the appropriate gravity, to make an array at a fixed gravity
    ! with a range of temperatures.
    k_temp=0
    do j_temp=1, n_teff
      col_temp(k_temp+1)=quadint(grav(1:n_grav(j_temp)), & 
           col%teff(j_temp)%logg(1:n_grav(j_temp)), grav_in, flag)
      !work_temp(j_temp)=grav(j_temp)
      if (flag == 'OK') then
        k_temp=k_temp+1
        work_temp(k_temp)=temp(j_temp)
        !print*, work_temp(k_temp), col_temp(k_temp)
      end if
    end do
    ! print*, grav_in, temp_in
    if (k_temp > 1) then
      interp=quadint(work_temp(1:k_temp), col_temp(1:k_temp), temp_in, flag)
      if (flag == '1') flg='1'
    else
      flg='1'
      interp=0.0
    end if

  end function interp

end module phoenix_bessell
