#-------------------------------------------------------- #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)