module bessell

    ! Given the gravity and colour, this interpolates
    ! the Bessell models to get the bolometric correction.

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

  use quad

  implicit none

  integer, parameter :: mgrav=11
  real, private, save, dimension(mgrav) :: gravity
  real, private, save, dimension(100,mgrav) :: teff=0.0
  integer, private, save, dimension(mgrav) :: nteff
  real, private, save, dimension(100,mgrav) :: bc_v=0.0, u_b=0.0, b_v=0.0, &
  v_i=0.0
  logical, private, save :: first=.true.
  integer, private, save :: ntycho
  real, private, save, dimension(:), allocatable :: delta_bc, delta_col, &
  col_bess, col_tycho
  logical, private, save :: got_tycho=.false.

  ! Do you want the colours of Vega to be zero (the V magnitude is always
  ! 0.03), or do you want it to have the colours implied from the Kurucz
  ! model atmospheres?  See the discussion in Bessell, Castelli & Plez (1998)
  ! vs that in Giradi et al (2002).
  logical :: vega_col_zero=.true.

  contains

  subroutine readfiles()

    integer :: jgrav
    logical :: table2
    integer :: iostat
    real :: teff_in, logg, bc_k, v_r, v_k, j_h, j_k, k_l

    character(len=30) :: data_dir

    first=.false.
    each_grav: do jgrav=1, mgrav
      nteff(jgrav)=0
      gravity(jgrav) = 0.5*real(jgrav) - 0.5
      ! Open up table 2.
      call get_environment_variable('CMDDATA', data_dir)
      open(unit=1, file=trim(data_dir)//'/Bessell/table2.dat', &
      action='read')
      ! And make it clear we have not yet got to table3.
      table2=.true.
      read: do 
        read(1,*, iostat=iostat) teff_in, logg, &
        bc_k, bc_v(nteff(jgrav)+1,jgrav), u_b(nteff(jgrav)+1,jgrav), &
        b_v(nteff(jgrav)+1,jgrav), v_r, v_i(nteff(jgrav)+1,jgrav)
        if (iostat /= 0) then
          if (table2) then
            close(1)
            open(unit=1, file=trim(data_dir)//'/Bessell/table3.dat', &
            action='read')
            read(1,*)
            read(1,*)
            table2=.false.
            cycle read
          else
            close (1)
            exit read
          end if
        end if
        if (abs(logg - (0.5*real(jgrav) - 0.5)) < tiny(gravity(jgrav))) then
          nteff(jgrav)=nteff(jgrav)+1
          teff(nteff(jgrav),jgrav)=log10(teff_in)
        end if
      end do read
      !print*, gravity(jgrav)
      !print*, teff(1:nteff(jgrav),jgrav)
    end do each_grav
    close(1)

  end subroutine readfiles


  real function bessell_col_ub(grav_in, teff_in, iflag)

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

    work=0.0
    if (vega_col_zero) work=0.004

    if (first) call readfiles()
    bessell_col_ub=bessell_interp(grav_in, teff_in, u_b, iflag)+work

  end function bessell_col_ub


  real function bessell_col_bv(grav_in, teff_in, iflag)

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

    work=0.0
    if (vega_col_zero) work=0.002

    if (first) call readfiles()
    bessell_col_bv=bessell_interp(grav_in, teff_in, b_v, iflag)+work

  end function bessell_col_bv


  real function bessell_col_vi(grav_in, teff_in, iflag)

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

    work=0.0
    if (vega_col_zero) work=0.003

    if (first) call readfiles()
    bessell_col_vi=bessell_interp(grav_in, teff_in, v_i, iflag)+work

  end function bessell_col_vi


  real function bessell_bc_v(grav_in, teff_in, iflag)

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

    if (first) call readfiles()
    bessell_bc_v=bessell_interp(grav_in, teff_in, bc_v, iflag)

  end function bessell_bc_v


  real function bessell_interp(grav_in, teff_in, col, iflag)

    real, intent(in) :: grav_in, teff_in
    real, dimension(:,:), intent(in) :: col
    integer, intent(out) :: iflag

    real, dimension(mgrav) :: col_grav, work_grav
    character(len=50) :: flag

    integer :: jgrav, kgrav

    iflag=0

    kgrav=0
    do jgrav=1, mgrav
      ! At each gravity interpolate at the requested temperature.
      col_grav(kgrav+1)=quadint(teff(1:nteff(jgrav),jgrav), &
      col(1:nteff(jgrav),jgrav), teff_in, flag)
      work_grav(kgrav+1)=gravity(jgrav)
      if (flag == 'OK') kgrav=kgrav+1
    end do
    if (kgrav > 1) then
      ! Now interpolate to the appropriate gravity.
      bessell_interp=quadint(work_grav(1:kgrav), col_grav(1:kgrav), &
      grav_in, flag)
      if (trim(flag) == 'OK') then
        !write(57,*) quadint(work_grav(1:kgrav), col_grav(1:kgrav), &
        !grav_in, flag), & 
        !work_grav(1:kgrav), grav_in, col_grav(1:kgrav)
      end if
      if (flag == '1') then
        iflag=11
        !print*, 'Range of gravities ', minval(work_grav(1:kgrav)), &
        !maxval(work_grav(1:kgrav)), ' at Teff ', 10.0**teff_in
        !print*, 'Required gravity ', grav_in
        !stop
      end if
    else
      iflag=12
      !print*, 'Only ', kgrav, ' gravities available at teff ', teff_in
      !stop
    end if

  end function bessell_interp

  subroutine read_tycho()

    character(len=50) :: data_dir
    integer :: i, iostat
    real :: dummy

    if (.not. got_tycho) then

      call get_environment_variable('CMDDATA', data_dir)
      open(unit=1, file=trim(data_dir)//'/Bessell/btvt.dat', &
      action='read')
      read(1,'(/,/)')

      ! Count the number of data points.
      ntycho=0
      count: do
        read(1,*,iostat=iostat) dummy
        if (iostat < 0) exit
        ntycho=ntycho+1
      end do count
      allocate(delta_bc(ntycho), delta_col(ntycho), col_bess(ntycho), &
      col_tycho(ntycho))
      rewind(1)
      read(1,'(/,/)')
      do i=1, ntycho
        read(1,*) col_tycho(i), delta_bc(i), delta_col(i)
      end do
      close(1)
      col_bess=col_tycho+delta_col
      ! print*, 'Read ', ntycho, ' tycho datapoints.'
      got_tycho=.true.

    end if

  end subroutine read_tycho

  real function tycho_bc(colour, bcflag)

    ! Given a Tycho B-V, this gives the change required to the V
    ! band bolometric correction.

    real, intent(in) :: colour
    integer, intent(inout) :: bcflag

    character(len=20) :: flag

    call read_tycho()

    ! Since the bc is subtracted from the data, sign is right.
    tycho_bc=quadint(col_tycho, delta_bc, colour, flag)

    bcflag=0
    if (flag == '1') bcflag=1

  end function tycho_bc

  subroutine tycho_col(colour, flag)

    ! Given a Bessell B-V, this gives the change required to the V
    ! band bolometric correction

    real, intent(inout) :: colour
    character(len=*), intent(inout) :: flag
    real :: work

    call read_tycho()

    ! Since the colour difference is subtracted from the bessell colour.
    colour=colour-quadint(col_bess, delta_col, colour, flag)

  end subroutine tycho_col
  

end module bessell
