In curvefit: TypeError: only size-1 arrays can be converted to Python scalars

Question

This error is recieved in jupyter, which is basicly python for as far I'm aware. I think I'm aware what this error means, however I don't see the problem with the code. Thanks in advance! 

import matplotlib.pyplot as plt

import numpy as np

import scipy.optimize 

from math import *

import matplotlib.patches as mpatches

X = [0,-0.0000959936,-0.000194899,-0.000273446,-0.000430558,-0.000366546,-0.000573153,-0.000657563,-0.000768193,0.000066904,0.000148355,0.000232717,

0.000285083,0.000427648,0.000401461,0.000526588,0.0000465422,-0.0000250891,0.0000436333,-0.0000261799,-0.0000610866,0.0000930847,-0.000101811]

Y = [0.96,0.47,0.11,0.25,0.15,0.12,0.13,0.11,0.11,0.78,0.16,0.19,0.13,0.14,0.16,0.13,0.6,2.02,1.91,1.4,0.94,0.68,0.97]

onzekerheidY = [0.0001,0.000001,-0.000101,-0.000182,-0.000344,-0.000278,-0.000491,-0.000578,-0.000692,0.000169,0.000253,0.00034,

0.000394,0.000541,0.000514,0.000643,0.000148,0.000074125,0.000145,0.000073,0.000037,0.000196,-0.000005]

def f(x, I0, theta, a, d):

    k = 9666438.934

    N = 1

    return I0 * (  ( (np.sin(k*(d/2)*np.sin(x)))/(k*(d/2)*np.sin(x)) * (np.sin(k*a*np.sin(x) *(N/2)))/(N*np.sin((k/2)*a*sin(x))) )   )^2

waarde,matrix = scipy.optimize.curve_fit(f,X,Y)

plt.plot(X,Y,'p',label='data')

a = np.linspace(-0.0300,0.0300,1000)

b = b = f(a, waarde[0], waarde[1],waarde[2],waarde[3])

plt.plot(a,b,'r',label= 'fit')

plt.errorbar(X,Y,onzekerheidY,linestyle='None',label = 'onzekerheden')

plt.grid()

plt.legend(loc = 'lower left')

plt.title('Fit van $R_{TM}$ in functie van $\Theta$ $_i$')

plt.xlabel('$\Theta$ $_i$')

plt.ylabel('$R_{TM}$')

And this is the error code:

TypeError                                 Traceback (most recent call last)
<ipython-input-9-b80cfce3b4ff> in <module>()
      4     return I0*( (np.sin(k*(d/2)*np.sin(x)))/(k*(d/2)*np.sin(x)) * (np.sin(k*a*np.sin(x) *(N/2)))/(N*np.sin((k/2)*a*sin(x)))    )^2
      5 
----> 6 waarde,matrix = scipy.optimize.curve_fit(f,X,Y)
      7 
      8 plt.plot(X,Y,'p',label='data')

C:\Users\boent\Anaconda2\lib\site-packages\scipy\optimize\minpack.pyc in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    749         # Remove full_output from kwargs, otherwise we're passing it in twice.
    750         return_full = kwargs.pop('full_output', False)
--> 751         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    752         popt, pcov, infodict, errmsg, ier = res
    753         cost = np.sum(infodict['fvec'] ** 2)

C:\Users\boent\Anaconda2\lib\site-packages\scipy\optimize\minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    381     if not isinstance(args, tuple):
    382         args = (args,)
--> 383     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    384     m = shape[0]
    385     if n > m:

C:\Users\boent\Anaconda2\lib\site-packages\scipy\optimize\minpack.pyc in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     25 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     26                 output_shape=None):
---> 27     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     28     if (output_shape is not None) and (shape(res) != output_shape):
     29         if (output_shape[0] != 1):

C:\Users\boent\Anaconda2\lib\site-packages\scipy\optimize\minpack.pyc in func_wrapped(params)
    461     if transform is None:
    462         def func_wrapped(params):
--> 463             return func(xdata, *params) - ydata
    464     elif transform.ndim == 1:
    465         def func_wrapped(params):

<ipython-input-9-b80cfce3b4ff> in f(x, I0, theta, a, d)
      2     k = 9666438.934
      3     N = 1
----> 4     return I0*( (np.sin(k*(d/2)*np.sin(x)))/(k*(d/2)*np.sin(x)) * (np.sin(k*a*np.sin(x) *(N/2)))/(N*np.sin((k/2)*a*sin(x)))    )^2
      5 
      6 waarde,matrix = scipy.optimize.curve_fit(f,X,Y)

TypeError: only size-1 arrays can be converted to Python scalars