Gaussian line profile using numeric python
Geoff Low
geoff at ou043085.otago.ac.nz
Thu Mar 29 06:03:58 EST 2001
Hi.
I have an extensive set of xy data that I want to convolute with a gaussian
profile. Below is my algorithm, which doesn't work. Can anyone see where
I'm being dumb?
import Numeric
from MLab import sum
from Numeric import arange,zeros,transpose,shape,exp,power,pi
from TableIO import readTableAsArray,TableFromColumns,writeTable
from math import sqrt
def get_data(filename):
d = readTableAsArray(filename,'#')
x = d[:,0]
y = d[:,1]
del d
return x,y
def main(inpfile,outfile='res.out',sfreq=350,efreq=5050,stepf=2,fwhm=2.5):
fin,iin = get_data(inpfile)
xpt = arange(sfreq,efreq,stepf,'d')
ypt = calcspec(fin,iin,xpt,fwhm)
del fin,iin
put_data(outfile,xpt,ypt)
def put_data(filename,x,y):
tmp = TableFromColumns([x,y])
tmp['filename'] = filename
writeTable(tmp)
del tmp
def calcspec(fin,iin,x,fw):
norm,expterm = 0.0,0.0
inten = zeros((len(x)),'d')
norm = float(1/(sqrt(2.0*pi)*fw))
expterm = float(1/(2*(fw**2)))
for i in range(len(x)):
print "Starting step ",i+1,"of ",len(x)
print "sum fin: ",sum(fin)," sum iin: ",sum(iin)
print "x[",i,"] = ",x[i],"; SUM(x[",i,"] - fin) = ",sum((x[i]-fin[:]),0)
#inten[i] = sum((norm*iin[:]*exp(-1.0*(power((expterm*(xpts[i]-fin[:])),2)))),0)
inten[i] = norm*sum(iin[:]*exp(-1.0*(power(((x[i]-fin[:])*expterm),2))))
#for j in range(len(fin)):
#inten[i] = inten[i] + float(norm*iin[j]*exp(-1.0*((expterm*(xpts[i]-fin[j]))**2)))
return inten
P.S. I'm no math wizard but this is my interpretation of the gaussian lineshape
function
G(v) = 1/(sqrt(2*PI)*FWHM)*exp(-(v-vo)^2/(2*FWHM^2))
Cheers
Geoff
http://www.zfree.co.nz
######
Fortune favours the brave....
but keeps an axe aside for the stupid!
######
More information about the Python-list
mailing list