# [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.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                                                         #
>
> 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 = []
>
>
> 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 =
>
>     282     else:
>
>     283         retval =
>
>
> in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
>
>     345             return
>
>     346         else:
>
> --> 347             return
>
>     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
>> >
>> > 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.
>>
>
>