#Data on section of text which this script is associated with
Chapter = '4.**'
Figure = '**'
#
#ToDo: Change values of kappa's, and include cos(thetabar) factor

import math

import phys
from ClimateUtilities import *


g = 9.8 # Acceleration of gravity on Earth


params = Dummy() #Empty object, to pass parameters

#Define Planck function
def Planck(wavenum,T):
    return 100.*math.pi*phys.c*phys.B(100.*wavenum*phys.c,T)
def dPlanck(wavenum,T):
    nu = 100.*wavenum*phys.c
    u = min(phys.h*nu/(phys.k*T),500.) #To prevent overflow
    dB = (2.*phys.h*nu**3/phys.c**2)*(phys.h*nu/phys.k)*(1/T**2)*math.exp(u)/(math.exp(u)-1.)**2
    return math.pi*100.*phys.c*dB


#Define absorption envelope here.  
def kappa(wave,wave0,kappa0,delta):
    if abs(wave-wave0)<= delta/2.:
        return kappa0
    elif abs(wave-wave0)<= delta1/2.:
        return kappa1
    else:
        return 0.


#Function to give temperature on dry adiabat. Can
#be replaced by any other function giving temperature
#as a function of pressure.
Rcp = 2./7.		
def T(p):
	return max(Tstrat,Ts*(p/ps)**Rcp)
def dTdp(p):
    if T(p)<= Tstrat:
        return 0.
    else:
        return -Rcp*(1./ps)*Ts*(p/ps)**(Rcp -1.)

#Note:  This could be sped up by using the absorption form instead
#of the transmissionform when the atmosphere is optically thin. It
#could be made faster also by choosing pLim so that the limit of integration
#is cut off where the integrand becomes small, in the optically thick case. 
def f(pPrime,params):
    wave = params.wave
    q = params.q
    return dPlanck(wave,T(pPrime))*math.exp(-q*kappa(wave,wave0,kappa0,delta)*pPrime/g)*dTdp(pPrime)
#Find OLR by a numerical quadrature
relTolerance = .01
def OLR(wave,q):
    tolerance = relTolerance*Planck(wave,Ts)
    pLim = ps
    params.wave = wave
    params.q = q
    quad = romberg(f)
    surfEmission = math.exp(-q*kappa(wave,wave0,kappa0,delta)*ps/g)
    BddTerm =  (Planck(wave,Ts)-Planck(wave,T(ps)))*surfEmission
    return Planck(wave,T(1.e-9)) + BddTerm + quad([pLim,1.e-9],params,tolerance)

#Computation of frequency-integrated net OLR
def fnet(wave,q):
    return OLR(wave,q)
def netOLR(q):
    tolerance = 1.
    limit = 5000.
    quad = romberg(fnet,10)
    return  quad([1.,limit],q,tolerance)

def spectrum(q):
    waveList = []
    BList = []
    BstratList = []
    OLRList = []
    for wave in [1. + 2.5*i for i in range(1000)]:
        waveList.append(wave)
        BList.append(Planck(wave,Ts))
        BstratList.append(Planck(wave,max(Tstrat,1.e-10))) #avoid overflow
        OLRList.append(OLR(wave,q))
    c = Curve()
    c.addCurve(waveList,'wavenum')
    c.addCurve(BList,'B','Planck Function Ts')
    c.addCurve(BstratList,'Bstrat','Planck Function Tstrat')
    c.addCurve(OLRList,'OLR','OLR')
    return c

#----------Main script starts here ------------------------
#Define constants needed
ps = 1.e5 #Surface pressure, in Pascals
Ts = 280. #Surface temperature
Tstrat = 200. #Stratospheric temperature
wave0 = 600. # center of absorption at 10 micron band
delta = 200. #Width of absorption band
delta1 = 600.
q0 = 300.e-6*(44./29.)
q = q0
kappa0 = 2.
kappa1 = .25

#A sample call:
c = spectrum(q)
w = plot(c)
#print 'The net OLR is',netOLR(),'W/m**2
#