program plancklam
! ktk March.2007

       implicit none
       
       double precision T, xlam, xlam5, blam, flam, flux, peak
       double precision c, h, k, pi, pi4, c1, c2, sigma

* 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

       pi = 3.1415926535897932d0  ! you had better know what this is
       pi4 = pi*pi*pi*pi
       c1 = pi * 2.d0*h*c*c ! a useful constant, factor of pi from solid angle
       c2 = ( h * c ) / k      ! another useful constant
       sigma = (pi4/15.d0)*c1/(c2*c2*c2*c2)  ! Stefan-Boltzmann constant

       print *
       write(6,*)'computes pi*B(lambda) cgs flux per Angstrom interval'

       print *
       write(6,*)'enter T(eff) of blackbody in degrees K'
       read(5,*) T

       print *
       write(6,*)'enter lambda in Angstroms'
       read(5,*) xlam

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

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

       flux = sigma * T * T * T * T

* the Wien Displacement Law Peak
       peak =  c2 * 0.2014052d8/T  ! Angstroms, peak of F_lambda

       print *
       write(6,10) flam 
10     format('specific blackbody flux is ',1p,e10.4, ' ergs/s/cm^2/A')
       print *
       write(6,11) flux 
11     format('integrated blackbody flux is ',1p,e10.4, ' ergs/s/cm^2')
       print *
       write(6,12) peak 
12     format('the blackbody peak lies at ',1p,e10.4, ' Angstroms')
       print *

       end