module phoenix_readfiles

  use quad

  implicit none

  ! The maximum number of gravities and temperatures.
  integer, parameter, public :: m_grav=13, m_teff=100
  type a_teff
    real, dimension(m_grav) :: logg
  end type a_teff
  type a_corr
    type(a_teff), dimension(m_teff) :: teff
  end type a_corr

contains

  subroutine readfiles(file, bc1, bc2, bc3, bc4, bc5, temp, grav, &
  n_teff, n_grav)

    ! Allows table of bolometric corrections to be read in.

    ! We make some assumptions about the file.
    ! 1) The file has records which have, in order, temperature, gravity and 
    !    bolometric corrections.
    ! 2) The gravities for all temperatures run from 5.5 downwards in steps
    !    of 0.5 (don't low they go). 

    character(len=*), intent(in) :: file
    type(a_corr), intent(inout) :: bc1, bc2, bc3, bc4, bc5
    real, dimension(m_teff), intent(inout) :: temp
    real, dimension(m_grav), intent(inout) :: grav
    integer, intent(inout) :: n_teff
    integer, dimension(m_teff), intent(inout) :: n_grav

    character(len=80), save :: lastfile=' '
    
    integer :: i, j, int_grav, int_temp, iostat, j_grav
    real :: col1_in, col2_in, col3_in, col4_in, col5_in
    real :: grav_in, temp_in, temp_safe

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

    if (trim(lastfile) == trim(file)) return

    lastfile=file

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

    ! Find the number of different temperatures.
    call get_environment_variable('CMDDATA', data_dir)
    open(10, file=trim(data_dir)//file, action='read') 
    read(10,*)
    read(10,*)
    n_teff=1
    read(10,*) temp_safe
    do
      read(10,*,iostat=iostat) temp_in
      if (iostat < 0) exit
      if (nint(temp_in) /= nint(temp_safe) ) then
        n_teff=n_teff+1
        temp_safe=temp_in
      end if
    end do
    if (n_teff > m_teff) then
      print*, 'In s/r phoenix_bessell m_teff needs enlarging to ', n_teff
      stop
    end if
    rewind(10)
    read(10,*)
    read(10,*)
    read(10,*) temp(1)
    i=1
    do
      read(10,*) temp_in
      if (nint(temp_in) /= nint(temp(i))) then
        i=i+1
        temp(i)=temp_in
      end if
      if (i == n_teff) exit
    end do

    temp=log10(temp)

    rewind(10)
    read(10,*)
    read(10,*)
    n_grav=0
    do 
      read(10,*,iostat=iostat)  temp_in, grav_in, col1_in, col2_in, &
      col3_in, col4_in, col5_in
      if (iostat < 0) exit
      ! We work in log temperature.
      temp_in=log10(temp_in)
      ! Find the array indices most appropriate for the read temperature 
      ! and gravity.
      int_grav=minloc(abs(grav-grav_in),1)
      int_temp=minloc(abs(temp-temp_in),1)
      ! And fill those elements.
      bc1%teff(int_temp)%logg(int_grav)=col1_in
      bc2%teff(int_temp)%logg(int_grav)=col2_in
      bc3%teff(int_temp)%logg(int_grav)=col3_in
      bc4%teff(int_temp)%logg(int_grav)=col4_in
      bc5%teff(int_temp)%logg(int_grav)=col5_in
      ! But what about the unfilled elements?  Well, all temperatures have
      ! gravity 5.5, but not all go all the way to -0.5.  Fortunately the
      ! file puts them in the order 5.5, 5.0, 4.5 etc.  So if we just count
      ! the number of gravities at a given temperature, the total must be
      ! the number of elements filled.
      n_grav(int_temp)=n_grav(int_temp)+1
    end do

    close(10)

    !write(33,*) n_grav
    !do int_temp=1, n_teff
    !  do int_grav=1, m_grav
    !    write(33,*) temp(int_temp), grav(int_grav), &
    !         bc3%teff(int_temp)%logg(int_grav)
    !  end do
    !end do

  end subroutine readfiles

  type(a_corr) function subtract_corr(bc1, bc2, n_teff)

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

    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

  real function interp(grav_in, temp_in, col, grav, temp, n_teff, n_grav, flg)

    real, intent(in) :: grav_in, temp_in
    type(a_corr), intent(in) :: col
    real, dimension(:), intent(in) :: grav, temp
    integer, intent(in) :: n_teff
    integer, intent(in), dimension(:) :: n_grav
    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_grtav(k_grav+1)=grav(j_grav)
    !print*, gtrav(j_grav)
    !  k_grav=tk_grav+1
    !end do   t
    !if (k_gratv > 1) then
    !  interp=tquadint(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
      if (n_grav(j_temp) > 0) then
        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 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_readfiles


module phoenix

  use phoenix_readfiles

  implicit none

  ! Things read in from the data file.
  character(len=80), private, parameter :: bessell='/phoenix_ubv_Z_0.0.dat' 
  character(len=80), private, parameter :: sloan='/phoenix_sloan_bessell.dat'
  type(a_corr), save, private :: bc_u, bc_b, bc_g, bc_v, bc_r, bc_i, bc_z
  real, save, dimension(m_teff), private :: temp
  real, save, dimension(m_grav), private :: grav
  integer, save :: n_teff
  integer, save, dimension(m_teff), private :: n_grav

contains



  ! Functions to create colours from effective temperatures from 
  ! tabulated ones bolometric corections.

  real function phoenix_bess_bv(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_bv=interp(grav_in, temp_in, &
    subtract_corr(bc_v, bc_b, n_teff), grav, temp, n_teff, n_grav, flg)

  end function phoenix_bess_bv

  real function phoenix_bess_vi(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_vi=interp(grav_in, temp_in, &
    subtract_corr(bc_i, bc_v, n_teff), grav, temp, n_teff, n_grav, flg)

  end function phoenix_bess_vi
  
  real function phoenix_sloan_gi(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_g, bc_r, bc_i, bc_z, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_gi=interp(grav_in, temp_in, &
    subtract_corr(bc_i, bc_g, n_teff), grav, temp, n_teff, n_grav, flg)

  end function phoenix_sloan_gi

  real function phoenix_sloan_iz(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_g, bc_r, bc_i, bc_z, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_iz=interp(grav_in, temp_in, &
    subtract_corr(bc_z, bc_i, n_teff), grav, temp, n_teff, n_grav, flg)

  end function phoenix_sloan_iz



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

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_u=interp(grav_in, temp_in, bc_u, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_bess_u

  real function phoenix_bess_b(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_b=interp(grav_in, temp_in, bc_b, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_bess_b

  real function phoenix_bess_v(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_v=interp(grav_in, temp_in, bc_v, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_bess_v

  real function phoenix_bess_r(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_r=interp(grav_in, temp_in, bc_r, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_bess_r

  real function phoenix_bess_i(grav_in, temp_in, flg)

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

    call readfiles(bessell, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_bess_i=interp(grav_in, temp_in, bc_i, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_bess_i

  real function phoenix_sloan_u(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_u=interp(grav_in, temp_in, bc_u, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_sloan_u

  real function phoenix_sloan_g(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_g=interp(grav_in, temp_in, bc_g, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_sloan_g

  real function phoenix_sloan_r(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_r=interp(grav_in, temp_in, bc_r, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_sloan_r

  real function phoenix_sloan_i(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_i=interp(grav_in, temp_in, bc_i, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_sloan_i

  real function phoenix_sloan_z(grav_in, temp_in, flg)

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

    call readfiles(sloan, bc_u, bc_b, bc_v, bc_r, bc_i, temp, grav, &
    n_teff, n_grav)
    phoenix_sloan_z=interp(grav_in, temp_in, bc_z, grav, temp, n_teff, &
    n_grav, flg)

  end function phoenix_sloan_z

end module phoenix
