[Tutor] Fitting data to error function
Colin Ross
colin.ross.dal at gmail.com
Mon Mar 16 22:50:50 CET 2015
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 <dyoo at 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.
>
More information about the Tutor
mailing list