HI Danny, Here is a simplified version:
-------------------------------------------------------------------- import numpy as np import pylab from pylab import * import matplotlib.pyplot as plt import scipy from scipy.integrate import quad from scipy.fftpack import fft, ifft, fftfreq ############################## # Load data from .txt file ############################## data1 = np.loadtxt('TL_pulseC.txt',skiprows = 2 ,usecols = (0,1)) data2 = np.loadtxt('cos_pulseC.txt',skiprows = 2 ,usecols = (0,1)) ################################ # Create data arrays ################################ pos = np.array(data1[:,0]) # micrometers pos_m = pos*1.0*10**(-6) # meters amp = np.array(data1[:,1]) # V amp_auto = np.array(data2[:,1]) # V ############################################################# # Finding the index of the maximum amplitude where the pulse delay is zero. # Then calculating the mean to shift the curve accordingly. ############################################################# peak_id = np.where(amp == np.max(amp))[0][0] mean = pos_m[peak_id] dif = pos_m - mean t_d =(2.*dif)/(2.99*10**8)*10**(12.) # picoseconds (ps) t = np.arange(-0.5,0.5,0.0001) # picoseconds (ps) ################################ # Define constants ################################ c = 2.99*10**8 # m/s alpha = np.pi # rad gamma = 200.0*10**(-15.) lamb_o = 1160.00*10**(-9.) # m omega_o = ((2.*np.pi*c)/lamb_o)*1.*10**(-12.) # ps^-1 delta = np.pi/2. # rad FWHM = 32.53*10**(-6) # m t_p = ((FWHM)/c)*10**(12.) # ps E_norm = 1. # normalized ######################################################### # Define functions ######################################################### ################################ # Input electric field ################################ def E_o(x): return E_norm*(cosh(1.76*x/t_p))**(-1.) def E_(x): return (1./2.)*E_o(x)*np.exp(-1j*omega_o*x) ################################# # Shaped electric field ################################# alpha = np.pi delta = np.pi/2. gamma = 200.*10**(-15) # sec dt = (t[1]-t[0])*(1*10**(-12)) # sec def phi(x): return alpha*np.cos(gamma*(x - omega_o) - delta) def M(x): return np.exp(1j*phi(x)) ################################## # Write E(w) as E(t) using FFT ################################## def E_out(x): E_in_w = fft(E_(x)) omega = fftfreq(len(x),dt)*2*np.pi E_out_w = E_in_w*M(omega) return ifft(E_out_w) def integrand(x,y): return abs(E_out(x))**2.*abs(E_(x - y))**2. def G(y): return quad(integrand, -np.infty, np.infty, args=(y)) auto_corr = [] for tau in t_d: integral,error = G(tau) auto_corr.append(integral) # Plotting plt.plot(t_d,auto_corr) plt.show() ------------------------------------------------------------------------- --------------------------------------------------------------------------- IndexError Traceback (most recent call last) /Users/colinross/Documents/Optics/Programming/five.py in <module>() 107 108 for tau in t_d: --> 109 integral,error = G(tau) 110 auto_corr.append(integral) 111 /Users/colinross/Documents/Optics/Programming/five.py in G(y) 102 103 def G(y): --> 104 return quad(integrand, -np.infty, np.infty, args=(y)) 105 106 auto_corr = [] /Users/colinross/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst) 279 args = (args,) 280 if (weight is None): --> 281 retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points) 282 else: 283 retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts) /Users/colinross/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points) 345 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 346 else: --> 347 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit) 348 else: 349 if infbounds != 0: /Users/colinross/Documents/Optics/Programming/five.py in integrand(x, y) 99 100 def integrand(x,y): --> 101 return abs(E_out(x))**2.*abs(E_(x - y))**2. 102 103 def G(y): /Users/colinross/Documents/Optics/Programming/five.py in E_out(x) 93 94 def E_out(x): ---> 95 E_in_w = fft(E_(x)) 96 omega = fftfreq(len(x),dt)*2*np.pi 97 E_out_w = E_in_w*M(omega) /Users/colinross/anaconda/lib/python2.7/site-packages/scipy/fftpack/basic.pyc in fft(x, n, axis, overwrite_x) 253 254 if n is None: --> 255 n = tmp.shape[axis] 256 elif n != tmp.shape[axis]: 257 tmp, copy_made = _fix_shape(tmp,n,axis) IndexError: tuple index out of range On Sun, Mar 15, 2015 at 11:57 PM, Danny Yoo <d...@hashcollision.org> wrote: > > What does fft expect to receive as an argument? We can read the > following: > > > > > http://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft > > > > Since fft is erroring out: there's only one possibility: E_(x) is not > > providing a value that's appropriate to fft(). > > > > Why would it do that? Two possibilities: > > > > 1. E_ is buggy and needs investigation. > > > > 2. E_ is fine, but the inputted value of x is one that E_x is not > > defined to handle. > > > > > > And so we start the investigation by considering those two > > possibilities. Start with #1. > > > > Do you know if E_ is fine? What is it supposed to do? What is the > > shape of its input? What is the shape of its output? Is its output > > something that fft() is designed to handle? > > Just to add: this is not to say that E_() is buggy. We just have to > start investigation somewhere, and it seemed as good a place as any, > since I don't know if _any_ of the functions are behaving. :P > > This is very much the reason why programmers like to do unit testing: > we want to know what functions at least do something that we know is > useful. We know all too well that whole programs break, and we want > to have some confidence on what components of our program are likely > ok. > > If #2 is the actual issue, then the question becomes: why is the > program producing an invalid input 'x'? And that's where we need to > start reasoning backwards, to discover how that value was constructed. > _______________________________________________ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor