diff --git a/.gitignore b/.gitignore index b4d87e9..7d4259d 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,10 @@ *.so *.so.dSYM +# Packages +*.egg +*.egg-info + +.ipynb_checkpoints + *.DS_Store diff --git a/pyrads/Absorption_Crosssections_HITRAN2016.py b/pyrads/Absorption_Crosssections_HITRAN2016.py index 71946a9..9c85355 100644 --- a/pyrads/Absorption_Crosssections_HITRAN2016.py +++ b/pyrads/Absorption_Crosssections_HITRAN2016.py @@ -2,13 +2,13 @@ ### THIS SCRIPT IS BASED ON PyTran, WHICH IS PART OF THE COURSEWARE IN ### Pierrehumbert, 2010, Principles of Planetary Climate ### -### MODIFIED BY dkoll +### MODIFIED BY Daniel Koll and Andrew Williams #-------------------------------------------------------- #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 +#absorption spectrum on a regular wavenumber grid. +# It also can be used in simple #line-by-line transmission function calculations # #The functions of interest to most users are: @@ -21,27 +21,6 @@ # 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. @@ -79,8 +58,8 @@ # #--------------------------------------------------------- -import string,math -from .ClimateUtilities import * +import math +import numpy from . import phys import os @@ -92,7 +71,7 @@ #------------Constants and file data------------ # -#Hitran field locations +#Hitran field locations (for reading par files) fieldLengths = [2,1,12,10,10,5,5,10,4,8,15,15,15,15,6,12,1,7,7] Sum = 0 fieldStart = [0] @@ -106,189 +85,85 @@ 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 + + +def QGenericPowerLaw(T,power): + """ + 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) + + power = 1: approximately correct for CO2 + power = 1.5: approximately correct for nonlinear molecules like H2O + + ToDo: Provide actual partition functions for CO2, H2O and CH4 + """ + return T**(power) + + +""" +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. +""" molecules = {} #Start an empty dictionary -molecules['H2O'] = [1,18.,QGenericNonLin] #Entry is [molecule number,mol wt,partition fn] -molecules['CO2'] = [2,44.,QGenericLin] -molecules['O3'] = [3,48.,QGenericNonLin] -molecules['N2O'] = [4,44.,QGenericLin] +molecules['H2O'] = [1,18.,1.5] #Entry is [molecule number,mol wt,partition fn] +molecules['CO2'] = [2,44.,1.] +molecules['O3'] = [3,48.,1.5] +molecules['N2O'] = [4,44.,1.] -molecules['CH4'] = [6,16.,QGenericNonLin] -molecules['NH3'] = [11,17.,QGenericNonLin] # linear structure? +molecules['CH4'] = [6,16.,1.5] +molecules['NH3'] = [11,17.,1.5] # linear structure? -molecules['HCN'] = [23,27.,QGenericNonLin] -molecules['C2H6'] = [27,30.,QGenericNonLin] +molecules['HCN'] = [23,27.,1.5] +molecules['C2H6'] = [27,30.,1.5] -molecules['SF6'] = [30,146.,QGenericNonLin] # careful: old file! +molecules['SF6'] = [30,146.,1.5] # careful: old file! #----------------------------------------------- -#Planck function def Planck(wavenum,T): + """ Planck Function """ 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(waveGrid,getGamma,p,T,dWave,numWidths = 1000.): - absGrid = numpy.zeros(len(waveGrid),numpy.Float) +def computeAbsorption(waveGrid,waveStart,waveEnd, + hitran_data,getGamma, + p,T,dWave, + nWidths = 1000.): + """ + 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 + """ + absGrid = numpy.zeros(len(waveGrid)) + waveList,sList,gamAirList,gamSelfList,ElowList,TExpList,QExponent = hitran_data + for i in range(len(waveList)): n = waveList[i] # Wavenumber of the line - #gam = gamList[i]*(p/1.013e5)*(296./T)**TExpList[i] # DKOLL: old - gam = getGamma(i)*(296./T)**TExpList[i] # DKOLL: new. getGamma includes p-scaling + gam = getGamma(i)*(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 - # + + S = sList[i] * (QGenericPowerLaw(296,QExponent)/QGenericPowerLaw(T,QExponent)) * Tfact*Tfact1 + iw = int(len(waveGrid)*(n-waveStart)/(waveEnd-waveStart)) - nsum = int(numWidths*gam/dWave) + nsum = int(nWidths*gam/dWave) i1 = max(0,iw-nsum) i2 = min(len(waveGrid)-1,iw+nsum) if i2>0: @@ -298,41 +173,34 @@ def computeAbsorption(waveGrid,getGamma,p,T,dWave,numWidths = 1000.): return absGrid -### DKOLL: add option to have a fixed cutoff. -### i.e., truncate line at N cm^-1 away from center instead of N halfwidths -### For example, MT_CKD continuum is defined as everything beyond 25cm^-1. -### -### DKOLL: also allow for option to remove the Lorenz line 'plinth', -## cf. MTCKD continuum references -def computeAbsorption_fixedCutoff(waveGrid,getGamma,p,T,dWave,numWidths=25.,remove_plinth=False): - absGrid = numpy.zeros(len(waveGrid),numpy.Float) +def computeAbsorption_fixedCutoff(waveGrid,waveStart,waveEnd, + hitran_data,getGamma, + p,T,dWave, + nWidths=25.,remove_plinth=False): + """ + DKOLL: add option to have a fixed cutoff. + i.e., truncate line at N cm^-1 away from center instead of N halfwidths + For example, MT_CKD continuum is defined as everything beyond 25cm^-1. + + DKOLL: also allow for option to remove the Lorenz line 'plinth', + cf. MTCKD continuum references + """ + absGrid = numpy.zeros(len(waveGrid)) + waveList,sList,gamAirList,gamSelfList,ElowList,TExpList,QExponent = hitran_data for i in range(len(waveList)): n = waveList[i] # Wavenumber of the line - #gam = gamList[i]*(p/1.013e5)*(296./T)**TExpList[i] # DKOLL: old - gam = getGamma(i)*(296./T)**TExpList[i] # DKOLL: new. getGamma includes p-scaling + gam = getGamma(i)*(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 - # + + S = sList[i]*(QGenericPowerLaw(296,QExponent)/QGenericPowerLaw(T,QExponent))*Tfact*Tfact1 + iw = int(len(waveGrid)*(n-waveStart)/(waveEnd-waveStart)) - #nsum = int(numWidths*gam/dWave) # DKOLL: old - nsum = int( numWidths/dWave ) # DKOLL: new + nsum = int( nWidths/dWave ) i1 = max(0,iw-nsum) i2 = min(len(waveGrid)-1,iw+nsum) # DKOLL: old @@ -343,152 +211,42 @@ def computeAbsorption_fixedCutoff(waveGrid,getGamma,p,T,dWave,numWidths=25.,remo # DKOLL: new elif (i2>0) and (remove_plinth==True): dn = numpy.arange(i1-iw,i2-iw)*dWave - abs = S*gam/math.pi * (1./(dn**2 + gam**2) - 1./(numWidths**2 + gam**2) ) + abs = S*gam/math.pi * (1./(dn**2 + gam**2) - 1./(nWidths**2 + gam**2) ) absGrid[i1:i2] += abs return 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 = gamAirList[i1:i2] #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): + """ + Gets the fieldNum'th data item from a Hitran2016 record + """ 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. -#**DKOLL: -# modified to speed up long line lists. -# skip lines until you're close to your spectral region of interest. -# then lead lines -# if you're outside region, break, skip rest of file. -# for this to work, make use of fact that line lists are ordered + def loadSpectralLines(molName,minWave=None,maxWave=None): - global waveList,sList,gamAirList,gamSelfList,ElowList,TExpList,Q + """ + 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. + + DKOLL: + modified to speed up long line lists. + skip lines until you're close to your spectral region of interest. + then lead lines + if you're outside region, break, skip rest of file. + for this to work, make use of fact that line lists are ordered + """ + molNum = molecules[molName][0] - Q = molecules[molName][2] #Partition function for this molecule + QExponent = molecules[molName][2] # Exponent for this molecule's partition function if molNum < 10: file = hitranPath+'0%d_hit16.par'%molNum elif molNum == 30: @@ -522,9 +280,7 @@ def loadSpectralLines(molName,minWave=None,maxWave=None): line = f.readline() while len(line)>0: count += 1 - #isoIndex = string.atoi(get(line,iso)) isoIndex = int(get(line,iso)) - #n = string.atof(get(line,waveNum)) n = float(get(line,waveNum)) #DKOLL: if nmaxWave: break # ignore rest of file else: - #S = string.atof(get(line,lineStrength)) S = float(get(line,lineStrength)) - #El = string.atof(get(line,Elow)) El = float(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 @@ -544,11 +299,8 @@ def loadSpectralLines(molName,minWave=None,maxWave=None): #**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)) gamAir = float(get(line,airWidth)) - #gamSelf = string.atof(get(line,selfWidth)) gamSelf = float(get(line,selfWidth)) - #TemperatureExponent = string.atof(get(line,TExp)) TemperatureExponent = float(get(line,TExp)) if isoIndex == 1: waveList.append(n) @@ -557,8 +309,10 @@ def loadSpectralLines(molName,minWave=None,maxWave=None): 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) @@ -567,41 +321,34 @@ def loadSpectralLines(molName,minWave=None,maxWave=None): ElowList = numpy.array(ElowList) TExpList = numpy.array(TExpList) f.close() + + hitran_lists = [waveList,sList,gamAirList,gamSelfList,ElowList,TExpList,QExponent] + return hitran_lists #==================================================================== -# -# End of function definitions # Main script starts here #==================================================================== -#Program data and initialization - - -#Standard wavenumbers used for spectral survey - -def getKappa_HITRAN(waveGrid,wave0,wave1,delta_wave,molecule_name,\ - press=1e4,temp=300.,lineWid=1000.,broadening="mixed", \ - press_self=None, \ - cutoff_option="relative",remove_plinth=False): - - # ! not a great solution ! - global waveStart, waveEnd, dWave, p, T, molName - waveStart = wave0 - waveEnd = wave1 - dWave = delta_wave - molName = molecule_name +def getKappa_HITRAN(waveGrid,waveStart,waveEnd,dWave, \ + press=1e4,temp=300.,lineWid=1000.,broadening="mixed", \ + press_self=None, molName=None, hitran_data=None, \ + cutoff_option="relative",remove_plinth=False): p = float(press) # DKOLL: make sure (p,T) are floats! T = float(temp) - loadSpectralLines(molName,minWave=wave0,maxWave=wave1) - -#-->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 + # Only run loadSpectralLines if necessary + if (molName is not None) and (hitran_data is None): + print("Running loadSpectralLines...") + hitran_data = loadSpectralLines(molName,minWave=waveStart,maxWave=waveEnd) + + + #-->Decide whether you want to compute air-broadened + #or self-broadened absorption. + gamAirList,gamSelfList = hitran_data[2], hitran_data[3] + p_tot = press/1.013e5 # need [atm]! p_self = press_self/1.013e5 p_air = (press-press_self)/1.013e5 @@ -615,17 +362,17 @@ def getKappa_HITRAN(waveGrid,wave0,wave1,delta_wave,molecule_name,\ if press_self==None: press_self = press getGamma = lambda i: gamAirList[i]*p_air + gamSelfList[i]*p_self + -#-->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 = lineWid + #-->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. if cutoff_option=="relative": - absGrid = computeAbsorption(waveGrid,getGamma,p,T,dWave,nWidths) + absGrid = computeAbsorption(waveGrid,waveStart,waveEnd,hitran_data,getGamma,p,T,dWave,nWidths= lineWid) elif cutoff_option=="fixed": - absGrid = computeAbsorption_fixedCutoff(waveGrid,getGamma,p,T,dWave,nWidths,remove_plinth=remove_plinth) + absGrid = computeAbsorption_fixedCutoff(waveGrid,waveStart,waveEnd,hitran_data,getGamma,p,T,dWave,nWidths= lineWid,remove_plinth=remove_plinth) return absGrid diff --git a/pyrads/ClimateGraphics.py b/pyrads/ClimateGraphics.py deleted file mode 100644 index c910627..0000000 --- a/pyrads/ClimateGraphics.py +++ /dev/null @@ -1,203 +0,0 @@ -from __future__ import division, print_function, absolute_import -#----------Section 3: Plotting utilities------------------------------- -#These need to be localized, for systems that don't support Ngl -# -#This is imported into ClimateUtilities. If you want to use a -#different graphics package (e.g. MatPlotLib) as a graphics driver -#in place of Ngl, you only need to rewrite this module. The most -#important thing to rewrite is the plot(...) function, which takes -#a Curve object as input and produces a line plot. The plot(...) -#function returns a plotObj object to allow further manipulation -#(mostly saving the plot in various formats), but if you have a plotting -#package that provides other ways of saving the plots, or if you just -#want to save the plots by doing screen dumps, you can just have -#the plot(...) function return None, or some dummy object. -# -#The other main function to customize is contour(...) which produces -#contour plots of an array. This function is not very extensively -#used in the courseware, so it is a lower priority for modification. -#----------------------------------------------------------------------- - - -#Note: On some older installations, Ngl has to -#be imported as PyNGL instead of Ngl. If there -#is a trouble with the import, fiddle with that. -try: - import Ngl -except: - try: - import PyNGL as Ngl - except: - print( "PyNGL was not found on your system.") - print( "You will not be able to use graphics functions.") - print( "PyNGL is available for Mac OSX,Linux and Windows/CygWin") - print( "Alternately, modify the module ClimateGraphics.py") - print( "to use some other graphics driver") - raise('Graphics Import Error') - -#A dummy class useful for passing parameters. -#Just make an instance, then add new members -#at will. Not sure why it also has to be defined -#here, since it is also defined in ClimateUtilities, -#but it seems to be necessary -class Dummy: - pass - -# A little class to make resources for Ngl plot options -class resource: - def __init__(self): - self.trYReverse = False - self.trXReverse = False - self.trXLog = False - self.trYLog = False - self.nglFrame = False - self.nglDraw = False - #ToDo: Missing data code resource, line styles, line colors - -# A little class for use as a return object from plot(), so -# the user has an easy way to delete a window or save a plot -# to an eps file. -class plotObj: - def __init__(self,workstation,plot,WorkstationResources=None): - self.workstation = workstation - self.plot = plot - self.WorkstationResources = WorkstationResources - #Deletes a plot window - def delete(self): - Ngl.destroy(self.workstation) - def save(self,filename = 'plot'): - #'eps' workstation doesn't work reliably, but 'ps' is OK - weps = Ngl.open_wks('ps',filename,self.WorkstationResources) - Ngl.change_workstation(self.plot,weps) - Ngl.draw(self.plot) - Ngl.change_workstation(self.plot,self.workstation) - Ngl.destroy(weps) - -#ToDo: -# *Implement use of missing data coding. -# *Provide some way to re-use window (e.g. by -# returning w, and having it be an optional -# argument. How to clear a window for re-use? -# *Make some use of dashed lines in line styles -# *The behavior of axis options like reverseX and -# XlogAxis is confusing when switchXY = True. We -# need to think about the semantics of what we -# mean by the "X" axis, and stick to it. As -# currently implemented the "X" axis means the -# horizontal axis for these options. (However, -# the semantics is different for the labeling options!) -# If we change this (and we probably should), the -# scripts plotting soundings will also need to be changed. -# as well as the Workbook intructions and problem sets. -# -lineColors =range(3,203,20)+range(3,203,20) -#Line colors for old versions of Ngl -#lineColors =range(3,16)+range(3,16) -lineThickness = [2,3,4,2,3,4,2,3,4,2,3,4,2,3,4,2,3,4,2,3,4] -plotSymbols = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5] -def plot(c): - r = resource() - r.trYReverse = c.reverseY - r.trXReverse = c.reverseX - r.trXLog = c.XlogAxis - r.trYLog = c.YlogAxis - # - #Line styles - r.xyLineColors = lineColors - r.xyLineThicknesses = lineThickness - #Plot title - r.tiMainString = c.PlotTitle - #Axis labels (ToDo: add defaults) - #X and Y axis labels - r.tiXAxisString = c.Xlabel - r.tiYAxisString = c.Ylabel - if c.switchXY: - r.tiXAxisString,r.tiYAxisString = r.tiYAxisString,r.tiXAxisString - - # Legends, for multicurve plot - legends = [] - for id in c.listVariables(): - if not id == c.Xid: - if len(c.label[id]) > 0: - legends.append(c.label[id]) - else: - legends.append(id) - #ToDo: Add option to skip legends - #Legends are redundant if there's only one curve - if len(legends) > 1: - r.pmLegendDisplayMode = "Always" - r.xyExplicitLegendLabels = legends - # - #Suppress line drawing and just plot symbol for scatter plot curves - r.xyMarkers = plotSymbols - r.xyMarkerColors = lineColors - r.xyMarkerSizeF = .01 - r.xyMarkLineModes = [] - for id in c.listVariables(): - if not id == c.Xid: - if c.scatter[id]: - r.xyMarkLineModes.append('Markers') - else: - r.xyMarkLineModes.append('Lines') - # - w = Ngl.open_wks('x11','Climate Workbook') - if c.switchXY: - plot = Ngl.xy(w,c.Y(),c.X(),r) - else: - plot = Ngl.xy(w,c.X(),c.Y(),r) - # - #Now draw the plot - r.nglDraw = True - Ngl.panel(w,[plot],[1,1],r) - return plotObj(w,plot) #So that user can delete window or save plot - -# A basic contour plotter, which will plot a contour plot -# of a Numeric array. The x and y scales can optionally -# be specified using keyword arguments x and y. For example, -# if we want the x scale to be the array (or list) lat, and -# the y scale to be the array (or list) lon, we would call -# contour as contour(A,x=lat,y=lon). -def contour(A,**kwargs): - #The following allows an expert user to pass - #Ngl options directly to the plotter. - #ToDo: Note that - #options explicitly specified later will over-ride - #what the user specified. This should be fixed, - #by checking for things already specified. We should - #also allow for using this resource to specify a color - #map. - if 'resource' in kwargs.keys(): - r = kwargs['resource'] - else: - r = Dummy() - # - r.cnFillOn = True #Uses color fill - # Set axes if they have been specified - # as keyword arguments - if 'x' in kwargs.keys(): - r.sfXArray = kwargs['x'] - if 'y' in kwargs.keys(): - r.sfYArray = kwargs['y'] - # - # Now create the plot - - rw = Dummy() - #Set the color map - if 'colors' in kwargs.keys(): - if (kwargs['colors'] == 'gray') or (kwargs['colors'] == 'grey') : - #Set the default greyscale - rw.wkColorMap = 'gsdtol' - else: - rw.wkColorMap = kwargs['colors'] - else: - #Default rainbow color table - rw.wkColorMap = "temp1" - - w = Ngl.open_wks('x11','Climate Workbook',rw) - r.nglDraw = False - r.nglFrame = False - plot = Ngl.contour(w,A,r) - #Now draw the plot - r.nglDraw = True - Ngl.panel(w,[plot],[1,1],r) - return plotObj(w,plot,rw) #So user can delete or save plot diff --git a/pyrads/ClimateGraphicsMPL.py b/pyrads/ClimateGraphicsMPL.py deleted file mode 100644 index 06f19b9..0000000 --- a/pyrads/ClimateGraphicsMPL.py +++ /dev/null @@ -1,241 +0,0 @@ -from __future__ import division, print_function, absolute_import -#----------Section 3: Plotting utilities------------------------------- -#Graphics interface customized for MatPlotLib -# --------------->UNDER DEVELOPMENT -# -#This is imported into ClimateUtilities. If you want to use a -#different graphics package (e.g. MatPlotLib) as a graphics driver -#in place of Ngl, you only need to rewrite this module. The most -#important thing to rewrite is the plot(...) function, which takes -#a Curve object as input and produces a line plot. The plot(...) -#function returns a plotObj object to allow further manipulation -#(mostly saving the plot in various formats), but if you have a plotting -#package that provides other ways of saving the plots, or if you just -#want to save the plots by doing screen dumps, you can just have -#the plot(...) function return None, or some dummy object. -# -#The other main function to customize is contour(...) which produces -#contour plots of an array. This function is not very extensively -#used in the courseware, so it is a lower priority for modification. -#This routine has not yet been implemented for MatPlotLib -#----------------------------------------------------------------------- - - -#Try to import matplotlib plotting routines -try: - import matplotlib as mpl - #Configure the backend so things work properly - #with the Python intepreter and idle. Not needed if you are - #using ipython. To work properly with idle, idle needs - #to be started up with the command "idle =n". Note that - #this trick limits the options for saving a file; only png - #format is supported. If you want higher resolution formats - #(notably eps) you should run using ipython -pylab , and eliminate - #the following two commands. - mpl.rcParams['backend'] = 'TkAgg' - mpl.rcParams['interactive'] = True - # - import pylab as pl -except: - print( 'matplotlib not found on your system') - print( 'You can still run the courseware, but') - print( 'will not be able to plot from inside Python') - -#A dummy class useful for passing parameters. -#Just make an instance, then add new members -#at will. Not sure why it also has to be defined -#here, since it is also defined in ClimateUtilities, -#but it seems to be necessary -class Dummy: - pass - -# A little class to make resources for Ngl plot options -class resource: - def __init__(self): - self.trYReverse = False - self.trXReverse = False - self.trXLog = False - self.trYLog = False - self.nglFrame = False - self.nglDraw = False - #ToDo: Missing data code resource, line styles, line colors - -# A little class for use as a return object from plot(), so -# the user has an easy way to delete a window or save a plot -# to an eps file. Plot objects not implemented yet for MPL -# Is there a way to save plots non-interactively in MPL? -class plotObj: - def __init__(self,workstation,plot,WorkstationResources=None): - self.workstation = workstation - self.plot = plot - self.WorkstationResources = WorkstationResources - #Deletes a plot window - def delete(self): - print( "In MatPlotLib just click the goaway button to dispose a plot") - def save(self,filename = 'plot'): - #**ToDo: Can we implement non-interactive save in MPL? - print( "In MatPlotLib save plots interactively using the file button") - print( " in the plot window") - -#ToDo: -# *Implement use of missing data coding. -# *Provide some way to re-use window (e.g. by -# returning w, and having it be an optional -# argument. How to clear a window for re-use? -# *Make some use of dashed lines in line styles -# *The behavior of axis options like reverseX and -# XlogAxis is confusing when switchXY = True. We -# need to think about the semantics of what we -# mean by the "X" axis, and stick to it. As -# currently implemented the "X" axis means the -# horizontal axis for these options. (However, -# the semantics is different for the labeling options!) -# If we change this (and we probably should), the -# scripts plotting soundings will also need to be changed. -# as well as the Workbook intructions and problem sets. -# - -#List of line colors and line styles to use -lineColors = ['b','g','r','c','m','y','k'] -lineStyles = ['-','--','-.',':'] -lineThickness = [2,3,4] -plotSymbols = ['.','o','v','<','s','*','+','x'] -def plot(c): - fig = pl.figure() #Always start a new figure - #**ToDo: Make sure log-axis options so they work consistently - # with Ngl implementation when x/y axes are switched. - # (Looks OK, but is that implementation the logical way?) - #r.trXLog = c.XlogAxis - #r.trYLog = c.YlogAxis - if c.XlogAxis: - pl.semilogx() - if c.YlogAxis: - pl.semilogy() - if c.XlogAxis & c.YlogAxis: - pl.loglog() - # - #Line styles (Not needed in MPL, handled automatically) - #r.xyLineColors = lineColors - #r.xyLineThicknesses = lineThickness - #Plot title - #r.tiMainString = c.PlotTitle - pl.title(c.PlotTitle) - #Axis labels (ToDo: add defaults) - #X and Y axis labels - #r.tiXAxisString = c.Xlabel - #r.tiYAxisString = c.Ylabel - if c.switchXY: - pl.ylabel(c.Xlabel) - pl.xlabel(c.Ylabel) - else: - pl.xlabel(c.Xlabel) - pl.ylabel(c.Ylabel) - - # Legends, for multicurve plot - legends = [] - for id in c.listVariables(): - if not id == c.Xid: - if len(c.label[id]) > 0: - legends.append(c.label[id]) - else: - legends.append(id) - #ToDo: Add option to skip legends - # - #Suppress line drawing and just plot symbol for scatter plot curves - #**ToDo: Implement for MatPlotLib - #r.xyMarkers = plotSymbols - #r.xyMarkerColors = lineColors - #r.xyMarkerSizeF = .01 - #r.xyMarkLineModes = [] - formatList = [] - count = 0 - for id in c.listVariables(): - if not id == c.Xid: - if c.scatter[id]: - #r.xyMarkLineModes.append('Markers') - color = lineColors[count%len(lineColors)] - symbol = plotSymbols[count%len(plotSymbols)] - formatList.append(color+symbol) - else: - color = lineColors[count%len(lineColors)] - style = lineStyles[count%len(lineStyles)] - formatList.append(color+style) - count += 1 - # - plotList = [] #Mainly so we can add legends. Could be done with label = ... - if c.switchXY: - count = 0 - for data in c.Y(): - plotList.append(pl.plot(data,c.X(),formatList[count])) - count += 1 - else: - count = 0 - for data in c.Y(): - plotList.append(pl.plot(c.X(),data,formatList[count])) - count += 1 - #Do the legends - pl.legend(plotList,legends) - # - #Do the axis reversal - #We do it here, since we don't know the axis limits until - #plotting is done - #r.trYReverse = c.reverseY - #r.trXReverse = c.reverseX - axes = pl.gca() #Gets current axis - if c.reverseX: - axes.set_xlim(axes.get_xlim()[::-1]) #::-1 reverses the array - if c.reverseY: - axes.set_ylim(axes.get_ylim()[::-1]) - #Now re-draw the plot - pl.draw() - #(Insert commands needed to show plot, if necessary) - return plotObj(None,fig) #Eventually we will use this to make subplots and do save option - -# A basic contour plotter, which will plot a contour plot -# of a Numeric array. The x and y scales can optionally -# be specified using keyword arguments x and y. For example, -# if we want the x scale to be the array (or list) lat, and -# the y scale to be the array (or list) lon, we would call -# contour as contour(A,x=lat,y=lon). -def contour(A,**kwargs): - #**ToDo: Add labeled contour lines, option to change contour levels, - # and to change palette. - # - fig = pl.figure() #Always start a new figure - # Set axes if they have been specified - # as keyword arguments - if 'x' in kwargs.keys(): - x = kwargs['x'] - else: - x = range(A.shape[1]) - if 'y' in kwargs.keys(): - y = kwargs['y'] - else: - y = range(A.shape[0]) - cs = pl.contourf(x,y,A) - cbar = pl.colorbar(cs) - return plotObj(None,fig) -## #The following allows an expert user to pass -## #Ngl options directly to the plotter. -## #ToDo: Note that -## #options explicitly specified later will over-ride -## #what the user specified. This should be fixed, -## #by checking for things already specified. We should -## #also allow for using this resource to specify a color -## #map. -## if 'resource' in kwargs.keys(): -## r = kwargs['resource'] -## else: -## r = Dummy() -## -## rw = Dummy() -## #Set the color map -## if 'colors' in kwargs.keys(): -## if (kwargs['colors'] == 'gray') or (kwargs['colors'] == 'grey') : -## #Set the default greyscale -## rw.wkColorMap = 'gsdtol' -## else: -## rw.wkColorMap = kwargs['colors'] -## else: -## #Default rainbow color table -## rw.wkColorMap = "temp1" diff --git a/pyrads/ClimateUtilities.py b/pyrads/ClimateUtilities.py deleted file mode 100644 index 8a3663c..0000000 --- a/pyrads/ClimateUtilities.py +++ /dev/null @@ -1,908 +0,0 @@ -from __future__ import division, print_function, absolute_import -### THIS SCRIPT IS TAKEN FROM THE COURSEWARE OF -### Pierrehumbert, 2010, Principles of Planetary Climate -### -### Downloaded from: -### https://geosci.uchicago.edu/~rtp1/PrinciplesPlanetaryClimate/Courseware/coursewarePortal.html - -#ToDo: *Check for right column length in setitem and addCurve -# -# *Implement show,hide for curves (in X() and Y()) -# *Implement missing data coding . -# self.setMissingDataCode(code) (char or numeric) -# possibly translate, check and force consistency -# *Possibly handle missing data for plotting using a -# fill() or interp() method to fill in missing data using -# the interpolation routine. The best way to handle -# missing data in computations is by using Masked Arrays -# -# *Make and import a PathFinder module, which the instructor -# will customize for the site. This module will help the -# students find locations of datasets and chapter scripts. -# Could provide "links" to find script that produced a given -# figure in the text. Alternately, we could just -# define directory strings like WorkbookDatasets here in -# ClimateUtilities. -# -# *Note: Other courseware modules should avoid -# importing numpy or Numeric directly. Get it -# by using "from ClimateUtilities import *" so -# that the preferred version is imported automatically -# -# *Looking ahead to Python 3, all print statements should -# be changed to some kind of "message" function, which -# can invoke either python 2.xx "print" or python 3 print(...) -# -# - -#----------------------------------------------- -# -#Import array package -# -#First try to import numpy, then fall back -#on Numeric. Arrange things so that the array -#package can be interchangeably referenced as -#either numpy.* or Numeric.* , though it is really -#numpy that is fully supported. (Almost everything -#will continue to work with Numeric, though, which -#may be useful for older installations) -#------------------------------------------------ - -try: - import numpy - import numpy as Numeric #For backwards compatibility - numpy.Float = numpy.float #For backwards compatibility -except: - try: - import Numeric - import Numeric as numpy #For frontwards compatibility - numpy.float = numpy.Float #For frontwards compatibility - print( "numpy not found. Using Numeric instead") - print( "Everything should still work, but consider upgrading to numpy") - except: - print( "Neither numpy nor Numeric found.") - print( "Please install numpy (preferred) or Numeric.") - -#----------------------------------------------------- -#Import graphics utilities -# -#First try to import ClimateGraphics, which contains the -#implementation of the plotting commands. The default version -#distributed with the courseware uses PyNgl as a driver, but -#ClimateGraphics can be fairly easily localized to use other -#graphics drivers (e.g. MatPlotLib). If the import fails, a -#dummy graphics stub module is imported, which allows the courseware -#to be run without the user needing to explicitly comment out -#the plot calls in the Chapter Scripts. -#------------------------------------------------------ - -try: - from .ClimateGraphicsMPL import * #Try importing MatPlotLib -# from ClimateGraphics import * #Try importing Ngl driver -except: - # If Ngl not found, try importing the graphics interface - # that uses MatPlotLib. If you have both installed and - #prefer MatPlotLib you can change the order of imports here - #or just do "from ClimateGraphicsMPL import *" after you - #import ClimateUtilities . - try: - from .ClimateGraphicsMPL import * #Try importing MatPlotLib - except: - print( " ") - print( "Graphics not implemented.") - print( "Plot routines will not produce graphs.") - print( "Instead, you can save results from a Curve") - print( "object c into a text file using c.dump()") - print( "and then plot the data using the graphics program") - print( "of your choice.") - from .DummyGraphics import * - -#Section 1: -----Data handling utilities------------------------------- - -# ToDo: Put documentation on use of Curve object here! -# -# Add keyword arguments for axes, etc, -# -# How to handle missing data codes? Note that -# since lists are converted to Numeric arrays, -# text codes can't be used. -# Provide a method to set which data is X. -# -# Add output option to put out a LaTeX formatted table -# -# Provide an easy way to add data column long-names -# -class Curve: - def __init__(self): - self.Xid = None # id of data to be considered as X - self.data = {} # Dictionary of data columns - self.label = {} #Data column label dictionary - self.scatter = {} #Marker for making a curve a scatter plot - #i.e. suppress line drawing and plot symbol - self.idList = [] #Keeps track of original order of data - self.NumCurves = 0 - self.description = None # A string providing general information - self.PlotTitle = '' # Title of the plot - self.switchXY = 0 # Switches X and Y axes for plotting - self.reverseX = 0 # Reverses X axis for plotting - self.reverseY = 0 # Reverses Y axis for plotting - self.XlogAxis = 0 # Use logarithmic axis for X - self.YlogAxis = 0 # Use logarithmic axis for Y - self.Xlabel = '' #X axis label - self.Ylabel = '' #Y axis label - # - # Colors and line titles - # - # - #Install a new curve in the data set, optionally with a variable name and label - #Any 1D indexable object can be installed here:a Numeric array, a Masked Array, - #a Masked Variable (as in cdms) or an ordinary list. - def addCurve(self,data,id = '',label = ''): - self.NumCurves += 1 - if len(id) == 0: - id = 'v%d'%(self.NumCurves-1) - #Transform data from list to a Numeric array here - if type(data) == type([]): - data = Numeric.array(data) - self.data[id] = data - if self.Xid == None: - self.Xid = id #Sets the default id for the X variable - self.idList.append(id) #Keep track of order in which columns added - self.label[id] = label - self.scatter[id] = False - #ToDo: Add checking for consistent lengths, type, etc. of what's being added - def listVariables(self): - return self.idList - def __getitem__(self,id): - return self.data[id] - def __setitem__(self,id,data): - try: - n = len(data[:]) - except: - print( "Object on RHS is not indexable") - return None - #Transform data from list to a Numeric array here - if type(data) == type([]): - data = Numeric.array(data) - if id in self.data.keys(): - self.data[id] = data - else: - self.addCurve(data,id) - - #Method to return Numeric abcissa array for plotting. - def X(self): - #Use of cross section lets us get data from any indexed object - #However, since Masked variables and Masked Arrays yield - #their same types as cross sections, we have to check - #explicitly for a _data component. - # - #This method is used mainly for generating arrays used in plotting, - #and should be streamlined at some point. - temp = self.data[self.Xid][:] - if hasattr(temp,'_data'): - temp = temp._data[:] - return Numeric.array(temp,Numeric.Float) - - #Method to return Numeric ordinate array for plotting - def Y(self): - #Use of cross section lets us get data from any indexed object - outArray = [] - for id in self.idList: - if not (id == self.Xid): - column = self.data[id] - if hasattr(column,'_data'): - outArray.append(column._data[:]) #Deals with masked arrays and variables - else: - outArray.append(column[:]) - return Numeric.array(outArray,Numeric.Float) - - #Dumps curve to a tab-delimited ascii file with column header - def dump(self,fileName = 'out.txt'): - outfile = open(fileName,'w') - # Write out the data description if it is available. - if not (self.description == None): - if not self.description[-1] == '\n': - self.description += '\n' #Put in a newline if needed - outfile.write(self.description) - header = "" - fmt = "" - ids = self.idList - for id in ids: - header += id+'\t' - fmt += '%e\t' - header = header[0:-1]+'\n' # Replaces last tab with a newline. - fmt = fmt[0:-1] + '\n' - outfile.write(header) - for i in range(len(self.data[ids[0]])): - out = tuple([(self.data[id])[i] for id in ids]) - outfile.write(fmt%out) - outfile.close() - - #Extracts a subset of the data and returns it as a new Curve. - #This is useful if you only want to plot some of the columns. - #The input argument dataList is a list of column names - def extract(self,dataList): - c = Curve() - for dataName in dataList: - c.addCurve(self[dataName],dataName) - return c - - - - - -#from string import atof -# deprecated in Python 2.7, gone in Python 3... just use float() - -# Scans a list of lines, locates data lines -# and size, splits of column headers and -# splits off general information text. -# -#A header can optionally be specified as input. -#The delimiter need not be specified if the file -#uses any whitespace character (including tabs) as -#column delimiters. Commas are not whitespace characters, -#and so need to be specified if they are used. -# -#ToDo: -# *Implement missing data coding (with default '-'). -# self.setMissingDataCode(code) (char or numeric) -# possibly translate, check and force consistency -# *Replace optional positional arguments with keyword arguments -def scan(buff,inHeader=None,delimiter = None): - if inHeader == None: - inHeader = [] - #First delete blank lines - buff = clean(buff) - # - #Now look for patterns that indicate data lines - # - startDataLine,endDataLine = findData(buff) - # Found number of items. Now read in the data - # - #Read in the first line. Is it a header? - header = [] - if delimiter == None: - line = buff[startDataLine].split() - else: - line = buff[startDataLine].split(delimiter) - # - try: - #atof(line[0]) - float(line[0]) - except: - header = line - if len(header) == 0: - header = ['V%d'%i for i in range(len(line))] - istart = 0 - else: - istart = 1 - # Replace with inHeader if inHeader has been specified - # (Only use the input header if its length is consistent) - if len(inHeader) == len(line): - header = inHeader - # - # Read in the rest of the lines - # - varlist = [[] for i in range(len(header))] - for line in buff[(startDataLine+istart):endDataLine]: - if delimiter == None: - items = line.split() - else: - items = line.split(delimiter) - try: - for i in range(len(varlist)): - #varlist[i].append(atof(items[i])) - varlist[i].append(float(items[i])) - except: - print( items) - vardict = {} - for name in header: - vardict[name] = varlist[header.index(name)] - return vardict,header #header is returned so we can keep cols in orig order - -#Eliminates blank lines -def clean(buff): - buff = [line.strip() for line in buff] - while(1): - try: - buff.remove('') - except: - break - return buff - -def findData(buff): - runStarts = [] - for i in range(len(buff)-1): - dn = abs(len(buff[i].split())-len(buff[i+1].split())) - if not dn==0: - runStarts.append(i+1) - runStarts.append(len(buff)) - #Find index of run with max length - nmax = -1 - for i in range(len(runStarts)-1): - n = runStarts[i+1]-runStarts[i] - if n > nmax: - nmax = n - imax = runStarts[i] - #Deal with case where entire file is one run - if nmax == -1: - nmax = len(buff) - imax = 0 - return imax,imax+nmax - - -# -# Function to read space or tab-delimited file into a curve object -# The input header is a list of names of the variables in each -# columns, which can be input optionally, mainly to deal with -# the case in which this information is not in the file being read -import string -def readTable(filename,inHeader = None,delimiter = None): - f = open(filename) - buff = f.readlines() - data,header = scan(buff,inHeader,delimiter) - c = Curve() - for key in header: - c.addCurve(data[key],key) - return c - - - - -#============================================== -#---Section 2: Math utilities------------------------------------------- -#============================================== - -#A dummy class useful for passing parameters. -#Just make an instance, then add new members -#at will. -class Dummy: - pass - -#---Polynomial interpolation and extrapolation (adapted -#from Numerical Recipes. -# -#It's used in Romberg extrapolation, but could be useful -#for polynomial OLR fits and so forth as well. Also -#needs online documentation -def polint(xa,ya,x): - n = len(xa) - if not (len(xa) == len(ya)): - print( "Input x and y arrays must be same length") - return "Error" - #Set up auxiliary arrays - c = Numeric.zeros(n,Numeric.Float) - d = Numeric.zeros(n,Numeric.Float) - c[:] = ya[:] - d[:] = ya[:] - #Find closest table entry - ns = 0 - diff = abs(xa[0]-x) - for i in range(n): - difft = abs(xa[i]-x) - if difft < diff: - diff = difft - ns = i - y=ya[ns] - for m in range(1,n): - for i in range(n-m): - ho=xa[i]-x - hp=xa[i+m]-x - w=c[i+1]-d[i] - c[i] = ho*w/(ho-hp) - d[i] = hp*w/(ho-hp) - if 2*ns < (n-m): - dy = c[ns] - else: - ns -= 1 - dy = d[ns] - y += dy - #You can also return dy as an error estimate. Here - #to keep things simple, we just return y. - return y - -#--------------------------------------------------------- -# Class for doing polynomial interpolation -# from a table, using polint -#--------------------------------------------------------- -# -#Usage: -# Let xa be a list of independent variable -# values and ya be a list of the corresponding -# dependent variable values. Then, to create a function -# f (actually a callable object, techically) that interpolates -# or extrapolates to any value x, create f using -# f = interp(xa,ya) -# Then you can get the value you want by invoking f(x) -# for your desired x. -# -# By default, the interpolator does fourth-order interpolation -# using the four nearest neighbors. You can change this by -# using an optional third argument to the creator. For -# example -# -# f = interp(xa,ya,8) -# will use the 8 nearest neighbors (if they are available) -#------------------------------------------------------------ -class interp: - def __init__(self,xa,ya,n=4): - self.xa = Numeric.array(xa) - self.ya = Numeric.array(ya) - self.n = n - def __call__(self,x): - #Find the closes index to x - if self.xa[0] < self.xa[-1]: - i = Numeric.searchsorted(self.xa,x) - else: - i = Numeric.searchsorted(-self.xa,-x) - i1 = max(i-self.n,0) - i2 = min(i+self.n,len(self.xa)) - return polint(self.xa[i1:i2],self.ya[i1:i2],x) - - -#----Quadrature (definite integral) by Romberg extrapolation. -#**ToDo: Add documentation and help string - -#Before developing a general quadrature class, we'll -#implement a class which efficiently carries out trapezoidal rule -#integration with iterative refinement -class BetterTrap: - def __init__(self,f,params,interval,nstart): - self.f = f - self.n = nstart - self.interval = interval - self.params = params - self.integral = self.dumbTrap(nstart) - def dumbTrap(self,n): - a = self.interval[0] - b = self.interval[1] - dx = (b-a)/n - sum = dx*(self.f(a,self.params)+self.f(b,self.params))/2. - for i in range(1,n): - x = a+i*dx - sum = sum + self.f(x,self.params)*dx - return sum - def refine(self): - #Compute the sum of f(x) at the - #midpoints between the existing intervals. - #To get the refinement of the trapezoidal - #rule sum we just add this to half the - #previous result - sum = 0. - a = self.interval[0] - b = self.interval[1] - dx = (b-a)/self.n - #Remember: n is the number of subintervals, - #not the number of endpoints. Therefore we - #have one midpoint per subinterval. Keeping that - #in mind helps us get the range of i right in - #the following loop - for i in range(self.n): - sum = sum + self.f(a+(i+.5)*dx,self.params)*(dx/2.) - #The old trapezoidal sum was multiplied by - #the old dx. To get its correct contribution - #to the refined sum, we must multiply it by .5, - #because the new dx is half the old dx - self.integral = .5*self.integral + sum - # - #Update the number of intervals - self.n = 2*self.n - - -#Here I define a class called -#romberg, which assists in carrying out evaluation of -#integrals using romberg extrapolation. It assumes polint has -#been imported - -class romberg: - def __init__(self,f,nstart=4): - self.nstart = nstart - self.trap = None - # - #------------------------------------------------- - #This snippit of code allows the user to leave the - #parameter argument out of the definition of f - #if it isn't needed - # - self.fin = f - #Find the number of arguments of f and append a - #parameter argument if there isn't any. - nargs = f.func_code.co_argcount - if nargs == 2: - self.f = f - elif nargs ==1: - def f1(x,param): - return self.fin(x) - self.f = f1 - else: - name = f.func_name - print( 'Error: %s has wrong number of arguments'%name) - #----------------------------------------------------- - # - #We keep lists of all our results, for doing - #Romberg extrapolation. These are re-initialized - #after each call - self.nList = [] - self.integralList = [] - def refine(self): - self.trap.refine() - self.integralList.append(self.trap.integral) - self.nList.append(self.trap.n) - dx = [1./(n*n) for n in self.nList] - return polint(dx,self.integralList,0.) - # - #Use a __call__ method to return the result. The - #__call__ method takes the interval of integration - #as its mandatory first argument,takes an optional - #parameter argument as its second argument, and - #an optional keyword argument specifying the accuracy - #desired. - #**ToDo: Introduce trick to allow parameter argument of - #integrand to be optional, as in Integrator. Also, make - #tolerance into a keyword argument - # - def __call__(self,interval,params=None,tolerance=1.e-6): - self.nList = [] - self.integralList = [] - #Make a trapezoidal rule integrator - self.trap = BetterTrap(self.f,params,interval,self.nstart) - self.nList.append(self.nstart) - self.integralList.append(self.trap.integral) - # - #Refine initial evaluation until - oldval = self.refine() - newval = self.refine() - while abs(oldval-newval)>tolerance: - oldval,newval = newval,self.refine() - return newval - - - -#-------Runge-Kutta ODE integrator -#**ToDo: -# * Implement a reset() method which resets to initial conditions. -# Useful for doing problem over multiple times with different -# parameters. -# -# * Referring to the independent variable as 'x' is awful, and -# confusing in many contexts. Introduce a variable name -# dictonary with default names like 'independent' and 'dependent' -# (and short synonyms) so if fi is the integrator object -# fi['indep'] is the current value of the independent variable -# and fi['dep'] is the current (vector) value of the dependent -# variable. Then allow user to rename or synonym these to -# the actual user-supplied names. This is an alternative -# to using the list returned by fi.next(). Then expunge -# all references to things like fi.x from the examples and -# chapter scripts. They are too confusing. Similarly, -# change the name of the increment from dx to something else -# -# * Similarly, we could introduce a dictionary of some sort -# to make it easier to set up multidimensional systems and -# refer to the different vector components by name -# (e.g. refer to v[0] at T, v[1] as dTdy , etc. ) -# -# * Make the integrator object callable. The call can return -# a list of results for all the intermediate steps, or optionally -# just the final value. -class integrator: - ''' - Runge-Kutta integrator, for 1D or multidimensional problems - - Usage: - - First you define a function that returns the - derivative(s), given the independent and dependent - variables as arguments. The independent variable (think - of time) is always a scalar. The dependent variable (think - of the x and y coordinates of a planet) can be either a scalar - or a 1D array, according to the problem. For the - multidimensional case, this integrator will work with any - kind of array that behaves like a Numeric array, i.e. supports - arithmetic operations. It will not work with plain Python lists. - The derivative function should return an array of derivatives, in - the multidimensional case. The derivative function can have any name. - - The derivative function can optionally have a third argument, to provide - for passing parameters (e.g. the radius of the Earth) to the - function. The "parameter" argument, if present, can be any Python - entity whatsoever. If you need to pass multiple constants, or - even tables or functions, to your derivative function, you can - stuff them all into a list or a Python object. - - - Example: - In the 1D case, to solve the equation - dz/dt = -a*t*t/(1.+z*z) - in which z is the dependent variable and t is the - independent variable, your derivative function would be - def g(t,z): - return -a*t*t/(1.+z*z) - - treating the parameter a as a global, or perhaps - def g(t,z,params): - return -params.a*t*t/(params.b+z*z) - - while in a 2D case, your function might look like: - def g(t,z): - return Numeric.array([-z[1],z[0]]) - - or perhaps something like: - def g(t,z): - return t*Numeric.sin(z) - - or even - def g(t,z,params): - return Numeric.matrixmultiply(params(t),z) - - where params(t) in this case is a function returning - a Numeric square matrix of the right dimension to multiply z. - - BIG WARNING: Note that all the examples which return a - Numeric array return a NEW INSTANCE (i.e. copy) of an - array. If you try to set up a global array and re-use - it to return your results from g, you will really be - just returning a REFERENCE to the same array each time, - and each call will change the value of all the previous - results. This will mess up the computation of intermediate - results in the Runge-Kutta step. An example of the sort of thing - that will NOT work is: - zprime = Numeric.zeros(2,Numeric.Float) - def g(t,z,params): - zprime[0] = z[1] - zprime[1] = -z[0] - return zprime - Try it out. This defines the harmonic oscillator, and a plot - of the orbit should give a circle. However, it doesn't The problem - reference/value distinction. The right way to define the function - would be - def g(t,z): - return Numeric.array([z[1],-z[0]]) - Try this one. It should work properly now. Note that any arithmetic - performed on Numeric array objects returns a new instance of an Array - object. Hence, a function definition like - def g(t,z): - return t*z*z+1. - will work fine. - - Once you have defined the derivitave function, - you then proceed as follows. - - First c reate an integrator instance: - int_g = integrator(g,0.,start,.01) - - where "0." in the argument list is the initial value - for the independent variable, "start" is the initial - value for the dependent variable, and ".01" is the - step size. You then use the integrator as follows: - - int_g.setParams(myParams) - while int_g.x < 500: - print( int_g.next()) - - The call to setParams is optional. Just use it if your - function makes use of a parameter object. The next() method - accepts the integration increment (e.g. dx) as an optional - argument. This is in case you want to change the step size, - which you can do at any time. The integrator continues - using the most recent step size it knows. - - Each call to int_g.next returns a list, the first of whose - elements is the new value of the independent variable, and - the second of whose elements is a scalar or array giving - the value of the dependent variable(s) at the incremented - independent variable. - - ''' - def __init__(self, derivs,xstart,ystart,dx=None): - self.derivsin = derivs - # - #The following block checks to see if the derivs - #function has a parameter argument specified, and - #writes a new function with a dummy parameter argument - #appended if necessary. This allows the user to leave - #out the parameter argument from the function definition, - #if it isn't needed. - nargs = derivs.func_code.co_argcount - if nargs == 3: - self.derivs = derivs - elif nargs == 2: - def derivs1(x,y,param): - return self.derivsin(x,y) - self.derivs = derivs1 - else: - name = derivs.func_name - print('Error: %s has wrong number of arguments'%name) - # - # - self.x = xstart - #The next statement is a cheap trick to initialize - #y with a copy of ystart, which works whether y is - #a regular scalar or a Numeric array. - self.y = 0.+ ystart - self.dx = dx #Can instead be set with the first call to next() - self.params = None - #Sets the parameters for the integrator (optional). - #The argument can be any Python entity at all. It is - #up to the user to make sure the derivative function can - #make use of it. - def setParams(self,params): - self.params = params - #Computes next step. Optionally, takes the increment - #in the independent variable as an argument. The - #increment can be changed at any time, and the most - #recently used value is remembered, as a default - def next(self,dx = None): - if not (dx == None): - self.dx = dx - h = self.dx - hh=h*0.5; - h6=h/6.0; - xh=self.x+hh; - dydx = self.derivs(self.x,self.y,self.params) - yt = self.y+hh*dydx - dyt = self.derivs(xh,yt,self.params) - yt =self.y+hh*dyt - dym = self.derivs(xh,yt,self.params) - yt =self.y+h*dym - dym += dyt - dyt = self.derivs(self.x+h,yt,self.params) - self.y += h6*(dydx+dyt+2.0*dym) - self.x += h - return self.x,self.y - - - -#**ToDo: -# -# Store the previous solution for use as the next guess(?) -# -# Handle arithmetic exceptions in the iteration loop -# -class newtSolve: - ''' - Newton method solver for function of 1 variable - A class implementing Newton's method for solving f(x) = 0. - - Usage: solver = newtSolve(f), where f is a function with - calling sequence f(x,params). Values of x such that - f(x,params) = 0 are - then found by invoking solver(guess), where guess - is the initial guess. The solver returns the string - 'No Convergence' if convergence fails. The argument - params allows parameters to be passed to the function. - It can be left out of the function definition if you don't - need it. Note that params can be any Python object at all - (including,e.g.,lists, functions or more complex user-defined objects) - - Optionally, one can specify the derivative function - in the creator,e.g. solver = newtSolve(f,fp). - If the derivative function isn't specified, the solver - computes the derivative approximately using a centered - difference. Note that in either case you can access - the derivative function by invoking solver.deriv(x) - As for f, fp can be optionally defined with a parameter - argument if you need it. The same parameter object is - passed to f and fp. - - Use solver.setParams(value) to set the parameter object - Alternately, the parameter argument can be passed as - an optional second argument in the solver call. (see - example below). - - Adjustable constants: - eps Increment for computing numerical approximation to - the derivative of f - tolerance Accuracy criterion for ending the iteration - (an approximation to the error in the root) - nmax maximum number of iterations - - e.g. to change the maximum number of iterations for an instance - of the class, set solver.nmax = 10 . - ----------------Usage Examples----------------------------- - - Example 1: Function without parameters: - def g(x): - return x*x - 1. - roots = newtSolve(g) - roots(2.) - - Example 2, Function with parameters: - def g(x,constants): - return constants.a*x*x - constants.b - roots = newtSolve(g) - constants = Dummy() - constants.a = 1. - constants.b = 2. - roots.setParam(constants) - roots(2.) - roots(1.) - - Example 2a: - Instead of using roots.setParam(...) we could do - roots(2.,constants) - roots(1.) the parameters are remembered - constants.a = 3. - roots(1.,constants) We changed the constants - - Example 3, using scan to find initial guesses: - def g(x): - return x*x - 1. - roots = newtSolve(g) - guesses = roots.scan([-2.,2.],100) - for guess in guesses: - print( roots(guess)) - ''' - def __init__(self,f,fprime = None): - self.fin = f - #Find the number of arguments of f and append a - #parameter argument if there isn't any. - nargs = f.func_code.co_argcount - if nargs == 2: - self.f = f - elif nargs ==1: - def f1(x,param): - return self.fin(x) - self.f = f1 - else: - name = f.func_name - print( 'Error: %s has wrong number of arguments'%name) - self.eps = 1.e-6 - def deriv(x,params): - return (self.f(x+self.eps,params)- self.f(x-self.eps,params))/(2.*self.eps) - if fprime == None: - self.deriv = deriv - else: - #A derivative function was explicitly specified - #Check if it has a parameter argument - nargs = fprime.func_code.co_argcount - if nargs == 2: - self.deriv = fprime #Has a parameter argument - elif nargs == 1: - self.fprimein = fprime - def fprime1(x,param): - return self.fprimein(x) - self.deriv = fprime1 - else: - name = fprime.func_name - print( 'Error: %s has wrong number of arguments'%name) - self.tolerance = 1.e-6 - self.nmax = 100 - self.params = None - def __call__(self,xGuess,params = None): - if not (params == None): - self.setParams(params) - x = xGuess - for i in range(self.nmax): - dx = (self.f(x,self.params)/self.deriv(x,self.params)) - x = x - dx - if abs(dx) < self.tolerance: - return x - return 'No Convergence' - def setParams(self,params): - #**ToDo: Check if f1 has a parameter argument - #defined, and complain if it doesn't - self.params = params - def scan(self,interval,n=10): - #Finds initial guesses to roots in a specified - #interval, subdivided into n subintervals. - #e.g. if the instance is called "solver" - #solver.scan([0.,1.],100) generates a list - #of guesses between 0. and 1., with a resolution - #of .01. The larger n is, the less is the chance that - #a root will be missed, but the longer the search - #will take. If n isn't specified, the default value is 10 - # - #ToDo: Replace this with a bisection search, allowing user - #to specify the maximum number of distinct guesses that - #need to be found. - guessList = [] - dx = (interval[1]-interval[0])/(n-1) - flast = self.f(interval[0],self.params) - for x in [interval[0]+ i*dx for i in range(1,n)]: - fnow = self.f(x,self.params) - if ((fnow >= 0.)&(flast <=0.)) or ((fnow <= 0.)&(flast >=0.)): - guessList.append(x) - flast = fnow - return guessList diff --git a/pyrads/DummyGraphics.py b/pyrads/DummyGraphics.py deleted file mode 100644 index 0486aec..0000000 --- a/pyrads/DummyGraphics.py +++ /dev/null @@ -1,177 +0,0 @@ -from __future__ import division, print_function, absolute_import -#----------Section 3: Plotting utilities------------------------------- -#This is a dummy graphics routine, to import if a graphics driver -#is not found. It is the fallback import if the import of ClimateGraphics -#fails in ClimateUtilities. This dummy routine allows courseware -#scripts to be run without the user needing to comment out plot commands. -#It also can be used as a template for customizing the plotter interface -#to work with you locally favored Python graphics driver (e.g. MatPlotLib). -#----------------------------------------------------------------------= - -#A dummy class useful for passing parameters. -#Just make an instance, then add new members -#at will. Not sure why this also has to be defined -#here, since it's defined in ClimateUtilities, but -#it seems to be necessary -class Dummy: - pass - -# A little class to make resources for Ngl plot options -class resource: - def __init__(self): - self.trYReverse = False - self.trXReverse = False - self.trXLog = False - self.trYLog = False - self.nglFrame = False - self.nglDraw = False - #ToDo: Missing data code resource, line styles, line colors - -# A little class for use as a return object from plot(), so -# the user has an easy way to delete a window or save a plot -# to an eps file. -class plotObj: - def __init__(self,workstation,plot,WorkstationResources = None): - self.workstation = workstation - self.plot = plot - self.WorkstationResources = WorkstationResources - #Deletes a plot window - def delete(self): - #Deletes a plot window, cleans up - pass - def save(self,filename = 'plot'): - #Saves a plot to a file - pass - -#Makes a line pot from a Curve object -def plot(c): - print( "Plotting not implemented") - #Set axis options according to information - #in the curve object c. - # - #c.reverseY: - # if True, reverse the Y axis - #c.reverseX: - # if True, reverse the X axis - #c.XlogAxis: - # if True, use a logarithmic X axis - #c.YlogAxis: - # if True, use a logarithmic Y axis - # - #Customize Line styles and colors here - # - #Set thePlot title - #c.PlotTitle: - # String containing the plot title - #Axis labels - #X and Y axis labels - #c.Xlabel: - # String containing the X axis label - #c.Ylabel: - # String containing the Y axis label - # - #Interchange the X and Y axes - if c.switchXY: - pass - #If True, exchang the axes - - # Legends, for multicurve plot - legends = [] - for id in c.listVariables(): - if not id == c.Xid: - if len(c.label[id]) > 0: - legends.append(c.label[id]) - else: - legends.append(id) - - # - #Suppress line drawing and just plot symbol for scatter plot curves - #Customize plotting symbols and marker sizes if desired - for id in c.listVariables(): - if not id == c.Xid: - if c.scatter[id]: - #If True, treat this data column as a scatter - #plot and don't draw lines - pass - else: - #If False, draw lines (default) - pass - # - #Initialize the plot window here, if necessary. - #w is the handle to the plot window - w = None - if c.switchXY: - #Put in the command for doing the line plot here - #Do the plot with the X and Y axes in the usual - #order - pass - else: - #Put in the command for doing the line plot here, - #but switch the order of the axes. - pass - # - #Now draw the plot - #Depending on your graphics software, the preceding - #command may already have drawn the plot. In some graphics - #packages, after the plot is created, another command needs - #to be executed in orter to display it. Either way, the - #variable plotHandle below is a handle referring to the plot. - #It is used to build a plotObj that can be used for further - #manipulation of the plot such as saving or deleting. In some - #graphics packages, which give control over such things from - #menus in the plot window, the use of a plotObj for this may - #be unnecessary. - plotHandle = None - return plotObj(w,plotHandle) #So that user can delete window or save plot - -# A basic contour plotter, which will plot a contour plot -# of a numpy array. The x and y scales can optionally -# be specified using keyword arguments x and y. For example, -# if we want the x scale to be the array (or list) lat, and -# the y scale to be the array (or list) lon, we would call -# contour as contour(A,x=lat,y=lon). -def contour(A,**kwargs): - print( "Plotting not implemented") - #The following allows an expert user to pass - # options directly to the plotter. - if 'resource' in kwargs.keys(): - r = kwargs['resource'] - else: - r = Dummy() - # - r.cnFillOn = True #Use color fill - - if 'x' in kwargs.keys(): - #Set the X array for the contour plot - XArray = kwargs['x'] - if 'y' in kwargs.keys(): - #Set the Y array for the contour plot - YArray = kwargs['y'] - # - # Now create the plot - rw = Dummy() - #Set the color map - if 'colors' in kwargs.keys(): - if (kwargs['colors'] == 'gray') or (kwargs['colors'] == 'grey') : - #Set the default greyscale - #(Substitute the appropriate command for your driver) - rw.wkColorMap = 'gsdtol' - else: - rw.wkColorMap = kwargs['colors'] - else: - #Default rainbow color table - rw.wkColorMap = "temp1" - - #Open/initialize a plot window - w = None - #Make the plot. plotHandle is the handle returned - #by the plotter, used for further manipulation of the - #plot. (Redundant for some kinds of plotting packages) - plotHandle = None - #Now draw the plot, if your driver needs this as a separate step - #(Insert command for drawing the plot here, e.g. - #ShowPlot(plotHandle). - # - #Return a plotObj with the necessary data for further - #manipulation. - return plotObj(w,plotHandle,rw) #So user can delete or save plot diff --git a/pyrads/Get_Fluxes.py b/pyrads/Get_Fluxes.py index 66bfa36..575ecc4 100644 --- a/pyrads/Get_Fluxes.py +++ b/pyrads/Get_Fluxes.py @@ -1,6 +1,5 @@ from __future__ import division, print_function, absolute_import import numpy as np -from . import phys from .Planck import Planck_n diff --git a/pyrads/OpticalThickness.py b/pyrads/OpticalThickness.py index 9995628..7853399 100644 --- a/pyrads/OpticalThickness.py +++ b/pyrads/OpticalThickness.py @@ -6,7 +6,7 @@ ''' from __future__ import division, print_function, absolute_import import numpy as np -from .Absorption_Crosssections_HITRAN2016 import getKappa_HITRAN +from .Absorption_Crosssections_HITRAN2016 import getKappa_HITRAN, loadSpectralLines from . import Absorption_Continuum_MTCKD from .Absorption_Continuum_MTCKD import get_H2OContinuum from scipy.integrate import cumtrapz @@ -19,21 +19,21 @@ def compute_tau_H2ON2(p,T,q,grid,params,RH=1.,use_numba=False): kappa = np.zeros( (grid.Np,grid.Nn) ) + + hitran_data_H2O = loadSpectralLines(molName="H2O",minWave=grid.n0,maxWave=grid.n1) + for pres,temp,q_H2O in zip(p,T,q): p_H2O = RH * params.esat(temp) # ... print( "compute kappa at p,T = ",pres,temp) if use_numba: - - from .Absorption_Crosssections_HITRAN2016_numba import getKappa_HITRAN_numba - kappaH2O = getKappa_HITRAN_numba(grid.n,grid.n0,grid.n1,grid.dn, \ - "H2O",press=pres,press_self=p_H2O, \ - temp=temp,broadening="mixed", lineWid=25., \ - cutoff_option="fixed",remove_plinth=True) + print("Numba functionality is not supported in this release! Coming soon...") + continue + else: kappaH2O = getKappa_HITRAN(grid.n,grid.n0,grid.n1,grid.dn, \ - "H2O",press=pres,press_self=p_H2O, \ + hitran_data=hitran_data_H2O,press=pres,press_self=p_H2O, \ temp=temp,broadening="mixed", lineWid=25., \ cutoff_option="fixed",remove_plinth=True) @@ -61,6 +61,10 @@ def compute_tau_H2ON2_CO2dilute(p,T,q,ppv_CO2,grid,params,RH=1.,use_numba=False) kappa = np.zeros( (grid.Np,grid.Nn) ) kappa_h2o = np.zeros( (grid.Np,grid.Nn) ) kappa_co2 = np.zeros( (grid.Np,grid.Nn) ) + + hitran_data_H2O = loadSpectralLines(molName="H2O",minWave=grid.n0,maxWave=grid.n1) + hitran_data_CO2 = loadSpectralLines(molName="CO2",minWave=grid.n0,maxWave=grid.n1) + for pres,temp,q_H2O in zip(p,T,q): p_H2O = RH * params.esat(temp) # ... R_mean = q_H2O*params.Rv + (1.-q_H2O)*params.R @@ -68,29 +72,21 @@ def compute_tau_H2ON2_CO2dilute(p,T,q,ppv_CO2,grid,params,RH=1.,use_numba=False) print( "compute kappa at p,T = ",pres,temp) + if use_numba: - from .Absorption_Crosssections_HITRAN2016_numba import getKappa_HITRAN_numba - - kappaH2O = getKappa_HITRAN_numba(grid.n,grid.n0,grid.n1,grid.dn, \ - "H2O",press=pres,press_self=p_H2O, \ - temp=temp,broadening="mixed", lineWid=25., \ - cutoff_option="fixed",remove_plinth=True) - - kappaCO2 = getKappa_HITRAN_numba(grid.n,grid.n0,grid.n1,grid.dn, \ - "CO2",press=pres,press_self=0., \ - temp=temp,broadening="air", lineWid=25., \ - cutoff_option="fixed",remove_plinth=False) # use of "air" broadening consistent with trace gas assumption - + print("Numba functionality is not supported in this release! Coming soon...") + continue + else: kappaH2O = getKappa_HITRAN(grid.n,grid.n0,grid.n1,grid.dn, \ - "H2O",press=pres,press_self=p_H2O, \ + hitran_data=hitran_data_H2O,press=pres,press_self=p_H2O, \ temp=temp,broadening="mixed", lineWid=25., \ cutoff_option="fixed",remove_plinth=True) - + kappaCO2 = getKappa_HITRAN(grid.n,grid.n0,grid.n1,grid.dn, \ - "CO2",press=pres,press_self=0., \ + hitran_data=hitran_data_CO2,press=pres,press_self=0, \ temp=temp,broadening="air", lineWid=25., \ - cutoff_option="fixed",remove_plinth=False) # use of "air" broadening consistent with trace gas assumption + cutoff_option="fixed",remove_plinth=True) # use of "air" broadening consistent with trace gas assumption # add continuum: # here I'm only using kappa from mtckd crosssection file, @@ -116,25 +112,25 @@ def compute_tau_H2ON2_CO2dilute(p,T,q,ppv_CO2,grid,params,RH=1.,use_numba=False) def compute_tau_dryCO2(p,T,q,ppv_CO2,grid,params,use_numba=False): kappa = np.zeros( (grid.Np,grid.Nn) ) + + hitran_data_CO2 = loadSpectralLines(molName="CO2",minWave=grid.n0,maxWave=grid.n1) + for pres,temp,q_H2O in zip(p,T,q): p_CO2 = pres * ppv_CO2 R_mean = params.R q_CO2 = convert_molar_to_mass_ratio(ppv_CO2,params.R_CO2,R_mean) print( "compute kappa at p,T = ",pres,temp) + if use_numba: - from .Absorption_Crosssections_HITRAN2016_numba import getKappa_HITRAN_numba - - kappaCO2 = getKappa_HITRAN_numba(grid.n,grid.n0,grid.n1,grid.dn, \ - "CO2",press=pres,press_self=p_CO2, \ - temp=temp,broadening="mixed", lineWid=25., \ - cutoff_option="fixed",remove_plinth=False) # don't take out plinth! + print("Numba functionality is not supported in this release! Coming soon...") + continue else: kappaCO2 = getKappa_HITRAN(grid.n,grid.n0,grid.n1,grid.dn, \ - "CO2",press=pres,press_self=p_CO2, \ + hitran_data=hitran_data_CO2,press=pres,press_self=p_CO2, \ temp=temp,broadening="mixed", lineWid=25., \ - cutoff_option="fixed",remove_plinth=False) # don't take out plinth! + cutoff_option="fixed",remove_plinth=False) # don't take out plinth! kappa[ p==pres,: ] = kappaCO2*q_CO2 # save print( "done! \n") diff --git a/pyrads/__init__.py b/pyrads/__init__.py index 959bde7..438076d 100644 --- a/pyrads/__init__.py +++ b/pyrads/__init__.py @@ -3,12 +3,10 @@ from . import Absorption_Continuum from . import Absorption_Continuum_MTCKD from . import Absorption_Crosssections_HITRAN2016 -from . import ClimateUtilities from . import Get_Fluxes from . import OpticalThickness from . import Planck from . import SetupGrids from . import Thermodynamics from . import VerticalStructure -from . import object_helpers from . import phys diff --git a/pyrads/object_helpers.py b/pyrads/object_helpers.py deleted file mode 100644 index 06b5bd6..0000000 --- a/pyrads/object_helpers.py +++ /dev/null @@ -1,21 +0,0 @@ -from __future__ import division, print_function, absolute_import -from copy import copy # to copy lists of objects... - -# ================================================================================== -# Retrieve object(s) from obj list that match obj.attribute == value -# or function(obj.attribute,value). -# Either return object or list. -# Usages: -# sublist = get_objects_from_list(masterlist,"name","Dry_50") -# sublist = get_objects_from_list(masterlist,"omega",5,comparefn=largerthan) -# (note: comparefn(x,y) must return a boolean!) - -def get_objects_from_list(obj_list,attribute,value,comparefn=None): - if comparefn is None: - x = [obj for obj in obj_list if getattr(obj,attribute)==value] - else: - x = [obj for obj in obj_list if comparefn(getattr(obj,attribute),value)] - - if len(x)==1: - x = x[0] - return x diff --git a/pyrads/phys.py b/pyrads/phys.py index c05b959..d4710ea 100644 --- a/pyrads/phys.py +++ b/pyrads/phys.py @@ -6,7 +6,7 @@ ### https://geosci.uchicago.edu/~rtp1/PrinciplesPlanetaryClimate/Courseware/coursewarePortal.html import math -from .ClimateUtilities import * #To get the math methods routines +#from .ClimateUtilities import * #To get the math methods routines # #All units are mks units # @@ -387,271 +387,3 @@ def B(nu,T): # # - -#----------Saturation Vapor Pressure functions--------------- -# - -#Saturation vapor pressure over ice (Smithsonian formula) -# Input: Kelvin. Output: Pascal -def satvpi(T): - # - # Compute es over ice (valid between -153 c and 0 c) - # see smithsonian meteorological tables page 350 - # - # Original source: GFDL climate model, circa 1995 - esbasi = 6107.1 - tbasi = 273.16 - # - aa = -9.09718 *(tbasi/T-1.0) - b = -3.56654 *math.log10(tbasi/T) - c = 0.876793*(1.0-T/tbasi) - e = math.log10(esbasi) - esice = 10.**(aa+b+c+e) - return .1*esice #Convert to Pascals - -#Saturation vapor pressure over liquid water (Smithsonian formula) -# Input: Kelvin. Output: Pascal -def satvpw(T): - # compute es over liquid water between -20c and freezing. - # see smithsonian meteorological tables page 350. - # - # Original source: GFDL climate model, circa 1995 - esbasw = 1013246.0 - tbasw = 373.16 - # - aa = -7.90298*(tbasw/T-1) - b = 5.02808*math.log10(tbasw/T) - c = -1.3816e-07*( 10.**( ((1-T/tbasw)*11.344)-1 ) ) - d = 8.1328e-03*( 10.**( ((tbasw/T-1)*(-3.49149))-1) ) - e = math.log10(esbasw) - esh2O = 10.**(aa+b+c+d+e) - return .1*esh2O #Convert to Pascals - -# An alternate formula for saturation vapor pressure over liquid water -def satvpw_Heymsfield(T): - ts=373.16 - sr=3.0057166 -# Vapor pressure over water. Heymsfield formula - ar = ts/T - br = 7.90298*(ar-1.) - cr = 5.02808*math.log10(ar); - dw = (1.3816E-07)*(10.**(11.344*(1.-1./ar))-1.) - er = 8.1328E-03*((10.**(-(3.49149*(ar-1.))) )-1.) - vp = 10.**(cr-dw+er+sr-br) - vp=vp*1.0e02 - return(vp) - -def satvpg(T): -#This is the saturation vapor pressure computation used in the -#GFDL climate model. It blends over from water saturation to -#ice saturation as the temperature falls below 0C. - if ((T-273.16) < -20.): - return satvpi(T) - if ( ((T-273.16) >= -20.)&((T-273.16)<=0.)): - return 0.05*(273.16-T)*satvpi(T) + 0.05*(T-253.16)*satvpw(T) - if ((T-273.16)>0.): - return satvpw(T) - -#Saturation vapor pressure for any substance, computed using -#the simplified form of Clausius-Clapeyron assuming the perfect -#gas law and constant latent heat -def satvps(T,T0,e0,MolecularWeight,LatentHeat): - Rv=Rstar/MolecularWeight - return e0*math.exp(-(LatentHeat/Rv)*(1./T - 1./T0)) - -#This example shows how to simplify the use of the simplified -#saturation vapor pressure function, by setting up an object -#that stores the thermodynamic data needed, so it doesn't have -#to be re-entered each time. Because of the __call__ method, -#once the object is created, it can be invoked like a regular -#function. -# -#Usage example: -# To set up a function e(T) that approximates the saturation -# vapor presure for a substance which has a latent heat of -# 2.5e6 J/kg, a molecular weight of 18 and has vapor pressure -# 3589. Pa at a temperature of 300K, create the function using: -# -# e = satvps_function(300.,3589.,18.,2.5e6) -# -# and afterward you can invoke it simply as e(T), where T -# is whatever temperature you want to evaluate it for. -# -#Alternately, satvps_function can be called with a gas object -#as the first argument, e.g. -# e = satvps_function(phys.CO2) -# -#If no other arguments are given, the latent heat of sublimation -#will be used when e(T) is called for temperatures below the triple -#point, and the latent heat of vaporization will be used for -#temperatures above the triple point. To allow you to force -#one or the other latent heats to be used, satvps_function takes -#an optional second argument when the first argument is a gas -#object. Thus, -# e = satvps_function(phys.CO2,'ice') -#will always use the latent heat of sublimation, regardless of T, -#while e = satvps_function(phys.CO2,'liquid') will always use -#the latent heat of vaporization. -class satvps_function: - def __init__(self,Gas_or_T0,e0_or_iceFlag=None,MolecularWeight=None,LatentHeat=None): - #Check if the first argument is a gas object. If not, assume - #that the arguments give T0, e0, etc. as numbers - self.iceFlag = e0_or_iceFlag - if isinstance(Gas_or_T0,gas): - self.gas = Gas_or_T0 - self.M = Gas_or_T0.MolecularWeight - self.T0 = Gas_or_T0.TriplePointT - self.e0 = Gas_or_T0.TriplePointP - if self.iceFlag == 'ice': - self.L = Gas_or_T0.L_sublimation - elif self.iceFlag == 'liquid': - self.L = Gas_or_T0.L_vaporization - else: - self.iceFlag = 'switch' - self.M = Gas_or_T0.MolecularWeight - else: - self.L = LatentHeat - self.M = MolecularWeight - self.T0 = Gas_or_T0 - self.e0 = e0_or_iceFlag - def __call__(self,T): - #Decide which latent heat to use - if self.iceFlag == 'switch': - if T ptop: - ans = ad.next() - pa = math.exp(ans[0]) - T = math.exp(ans[1]) - p = pa+self.satvp(T) - pL.append(p) - molarConL.append(self.satvp(T)/p) - TL.append(T) - #Numeric.array turns lists into arrays that one - #can do arithmetic on. - pL = Numeric.array(pL) - TL = Numeric.array(TL) - molarConL = Numeric.array(molarConL) - #Now compute mass specific concentration - Mc = self.condensible.MolecularWeight - Mnc = self.noncon.MolecularWeight - Mbar = molarConL*Mc +(1.-molarConL)*Mnc - qL = (Mc/Mbar)*molarConL - # - #The else clause below interpolates to a - #specified pressure array pgrid, if desired. - # interp is a class defined in ClimateUtilities - #which creates a callable object which acts like - #an interpolation function for the listed data give - #as arguments. - if pgrid == None: - return pL,TL,molarConL,qL - else: - T1 = interp(pL,TL) - mc1 = interp(pL,molarConL) - q1 = interp(pL,qL) - T = Numeric.array([T1(pp) for pp in pgrid]) - mc = Numeric.array([mc1(pp) for pp in pgrid]) - q = Numeric.array([q1(pp) for pp in pgrid]) - return Numeric.array(pgrid),T, mc, q