0

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'])
James Draper
  • 5,110
  • 5
  • 40
  • 59
Tyler Ram
  • 13
  • 5
  • Working code belongs on [CodeReview](https://codereview.stackexchange.com/) I suggest you migrate your question there. – James Draper Aug 21 '17 at 22:36
  • I think your code is not complete: `TempData` is not defined, and `FitDataup` is defined only if `MaxIndexup > MaxIndexdown`, neither of which are defined. Please create a self-contained example. – M Newville Aug 23 '17 at 03:32

0 Answers0