def calcChiSq(p, x, y, yerr):
    '''
    Error function for fit
    '''
    e = (y - fitFunc(x, *p))/yerr
    return e

# First do the plot of Refractive index versus wavelength
fig = plt.figure(figsize = (8, 6))


xdata = wavelengths
ydata = refractive_indices

xerror = 0
yerror = d_refractive_indices

fig = plt.figure(figsize = (8, 6))
plt.title('Refractive index as a function of wavelength')
plt.xlabel('Wavelength (m)')
plt.ylabel('Refractive Index')
plt.grid(color = 'g')
plt.errorbar(xdata, ydata, yerror, xerror, marker=".", linestyle = '')
plt.show



# Now we need to do a fit to the data. The fitting function is defined below, and is non-linear. You should update the variable
# mu_0 to reflect _your_ own value for mu_0, the refractive index at lambda_0
lambda_0 = 4.7999e-7
mu_0 = 1.747

def fitFunc(x, alpha, beta):
    return mu_0 + alpha*x + beta*x**2

# Do the fit with our variables. You might need to edit the fit bounds to get an optimised fit

popt, pcov = curve_fit(fitFunc, (wavelengths-lambda_0), refractive_indices, bounds=([-2.5E5, -2.5E11], [-2.0E5, -2.0E11]))
print('Alpha: {0:.3e}, Beta: {1:.3e}'.format(*popt))

# Plot the fit
plt.plot(wavelengths, fitFunc( (wavelengths-lambda_0), *popt))

# Calculate the ChiSq/NDF of the fit, and print
chiarr = calcChiSq(popt, xdata, ydata, yerror)**2
chisq = np.sum(chiarr)
NDF = nPoints - nPars
chisqndf = chisq/NDF
print('ChiSq/NDF: {0:.2f}'.format(chisqndf))