#Script for making various plots illustrating the Planck function
#
#Data on section of text which this script is associated with
Chapter = '3.**'
Figure = '**'
#

import phys,math
from ClimateUtilities import *

#Under construction.

#Functions for computing integrated emission
def fB(nu,Junk,param):
    #Note that we include the factor of pi here, so that B
    #will integrate up to the net emissions over all angles
    return math.pi*phys.B(nu,param.T)

#This function computes the cumulative emission for a range of
#N frequenices between specified frequencies nu0 and nu1, by integrating
#the Planck function over frequency.  It uses the Runge-Kutta
#integrator to do the integration, since that provides a
#simple way of getting the integral at a list of consecutive values.
def cumEmission(nu0,nu1,N,T):
    param = Dummy()
    param.T = T
    dnu = (nu1-nu0)/N
    Bcum = integrator(fB,nu0,0.,dnu)
    Bcum.setParams(param)
    BcumList = [0.]
    nuList = [nu0]
    nu = nu0
    while nu <= nu1:
        (nu,val) = Bcum.next()
        nuList.append(nu)
        BcumList.append(val)
    return nuList,BcumList

#Specify list of temperatures.
#Plot vs. wavenumber, and vs. wavelength.

Tlist = [250.,300.,350.,400.]

#Generate a good range of frequencies for plotting
#This assumes that the first temperature in the list
#is the minimum, and the last is the maximum
nuA = phys.k*Tlist[0]/phys.h
dnu = nuA/20.
nuB = phys.k*Tlist[-1]/phys.h
N = int(10.*nuB/dnu)
nuMax = (N-1)*dnu

nuList = [i*dnu for i in range(1,N)]
c = Curve()
wavenum = [.01*nu/phys.c for nu in nuList] #Wavenumber in 1/cm
c.addCurve(wavenum,'n')

#Note: , since we are plotting as a function
#of wavenumber, we rescale B, since it's a density.
for T in Tlist:
  Blist = [(phys.c/.01)*phys.B(nu,T) for nu in nuList]
  c.addCurve(Blist,'B%f'%T,'B(n,%f)'%T)

w= plot(c)

#Plot cumulative emission vs. wavenumber
c1 = Curve()
nu,P = cumEmission(dnu,nuMax,500,300.) #Just to get freq list
wavenum = [.01*x/phys.c for x in nu]
c1.addCurve(wavenum,'wavenum','Wavenumber (1/cm)')

for T in Tlist:
    nu,P = cumEmission(dnu,nuMax,500,T)
    c1.addCurve(P,'Cum%f'%T,'T = %f'%T)
plot(c1)