From 2b90aafebb6885502da7fbbd90a10ec3bfcb882c Mon Sep 17 00:00:00 2001
From: kamischi <karl-michael.schindler@web.de>
Date: Tue, 19 Mar 2024 09:24:20 +0100
Subject: [PATCH] Add the python code from Hannes

---
 Python code/README.md            |   3 +
 Python code/nice_fitting_code.py | 592 +++++++++++++++++++++++++++++++
 2 files changed, 595 insertions(+)
 create mode 100644 Python code/README.md
 create mode 100755 Python code/nice_fitting_code.py

diff --git a/Python code/README.md b/Python code/README.md
new file mode 100644
index 0000000..a966f43
--- /dev/null
+++ b/Python code/README.md	
@@ -0,0 +1,3 @@
+# Python code for fitting EELS spectra using the modified Lambin code.
+
+From the thesis of Hannes Hermman.
diff --git a/Python code/nice_fitting_code.py b/Python code/nice_fitting_code.py
new file mode 100755
index 0000000..c52ca8e
--- /dev/null
+++ b/Python code/nice_fitting_code.py	
@@ -0,0 +1,592 @@
+#################### Load packages for fitting
+from os import truncate
+import matplotlib.pyplot as plt
+import numpy as np
+from libhreels import calcHREELS4 as ch4
+from libhreels.HREELS import HREELS
+from copy import deepcopy
+import lmfit 
+import os
+here = os.path.dirname(os.path.realpath(__file__)) 
+import re, json
+
+#################### Method to add a Drude term to phonon contribution
+def add3ParameterDrude(wTOPlasma, gTOPlasma, wLOPlasma, gLOPlasma, material):
+    newMaterial = deepcopy(material)
+    try:
+        if len(newMaterial['wTO']) > 0:
+            newMaterial['wTO'] += [0.]#postive -> flag for choosing kurosawa with wTO->0
+            newMaterial['gTO'] += [gTOPlasma]# positive!!!
+            newMaterial['wLO'] += [wLOPlasma]#-> wLO of plasma (sign not relevant!!!)
+            newMaterial['gLO'] += [gLOPlasma]# positive!!!
+            return newMaterial
+    except:
+        print('Cannot add Drude to material',material)
+    return material
+
+#################### Defining the fit function for HREELS
+def lossFunction(x, **par):
+    #################### Phonon contribution without any drude
+    stack_without_drude = ({'eps': par['epsi'],
+        'wTO': [par['wTO1'], par['wTO2'], par['wTO3']],
+        'gTO': [par['gTO1'], par['gTO2'], par['gTO3']],
+        'wLO': [par['wLO1'], par['wLO2'], par['wLO3']],
+        'gLO': [par['gLO1'], par['gLO2'], par['gLO3']]
+        })
+    
+    if FlagHighFreqPlasmonSetZero==True:
+        par['gTOP'] = 0.0
+        par['wLOP'] = 0.0
+        par['gLOP'] = 0.0
+    if FlagHighFreqPlasmonSetSymmetric==True:
+        par['gTOP'] = par['gLOP']
+
+    #################### Adding first Drude term to Phonon contribution
+    stack_with_drude = add3ParameterDrude(par['wTOP'],par['gTOP'],par['wLOP'],par['gLOP'],stack_without_drude)
+
+    if FlagLowFreqPlasmonSetZero==True:
+        par['gTOPP'] = 0.0
+        par['wLOPP'] = 0.0
+        par['gLOPP'] = 0.0
+    if FlagLowFreqPlasmonSetSymmetric==True:
+        par['gTOPP'] = par['gLOPP']
+
+    #################### Adding second Drude term to Phonon contribution + first drude term
+    stack_with_2nd_drude = add3ParameterDrude(par['wTOPP'],par['gTOPP'],par['wLOPP'],par['gLOPP'],stack_with_drude)
+
+    #################### Combining Layers with and without Drude part
+    extended_stack = [[stack_with_2nd_drude,par['DrudeThickness']],[stack_without_drude,par['STOthickness']]]  #STO+Drude+drude
+    
+    #################### Calculate HREELS sepctra with Lambin-Code
+    film1 = ch4.lambin(film=extended_stack, setup=setup)
+    film1.gauss = par['gauss'] 
+    film1.width = par['width']         
+    if FlagAsymSetOneMinusGauss==True:
+        film1.asym = 1 - par['gauss'] 
+    else:
+        film1.asym = par['asym']
+    xs, spectra_drude = film1.calcHREELS(x)
+    xReturn, model = xs, spectra_drude
+
+    return model*par['amp']
+
+
+
+
+#################### Get then data ################################
+d = "/mnt/c/Users/user/ownCloud/Phd_Halle/_code/_GPH-Data"
+file = "EBE02";    T=328.1     ;dope=0.05 # 0.05%
+T=300 ;dope=0.1
+
+d1 = HREELS(file,d)  
+
+
+##########################################################
+#################### Experimental geometry and parameters
+setup = {
+    "e0": 4.0,
+    "theta": 60.,
+    "phia": 0.35,        
+    "phib": 2.0,
+    "temperature": T,
+    "debug": False}
+
+
+#################### Initital Options for the Fit Model
+FlagWriteParameters=        True 
+FlagReadSpecial=            False
+FlagReadOldInit=            True
+FlagReadOldInitOnFirstCycle=False 
+FlagAdd3ParameterDrude=     True
+FlagPrintResult=            True
+FlagPrintEveryResult=       False
+
+
+initial_thickness=20.  
+FlagVaryThickness=True
+initial_wlop = 1200
+
+FlagAsymSetOneMinusGauss =          True
+FlagHighFreqPlasmonSetZero =        False
+FlagHighFreqPlasmonSetSymmetric =   False
+FlagLowFreqPlasmonSetZero =         True
+FlagLowFreqPlasmonSetSymmetric =    True
+
+
+#################### Define order fo Fitting steps as list
+FittingSteps=[0,4,0,9]
+
+m=-1
+for FittingStep in FittingSteps:
+    m=m+1
+
+    #################### Should the first fitting step use the last saved parameterset?
+    if m==0:
+        FlagReadOldInit=FlagReadOldInitOnFirstCycle
+    if m>0:
+        FlagReadSpecial=False
+        FlagReadOldInit=True  
+    
+    print('--------------------Start to fit with fittingstep:', FittingStep, ", step ",m, "/", len(FittingSteps))
+
+    #################### Assign flags to the current fitting step
+    if FittingStep<0: #just read
+        FlagReadOldInit=True
+        FlagWriteParameters=False
+        FlagFitElastic=False
+        FlagFitRest=False
+        FlagVaryPlasma=False
+        FlagVaryPlasmaFrequency=False
+        FlagWeightsForBothAreas=False
+    elif FittingStep==0: # fit elastic only
+        FlagFitElastic=True
+        FlagFitRest=False
+        FlagVaryPlasma=False
+        FlagVaryPlasmaFrequency=False
+        FlagWeightsForBothAreas=False
+    elif FittingStep==1: # fit phonons except elastic
+        FlagFitElastic=False
+        FlagFitRest=True
+        FlagVaryPlasma=False
+        FlagVaryPlasmaFrequency=False
+        FlagWeightsForBothAreas=False
+    elif FittingStep==2: # fit elastic + plasma
+        FlagFitElastic=True
+        FlagFitRest=False
+        FlagVaryPlasma=True
+        FlagVaryPlasmaFrequency=True
+        FlagWeightsForBothAreas=False
+    elif FittingStep==3: # fit phonons + plasma
+        FlagFitElastic=False
+        FlagFitRest=True
+        FlagVaryPlasma=True
+        FlagVaryPlasmaFrequency=True
+        FlagWeightsForBothAreas=True#False
+    elif FittingStep==4: # fit plasma only
+        FlagFitElastic=False
+        FlagFitRest=False
+        FlagVaryPlasma=True
+        FlagVaryPlasmaFrequency=True
+        FlagWeightsForBothAreas=True
+        #FlagVaryThickness=False
+    elif FittingStep==5: #only frequency of plasma
+        FlagFitElastic=False
+        FlagFitRest=True 
+        FlagVaryPlasma=True#False
+        FlagWeightsForBothAreas=True
+        FlagVaryPlasmaFrequency=False
+        #FlagVaryThickness=False
+    elif FittingStep==6:#vary thickness of layer on sto
+        FlagFitElastic=False
+        FlagFitRest=False
+        FlagVaryPlasma=False
+        FlagVaryPlasmaFrequency=False
+        FlagWeightsForBothAreas=True 
+    elif FittingStep==7:#vary thickness and plasma of layer on sto
+        FlagFitElastic=False
+        FlagFitRest=False
+        FlagVaryPlasma=True
+        FlagVaryPlasmaFrequency=True
+        FlagWeightsForBothAreas=True 
+    elif FittingStep==8: #fit only gLOP
+        FlagFitElastic=False
+        FlagFitRest=True 
+        FlagVaryPlasma=True#False
+        FlagWeightsForBothAreas=True
+        FlagVaryPlasmaFrequency=False
+    elif FittingStep==9: #vary everything
+        FlagFitElastic=True
+        FlagFitRest=True 
+        FlagVaryPlasma=True#False
+        FlagWeightsForBothAreas=True
+        FlagVaryPlasmaFrequency=True
+      
+    ############### Create a new figure for each fitting step
+    fig = plt.figure(num=None, figsize=(16, 9), dpi=80, facecolor='w', edgecolor='k')
+    ax =  fig.add_subplot(337)               # surface loss with dampings
+    ax2 = fig.add_subplot(331,sharex=ax)    # HREELS intensity with overlayer
+    ax3 = fig.add_subplot(334,sharex=ax)    # fit of rest
+    ax4 = fig.add_subplot(338)              # chi     
+    #ax5 = fig.add_subplot(339)             # empty
+    ax3b = fig.add_subplot(335)#,sharex=ax) #elastic peak fit
+
+    #################### Crop measured data
+    xd= d1.xdata[d1.findIndex(-80.1):d1.findIndex(1000)]
+    xd= d1.xdata[d1.findIndex(-80.1):]
+    yd= d1.ydata[d1.findIndex(-80.1):d1.findIndex(1000)]
+    yd= d1.ydata[d1.findIndex(-80.1):]
+    
+    #################### Height normalization and plot of data spectra
+    yd_save=deepcopy(yd)
+    maxCountRate = max(yd)
+    yd = yd_save/max(yd_save)
+    ax2.plot(xd,yd,'k-', label=file)
+    ax3.plot(xd,yd,'k-', label=file)
+    ax3b.plot(xd,yd,'k-', label=file)
+
+    #################### Flag management for ranges and weights
+    weights_elastic=1.5 
+    weights_rest=1
+    eps_parameter=5.2          
+
+    weights = [1. if (aa<40 and aa>-40) else 1. for aa in xd] #default weights
+
+    if FlagFitElastic==True:
+        print("fitting elastic peak")
+        FlagVaryElastic=True
+        weights = [1 if (aa<40 and aa>-40) else 0.0005 for aa in xd]   # Weights to fix elastic peak 
+    else:
+        FlagVaryElastic=False
+
+    if FlagFitRest==True:
+        print("fitting everything except elastic peak")
+        FlagVaryRest=True
+        weights = [1 if (aa>40 and aa<2000) else 0.0005 for aa in xd]   # No weight on elastic peak
+    else:
+        FlagVaryRest=False  
+
+    if FlagWeightsForBothAreas==True:
+        print("both regions are fitted")
+        weights = [weights_elastic if (aa<40 and aa>-40) else weights_rest for aa in xd]
+
+    #################### Define additional regions with heigher weights to get better fits of lineshpaes
+    ii=-1
+    for xx in xd:
+        ii=ii+1
+        if (xx > 10) and (xx <140):
+           weights[ii] =weights[ii]*1.0 
+        if (xx > 140) and (xx <200):
+           weights[ii] =weights[ii]*15.0 
+        if (xx > 350) and (xx <1530):
+           weights[ii] =weights[ii]*5.0
+        if (xx>450) and (xx<500):
+            weights[ii] =weights[ii]*5.0 
+        if (xx>725) and (xx<785):
+            weights[ii] =weights[ii]*5.0 
+        if (xx>1100) and (xx<1300):
+            weights[ii] =weights[ii]*3.0 
+        if (xx>1400) and (xx<1550):
+            weights[ii] =weights[ii]*3.0 
+
+
+    #################### Define fit function 
+    mod = lmfit.Model(lossFunction)
+    params = lmfit.Parameters()
+
+    #################### Define fitparamter for elastic peak
+    params.add('amp',   value= 1.,         min =0.001) 
+    params.add('gauss', value= 0.722,      min=0.6,    max=0.99)
+    params.add('width', value= 13.7189,    min=10.0,   max=35.0)
+    if FlagAsymSetOneMinusGauss == True:
+        params.add('asym',  value= -1.0,    min=-1.0,   max=0.99)
+    else:
+        params.add('asym',  value= 1-params['gauss'],   min=0.0, max=0.99)
+
+    #################### Define fitparamter for phonons
+    #Phonon parameter from LN2 STO fit
+    params.add('epsi',  value=5.6,          min=2.0,    max=8.0)
+
+    params.add('wTO1', value= 74.5533,      min=30.,    max=150.) 
+    params.add('gTO1', value=27.9999,       min=0.0001, max=12 )
+    params.add('wLO1', value=171.967,       min=80,     max=180.) 
+    params.add('gLO1', value=3.09999,       min=0.0001, max=10.)
+
+    params.add('wTO2', value=173.999 ,      min=160,    max=540.)
+    params.add('gTO2', value= 5.99999,      min=1.,     max=150.)
+    params.add('wLO2', value=442.999,       min=200.,   max=800.)
+    params.add('gLO2', value=4.62163,       min=0.2,    max=150.)
+
+    params.add('wTO3', value=542.0,         min=500,    max=800.)
+    params.add('gTO3', value=18.7496 ,      min=5.,     max=150.)
+    params.add('wLO3', value=661.556,       min=550.,   max=1000.) 
+    params.add('gLO3', value=25.0,          min=5.,     max=150.)
+
+    #################### Define auxillary fitparamter 
+    params.add('T',    value= T,    vary=False)
+    params.add('dope', value= dope, vary=False)
+    params.add('MaxTO3Data',  value=max(yd[100:]))
+    params.add('ratio',value=0.00,          min=0.01,   max=0.99, vary=True)
+    try:
+        params.add('MaxTO3NoDamp',   value= max(p2[100:]))
+    except:
+        params.add('MaxTO3NoDamp',   value= 0.)
+
+    #################### Define fitparamter for thickness of film on bulk
+    params.add('DrudeThickness', value=initial_thickness, min=4., max=30000., vary=True)
+    params.add('STOthickness',   value=30000,             min=0., max=30000., vary=False)
+
+    #################### Define fitparamter for drude term
+    params.add('wTOP',  value=0.,           min=0.,     max=0.1) 
+    params.add('gTOP',  value=20,           min=0.1 ,   max=10000000.)
+    params.add('wLOP',  value=initial_wlop, min=0.,     max=4000000.)
+    params.add('gLOP',  value=20,           min=0.1,    max=40000000.)
+
+    params.add('wTOPP', value=0.,           min=0.,     max=0.1) 
+    params.add('gTOPP', value=1.0,          min=0.)
+    params.add('wLOPP', value=5300,         min=100.,   max=20000.)
+    params.add('gLOPP', value=1.0,          min=0.)
+
+    if FlagHighFreqPlasmonSetSymmetric == True:
+        params.add('gTOP',    value=-1, min=-1. , max=-1.01)
+    if FlagLowFreqPlasmonSetSymmetric == True:
+        params.add('gTOPP',    value=-1, min=-1. , max=-1.01)
+
+    if FlagHighFreqPlasmonSetZero==True:
+        params.add('wTOP',    value=0., min=0.,  max=0.1) 
+        params.add('gTOP',    value=0.1, min=0.,  max=1000.)  
+        params.add('wLOP',    value=0., min=0.,  max=1000.)
+        params.add('gLOP',    value=0.1, min=0.,  max=1000.)
+        FlagVaryPlasma=False 
+        FlagVaryPlasmaFrequency=False 
+    if FlagLowFreqPlasmonSetZero==True:
+        params.add('wTOPP',    value=0., min=0.,  max=0.1) 
+        params.add('gTOPP',    value=0.1, min=0.,  max=1000.)  
+        params.add('wLOPP',    value=0., min=0.,  max=1000.)
+        params.add('gLOPP',    value=0.1, min=0.,  max=1000.)
+        FlagVaryLowFreqPlasma=False 
+        FlagVaryLowFreqPlasmaFrequency=False 
+    
+        
+
+    #################### Read old paramter from parameterfile
+    if FlagReadOldInit==True:
+        try:
+            oldInit_params = lmfit.Parameters()
+            if FlagAdd3ParameterDrude==False:
+                with open(d+"/"+file+"_init_parameters.txt","r") as h:
+                    init_params=h.readline()
+                filetext=d+"/"+file+"_init_parameters.txt"
+            if FlagAdd3ParameterDrude==True:
+                with open(d+"/"+file+"_init_3parameters.txt","r") as h:
+                    init_params=h.readline()
+                filetext=d+"/"+file+"_init_3parameters.txt"
+            if FlagReadSpecial==True:
+                with open(d+"/"+specialfile+"_init_3parameters.txt","r") as h:
+                    init_params=h.readline() 
+                filetext=d+"/"+specialfile+"_init_3parameters.txt"
+            #print(init_params)
+            init_params = init_params.replace(" (fixed)", ", vary=false")
+            for i in range(len(init_params)):
+                if init_params[i]=='+':
+                    for j in range(15):
+                        if init_params[i+j]==',':
+                            init_params=init_params[:i-1] + init_params[i+j:]
+                            break
+                if init_params[i:i+5]==">)])":
+                    break
+            matches = re.findall('<Parameter ([^>]*)>', init_params)
+            for e in matches:
+                e = json.loads("{\"name\":"+e.replace("=","\":").replace(", ",", \"").replace("'","\"").replace("[","\"").replace("]","\"")+"}")
+                oldInit_params.add(e["name"], value=e["value"], min=float((e["bounds"].split(":"))[0]),max=float((e["bounds"].split(":"))[1]), vary=(False if "vary" in e else True ))
+            oldInit_params.add('T',   value= T)
+            oldInit_params.add('dope',   value= dope)
+            params = oldInit_params
+            print("Old initial parameters loaded from ", filetext)    
+        except:
+            print('no initfile found')    
+
+
+    #################### Manage all parameter which should be varried this fitting step
+    params['T'].vary =              False
+    params['dope'].vary =           False
+    params['MaxTO3Data'].vary =     False
+    params['MaxTO3NoDamp'].vary =   False
+
+    params['amp'].vary =            FlagVaryElastic
+    params['gauss'].vary =          FlagVaryElastic
+    params['width'].vary =          FlagVaryElastic
+    params['asym'].vary =           FlagVaryElastic
+    params['DrudeThickness'].vary = FlagVaryThickness
+    params['STOthickness'].vary =   False
+
+    params['epsi'].vary =           True #FlagVaryPlasma
+
+    params['wTO1'].vary = FlagVaryRest
+    params['gTO1'].vary = FlagVaryRest
+    params['wLO1'].vary = FlagVaryRest
+    params['gLO1'].vary = FlagVaryRest
+    params['wTO2'].vary = FlagVaryRest
+    params['gTO2'].vary = FlagVaryRest
+    params['wLO2'].vary = FlagVaryRest
+    params['gLO2'].vary = FlagVaryRest
+    params['wTO3'].vary = FlagVaryRest
+    params['gTO3'].vary = FlagVaryRest
+    params['wLO3'].vary = FlagVaryRest
+    params['gLO3'].vary = FlagVaryRest
+
+    params['wTOP'].vary = False
+    params['gTOP'].vary = FlagVaryPlasma
+    params['wLOP'].vary = FlagVaryPlasma 
+    params['gLOP'].vary = FlagVaryPlasma 
+
+    params['wTOPP'].vary = False
+    params['gTOPP'].vary = FlagVaryPlasma 
+    params['wLOPP'].vary = FlagVaryPlasma
+    params['gLOPP'].vary = FlagVaryPlasma 
+
+    params.add('MaxTO3Data',   value= max(yd[100:]))
+    params.add('T',   value= T)
+    params.add('dope',   value= dope)
+
+    #################### Do the actual fitting
+    initial = mod.eval(params, x=xd)
+    result = mod.fit(yd, params=params, x=xd, weights=weights)
+
+
+    #################### Plot fit results
+    ax2.plot(xd,result.best_fit, 'r--', label='Fit')
+    ax3.plot(xd,result.best_fit, 'r--', label='Fit')
+    ax3b.plot(xd,result.best_fit, 'r--', label='Fit')
+
+    #################### Create different scenarios f.e. only phonons
+    para2 = deepcopy(result.params) #no damping case
+    if FlagAdd3ParameterDrude==True:
+        para2.add('wTOP', value=0., vary=False)
+        para2.add('gTOP', value=10., vary=False)
+        para2.add('wLOP', value=0., vary=False)
+        para2.add('gLOP', value=10., vary=False)
+    p2 = mod.eval(para2, x=xd)
+
+    para2b = deepcopy(result.params) #no damping case
+    para2b.add('wTOP', value=0., vary=False)
+    para2b.add('gTOP', value=10., vary=False)
+    para2b.add('wLOP', value=0., vary=False)
+    para2b.add('gLOP', value=10., vary=False)
+    para2b.add('wTOPP', value=0., vary=False)
+    para2b.add('gTOPP', value=10., vary=False)
+    para2b.add('wLOPP', value=0., vary=False)
+    para2b.add('gLOPP', value=10., vary=False)
+    p2b = mod.eval(para2b, x=xd)
+
+    para3 = deepcopy(result.params)
+    para3.add('wTOP', value=0., vary=False)
+    para3.add('wLOP', value=0., vary=False)
+    para3.add('T',value=300., vary=False)
+    p3 = mod.eval(para3, x=xd)
+
+    para4 = deepcopy(result.params)
+    para5 = deepcopy(result.params)
+    para6 = deepcopy(result.params)
+    para7 = deepcopy(result.params)
+    para4.add('wTOP', value=0.,   vary=False)
+    para4.add('wLOP', value=0.,   vary=False)
+    para4.add('T',    value=100., vary=False)
+    para5.add('DrudeThickness', value=0.3, vary=False)
+    para6.add('DrudeThickness', value=30, vary=False)
+    para7.add('DrudeThickness', value=300, vary=False)
+    p4 = mod.eval(para4, x=xd)
+    p5 = mod.eval(para5, x=xd)
+    p6 = mod.eval(para6, x=xd)
+    p7 = mod.eval(para7, x=xd)
+
+    #################### Plot different models
+    ax2.plot(xd, p2,'b-', label='No damping second plasmon, TO3 at '+str(max(p2[100:]).round(4)))
+    ax2.plot(xd, p2b,'g-', label='No damping both plasmons, TO3 at '+str(max(p2b[100:]).round(4)))
+    ax2.plot(xd, p3,'c-', label='T=300 K')
+    ax2.plot(xd, p4,'-', label='T=100 K, Thick='+str(para4['DrudeThickness'].value))
+    ax2.plot(xd, p5,'-', label='T=100 K, Thick='+str(para5['DrudeThickness'].value))
+    ax2.plot(xd, p6,'-', label='T=100 K, Thick='+str(para6['DrudeThickness'].value))
+    ax2.plot(xd, p7,'-', label='T=100 K, Thick='+str(para7['DrudeThickness'].value))
+
+
+    #################### Plot fit weigths
+    weights_scaled = weights
+    for ii in range(len(weights)):
+        weights_scaled[ii-1] = weights[ii-1]/100
+    ax3.plot(xd, weights_scaled,'g--', label='weights/100')
+    ax3b.plot(xd, weights_scaled,'g--', label='weights/100')
+
+    '''#################### Plot surface loss
+    xSL, dataSurfaceLoss = modSurfaceLoss.eval(para2, x=xd)
+    ax.plot(xSL, dataSurfaceLoss,'b-', label="no damping from high freq plamon")
+    xSL, dataSurfaceLoss = modSurfaceLoss.eval(para2b, x=xd)
+    ax.plot(xSL, dataSurfaceLoss,'g-', label="no damping, both plasmons")
+    xSL, dataSurfaceLoss = modSurfaceLoss.eval(result.params, x=xd)
+    ax.plot(xSL, dataSurfaceLoss,'r-', label="with damping")
+    ax4.plot(xd,yd-result.best_fit,'k-', label='Data - Fit, Chi^2*1E6: '+str(1000000*result.chisqr.round(10)))'''
+    
+    #################### Plot legends
+    ax.legend()
+    ax2.legend()
+    ax3.legend()
+    ax3b.legend()
+    ax4.legend()
+
+    #################### Setting axis ranges
+    ax2.set(ylabel='HREELS Intensity',
+        ylim=[-0.001, 1.1*max(p2[300:])]
+        )
+    ax3.set(ylabel='HREELS Intensity',
+        ylim=[0.0, 1.1*max(p2[300:])]
+        )
+    ax3b.set(ylabel='HREELS Intensity', #elastic peak ist highlighted
+        xlim=[-60.0, 60.0], ylim=[-0.001, 1.1*max(yd)]
+        )    
+    ax4.set(ylabel='Fitting error')#,ylim=[-1.0, 1.0])    
+    ax.set(ylabel='Surface loss',
+        xlabel='Energy Loss (cm$^{-1}$)',
+        xlim=[-75,1600],ylim=[0,0.009])
+
+    #################### Write fit summary in plot figure
+    plt.text(0.75, 0.99, result.fit_report(), fontsize=10, ha='left', va='top', transform=plt.gcf().transFigure)
+    plt.text(0.99, 0.03,'Iterationstep '+str(m)+'/'+str(len(FittingSteps))+', Fittingstep: '+str(FittingStep), fontsize=10, ha='right', va='bottom', transform=plt.gcf().transFigure)
+    fig.tight_layout(pad=3)
+
+
+    #################### Write parameters in auxillary file
+    if FlagWriteParameters==True:
+        if FlagAdd3ParameterDrude==False:
+                hilf=d+"/"+file+"_fit_parameters.txt"
+        if FlagAdd3ParameterDrude==True:
+                hilf=d+"/"+file+"_fit_3parameters.txt"    
+        print("Parameters exported to: ", hilf)
+        f = open(hilf, "w")
+        f.write(result.fit_report())
+        f.close()
+
+        if FlagAdd3ParameterDrude==False:
+                hilf=d+"/"+file+"_init_parameters.txt"
+        if FlagAdd3ParameterDrude==True:
+                hilf=d+"/"+file+"_init_3parameters.txt"
+
+        g = open(hilf, "w")
+        init_text=str(deepcopy(result.params))
+        init_text=init_text.replace("Parameter 'MaxTO3NoDamp', value=0.0", "Parameter 'MaxTO3NoDamp', value="+str(max(p2[100:])))
+        g.write(init_text)
+        g.close()
+
+    #################### Print fit result of every fitting step in console
+    if FlagPrintEveryResult==True:
+        print(result.fit_report())
+
+#################### Set the dataname for saved figure
+figurelabel = "_O"
+if FlagAsymSetOneMinusGauss == True:
+    figurelabel = figurelabel+"_AsymDepOnGauss"
+if FlagHighFreqPlasmonSetZero == True:
+    figurelabel = figurelabel+"_HFreqPlas0"
+if FlagHighFreqPlasmonSetSymmetric == True:
+    figurelabel = figurelabel+"_HFreqPlasSymm"
+if FlagLowFreqPlasmonSetZero == True:
+    figurelabel = figurelabel+"_LFreqPlas0"
+if FlagLowFreqPlasmonSetSymmetric == True:
+    figurelabel = figurelabel+"_HFreqPlasSymm"
+
+#################### Write dataname on plot figure and save it
+plt.text(0.99, 0.05, "Selected options: "+figurelabel, fontsize=10, ha='right', va='bottom', transform=plt.gcf().transFigure)
+plt.text(0.99, 0.01,os.path.basename(__file__), fontsize=10, ha='right', va='bottom', transform=plt.gcf().transFigure)
+output_filename = os.path.splitext(__file__)[0] + '_'+file+'_'+figurelabel+'.png'
+plt.savefig(output_filename)
+
+#################### Print fit result of final fitting step in console
+if FlagPrintResult==True:
+    print(result.fit_report())
+print("Chi squared: ", result.chisqr)
+print("nw. Chi squared: ", sum(yd-result.best_fit))#result.chisqr)
+print("TO3 without damping at: ", str(max(p2[100:]).round(6)))
+
+#################### Plot everything
+plt.show()
+
+
+
-- 
GitLab