#---------------------------------------------------------------
#Description:
#
#Plots spectra of outgoing infrared and back-radiation to
#the surface, for a water-vapor/air atmosphere.  Compares
#results with and without the water vapor continuum
#---------------------------------------------------------------

#Data on section of text which this script is associated with
Chapter = '4'
Section = 'section:continuum'
Figure = 'fig:WVradContinEffects'

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
#Load air-broadened exponential sum table
rad.bandData = rad.loadExpSumTable(EsumPath+'H2OTable.260K.100mb.air.data')
rad.SelfRat = 5. #set mean self-broadening ratio. 7 might be a better
                 #value if the water vapor amount isn't too high
#If you are doing a pure water atmosphere, it's better to load
#in the self-broadened data and use that. This is useful for
#runaway greenhouse calculations.
#rad.bandData = rad.loadExpSumTable(EsumPath+'H2OTable.260K.100mb.self.data')
#rad.SelfRat = 1.
#
rad.ForeignContinuum = rad.NoContinuum #No w.v. foreign continuum
rad.BackgroundGas = phys.air #The transparent background gas
rad.GHG =phys.H2O #The greenhouse gas
rad.Tstar = 0. #Ignore temperature scaling

def fplot(Tg,n=20,psAir = 1.e5):
    #-------------------Initialize arrays
    psat = phys.satvps_function(phys.water)
    ps = psAir + psat(Tg) #Total surface pressure for this temperature
    ptop = 100. #For high temperatures (350K or more) might need
                #to reduce ptop to perhaps .1 Pa
    p = rad.setpLin(ps,ptop,n)
    #p = rad.setpLog(ps,ptop,n) #Logarithmic option. Useful for high Tg
                                #Gives more resolution in upper atmosphere
    #
    #Set temperature to moist adiabat and mass mixing ratio
    #to saturation
    m = phys.MoistAdiabat()
    m.ptop = ptop #Set top pressure consistently for adiabat computation
    p,T,molarCon,q = m(psAir,Tg,p)
    #
    #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
#Do the plot with the continuum
rad.SelfContinuum = rad.H2OSelfContinuum
c = fplot(Tg,20) #Use the optional second argument to increase
                  #the resolution if you need to, e.g. fplot(Tg,100)
                  #to do the calculation with 100 points in vertical
                  #instead of default 20
c.PlotTitle = 'With Continuum'
plot(c)
#Do it without the continuum
rad.SelfContinuum = rad.NoContinuum
c1 = fplot(Tg,20)
c1.PlotTitle = 'Without Continuum'
plot(c1)