I am trying to fit two sweeps that take on a Fraunhofer distribution. Please help me to understand a snippet in my code to make it run better, because it changes the fitting, but I don't know why! This is the line I mean specifically:
if MaxIndexup>MaxIndexdown: FitDataup = TempData.loc[TempData["Field (Oe)"] <= MaxIndexup]
When this runs, it puts a bunch of files with a lot of data that is irrelevant to the question at hand. Which is how do you identify the peak of the (almost) Gaussian distribution, and then how do you fit it so the lobes work better (if you understand lobes given by a Fraunhofer diffraction pattern). I have my code working, I just want/need to understand more so how to parse through and fit the data better:
from os import walk
import numpy as np
import pandas as pd
import re
import math
from scipy import stats
import scipy.special as sp
from lmfit import Model
#create function that computes the Fruanhofer diffraction pattern: reference--> https://en.wikipedia.org/wiki/Fraunhofer_diffraction; also see the Guassian section
def Fraunup(x, Imaxup, Lambda, M, CorrectedBitSize):
y1 = Imaxup * abs(2*sp.jv(1,((float(CorrectedBitSize)*(2*abs(Lambda) + 2*dN + dF)*(x+M))/(phi0*math.pow(10,10))))/((float(CorrectedBitSize)*(2*Lambda + 2*dN + dF)*(x+M))/(phi0*math.pow(10,10))))
return y1
def Fraundown(x, Imaxdown, Lambda, M, CorrectedBitSize):
y2 = Imaxdown * abs(2 * sp.jv(1, ((float(CorrectedBitSize) * (2 * abs(Lambda) + 2 * dN + dF) * (x + M)) / (phi0 * math.pow(10, 10)))) / ( (float(CorrectedBitSize) * (2 * Lambda + 2 * dN + dF) * (x + M)) / (phi0 * math.pow(10, 10))))
return y2
Iparamup = 1.5 * TempData["Ic+ (mA)"].max()
Iparamdown = 1.5 * TempData2["Ic+ (mA)"].max()
Lambparam = (0.00000009)
Mparam = (200 / CorrectedBitSize_Squared) #gives you H_shift
#Create the Model and add specfic parameters. M deals with the H_Shift
modelup = Model(Fraunup)
params = modelup.make_params()
params.add('Imaxup',value=Iparamup)
params.add('Lambda',value=Lambparam)
params.add('M',value=np.mean(Mparam),min=0,max=500)
modeldown = Model(Fraundown)
params = modeldown.make_params()
params.add('Lambda',value=Lambparam)
params.add('M',value=np.mean(Mparam),min=0,max=500)
params.add('Imaxdown', value=Iparamdown)
if MaxIndexup>MaxIndexdown:
FitDataup = TempData.loc[TempData["Field (Oe)"] <= MaxIndexup]
FitDataup = FitDataup[~np.isnan(FitDataup['Ic+ (mA)'])].reset_index()
FitDatadown = FitDatadown[~np.isnan(FitDatadown['Ic+ (mA)'])].reset_index()
for x in range(len(Mparam)):
resultsup = modelup.fit(FitDataup['Ic+ (mA)'],x=-FitDataup['Field (Oe)'], CorrectedBitSize=CorrectedBitSize[x], Imaxup=Iparamup, Lambda=Lambparam, M= Mparam[x], method='Powell')
resultsdown = modeldown.fit(FitDatadown['Ic+ (mA)'], x=-FitDatadown['Field (Oe)'], CorrectedBitSize=CorrectedBitSize[x],Imaxdown=Iparamdown, Lambda=Lambparam, M=Mparam[x], method='Powell')
#fit the value using the best results
YFitValueup = Fraunup(TempData["Field (Oe)"], resultsup.best_values['Imaxup'], resultsup.best_values['Lambda'], resultsup.best_values['M'], resultsup.best_values['CorrectedBitSize'])
YFitValuedown = Fraundown(TempData2["Field (Oe)"], resultsdown.best_values['Imaxdown'], resultsdown.best_values['Lambda'], resultsdown.best_values['M'], resultsdown.best_values['CorrectedBitSize'])