program fakephot
! ktk Jan.2009

       implicit none
       integer i, npts, japprox
       parameter(npts = 3001)
       double precision wav, Teff, dtau, q, a, h, c, k, c2, peak
       double precision tau(npts), T(npts)
       double precision plancklam, inten(npts), sum

* some constants CODATA (2006)

       c = 2.99792458d10    ! speed of light in the vacuum
       h = 6.62606896d-27    ! Planck's constant
       k = 1.3806504d-16     ! Boltzmann's constant

       c2 = ( h * c ) / k   ! another useful constant

       print *
       write(6,*)'Computes emergent continuum intensity I_lambda'  
       write(6,*)'at a specified wavelength vs. optical depth for'  
       write(6,*)'a toy model stellar photosphere.'
       write(6,*)'' 
       print *

       write(6,*)'enter effective temperature of photosphere'
       read(5,*) Teff

* the Wien Displacement Law Peak
       peak =  c2 * 0.2014052d8/Teff  ! Angstroms, peak of B_lambda
       write(6,12) peak 
12     format('the blackbody peak lies at ',1p,e10.4, ' Angstroms')
       print *

       write(6,*)'continuum wavelength (Ang.) for B_lambda computation'
       read(5,*) wav
       print *
*
69     write(6,*)'choose a simple analytic approximation for T(tau)'
       write(6,*)'(1) Eddington Approximation'
       write(6,*)'(2) Modified Schuster-Schwarzschild Approximation'
       write(6,*)'(3) Milne-Eddington Approximation'
       write(6,*)'(4) yet another simple approximation'
       write(6,*)'please choose either 1, 2, 3, or 4; enter number here'
       read(5,*) japprox
       print *

       if (japprox .eq. 1) then

* Eddington Approximation
         q = (2.d+00)/(3.d+00)
         a = 1./((2.d+00)/(3.d+00) + q) ! reach T_eff at optical depth 2/3

       elseif (japprox .eq. 2) then
* modified Schuster-Schwarzschild Approximation
         q = 1./sqrt(3.d+00)
         a = 1./((2.d+00)/(3.d+00) + q) ! reach T_eff at optical depth 2/3

       elseif (japprox .eq. 3) then 
* Milne-Eddington Approximation
         q = (2.d+00)/(3.d+00)  ! this is just a place holder, q = q(tau) below
         a = (1.d+00)/((2.d+00)/(3.d+00) + q) ! reach T_eff at optical depth 2/3
         ! these are placeholders for q(tau), a, computed below.

       elseif (japprox .eq. 4) then
* another approximation
         q = (1.d+00)/(3.d+00)
         a = (1.d+00)/((2.d+00)/(3.d+00) + q) ! reach T_eff at optical depth 2/3

       else 
         goto 69

       endif
*
       dtau = 0.01d+00  ! increment for optical depth

* zero out the arrays
       do i = 1, npts
         T(i) = 0.
         tau(i) = 0.
         inten(i) = 0.
       end do

       open(unit=16,file='fakephot.dat')
        write(16,90) Teff
90      format('# emergent I from toy photosphere, Teff = ',1pe10.4)
        write(16,91) wav
91      format('# intensity determined at ', 1pe10.4, ' Angstroms')
        write(16,92)
92      format('# opt. depth   T(tau)  I_lambda cgs/Angstrom   Sum')

* do the loop on optical depth
       sum = 0.
       tau(0) = -dtau  ! initialize the loop
       do i = 1, npts
         tau(i) = tau(i-1) + dtau  !the optical depth

         if (japprox .eq. 3) then  ! Milne-Eddington Approximation
           q = 0.7104d+00 - ( 0.1331d+00 * dexp((-3.449d+00)*tau(i)) )
           a = (1.d+00)/((2.d+00)/(3.d+00) + q) ! reach T_eff at tau = 2/3
         endif

* now compute the run of T(tau) for chosen approximation

       T(i) = (a * Teff*Teff*Teff*Teff*(tau(i) + q))**(0.25d+00)

* plancklam: a function to return BB intensity B_lambda in erg/s/cm^2/A

       inten(i) = plancklam(T(i), wav) * dexp(-tau(i))
       sum = sum + inten(i)*dtau

       write(16,100) tau(i), T(i), inten(i), sum
100    format(1x, 1p, 3(e10.4,2x), 6x, e11.5)
       enddo

       close(unit=16)

       write(6,93) sum
93     format('the summed intensity (cgs/A) is ', 1pe11.5)

       write(6,*) ' output sent to fakephot.dat'
       print*

       stop
       end

**************************************************
       double precision function plancklam(Temp, wav)

       double precision blam
       double precision c, h, k, pi, pi4
       double precision Temp, wav, xlam, xlam5, c1, c2

       c = 2.99792458d10    ! speed of light in the vacuum
       h = 6.62606896d-27    ! Planck's constant
       k = 1.3806504d-16     ! Boltzmann's constant

       pi = 3.1415926535897932d0  ! you had better know what this is
       pi4 = pi*pi*pi*pi
       c1 = 2.d0*h*c*c ! a useful constant
       c2 = ( h * c ) / k      ! another useful constant

* convert wavelength in Angstroms to cm
       xlam = wav / 1.0d8
       xlam5 = xlam * xlam * xlam * xlam * xlam

       blam = (c1/xlam5) / (dexp(c2/(xlam*Temp)) - 1.d0)
       plancklam = blam/1.0d+08   ! output in per unit Angstrom interval
       
       end
**************************************************