#--------------------------------------------------------
#Description:
#Utilities for reading the HITRAN line database, synthesizing
#absorption spectrum on a regular wavenumber grid, and computing
#exponential sum coefficients. It also can be used in simple
#line-by-line transmission function calculations
#
#The functions of interest to most users are:
# loadSpectralLines(molName)
# Loads spectral line data for molecule molName.
#
# computeAbsorption(p,T,dWave,LineTailCutoff)
# Synthesizes the absorption coefficient vs. wavenumber
# from the line data, using an assumed pressure-broadened
# Lorentz line shape with tails cut off LineTailCutoff
# cm**-1 from the line centers.
#
# plotLogQuartiles(bandWidth)
# Makes the plots of absorption statistics in
# bands of width bandWidth, as discussed in
# the greenhouse gas spectral survey in Chapter 4.
#
# makeEsumTable(waveStart,waveEnd,Delta=50.,N=21)
# Constructs the exponential sum tables used in
# the homebrew radiation code.
#
#
#Scroll down to "Main script starts here" to see a basic example
#of use of the routines to plot absorption coefficients and make
#the plot of the quartile summaries in the spectral survey in the book.
#
#
#There are a few other things in PyTran for computing line-by-line
#band averaged transmission functions, and transmission along
#a path computed exactly at a fixed wavenumber without making
#pressure scaling assumptions. These are used in Chapter 4 to
#do some of the LBL transmission function calculations, but are
#probably of less general interest.
#
#Note that a limitation of PyTran is that it uses a cutoff
#Lorentz line shape to synthesize absorption from line data.
#This is mainly to keep things simple and easy to understand,
#and it covers most cases of interest adequately well. However,
#a "professional strength" code should use the Voigt line shape
#instead, to allow for the dominance of Doppler broadening near
#line centers when pressure is low. In addition, to get the line
#overlap right in high pressure situations (including Early Mars)
#one ought to consider a more careful treatment of the non-Lorentz
#far tails of collisionally broadened lines. The student who wishes
#to explore these effects (the latter of which is leading-edge research)
#will find it easy to modify the code to accomodate different assumptions
#on line shape.
#
#As currently written, Pytran loads in the lines for the dominant
#isotopologue for each molecule (based on abundance in Earth's
#atmosphere). If you want to modify the code to look at minor
#isotopologues, it is important to note that the HITRAN database
#downweights the line strengths for each isotopologue according
#to relative abundance in Earth's atmosphere.
#--------------------------------------------------------
#
#Change Log
# 3/13/2012: Corrected algebraic prefactor in temperature
# scaling of line strength, and put in a more general
# line-dependent scaling formula for the exponential factor
# In this release, a generic power-law dependence for the
# partition function Q(T) is used, but in the next release
# I will implement the exact Q(T) for selected molecules,
# based on routines provided as part of the HITRAN distribution.
#
# 6/10/2013: There was a bug in the corrected line strength scaling, which
# has now been fixed. See comments in computeAbsorption(...)
#
#---------------------------------------------------------
Chapter = '4'
#Path to the workbook datasets
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
#Path to the hitran by-molecule database
hitranPath = datapath+'Chapter4Data/hitran/ByMolecule/'
import string,math
from ClimateUtilities import *
import phys
#**ToDo: Add isotope handling option. Right now, the code just does
# the dominant isotope.
#------------Constants and file data------------
#
#Hitran field locations
fieldLengths = [2,1,12,10,10,5,5,10,4,8,15,15,15,15,6,12,1,7,7]
Sum = 0
fieldStart = [0]
for length in fieldLengths:
Sum += length
fieldStart.append(Sum)
iso = 1
waveNum = 2
lineStrength = 3
airWidth = 5
selfWidth = 6
Elow = 7
TExp = 8
#
#
#Total internal partition functions (or generic approximations).
#These are used in the temperature scaling of line strength.
#The generic partition functions are OK up to about 500K
#(somewhat less for CO2)
def QGenericLin(T): #Approx. for linear molecules like CO2
return T
def QGenericNonLin(T): #Approx for nonlinear molecules like H2O
return T**1.5
#**ToDo: Provide actual partition functions for CO2, H2O and CH4
#Molecule numbers and molecular weights
#Add more entries here if you want to do other
#molecules in the HITRAN database. These entries are
#for the major isotopologues, but by using different
#molecule numbers you can do other isotopologues.
#The dictionary below maps the molecule name to the HITRAN
#molecule number (see documentation) and corresponding molecular
#weight.
#
#**ToDo:*Correct this to allow for minor isotopomers.
# *Improve structure of the molecule dictionary,
# e.g. use objects instead of arrays, allow for isotopologues
# *Add in entries for the rest of the molecules
molecules = {} #Start an empty dictionary
molecules['H2O'] = [1,18.,QGenericNonLin] #Entry is [molecule number,mol wt,partition fn]
molecules['CO2'] = [2,44.,QGenericLin]
molecules['CO'] = [5,28.,QGenericLin]
molecules['CH4'] = [6,16.,QGenericNonLin]
molecules ['HF'] = [14,20.,QGenericLin]
molecules['HCl'] = [15,36.,QGenericLin]
molecules['SO2'] = [9,64.,QGenericNonLin]
#-----------------------------------------------
#Planck function
def Planck(wavenum,T):
return 100.*math.pi*phys.c*phys.B(100.*wavenum*phys.c,T)
#
#Function to make an exponential sum table
#
def makeEsumTable(waveStart,waveEnd,Delta=50.,N=21):
wave = waveStart
c = Curve()
while (wave+Delta ) <= waveEnd:
header1 = 'logKappa.%d.%d'%(wave,wave+Delta)
header2 = 'dH.%d.%d'%(wave,wave+Delta)
bins,hist,cum = logHisto([wave,wave+Delta],N)
c.addCurve(bins,header1)
c.addCurve(hist,header2)
wave += Delta
return c
#Note: Log quartile values are exponentiated for plotting
def plotLogQuartiles(bandWidth):
waveList = []
statsList = []
wave = waveStart
while wave + bandWidth <= waveEnd:
waveList.append(wave+bandWidth/2.)
statsList.append(logQuartiles([wave,wave+bandWidth])[0])
wave = wave+bandWidth
c = Curve()
c.addCurve(waveList,'Wavenumber')
c.addCurve([math.exp(stat[0]) for stat in statsList],'Min')
c.addCurve([math.exp(stat[1]) for stat in statsList],'q25')
c.addCurve([math.exp(stat[2]) for stat in statsList],'Median')
c.addCurve([math.exp(stat[3]) for stat in statsList],'q75')
c.addCurve([math.exp(stat[4]) for stat in statsList],'Max')
return c
def logQuartiles(waveRange):
i1 = numpy.searchsorted(waveGrid,waveRange[0])
i2 = numpy.searchsorted(waveGrid,waveRange[1])
logAbs = [math.log(kappa) for kappa in absGrid[i1:i2] if kappa>0.]
numZeros = (i2-i1) - len(logAbs)
return quartiles(logAbs),numZeros
def quartiles(data):
vs = numpy.sort(data)
N = len(data)
if N > 0:
return [vs[0],vs[N/4],vs[N/2],vs[3*N/4],vs[-1]]
else:
return [-1.e20,-1.e20,-1.e20,-1.e20,-1.e20]
#Utility to compute histogram of absorption data
#Returns bins, histogram and cumulative PDF
def histo(data,nBins=10,binMin=None,binMax=None):
#Flag bands with no lines in them
if len(data) == 0:
binMid = numpy.zeros(nBins-1,numpy.Float)
hist = numpy.zeros(nBins-1,numpy.Float)
hist[0] = 1.
binMid[0] = -1.e30
cumHist = numpy.zeros(nBins-1,numpy.Float)+1.
return binMid,hist,cumHist
#
if binMin == None:
binMin = min(data)
if binMax == None:
binMax = max(data)
d = (binMax-binMin)/(nBins-1)
bins = [binMin + d*i for i in range(nBins)]
histo = [0. for i in range(nBins)]
for value in data:
i = int((value-binMin)/d - 1.e-5)
if (i>=0) and (i0.]
return histo(logAbs,nBins)
def plotLogHisto(waveRange,data,nBins=10):
bins,hist,cum = logHisto(waveRange,data,nBins)
c = Curve()
c.addCurve(bins)
c.addCurve(hist)
return c
#Computes the absorption spectrum on a wave grid, by summing up
#contributions from each line. numWidths is the number
#of line widths after which the line contribution is cut off.
#Typically we use 100-1000 for Earth troposphere, but in low pressure
#(like Mars or upper strat) values as high as 10000 might be needed.
#The validity of the Lorentz shape at such large cutoffs is dubious.
#At The cutoff can affect the absorption significantly in the
#water vapor or CO2 window, or "continuum" regions
def computeAbsorption(p,T,dWave,numWidths = 1000.):
waveGrid = [waveStart + dWave*i for i in range(int((waveEnd-waveStart)/dWave))]
waveGrid = numpy.array(waveGrid)
absGrid = numpy.zeros(len(waveGrid),numpy.Float)
for i in range(len(waveList)):
n = waveList[i] # Wavenumber of the line
gam = gamList[i]*(p/1.013e5)*(296./T)**TExpList[i]
#Temperature scaling of line strength
Tfact = math.exp(-100.*(phys.h*phys.c/phys.k)*ElowList[i]*(1/T-1/296.))
#The following factor is usually pretty close to unity
#for lines that aren't far from the peak of the Planck spectrum
#for temperature T, but it can become important on the low frequency
#side, and is easy to incorporate.
Tfact1 = (1.- math.exp(-100.*(phys.h*phys.c/phys.k)*n/T))/ \
(1.- math.exp(-100.*(phys.h*phys.c/phys.k)*n/296.))
#The following has the incorrect algebraic prefactor used in
# the original version of PyTran (see Errata/Updates document)
#S = sList[i]*(T/296.)**TExpList[i]*Tfact
#The following is the corrected version, including also the
# low frequency factor Tfact1
#S = sList[i]*(Q(296.)/Q(T))*TExpList[i]*Tfact*Tfact1
#Preceding line didn't delete "*TExpList" factor. Results now
#checked against LMD kspectrum code, for Lorentz line case
#-->Corrected on 6/10/2013
S = sList[i]*(Q(296.)/Q(T))*Tfact*Tfact1
#
iw = int(len(waveGrid)*(n-waveStart)/(waveEnd-waveStart))
nsum = int(numWidths*gam/dWave)
i1 = max(0,iw-nsum)
i2 = min(len(waveGrid)-1,iw+nsum)
if i2>0:
dn = numpy.arange(i1-iw,i2-iw)*dWave
abs = S*gam/(math.pi*( dn**2 + gam**2))
absGrid[i1:i2] += abs
return waveGrid,absGrid
#Function to to compute absorption coefficient
#vs. pressure at a given wavenumber wavenum. and temperature T. Line parameter data
#is global and must be read in first.
def absPath(p,wavenum,T=296.,numWidths = 1000.):
width = .1*p/1.e5
i1 = numpy.searchsorted(waveList,wavenum-numWidths*width)-1
i2 = numpy.searchsorted(waveList,wavenum+numWidths*width)+1
i1 = max(i1,0)
i2 = min(i2,len(waveList))
lineCenters = waveList[i1:i2]
lineStrengths = sList[i1:i2]
lineWidths = gamList[i1:i2]
TExpList1 = TExpList[i1:i2]
ElowList1 = ElowList[i1:i2]
gam = lineWidths*(p/1.013e5)*(296./T)**TExpList1
#Temperature scaling of line strength
Tfact = numpy.exp(-100.*(phys.h*phys.c/phys.k)*ElowList1*(1/T-1/296.))
S = lineStrengths*(T/296.)**TExpList1*Tfact
terms = S*gam/ \
(math.pi*( (lineCenters-wavenum)**2 + gam**2))
return sum(terms)
#Function to compute a list of values of the transmission
#at a given frequency, for each of a set of pressures.
#This can be used to compute the band averaged transmission
#in a given band between wave1 and wave2 by calling transPath
#for a regular list of wavenumbers covering the interval, and summing
#the results. For the transmission computed in the book, we actually
#did this integration in a Monte-Carlo fashion, by picking random
#wavenumbers in the interval and accumulating the results, This
#has the advantage that you can watch the progress of the computation
#as it converges, and you get some information at once from the
#entire band. For a typical band of width 50 cm**-1, even 1000
#samples can give quite good results
#
#**ToDo: Modify this to take a list of pressure and temperature
#instead. Also, fix the generation of the pressure
#list, so it ends at p2 instead of running over. Perhaps also
#allow for changing T and numWidths in the absorption computation
def transPath(wave,p1,p2,q,ndiv):
def f(p,q):
return absPath(p,wave)*q/g
dtau1 = romberg(f)
dtau = 0.
dp = (p2-p1)/ndiv
p = p1
pList = [p1]
transList = [1.]
trans = 1.
while p <= p2:
dtau += dtau1([p,p+dp],q)
p += dp
pList.append(p)
transList.append(math.exp(-abs(dtau)))
return numpy.array(pList),numpy.array(transList)
#Returns a curve for making a plot of absorption vs. pressure
#at a given frequency
def plotP(n,T=296.,numWidths = 1000.):
press = [10.**(.25*i) for i in range(35)]
absp = [absPath(p,n,T,numWidths) for p in press]
c = Curve()
c.addCurve(press)
c.addCurve(absp)
c.YlogAxis = c.XlogAxis = True
return c
#Gets the fieldNum'th data item from a Hitran2004 record
def get(line,fieldNum):
return line[fieldStart[fieldNum]:fieldStart[fieldNum]+fieldLengths[fieldNum]]
#Computes the mean transmission between two specified
#wavenumbers, for a given absorber path. This
#version does not take into account pressure broadening.
#It should be modified to do so.
def TransBar(wave1,wave2,massPath):
i1 = numpy.searchsorted(waveGrid,wave1)
i2 = numpy.searchsorted(waveGrid,wave2)
trans = numpy.exp(-massPath*absGrid[i1:i2])
return numpy.average(trans)
def plotTransBar(wave1,wave2,massPathStart,massPathEnd):
nplot = 100
r = math.log(massPathEnd/massPathStart)/nplot
mPath = [massPathStart*math.exp(r*i) for i in range(nplot)]
trans = [TransBar(wave1,wave2,path) for path in mPath]
c = Curve()
c.addCurve(mPath)
c.addCurve(trans)
return c
#Function to compute OLR spectrum line-by-line This
#version is to demonstrate the concept. It doesn't take
#into account the pressure dependence of the absorption
def Tprof(pps):
return max(280.*pps**(2./7.),200.)
def OLRspec(pathMax):
ps = 1.e5
dp = .05*ps
transLast = 1.+ numpy.zeros(len(absGrid),numpy.Float)
OLR = numpy.zeros(len(absGrid),numpy.Float)
pp = 0.
while pp < ps:
pp += dp
path = pathMax*pp/ps
trans = numpy.exp(-path*absGrid)
TT = (Tprof(pp/ps)+ Tprof((pp-dp)/ps))/2.
B = numpy.array([Planck(w,TT) for w in waveGrid])
OLR += B*(transLast-trans)
transLast = trans
TT = Tprof(1.)
B = numpy.array([Planck(w,TT) for w in waveGrid])
OLR += B*trans
return OLR
#Computes a function smoothed over specified wavenumber band
def smooth(wAvg,data):
navg = int(wAvg/dWave)
nmax = len(waveGrid)-navg
return [numpy.average(data[i:i+navg]) for i in range(0,nmax,navg)]
#Open the Hitran file for the molecule in question,
#and read in the line data
#**ToDo: add isotope index to argument
# Note that HITRAN weights the line strengths
# according to relative abundance of each isotopologue
# in Earth's atmosphere. When loading minor isotopologues,
# it is important to undo this weighting, so that the
# correct abundance weighting for the actual mixture of
# interest can be computed.
def loadSpectralLines(molName):
global waveList,sList,gamList,gamAirList,gamSelfList,ElowList,TExpList,Q
molNum = molecules[molName][0]
Q = molecules[molName][2] #Partition function for this molecule
if molNum < 10:
file = hitranPath+'0%d_hit04.par'%molNum
else:
file = hitranPath+'%d_hit04.par'%molNum
f = open(file)
waveList = []
sList = []
gamList =[]
gamAirList = []
gamSelfList = []
TExpList = []
ElowList = [] #Lower state energy
#
count = 0
line = f.readline()
while len(line)>0:
count += 1
isoIndex = string.atoi(get(line,iso))
n = string.atof(get(line,waveNum))
S = string.atof(get(line,lineStrength))
El = string.atof(get(line,Elow))
#Convert line strength to (m**2/kg)(cm**-1) units
#The cm**-1 unit is there because we still use cm**-1
#as the unit of wavenumber, in accord with standard
#practice for IR
#
#**ToDo: Put in correct molecular weight for the
# isotope in question.
S = .1*(phys.N_avogadro/molecules[molName][1])*S
gamAir = string.atof(get(line,airWidth))
gamSelf = string.atof(get(line,selfWidth))
TemperatureExponent = string.atof(get(line,TExp))
if isoIndex == 1:
waveList.append(n)
sList.append(S)
#gamList.append(gamSelf) #Put in switch to control this
gamAirList.append(gamAir)
gamSelfList.append(gamSelf)
ElowList.append(El)
TExpList.append(TemperatureExponent)
#Read the next line, if there is one
line = f.readline()
#Convert the lists to numpy/Numeric arrays
waveList = numpy.array(waveList)
sList = numpy.array(sList)
#gamList = numpy.array(gamList)
gamAirList = numpy.array(gamAirList)
gamSelfList = numpy.array(gamSelfList)
ElowList = numpy.array(ElowList)
TExpList = numpy.array(TExpList)
f.close()
#====================================================================
#
# End of function definitions
# Main script starts here
#====================================================================
#Program data and initialization
#-->The following block gives some choices for
# what wavenumber range to cover. Uncomment the
# one you want, or put in another to suit your purposes
#Standard wavenumbers used for spectral survey
waveStart = 25.
waveEnd = 3000.
#Extended wavenumber range for runaway greenhouse
#waveStart = 25.
#waveEnd = 5000.
#Super extended range for magma ocean and near IR absorption
#waveStart = 25.
#waveEnd = 30000. #30000. for water vapor
#Water vapor window region
#waveStart = 400
#waveEnd = 1500
#
#Titan CH4 region
#waveStart = 10.
#waveEnd = 300.
#CO2 principal band (Relevent to Earth up to around 1% CO2)
#waveStart = 450.
#waveEnd = 1000.
#Venus near-IR emission
#waveStart = 2300.
#waveEnd = 20000.
#Zoom in to look at line structure
#waveStart = 500
#waveEnd = 520
#or:
#waveStart = 2250.
#waveEnd = 5000.
#-->Load spectral line data
molName = 'H2O' #Choose which molecule you want here
#Currently implemented choices are H2O,
#CO2, CO, CH4, HF, HCL and SO2
#but you can do anything HITRAN knows about
#by adding more entries to the molecule
#dictionary near the top of the script.
loadSpectralLines(molName)
#-->Now compute the absorption on the wave grid
#from the line data. Edit the following 2 lines
#for the desired base pressure and temperature.
#The plots of temperature scaling coefficients in
#the book were made by running the code for two different
#temperatures and fitting to an exponential dependence to
#get T*.
#
#p = 1.e4 #Standard value for book is 1.e4
#T = 260.#Standard value for book is 260.
p = 1.e4
T = 260.
#-->Decide whether you want to compute air-broadened
#or self-broadened absorption. The plots of self/air
#ratio in the book were done by running this script for
#each choice and computing the ratio of the absorption coefficients
gamList = gamAirList
#gamList = gamSelfList
#
#-->Set wavenumber interval for computing the absorption
#Smaller increments give more accuracy, but take longer and
#generate much bigger files of absorption coefficients
#dWave = .01*(p/1.e5) #Standard value
dWave = .1*(p/1.e5) #Coarse-resolution for quick-look
#-->Compute the absorption on the wavenumber grid
#Set nWidths to the number of line widths to go out to in
#superposing spectral lines to compute the absorption coefficient.
#This is a very crude implementation of the far-tail cutoff.
#I have been using 1000 widths as my standard value.
nWidths = 1000.
waveGrid,absGrid = computeAbsorption(p,T,dWave,nWidths)
#-->Once the absorption is computed, you can use other utilities to
#compute exponential sum coefficients, quartile plots, and so forth.
#Some of the utilities, like the computation of absorption vs. pressure
#and the computation of transmission at a fixed wavenumber, can be
#used without computing the absorption on the regular waveGrid
#
#The remainder of the script does some sample plots based on the
#absorption spectrum on waveGrid
c1 = Curve()
c1.addCurve(waveGrid)
c1.addCurve(absGrid+1.e-8) #Cheap but useful hack for log plot
c1.YlogAxis = True
c1.PlotTitle = 'Absorption cross section for '+molName
c1.Xlabel = 'Wavenumber (1/cm)'
c1.Ylabel = 'Cross section (m**2/kg)'
plot(c1)
#-->Make a spectral summary plot showing the quartiles of
#the absorption coefficient, plus max and min, in each
#band of width bandWidth (in 1/cm). These are the kind of
#plots shown in the spectral survey of greenhouse gases
#in Chapter 4.
bandWidth = 50.
c2 = plotLogQuartiles(bandWidth)
c2.Ylabel = 'Cross section (m**2/kg)'
c2.Xlabel = 'Wavenumber (1/cm)'
c2.PlotTitle = 'Logarithmic quartile summary for '+molName
c2.YlogAxis = True
#The following code is necessary to allow the
#results to be plotted on a logarithmic axis. Some
#of the absorption coefficients are zero, and because
#of limitations in Ngl, we need to reset those to an
#arbitrary small number in order to use a logarithmic
#vertical axis
variables = c2.listVariables()[1:] #List of dependent variables
for v in variables:
#The following line replaces zeros or negatives
c2[v] = numpy.where(c2[v]<= 0.,1.e-10,c2[v])
plot(c2)
#Compute the transmission for a typical path. This
#is just for display purposes, to give an idea of
#where things are optically thick or optically thin
#**ToDo: Compute this for an actual pressure-weighted
# path, ignoring temperature variation along path
# or find some better example plot.
#
path = 1. #A standard path in kg/m**2
trans = numpy.exp(-path*absGrid)
T=260.
PlanckRef = [Planck(w,T) for w in waveGrid]
PlanckRef = numpy.array(PlanckRef)/max(PlanckRef) #Normalize for plotting
c3 = Curve()
c3.addCurve(waveGrid)
c3.addCurve(trans,'Transmission')
c3.addCurve(PlanckRef,'Planck%f'%T)
c3.Xlabel = 'Wavenumber (1/cm)'
c3.Ylabel = 'Transmission or normalized Planck function'
c3.PlotTitle = 'Transmission through a path of %f kg/m**2 for %s '%(path,molName)
plot(c3)