[Tutor] Fitting data to error function
Colin Ross
colin.ross.dal at gmail.com
Mon Mar 16 22:55:46 CET 2015
What I am trying to do is calculate the non-colinear autocorrelation:
G(t_d) = \int_{-\infty}^{+\infty} |E(t)|^2 * |E(t - t_d)|^2 dt
So I need to loop through an array of t_d values (len = 376) and calculate
G(t_d) for as many t values as possible to eliminate sampling issues.
Colin
On Mon, Mar 16, 2015 at 6:50 PM, Colin Ross <colin.ross.dal at gmail.com>
wrote:
> 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