program read

  integer :: i
  double precision :: a0, a1, a2, I00L
  double precision :: ct, ii
  integer :: lambda
  integer :: it
  double precision :: s
  double precision :: c, h, k, L, T
  character (len=8) :: name
  double precision :: tmin, tmax, tau, tau_max, tau_min, target, tau6500, tau600

  open(unit=10, file="table")
  do i = 1, 8
!
!     -----loop over the 8 coefficients in the table
!
    read(10, *) lambda, a0, a1, a2, I00L
!
!     -----make a file name with the correct lambda value
!
    write(name,'(i5)')  lambda
    open(unit=11,file=name)
!
!     -----loop over possible tau values from 0 to 3.0
!
    do it = 0, 30
      tau = dble(it) / 10.0d0

!
!     -----find the sources function
!
      S = i00L * (a0 + a1 * tau + 0.5d0 * a2 * tau**2 )

      c = 3d10;      
      h = 6.57d-27;  
      k = 1.38d-16;  
      L = DBLE(lambda) / 1d8

!
!     ----- find temp by solving for T using SLO and Plancks law
!
      T = h * c / (k * L * log(1.0d0  + 2.0d0 *h*c**2/(L**5 * S)));
!
!     -----write the data into the file
!
      write(11,'(2f20.8, e20.8)') tau, T, s
    enddo
    close(11)

!
!     -----write data to a gnuplot file
!
    if (i == 1) write(27,'(a)', advance='no') "plot "
    write(27,'(5a)',advance='no') '"',trim(name),'" using 1:b ti "',trim(name),'" w li'
    if (i<8)  then
      write(27,'(a)', advance='no') ','
    else
      write(27,*)''
    endif

!
!     -----find the temp range between tau= 1 and tau = 0
!          for the current wavelength
    tau_min = 0.0
    tau_max = 3.0
    tau = 0.5d0 * (tau_min + tau_max)
    S = i00L * (a0 + a1 * tau_min + 0.5d0 * a2 * tau_min**2 )
    Tmin =  h * c / (k * L * log(1.0d0  + 2.0d0 *h*c**2/(L**5 * S)));
    S = i00L * (a0 + a1 * tau_max + 0.5d0 * a2 * tau_max**2 )
    Tmax= h * c / (k * L * log(1.0d0  + 2.0d0 *h*c**2/(L**5 * S)));
!
!     -----begin an iteration loop to find tau at 6000 K
!
    target = 6000
    do it = 0, 30
!
!     -----averge the last two guesses and recalculate temp
!
      tau = 0.5d0 * (tau_min + tau_max)
      S = i00L * (a0 + a1 * tau + 0.5d0 * a2 * tau**2 )
      T = h * c / (k * L * log(1.0d0  + 2.0d0 *h*c**2/(L**5 * S)));
!
!     -----adjust the lower or upper limit as appropriate to bracket
!          the root by bisection
!
      if (T < target ) then
        tau_min = tau
      else
        tau_max = tau
      endif
    enddo
    tau6000 = tau

!
!     -----repeat with temp = 6500K
!
    tau_min = 0
    tau_max = 3
    target = 6500
    do it = 0, 30
!
!     -----averge the last two guesses and recalculate temp
!
      tau = 0.5d0 * (tau_min + tau_max)
      S = i00L * (a0 + a1 * tau + 0.5d0 * a2 * tau**2 )
      T = h * c / (k * L * log(1.0d0  + 2.0d0 *h*c**2/(L**5 * S)));
!
!     -----adjust the lower or upper limit as appropriate to bracket
!          the root by bisection

      if (T < target ) then
        tau_min = tau
      else
        tau_max = tau
      endif
    enddo
    tau6500 = tau

!
!     -----write the result out
!
    write(*,'(i7,2e20.9)') lambda, tau6000, tau6500





  enddo
  close(10)
  stop



end program read