#This script computes the transmission profile of
#incoming stellar radiation.  It can either be run
#with a well-mixed gas ignoring temperature scaling (using miniClimt)
#or with an inhomogeneous profile including temperature scaling
#(using miniClimtFancy)
#

#ToDo: Move the exponential sum tables to the Workbook data directory
#and modify the load statement

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

from ClimateUtilities import *
import phys
import miniClimt as radmodel
import miniClimtFancy as radmodelF
import planets

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


#Computation of out-of range stellar flux
mPlanck = romberg(radmodel.Planck)
def CompOutBandStellar(r):
    nu2 = r.bandData[-1].nu2
    return .25*(r.Lstellar/(phys.sigma*radmodel.Tstellar**4))*mPlanck([nu2,10.*nu2],radmodel.Tstellar)

#Set the star's characteristics
radmodelF.Tstellar = radmodel.Tstellar = 3000.
radmodelF.Lstellar = radmodel.Lstellar = 1400.
StarType = 'M'

#Flux computation for pure CO2

#Load the appopriate band data
#REMEMBER: Load the self-broadened or foreign-broadened tables,
#          depending on what kind of calculation you are doing
#(The following is a self-broadened table)
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableExtendedWave.260K.100mb.self.data')
#Set the continua
radmodel.Continuum = radmodel.CO2Continuum
q = 1. #Concentration for a pure gas

Tg = 260. #Since we don't do temperature scaling, this doesn't matter
n = 200
ps,ptop = 2.e5,100.
T = Numeric.zeros(n,Numeric.Float) + Tg
p = radmodel.setpLin(ps,ptop,n)
flux,heat = flux,heat = radmodel.radCompStellar(p,T,Tg,q)
flux = flux - CompOutBandStellar(radmodel) #Incoming flux is negative!

c1 = Curve()
c1.addCurve(p/100.,'p')
c1.addCurve(flux,'flux')
c1.Xlabel = 'Pressure (mb)'
c1.Ylabel = 'Flux (W/m**2)'
c1.switchXY = c1.reverseY = True
c1.PlotTitle = StarType+' star, Pure CO2'
plot(c1)
c1.dump('FluxProfilePureCO2'+StarType+'.txt')


#Now do computation for pure water vapor
radmodel.bandData = \
        radmodel.loadExpSumTable(EsumPath+'H2OTableExtendedWave.260K.100mb.self.data')
radmodel.Continuum = radmodel.H2OSelfContinuum
q = 1.

Tg = 260. #Since we don't do temperature scaling, this doesn't matter
n = 200
ps,ptop = 2.e5,100.
T = Numeric.zeros(n,Numeric.Float) + Tg
p = radmodel.setpLin(ps,ptop,n)
flux,heat = flux,heat = radmodel.radCompStellar(p,T,Tg,q)
flux = flux - CompOutBandStellar(radmodel) #Incoming flux is negative!

c2 = Curve()
c2.addCurve(p/100.,'p')
c2.addCurve(flux,'flux')
c2.Xlabel = 'Pressure (mb)'
c2.Ylabel = 'Flux (W/m**2)'
c2.switchXY = c2.reverseY = True
c2.PlotTitle = StarType+' star, Pure WV'
plot(c2)
c2.dump('FluxProfilePureWV'+StarType+'.txt')


#Now do computation for CO2 in air
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableExtendedWave.260K.100mb.air.data')
radmodel.Continuum = radmodel.CO2Continuum

eta = .2
q =  eta*44./(eta*44. + (1-eta)*29.)

Tg = 260. #Since we don't do temperature scaling, this doesn't matter
n = 200
ps,ptop = 2.e5,100. 
T = Numeric.zeros(n,Numeric.Float) + Tg
p = radmodel.setpLin(ps,ptop,n)
flux,heat = flux,heat = radmodel.radCompStellar(p,T,Tg,q)
flux = flux - CompOutBandStellar(radmodel) #Incoming flux is negative!

c3 = Curve()
c3.addCurve(p/100.,'p')
c3.addCurve(flux,'flux')
c3.Xlabel = 'Pressure (mb)'
c3.Ylabel = 'Flux (W/m**2)'
c3.switchXY = c3.reverseY = True
c3.PlotTitle = StarType+' star, %f molar CO2 in air'%eta
plot(c3)
c3.dump('FluxProfileCO2inAir'+StarType+'.txt')


#Now do computation for water vapor in air, on moist adiabat


radmodelF.bandData = \
        radmodelF.loadExpSumTable(EsumPath+'H2OTableExtendedWave.260K.100mb.air.data')
radmodelF.Continuum = radmodelF.H2OSelfContinuum
radmodelF.Tstar = 0. #Temperature scaling; set a better value?
radmodelF.SelfRat = 7. #Self to Air broadened ratio

Tg = 300.
n = 200
ps,ptop = 1.e5,100.

#Compute the moist adiabat
m = phys.MoistAdiabat()
p,T,MolarCon,q = m(ps,Tg,p)

flux,heat = flux,heat = radmodelF.radCompStellar(p,T,Tg,q)
flux = flux - CompOutBandStellar(radmodelF) #Incoming flux is negative!

c4 = Curve()
c4.addCurve(p/100.,'p')
c4.addCurve(flux,'flux%f'%Tg)
c4.PlotTitle = StarType+' star, Saturated WV in air, moist adiabat'
c4.Xlabel = 'Pressure (mb)'
c4.Ylabel = 'Flux (W/m**2)'
c4.switchXY = c4.reverseY = True

#Now do it again for a different temperature
Tg = 250.
#Compute the moist adiabat
m = phys.MoistAdiabat()
p,T,MolarCon,q = m(ps,Tg,p)

flux,heat = flux,heat = radmodelF.radCompStellar(p,T,Tg,q)
flux = flux - CompOutBandStellar(radmodelF) #Incoming flux is negative!

c4.addCurve(flux,'flux%fK'%Tg)
plot(c4)
c4.dump('FluxProfileWVinAir'+StarType+'.txt')