#---------------------------------------------------------------
#This script is like WaterVaporRealGasFluxSpectrum.py, and computes
#the spectrum of OLR for a real gas. The difference is that it
#computes it for a well mixed gas (CO2 by default), and doesn't
#generate plots with and without the continuum. 
#---------------------------------------------------------------




#Data on section of text which this script is associated with
Chapter = '4'
Section = '**'
Figure = None
#

import miniClimtFancy as rad
import phys
from ClimateUtilities import *

#Path to the workbook datasets
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
#Path to the exponential sum tables
EsumPath = datapath + 'Chapter4Data/ExpSumTables/'

#Initialize the radiation calculation, and over-ride defaults
rad.bandData = rad.loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.air.data')
rad.ForeignContinuum = rad.NoContinuum #Do the calculation without contin
rad.SelfContinuum = rad.NoContinuum #Do the calculation without contin
rad.BackgroundGas = phys.air #The transparent background gas
rad.GHG =phys.CO2 #The greenhouse gas
rad.Tstar = 0. #Ignore temperature scaling

def fplot(Tg,qGHG,n=20):
    #qGHG is the mass mixing ratio of the greenhouse gas
    #-------------------Initialize arrays
    psAir = 1.e5 #Partial pressure of air at surface
    ps = psAir # Assumes partial pressure of GHG is small
    ptop = 100.
    p = rad.setpLin(ps,ptop,n)
    #
    #Set temperature to moist adiabat
    #to saturation
    m = phys.MoistAdiabat()
    p,T,molarCon,qw = m(psAir,Tg,p)
    #
    #Set constant gas mixing ratio
    q = Numeric.ones(len(p),Numeric.Float)*qGHG
    #Now compute spectrum of OLR and back radiation
    #
    wave = []
    surf = []
    OLR = []
    for band in rad.bandData:
        flux,heat = rad.radCompBand(p,T,Tg,q,band)
        wave1 = (band.nu2+band.nu1)/2.
        Delta = (band.nu2-band.nu1)
        wave.append(wave1)
        surf.append(rad.Planck(wave1,Tg)-flux[-1]/Delta) #Convert flux to per wavenumber
        OLR.append(flux[0]/Delta)#Convert flux to per wavenumber
    B = [rad.Planck(wave1,Tg) for wave1 in wave]
    c = Curve()
    c.addCurve(wave,'Wavenumber')
    c.addCurve(surf,'surf')
    c.addCurve(OLR,'OLR')
    c.addCurve(B,'Planck')
    return c

#Do the plots
Tg = 300. #Specify surface temperature
ppmv = 10000.#GHG molar concentration, in ppmv
c = fplot(Tg, ppmv*1.e-6 * (29./44.))
c.PlotTitle = 'OLR Spectrum, CO2 = %f ppmv'%ppmv
w = plot(c)