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