I am have experimental data for the double slit experiment that I wish to fit a numerical curve to. I have opted to using the lmfit library because I like their fit result output a lot more than other curve fitting packages available to python. To start, I will post my code, then go on to explain the errors I am seeing and what I have read to try and better explain my problem
So here is my code
from PIL import Image
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.signal import argrelextrema
from lmfit.models import Model
#Determine the observed intensity distribution
img_data = pd.read_csv('Image31 intensity profile.csv')
pixels = np.array(img_data['Distance_(_)'])
counts = np.array(img_data['Gray_Value'])
norm_intensity = counts / np.linalg.norm(counts)
# Double slit interferometry coherent beam
def keV2m(keV):
"""Converts Photon Energy [keV] to wavelength [m].
.. note:: Calculation in Vacuum.
"""
wl = 1./(keV*1000)*4.1356*(10**(-7))*2.9998
return wl
ener = 10#kev
I_0 = 1
a = 2.0e-6 # slit width in m
d = 10e-6 # slit seperation in m
L = 6.0 # slit to detec dist in m
detec_pixel_size = 3.1e-6 #pixel size
x = np.arange(-500,500)*detec_pixel_size + 1e-9
def coh_inter_patt(x):
lam = keV2m(ener)
sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2)
sinq_part[sinq_part == np.nan] = 1
inter_pat = I_0*(np.cos((np.pi*x*d)/(lam*L))**2)*sinq_part
return inter_pat
def coh_inter_patt_sft(x,detec_pixel_size,sft=10):
lam = keV2m(ener)
sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2)
sinq_part[sinq_part == np.nan] = 1
inter_pat = I_0*(np.cos((np.pi*(x+(sft*detec_pixel_size))*d)/(lam*L))**2)*sinq_part
return inter_pat
inter_pat = coh_inter_patt(x)
inter_pat_sft = coh_inter_patt_sft(x,detec_pixel_size,sft=10)
max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])
min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])
visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])
# 'Blurred' patterns
extension = 40 # 10*detec_pixel_size long extended source
blur_patt = np.zeros_like(coh_inter_patt(x))
for i in range(extension):
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)
blur_patt = blur_patt/np.linalg.norm(blur_patt)
max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])
min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])
visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])
blur_max_vals = np.sort(argrelextrema(blur_patt, np.greater)[0])
blur_min_vals = np.sort(argrelextrema(blur_patt, np.less)[0])
blur_visibility = (blur_patt[max_vals] - blur_patt[min_vals[1:]])/(blur_patt[max_vals] + blur_patt[min_vals[1:]])
#Look at plot output
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(111)
ax1.plot(pixels, norm_intensity, color='red', label='Observed')
ax1.plot(blur_patt,label='Analytical, vis %f'%blur_visibility.mean(), color='k')
ax1.legend()
#iterative curve fitting
def Blurred_fringes(N):
extension = N # 10*detec_pixel_size long extended source
blur_patt = np.zeros_like(coh_inter_patt(x))
for i in range(extension):
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)
#inter_pat = inter_pat/np.linalg.norm(inter_pat)
blur_patt = blur_patt/np.linalg.norm(blur_patt)
return blur_patt
model = Model(Blurred_fringes)
params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)
Any variable or calculation involving maxima, minima, and visibility are a seperate calculation and not involved in my problem as far I can tell. So my issue with the code that I posted is that I run the error, "Key Error: N". I know in general that the lmfit fit() function wants dictionary objects made from the make_params() function. And I know from this article here that the first argument to the model function is considered an independent variable by default
To try and work around that, I changed my curve fitting block to the following:
X = np.arange(0,1000,1)
def Blurred_fringes(x, N):
extension = N # 10*detec_pixel_size long extended source
blur_patt = np.zeros_like(coh_inter_patt(x))
for i in range(int(extension)):
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)
#inter_pat = inter_pat/np.linalg.norm(inter_pat)
blur_patt = blur_patt/np.linalg.norm(blur_patt)
return blur_patt
test_y = Blurred_fringes(x, 20)
plt.plot(X,test_y)
plt.show()
model = Model(Blurred_fringes)
params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)
That is, I defined a numpy arange that covers the same x range as before, and fed that into my function and fitting commands. And I am instead met with this error "ValueError: The model function generated NaN values and the fit aborted! Please check your model function and/or set boundaries on parameters where applicable.". I believe these NaN values are associated with these output warnings:
:30: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pixa)/(lamL))/((np.pixa)/(lamL)))2) :41: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pi(x+(sftdetec_pixel_size))a)/(lamL))/((np.pi*(x+(sft*detec_pixel_size))a)/(lamL)))**2)
In general, you may run into divide by 0 errors with sinc functions. But I sorted that out with sinq_part[sinq_part == np.nan] = 1.
All in all, my end goal is to get a fit working on this function before I add an iteration step where I look at the confidence levels for a range of N values and find the best fit that way. N should be the only fitting parameter in my case. All of the other variables defined much ealier are experimental parameters and will be fixed (may change slightly due to associated error within those values, for example the sampe-to-detector distance might be 6.2 instead of 6.0). If anyone has any insight on how to proceed, please let me know. Thanks!