From 34c7192ab321b5ca64c6677657f5bc16d6a62f7c Mon Sep 17 00:00:00 2001
From: kamischi <karl-michael.schindler@web.de>
Date: Tue, 19 Mar 2024 15:08:36 +0100
Subject: [PATCH] remove obsolete version numbers.
Update the versioned file name to the basic name and remove lower versions.
---
Python code/libhreels/calcHREELS.py | 116 ++++++----
Python code/libhreels/calcHREELS3.py | 250 ---------------------
Python code/libhreels/calcHREELS4.py | 250 ---------------------
Python code/libhreels/dielectrics.py | 243 +++++++++++++-------
Python code/libhreels/dielectrics2.py | 245 --------------------
Python code/libhreels/dielectrics2_save.py | 250 ---------------------
Python code/nice_fitting_code.py | 4 +-
7 files changed, 239 insertions(+), 1119 deletions(-)
delete mode 100755 Python code/libhreels/calcHREELS3.py
delete mode 100755 Python code/libhreels/calcHREELS4.py
delete mode 100755 Python code/libhreels/dielectrics2.py
delete mode 100755 Python code/libhreels/dielectrics2_save.py
diff --git a/Python code/libhreels/calcHREELS.py b/Python code/libhreels/calcHREELS.py
index ef71336..989dd5b 100755
--- a/Python code/libhreels/calcHREELS.py
+++ b/Python code/libhreels/calcHREELS.py
@@ -5,6 +5,8 @@ import sys, re, os
from libhreels.HREELS import myPath
from copy import deepcopy
+import scipy.integrate as integrate
+
libDir = os.path.dirname(os.path.realpath(__file__))
if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')):
@@ -12,11 +14,10 @@ if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')):
have been complied with the proper f2py3 for the right
python version!!''')
try:
- from libhreels import myEels2 as LambinEELS # wrapper for myEels2.f90
- from libhreels import myBoson as LambinBoson # wrapper for myBoson.f90
- from libhreels import myEels3 as LambinEELS3 # wrapper for myEels3.f90
+ from libhreels import DoEELS as LambinEELS # wrapper for Doeels.f90
+ from libhreels import DoBoson as LambinBoson # wrapper for Doboson.f90
except:
- print('myEels2 and MyBoson are not available here (Check your version)')
+ print('DoEELS and DoBoson are not available here (Check your version)')
# Experimental setup as dictionary:
setup = {
@@ -143,25 +144,6 @@ class lambin:
i += 1
return wn_array[i-1:], loss_array[i-1:]
- def calcSurfaceLoss3(self,x):
- ''' Calculate the surface loss spectrum for the array of x, which needs to be an equidistant array.
- All parameters are defined in the class __init__() call.'''
- wmin = min(x)
- wmax = max(x)-0.001
- dw = (wmax-wmin)/(len(x)-1) # assumes that x is an equidistant array
- wn_array_size = len(x) # size of array for x and epsilon (wn_array, loss_array)
- nper = 1.
- contrl = '{:<10}'.format('None'[:10]) # Can be 'image' to include image charge
- mode = '{:<10}'.format('kurosawa'[:10])
- wn_array,loss_array = LambinEELS3.mod_doeels.doeels(self.e0,self.theta,self.phia,self.phib,
- wmin,wmax,dw,self.layers,self.neps,nper,self.name,
- self.thick,self.epsinf,self.listNOsci,self.osc,contrl,mode,wn_array_size)
- i=0
- for item in wn_array:
- if item > 0: break
- i += 1
- return wn_array[i-1:], loss_array[i-1:]
-
def calcHREELS(self,x, normalized=True, areanormalized=False):
emin = min(x)
emax = max(x)-0.001
@@ -169,10 +151,11 @@ class lambin:
xLoss,loss_array = self.calcSurfaceLoss(x)
wmin = min(xLoss)
wmax = max(xLoss)
- xOut,spectrum,n = LambinBoson.doboson3(self.temperature,self.width,self.gauss,self.asym,
+ xOut,spectrum,n = LambinBoson.doboson(self.temperature,self.width,self.gauss,self.asym,
emin,emax,wmin,wmax,loss_array,self.debug,len(loss_array))
if normalized:
norm = max(spectrum[:n])
+ #print("peak normalized, norm=", norm)
if areanormalized: #edit by HHE
try:
areanormalize_xstart = np.argmin(abs(x+100.)) #seems to be oddly complicated, but is way more stable than x.index(-100.) or where()
@@ -182,8 +165,17 @@ class lambin:
areanormalize_xend = np.argmin(abs(x-1000.))
except:
areanormalize_xend = len(x)
- norm=spectrum[areanormalize_xstart:areanormalize_xend].sum()
+ #norm=spectrum[areanormalize_xstart:areanormalize_xend].sum()
+ cropped_spectra=spectrum[areanormalize_xstart:areanormalize_xend]
+ cropped_x=x[areanormalize_xstart:areanormalize_xend]
+ #norm=cropped_spectra.sum()/len(cropped_x)*(x[areanormalize_xend]-x[areanormalize_xstart])
+ #print(len(cropped_x), len(cropped_spectra), len(x), len(spectrum))
+ #print(cropped_spectra.sum(), len(cropped_x),x[areanormalize_xend]-x[areanormalize_xstart])
+ norm=integrate.simps(cropped_spectra, dx=x[areanormalize_xstart+1]-x[areanormalize_xstart])#/(1100/(max(cropped_x)-min(cropped_x)))
+ #print("areanormalized, norm=", norm, " dx=", x[areanormalize_xstart+1]-x[areanormalize_xstart])
+ else:
+ print("not normalized")
return xOut[:len(x)], spectrum[:len(x)]/norm
def calcEps(self, x):
@@ -194,29 +186,65 @@ class lambin:
epsArray.append(yn)
return np.transpose(np.array(epsArray))
- def calcEps3(self, x):
- epsArray = []
- nOsci = len(self.wOsc)
- for wn in x:
- yn = LambinEELS3.mod_doeels.seteps(self.listNOsci,nOsci,self.osc,self.epsinf,wn,self.layers)
- epsArray.append(yn)
- return np.transpose(np.array(epsArray))
-
####################################################################################
def myMain():
import matplotlib.pyplot as plt
- oxide = importMaterials('NiO')
- substrate = importMaterials('Ag')
- x = np.linspace(-100.,1000.,1000)
- for d in [4., 20., 50., 200., 10000.]:
- film1 = lambin(film=[[oxide,d],[substrate,1000.]])
- # xs, spectrum = film1.calcSurfaceLoss(x)
- xs, spectrum = film1.calcHREELS(x, normalized=True, areanormalized=True)
- plt.plot(xs,spectrum, label='{:.0f}'.format(d))
- plt.ylabel('HREELS Intensity')
+ import numpy as np
+ from copy import deepcopy
+ import os
+
+ x = np.linspace(-100.,1000,400)
+ x2 = np.linspace(-100.,1000,2400)
+
+
+ material = {'eps': 4.,
+ 'wTO': [0], 'gTO': [20], 'wLO': [598.7], 'gLO': [20],
+ }
+
+ film1 = lambin(film=[[material,10000.]])
+ film1.temperature = 100
+ #xs, spectrum = film1.calcSurfaceLoss(x)
+ xs, spectrum = film1.calcHREELS(x,normalized=False,areanormalized=False)
+ xs1, spectrum1 = film1.calcHREELS(x2,normalized=False,areanormalized=False)
+ xs2, spectrum2 = film1.calcHREELS(x,areanormalized=True)
+ xs3, spectrum3 = film1.calcHREELS(x2,areanormalized=True)
+
+ norm_test=integrate.simps(spectrum, dx=xs[2]-xs[1])/(max(xs)-min(xs))
+ norm_test1=integrate.simps(spectrum1, dx=xs1[2]-xs1[1])/(max(xs1)-min(xs1))
+ norm_test2=integrate.simps(spectrum2, dx=xs2[2]-xs2[1])/(max(xs2)-min(xs2))
+ norm_test3=integrate.simps(spectrum3, dx=xs3[2]-xs3[1])/(max(xs3)-min(xs3))
+ print(norm_test,norm_test1,norm_test2,norm_test3)
+ print(xs[2]-xs[1],xs1[2]-xs1[1],xs2[2]-xs2[1],xs3[2]-xs3[1])
+
+ plt.plot(xs[:-1],spectrum[:-1], label='normalized=Flase, '+str(len(x))+' points, E0 at '+str(max(spectrum)))
+ plt.plot(xs1[:-2],spectrum1[:-2], label='normalized=Flase, '+str(len(x2))+' points, E0 at '+str(max(spectrum1)))
+ plt.plot(xs2[:-10],spectrum2[:-10],label='area normalized=True, '+str(len(x))+' points, E0 at '+str(max(spectrum2)))
+ plt.plot(xs3[:-20],spectrum3[:-20],label='area normalized=True, '+str(len(x2))+' points, E0 at '+str(max(spectrum3)))
+
+ print('spec2/spec0=', max(spectrum2)/max(spectrum))
+ print('spec3/spec0=', max(spectrum3)/max(spectrum))
+ print('spec0/spec2=', max(spectrum)/max(spectrum2))
+ print('spec0/spec3=', max(spectrum)/max(spectrum3))
+
+ # pureDrude = {'eps': 4.,
+ # 'wTO': [-400], 'gTO': [-20], 'wLO': [0], 'gLO': [0],
+ # }
+ # film2 = lambin(film=[[pureDrude,10000.]])
+ # film2.temperature = 100
+ # xs, spectrum = film2.calcSurfaceLoss(x)
+ # plt.plot(xs,spectrum,'g-.', label='pure Drude 400')
+
+
+ plt.ylabel('Surface Loss')
plt.xlabel('Energy Loss (cm$^{-1}$)')
- plt.legend(title='Thickness in ($\AA$)')
+ plt.legend(title=r'Plasma frequency')
+
+ 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] + '.png'
+ plt.savefig(output_filename)
+
plt.show()
+
if __name__ == '__main__':
- myMain()
\ No newline at end of file
+ myMain()
diff --git a/Python code/libhreels/calcHREELS3.py b/Python code/libhreels/calcHREELS3.py
deleted file mode 100755
index c832bde..0000000
--- a/Python code/libhreels/calcHREELS3.py
+++ /dev/null
@@ -1,250 +0,0 @@
-#!/usr/bin/env python3
-import numpy as np
-import json
-import sys, re, os
-from libhreels.HREELS import myPath
-from copy import deepcopy
-
-import scipy.integrate as integrate
-
-libDir = os.path.dirname(os.path.realpath(__file__))
-
-if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')):
- print('''Make sure the Fortran routines 'myEels2' and 'myBoson'
- have been complied with the proper f2py3 for the right
- python version!!''')
-try:
- from libhreels import myEels3 as LambinEELS3 # wrapper for myEels3.f90
- from libhreels import myBoson as LambinBoson # wrapper for myBoson.f90
-except:
- print('myEels2 and MyBoson are not available here (Check your version)')
-
-# Experimental setup as dictionary:
-setup = {
- "e0": 4.0,
- "theta": 60.,
- "phia": 0.33,
- "phib": 2.0,
- "temperature": 298.,
- "debug": False
-}
-# Instrumental function describing elastic peak shape:
-instrument = {
- "width": 18.,
- "intensity": 100000.,
- "asym": 0.01,
- "gauss": 0.88
-}
-
-
-def importMaterials(string='', path=libDir):
- ''' Returns a dictionary with all phonon parameters for the material provided
- as string argument. If the string is empty or not matching, a list of all available
- materials is printed.
- '''
- file = os.path.join(myPath(path),'materials.json')
- with open(file) as json_file:
- materials = json.load(json_file)
- try:
- mat = materials[string]
- except:
- print('No data for material >>{}<< found in {} materials.json!!'.format(string, path))
- print('Available materials:\n{}\n'.format(materials.keys()))
- mat = 'None'
- return mat
-
-
-def addDrude(wPlasma, gPlasma, material):
- ''' Adds a simple Drude response to the materials properties (which are provided
- as 3rd argument) and returns a new materials dictionary with all phonon parameters. Note
- that at least the eps_infinity has to given before.
- '''
- newMaterial = deepcopy(material)
- # To select the Drude response, the Lambin Fortran code requires a negative first oscillator frequency
- # which is then interpreted as plasma frequency.
- # The values of the former LO parameter are irrelevant (but need to be provided).
- if not wPlasma > 0:
- print('Cannot add Drude to material',material)
- return material
- try:
- if len(newMaterial['wTO']) > 0:
- newMaterial['wTO'] += [-1*wPlasma]
- newMaterial['gTO'] += [-1*gPlasma]
- newMaterial['wLO'] += [0]
- newMaterial['gLO'] += [0]
- return newMaterial
- except:
- print('Cannot add Drude to material',material)
- return material
-
-################################################################################
-################################################################################
-class lambin:
- def __init__(self, film, setup=setup, instrument=instrument):
- self.e0 = setup['e0']
- self.theta = setup['theta']
- self.phia = setup['phia']
- self.phib = setup['phib']
- self.temperature = setup['temperature']
- self.debug = setup['debug']
- self.width = instrument['width']
- self.gauss = instrument['gauss']
- self.intensity = instrument['intensity']
- self.asym = instrument['asym']
- self.layers = len(film) # number of layers
- self.neps = self.layers
- name_size = self.layers
- self.name = []; self.thick=[]; self.listNOsci=[]; self.epsinf =[]; Q = []
- allTO=[]; allgTO=[]; allgLO=[]; nDrude=0
- name2 = []
- for layer in film:
- try:
- a = layer[0]['name']
- except:
- a = 'None'
- self.name.append('{:<10}'.format(a[:10])) # film name and material
- name2.append(a)
- try:
- a = layer[1]
- except:
- a = 10000.
- self.thick.append(a)
- self.epsinf.append(layer[0]['eps'])
- nTO = 2 * len(layer[0]['wTO'])
- allTO.extend(layer[0]['wTO'])
- allgTO.extend(layer[0]['gTO'])
- allTO.extend(layer[0]['wLO'])
- allgTO.extend(layer[0]['gLO'])
- Q.extend(2*len(layer[0]['wTO'])*[10.])
- self.listNOsci.append(nTO)
-
- if len(allTO)!=sum(self.listNOsci) or len(allgTO)!=sum(self.listNOsci):
- print('Error in materials: ', layer[0])
- self.wOsc = np.array(allTO)
- self.gOsc = np.array(allgTO)
- self.osc = np.array([self.wOsc, np.array(Q), self.gOsc])
- return
-
- def calcSurfaceLoss(self,x):
- ''' Calculate the surface loss spectrum for the array of x, which needs to be an equidistant array.
- All parameters are defined in the class __init__() call.'''
- wmin = min(x)
- wmax = max(x)-0.001
- dw = (wmax-wmin)/(len(x)-1) # assumes that x is an equidistant array
- wn_array_size = len(x) # size of array for x and epsilon (wn_array, loss_array)
- nper = 1.
- contrl = '{:<10}'.format('None'[:10]) # Can be 'image' to include image charge
- mode = '{:<10}'.format('kurosawa'[:10])
- wn_array,loss_array = LambinEELS3.mod_doeels.doeels(self.e0,self.theta,self.phia,self.phib,
- wmin,wmax,dw,self.layers,self.neps,nper,self.name,
- self.thick,self.epsinf,self.listNOsci,self.osc,contrl,mode,wn_array_size)
- i=0
- for item in wn_array:
- if item > 0: break
- i += 1
- return wn_array[i-1:], loss_array[i-1:]
-
- def calcHREELS(self,x, normalized=True, areanormalized=False):
- emin = min(x)
- emax = max(x)-0.001
- norm = 1
- xLoss,loss_array = self.calcSurfaceLoss(x)
- wmin = min(xLoss)
- wmax = max(xLoss)
- xOut,spectrum,n = LambinBoson.doboson3(self.temperature,self.width,self.gauss,self.asym,
- emin,emax,wmin,wmax,loss_array,self.debug,len(loss_array))
- if normalized:
- norm = max(spectrum[:n])
- #print("peak normalized, norm=", norm)
- if areanormalized: #edit by HHE
- try:
- areanormalize_xstart = np.argmin(abs(x+100.)) #seems to be oddly complicated, but is way more stable than x.index(-100.) or where()
- except:
- areanormalize_xstart = 0
- try:
- areanormalize_xend = np.argmin(abs(x-1000.))
- except:
- areanormalize_xend = len(x)
- #norm=spectrum[areanormalize_xstart:areanormalize_xend].sum()
- cropped_spectra=spectrum[areanormalize_xstart:areanormalize_xend]
- cropped_x=x[areanormalize_xstart:areanormalize_xend]
-
- #norm=cropped_spectra.sum()/len(cropped_x)*(x[areanormalize_xend]-x[areanormalize_xstart])
- #print(len(cropped_x), len(cropped_spectra), len(x), len(spectrum))
- #print(cropped_spectra.sum(), len(cropped_x),x[areanormalize_xend]-x[areanormalize_xstart])
- norm=integrate.simps(cropped_spectra, dx=x[areanormalize_xstart+1]-x[areanormalize_xstart])#/(1100/(max(cropped_x)-min(cropped_x)))
- #print("areanormalized, norm=", norm, " dx=", x[areanormalize_xstart+1]-x[areanormalize_xstart])
- else:
- print("not normalized")
- return xOut[:len(x)], spectrum[:len(x)]/norm
-
- def calcEps(self, x):
- epsArray = []
- nOsci = len(self.wOsc)
- for wn in x:
- yn = LambinEELS3.mod_doeels.seteps(self.listNOsci,nOsci,self.osc,self.epsinf,wn,self.layers)
- epsArray.append(yn)
- return np.transpose(np.array(epsArray))
-
-####################################################################################
-def myMain():
- import matplotlib.pyplot as plt
- import numpy as np
- from copy import deepcopy
- import os
-
- x = np.linspace(-100.,1000,400)
- x2 = np.linspace(-100.,1000,2400)
-
-
- material = {'eps': 4.,
- 'wTO': [0], 'gTO': [20], 'wLO': [598.7], 'gLO': [20],
- }
-
- film1 = lambin(film=[[material,10000.]])
- film1.temperature = 100
- #xs, spectrum = film1.calcSurfaceLoss(x)
- xs, spectrum = film1.calcHREELS(x,normalized=False,areanormalized=False)
- xs1, spectrum1 = film1.calcHREELS(x2,normalized=False,areanormalized=False)
- xs2, spectrum2 = film1.calcHREELS(x,areanormalized=True)
- xs3, spectrum3 = film1.calcHREELS(x2,areanormalized=True)
-
- norm_test=integrate.simps(spectrum, dx=xs[2]-xs[1])/(max(xs)-min(xs))
- norm_test1=integrate.simps(spectrum1, dx=xs1[2]-xs1[1])/(max(xs1)-min(xs1))
- norm_test2=integrate.simps(spectrum2, dx=xs2[2]-xs2[1])/(max(xs2)-min(xs2))
- norm_test3=integrate.simps(spectrum3, dx=xs3[2]-xs3[1])/(max(xs3)-min(xs3))
- print(norm_test,norm_test1,norm_test2,norm_test3)
- print(xs[2]-xs[1],xs1[2]-xs1[1],xs2[2]-xs2[1],xs3[2]-xs3[1])
-
- plt.plot(xs,spectrum, label='normalized=Flase, '+str(len(x))+' points, E0 at '+str(max(spectrum)))
- plt.plot(xs1,spectrum1, label='normalized=Flase, '+str(len(x2))+' points, E0 at '+str(max(spectrum1)))
- plt.plot(xs2,spectrum2,label='area normalized=True, '+str(len(x))+' points, E0 at '+str(max(spectrum2)))
- plt.plot(xs3,spectrum3,label='area normalized=True, '+str(len(x2))+' points, E0 at '+str(max(spectrum3)))
-
- print('spec2/spec0=', max(spectrum2)/max(spectrum))
- print('spec3/spec0=', max(spectrum3)/max(spectrum))
- print('spec0/spec2=', max(spectrum)/max(spectrum2))
- print('spec0/spec3=', max(spectrum)/max(spectrum3))
-
- # pureDrude = {'eps': 4.,
- # 'wTO': [-400], 'gTO': [-20], 'wLO': [0], 'gLO': [0],
- # }
- # film2 = lambin(film=[[pureDrude,10000.]])
- # film2.temperature = 100
- # xs, spectrum = film2.calcSurfaceLoss(x)
- # plt.plot(xs,spectrum,'g-.', label='pure Drude 400')
-
-
- plt.ylabel('Surface Loss')
- plt.xlabel('Energy Loss (cm$^{-1}$)')
- plt.legend(title=r'Plasma frequency')
-
- 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] + '.png'
- plt.savefig(output_filename)
-
- plt.show()
-
-
-if __name__ == '__main__':
- myMain()
\ No newline at end of file
diff --git a/Python code/libhreels/calcHREELS4.py b/Python code/libhreels/calcHREELS4.py
deleted file mode 100755
index ea3c233..0000000
--- a/Python code/libhreels/calcHREELS4.py
+++ /dev/null
@@ -1,250 +0,0 @@
-#!/usr/bin/env python3
-import numpy as np
-import json
-import sys, re, os
-from libhreels.HREELS import myPath
-from copy import deepcopy
-
-import scipy.integrate as integrate
-
-libDir = os.path.dirname(os.path.realpath(__file__))
-
-if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')):
- print('''Make sure the Fortran routines 'myEels2' and 'myBoson'
- have been complied with the proper f2py3 for the right
- python version!!''')
-try:
- from libhreels import myEels4 as LambinEELS4 # wrapper for myEels4.f90
- from libhreels import myBoson as LambinBoson # wrapper for myBoson.f90
-except:
- print('myEels2 and MyBoson are not available here (Check your version)')
-
-# Experimental setup as dictionary:
-setup = {
- "e0": 4.0,
- "theta": 60.,
- "phia": 0.33,
- "phib": 2.0,
- "temperature": 298.,
- "debug": False
-}
-# Instrumental function describing elastic peak shape:
-instrument = {
- "width": 18.,
- "intensity": 100000.,
- "asym": 0.01,
- "gauss": 0.88
-}
-
-
-def importMaterials(string='', path=libDir):
- ''' Returns a dictionary with all phonon parameters for the material provided
- as string argument. If the string is empty or not matching, a list of all available
- materials is printed.
- '''
- file = os.path.join(myPath(path),'materials.json')
- with open(file) as json_file:
- materials = json.load(json_file)
- try:
- mat = materials[string]
- except:
- print('No data for material >>{}<< found in {} materials.json!!'.format(string, path))
- print('Available materials:\n{}\n'.format(materials.keys()))
- mat = 'None'
- return mat
-
-
-def addDrude(wPlasma, gPlasma, material):
- ''' Adds a simple Drude response to the materials properties (which are provided
- as 3rd argument) and returns a new materials dictionary with all phonon parameters. Note
- that at least the eps_infinity has to given before.
- '''
- newMaterial = deepcopy(material)
- # To select the Drude response, the Lambin Fortran code requires a negative first oscillator frequency
- # which is then interpreted as plasma frequency.
- # The values of the former LO parameter are irrelevant (but need to be provided).
- if not wPlasma > 0:
- print('Cannot add Drude to material',material)
- return material
- try:
- if len(newMaterial['wTO']) > 0:
- newMaterial['wTO'] += [-1*wPlasma]
- newMaterial['gTO'] += [-1*gPlasma]
- newMaterial['wLO'] += [0]
- newMaterial['gLO'] += [0]
- return newMaterial
- except:
- print('Cannot add Drude to material',material)
- return material
-
-################################################################################
-################################################################################
-class lambin:
- def __init__(self, film, setup=setup, instrument=instrument):
- self.e0 = setup['e0']
- self.theta = setup['theta']
- self.phia = setup['phia']
- self.phib = setup['phib']
- self.temperature = setup['temperature']
- self.debug = setup['debug']
- self.width = instrument['width']
- self.gauss = instrument['gauss']
- self.intensity = instrument['intensity']
- self.asym = instrument['asym']
- self.layers = len(film) # number of layers
- self.neps = self.layers
- name_size = self.layers
- self.name = []; self.thick=[]; self.listNOsci=[]; self.epsinf =[]; Q = []
- allTO=[]; allgTO=[]; allgLO=[]; nDrude=0
- name2 = []
- for layer in film:
- try:
- a = layer[0]['name']
- except:
- a = 'None'
- self.name.append('{:<10}'.format(a[:10])) # film name and material
- name2.append(a)
- try:
- a = layer[1]
- except:
- a = 10000.
- self.thick.append(a)
- self.epsinf.append(layer[0]['eps'])
- nTO = 2 * len(layer[0]['wTO'])
- allTO.extend(layer[0]['wTO'])
- allgTO.extend(layer[0]['gTO'])
- allTO.extend(layer[0]['wLO'])
- allgTO.extend(layer[0]['gLO'])
- Q.extend(2*len(layer[0]['wTO'])*[10.])
- self.listNOsci.append(nTO)
-
- if len(allTO)!=sum(self.listNOsci) or len(allgTO)!=sum(self.listNOsci):
- print('Error in materials: ', layer[0])
- self.wOsc = np.array(allTO)
- self.gOsc = np.array(allgTO)
- self.osc = np.array([self.wOsc, np.array(Q), self.gOsc])
- return
-
- def calcSurfaceLoss(self,x):
- ''' Calculate the surface loss spectrum for the array of x, which needs to be an equidistant array.
- All parameters are defined in the class __init__() call.'''
- wmin = min(x)
- wmax = max(x)-0.001
- dw = (wmax-wmin)/(len(x)-1) # assumes that x is an equidistant array
- wn_array_size = len(x) # size of array for x and epsilon (wn_array, loss_array)
- nper = 1.
- contrl = '{:<10}'.format('None'[:10]) # Can be 'image' to include image charge
- mode = '{:<10}'.format('kurosawa'[:10])
- wn_array,loss_array = LambinEELS4.mod_doeels.doeels(self.e0,self.theta,self.phia,self.phib,
- wmin,wmax,dw,self.layers,self.neps,nper,self.name,
- self.thick,self.epsinf,self.listNOsci,self.osc,contrl,mode,wn_array_size)
- i=0
- for item in wn_array:
- if item > 0: break
- i += 1
- return wn_array[i-1:], loss_array[i-1:]
-
- def calcHREELS(self,x, normalized=True, areanormalized=False):
- emin = min(x)
- emax = max(x)-0.001
- norm = 1
- xLoss,loss_array = self.calcSurfaceLoss(x)
- wmin = min(xLoss)
- wmax = max(xLoss)
- xOut,spectrum,n = LambinBoson.doboson3(self.temperature,self.width,self.gauss,self.asym,
- emin,emax,wmin,wmax,loss_array,self.debug,len(loss_array))
- if normalized:
- norm = max(spectrum[:n])
- #print("peak normalized, norm=", norm)
- if areanormalized: #edit by HHE
- try:
- areanormalize_xstart = np.argmin(abs(x+100.)) #seems to be oddly complicated, but is way more stable than x.index(-100.) or where()
- except:
- areanormalize_xstart = 0
- try:
- areanormalize_xend = np.argmin(abs(x-1000.))
- except:
- areanormalize_xend = len(x)
- #norm=spectrum[areanormalize_xstart:areanormalize_xend].sum()
- cropped_spectra=spectrum[areanormalize_xstart:areanormalize_xend]
- cropped_x=x[areanormalize_xstart:areanormalize_xend]
-
- #norm=cropped_spectra.sum()/len(cropped_x)*(x[areanormalize_xend]-x[areanormalize_xstart])
- #print(len(cropped_x), len(cropped_spectra), len(x), len(spectrum))
- #print(cropped_spectra.sum(), len(cropped_x),x[areanormalize_xend]-x[areanormalize_xstart])
- norm=integrate.simps(cropped_spectra, dx=x[areanormalize_xstart+1]-x[areanormalize_xstart])#/(1100/(max(cropped_x)-min(cropped_x)))
- #print("areanormalized, norm=", norm, " dx=", x[areanormalize_xstart+1]-x[areanormalize_xstart])
- else:
- print("not normalized")
- return xOut[:len(x)], spectrum[:len(x)]/norm
-
- def calcEps(self, x):
- epsArray = []
- nOsci = len(self.wOsc)
- for wn in x:
- yn = LambinEELS4.mod_doeels.seteps(self.listNOsci,nOsci,self.osc,self.epsinf,wn,self.layers)
- epsArray.append(yn)
- return np.transpose(np.array(epsArray))
-
-####################################################################################
-def myMain():
- import matplotlib.pyplot as plt
- import numpy as np
- from copy import deepcopy
- import os
-
- x = np.linspace(-100.,1000,400)
- x2 = np.linspace(-100.,1000,2400)
-
-
- material = {'eps': 4.,
- 'wTO': [0], 'gTO': [20], 'wLO': [598.7], 'gLO': [20],
- }
-
- film1 = lambin(film=[[material,10000.]])
- film1.temperature = 100
- #xs, spectrum = film1.calcSurfaceLoss(x)
- xs, spectrum = film1.calcHREELS(x,normalized=False,areanormalized=False)
- xs1, spectrum1 = film1.calcHREELS(x2,normalized=False,areanormalized=False)
- xs2, spectrum2 = film1.calcHREELS(x,areanormalized=True)
- xs3, spectrum3 = film1.calcHREELS(x2,areanormalized=True)
-
- norm_test=integrate.simps(spectrum, dx=xs[2]-xs[1])/(max(xs)-min(xs))
- norm_test1=integrate.simps(spectrum1, dx=xs1[2]-xs1[1])/(max(xs1)-min(xs1))
- norm_test2=integrate.simps(spectrum2, dx=xs2[2]-xs2[1])/(max(xs2)-min(xs2))
- norm_test3=integrate.simps(spectrum3, dx=xs3[2]-xs3[1])/(max(xs3)-min(xs3))
- print(norm_test,norm_test1,norm_test2,norm_test3)
- print(xs[2]-xs[1],xs1[2]-xs1[1],xs2[2]-xs2[1],xs3[2]-xs3[1])
-
- plt.plot(xs[:-1],spectrum[:-1], label='normalized=Flase, '+str(len(x))+' points, E0 at '+str(max(spectrum)))
- plt.plot(xs1[:-2],spectrum1[:-2], label='normalized=Flase, '+str(len(x2))+' points, E0 at '+str(max(spectrum1)))
- plt.plot(xs2[:-10],spectrum2[:-10],label='area normalized=True, '+str(len(x))+' points, E0 at '+str(max(spectrum2)))
- plt.plot(xs3[:-20],spectrum3[:-20],label='area normalized=True, '+str(len(x2))+' points, E0 at '+str(max(spectrum3)))
-
- print('spec2/spec0=', max(spectrum2)/max(spectrum))
- print('spec3/spec0=', max(spectrum3)/max(spectrum))
- print('spec0/spec2=', max(spectrum)/max(spectrum2))
- print('spec0/spec3=', max(spectrum)/max(spectrum3))
-
- # pureDrude = {'eps': 4.,
- # 'wTO': [-400], 'gTO': [-20], 'wLO': [0], 'gLO': [0],
- # }
- # film2 = lambin(film=[[pureDrude,10000.]])
- # film2.temperature = 100
- # xs, spectrum = film2.calcSurfaceLoss(x)
- # plt.plot(xs,spectrum,'g-.', label='pure Drude 400')
-
-
- plt.ylabel('Surface Loss')
- plt.xlabel('Energy Loss (cm$^{-1}$)')
- plt.legend(title=r'Plasma frequency')
-
- 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] + '.png'
- plt.savefig(output_filename)
-
- plt.show()
-
-
-if __name__ == '__main__':
- myMain()
\ No newline at end of file
diff --git a/Python code/libhreels/dielectrics.py b/Python code/libhreels/dielectrics.py
index ab35bca..6438a9f 100755
--- a/Python code/libhreels/dielectrics.py
+++ b/Python code/libhreels/dielectrics.py
@@ -2,16 +2,26 @@
import matplotlib.pyplot as plt
import numpy as np
from numpy import real, imag, sqrt
+from scipy import constants
-#from scipy.constants import physical_constants
-# Thz = 100*physical_constants['speed of light in vacuum'][0] # conversion from cm^-1 to Hz
+def doping2chargecarrierdensity(doping): #use doping in percentage by mass, not volume
+ '''Returns the charge carrier density for a given doping in percentage of mass.'''
+ return float(doping*3.31664067935517E+020*1E6) # doping[%] * charge carrier/cm^3 * 1E6 = n_e [m^-3]
-def doping2plasmaFrequency(doping,epsInfinity=1.):
- '''Returns the bulk plasma frequency in cm-1 for a doping
- given by the argument in cm-3.
- Check if eps_Infinity handling is correct.'''
- return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(epsInfinity)
+def doping2plasmaFrequency(doping,epsInfinity=1.): #use doping in percentage by mass, not volume
+ '''Returns the plasma frequency for a given doping in percentage of mass.'''
+ eps_0 = constants.value("vacuum electric permittivity")
+ e = constants.value('elementary charge')
+ m_e = constants.value('electron mass')
+ c = constants.value('speed of light in vacuum')
+ eps_r = epsInfinity
+ n_e = doping2chargecarrierdensity(doping)
+
+ w_P = np.sqrt(n_e*e*e/(eps_r*eps_0*m_e)) #Hz
+ w_P = w_P /1E+12 * 33.35641 #Hz -> THz -> cm^-1
+ return w_P
+
def doping2surfacePlasma(doping,epsInfinity=1.):
'''Returns the surface plasma frequency in cm-1 for a doping
@@ -19,140 +29,217 @@ def doping2surfacePlasma(doping,epsInfinity=1.):
return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(1+epsInfinity)
def loss(eps):
- return np.imag(1/eps)
+ '''Returns the loss function for a given eps'''
+ return np.imag(eps) #HHE changed former sign
def surfaceLoss(eps):
- return np.imag(1/(1+eps))
+ '''Returns the surface loss function for a given eps'''
+ return np.imag(-1/(1+eps)) #HHE changed former sign, according to Phd of FSc
def reflectivity(eps): # IR reflectivity
- #print("eps in reflect: ",eps[1])
+ '''Returns the reflectivity for a given eps'''
sq = np.sqrt(eps)
- a =np.abs((sq-1)/(sq+1))
- #print("a in reflect: ",a[1])
+ a = np.abs((sq-1)/(sq+1))
return a*a
-def sigma(eps,w,eps_Infinity=1): # Complex conductivity
- return np.imag((eps-eps_Infinity)*w)
+def sigma(eps,w):
+ '''Returns the complex conductivity for a given eps'''
+ return (eps-1)*w/1j
-def plotDielectrics(x,eps,titlee=" ", flag_plot_show=True):
- fig, axs = plt.subplots(3, 1, sharex=True, figsize=(6,6))
- axs[0].plot(x, -np.imag(eps), label='-Im $\epsilon (\omega )$')
- axs[0].plot(x, np.real(eps), label='Re $\epsilon (\omega )$')
+def plotDielectrics(x,eps, title=" ", plot_show=True):
+ '''This method plots a given (x,eps)-dataset as Re/Im(eps), SurfaceLoss and Reflectivity
+ '''
+ fig, axs = plt.subplots(3, 1, sharex=True)
+ fig.suptitle(str(title), fontsize=16)
+ axs[0].plot(x, np.imag(eps), label='Im( $\epsilon (\omega )$ )')
+ axs[0].plot(x, np.real(eps), label='Re( $\epsilon (\omega )$ )')
axs[0].legend()
axs[0].set_ylabel('Dielectric Function')
- axs[1].plot(x, surfaceLoss(eps), linestyle='-')
+ axs[1].plot(x, surfaceLoss(eps), lineStyle='-', label='Loss function')
axs[1].set_ylabel('Surface Loss Function')
- axs[1].set_ylim([0,None])
+ axs[1].set_ylim([-0.10,None])
axs[2].plot(x, reflectivity(eps))
axs[2].set_ylabel('Reflectivity')
axs[2].set_xlabel('Frequency')
- axs[2].set_ylim([0,1])
+ axs[2].set_ylim([0.0,1])
axs[0].set_xlim(left=0, right=max(x))
- fig.suptitle(str(titlee), fontsize=16)
- if flag_plot_show == True:
+ if plot_show==True:
plt.show()
+
class oscillator:
- #print("class oscillator is called") #added by HHe
- def __init__(self, wTO, wLO, gammaTO=20, gammaLO=20):
+ '''This defines epsilon of an oscillator with a given (wTO, gTO, wLO, gLO)
+ '''
+ def __init__(self, wTO=177., gTO=20., wLO=184., gLO=20.):
self.wTO = wTO
self.wLO = wLO
- self.gammaTO = gammaTO
- self.gammaLO = gammaLO
-
- #print("wTO, wLO, gTO, gLO",wTO,wLO,gammaTO,gammaLO) #added by HHe
+ self.gammaTO = gTO
+ self.gammaLO = gLO
def eps(self,w):
- nom = self.wLO*self.wLO - w*w + 1j*self.gammaLO*w
- denom = self.wTO*self.wTO - w*w + 1j*self.gammaTO*w
+ nom = self.wLO*self.wLO - w*w - 1j*w*self.gammaLO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa
+ denom = self.wTO*self.wTO - w*w - 1j*w*self.gammaTO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa
return nom/denom
def __call__(self, w):
return self.eps(w)
-class simpleOsci:#(oscillator): #added by HHe
- #print("class simpleOsci is called") #added by HHe
- def __init__(self, wTO, Q, gammaTO=20):
- self.wTO = wTO
- self.Q = Q
- self.gammaTO = gammaTO
+class drude:
+ '''This defines an additive Drude component to a dielectric function:
+ eps_drude = drude2(omega_p, gamma_p, gamma_0)
+ ... eps(omega) + eps_drude(omega)
+ '''
+ def __init__(self, omega_p, gamma_p, gamma_0):
- print("wTO, Q, gTO",wTO,Q,gammaTO) #added by HHe
+ self.omega_p = omega_p
+ self.gamma_p = gamma_p
+ self.gamma_0 = gamma_0
+
+ print("called drude with w_P, g_P, g_0:", omega_p, gamma_p, gamma_0)
+
def eps(self,w):
- nom = self.wTO*self.wTO*self.Q
- denom = self.wTO*self.wTO - w*w + 1j*self.gammaTO*w
- return nom/denom
+ newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
+ w = newW
+ #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
- def __call__(self, w): #added by HHe
- return self.eps(w)
+ nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p)
+ denom = w*w + 1j*w*self.gamma_0 #sign changed according to Lambin notation
+ return -(nom/denom) #HHE changed
-class drude:
- '''Returns the additive Drude response with plasma frequency and damping as parameters. Note that for the
- full dielctric response, this Drude contribution has to be added to eps_infinity and any phonon
- contribution.
+ def __call__(self, w):
+ return self.eps(w)
+
+
+class cole_davidson:
+ '''This defines an additive Cole-Davidson component to a dielectric function:
+ eps_cole_davidson = cole_davidson(omega_p, gamma_p, gamma_0)
+ ... eps(omega) + eps_drude(omega)
'''
- def __init__(self, wPL, gamma):
- self.wPL = wPL
- self.gamma = gamma
+ def __init__(self, omega_p, gamma_p, gamma_0, beta):
- def eps(self,w):
- # Starting from libhreels version 1.1.4, it has been corrected:
- return -(self.wPL*self.wPL)/(w*w - 1j*self.gamma*w)
+ self.omega_p = omega_p
+ self.gamma_p = gamma_p
+ self.gamma_0 = gamma_0
+ self.beta = beta
+
+ print("called coledavidson with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
+ def eps(self,w):
+ newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
+ w = newW
+ #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
+ #return -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)) #HHE changed
+
+ #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
+ #denom = (w*w - 1j*w*self.gamma_0)**self.beta #sign changed according to gervais paper, formular 7
+
+ hilf = 1j * w * self.gamma_0 * (1-1j*w*(1/self.gamma_0))**(self.beta)
+
+ nom = self.omega_p**2 - 1j*w*(self.gamma_p) - w*w + hilf #unfortunately also the nominator changes due to cole-davidson
+ denom = hilf #sign changed according to gervais paper, formular 7 -> did not match, changed back to lambin formularism
+
+
+ return -(nom/denom)
+
def __call__(self, w):
return self.eps(w)
-class drude2:
- '''Returns the additive extended Drude response with three arguments: plasma frequency and two damping
- parameters. The first damping parameter is the plasma damping and the second the damping at frequency
- zero. Note that for the full dielectric response, this Drude contribution has to be added to eps_infinity
- and any phonon contribution.
-
- eps = drude2(omega_p, gamma_p, gamma_0)
+class cole_cole:
+ '''This defines an additive Cole-Davidson component to a dielectric function:
+ eps_cole_cole = cole_cole(omega_p, gamma_p, gamma_0, beta)
+ ...
'''
- def __init__(self, omega_p, gamma_p, gamma_0=None):
+ def __init__(self, omega_p, gamma_p, gamma_0, beta):
+
self.omega_p = omega_p
self.gamma_p = gamma_p
- if gamma_0:
- self.gamma_0 = gamma_0
- else:
- self.gamma_0 = gamma_p
+ self.gamma_0 = gamma_0
+ self.beta = beta
+
+ print("called colecole with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
def eps(self,w):
newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
w = newW
- return -((self.omega_p**2 + 1j*(self.gamma_p-self.gamma_0)*w)/(w*(w - 1j *self.gamma_0)))
+ #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
+ #denom = w*w + w*(1j*self.gamma_0)**self.beta
+
+ hilf = 1j**(self.beta+1) * w**(self.beta+1) * self.gamma_0**(1-self.beta)
+
+ nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) - w*w - hilf #unfortunately also the nominator changes due to cole-cole
+ denom = 1j*w*(self.gamma_0) - hilf #sign changed according to lambin formularism
+
+ return -(nom/denom)
def __call__(self, w):
return self.eps(w)
+
+class simpleOscillator(oscillator):
+ '''This defines epsilon of an oscillator with a given (wPl, Q and gTO)
+ '''
+ def __init__(self, wPl, Q, gTO=20):
+ self.wPl = wPl
+ self.Q = Q
+ self.gammaTO = gTO
+
+ def eps(self,w):
+ nom = self.wPl*self.wPl*self.Q
+ denom = self.wPl*self.wPl - w*w - 1j*self.gammaTO*w #HHE changed, wrong sign of Im
+ return nom/denom
+
+class simpleDrude:
+ '''This defines epsilon of an oscillator with a given (wPl, gamma). gamma = gammaTO = gammaLO
+ '''
+ def __init__(self, wPL, gamma):
+ self.wPL = wPL
+ self.gamma = gamma
+
+ def eps(self,w):
+ return -(self.wPL*self.wPL)/(w*w-1j*w*self.gamma) #HHE changed
+
+ def __call__(self, w):
+ return self.eps(w)
+
+
+
def myMain():
# BaTiO3
- data = [
+ simple_data = [
[178.6, 8.02, 3.09800E-02],
[270.6, 24.00, 10.0800E-02],
[522.9, 1.20, 6.97200E-02]]
- # Width (3rd parameter) is given in fraction of TO frequency
- oscis = [simpleOsci(TO, Q, f*TO) for (TO, Q, f) in data]
+ data = [
+ [177.4 , 1.9262 , 184.003 , 9.718],
+ [272.84 , 93.3091 , 470.972 , 14.31 ],
+ [506.379 , 42.60 , 739.651 , 33.044],
+ ]
- x = np.linspace(5,1000,num=1200)
+ x = np.linspace(0,1000,num=1200)
epsInfinity = 5.25
+
+ simpleoscis = [simpleOscillator(TO, Q, f*TO) for (TO, Q, f) in simple_data]
eps = epsInfinity
- for each in oscis:
- eps += each(x)
+ for each in simpleoscis:
+ eps += each(x) #sum
+ plotDielectrics(x, eps, title="Eps of BaTiO calced with simpleOscillator", plot_show=False)
- # Now add Drude contribution:
- eps += drude2(1200,60)(x)
- plotDielectrics(x, eps)
+ oscis = [oscillator(wTO=a, wLO=c, gTO=b, gLO=d) for (a, b, c, d) in data]
+ eps2 = epsInfinity
+ for each in oscis:
+ eps2 *= each(x) #kurosawa
+ plotDielectrics(x, eps2, title="Eps of BaTiO calced with oscillator", plot_show=False)
- print("hallo:")
- print(doping2plasmaFrequency(1E19,epsInfinity=5.2))
+ plt.show()
if __name__ == '__main__':
- myMain()
+ print(doping2plasmaFrequency(0.05,3.024434))
+ print(doping2plasmaFrequency(0.1,4.57655596))
+
+ myMain()
diff --git a/Python code/libhreels/dielectrics2.py b/Python code/libhreels/dielectrics2.py
deleted file mode 100755
index 6438a9f..0000000
--- a/Python code/libhreels/dielectrics2.py
+++ /dev/null
@@ -1,245 +0,0 @@
-# Testing dielectric models and their surface loss function
-import matplotlib.pyplot as plt
-import numpy as np
-from numpy import real, imag, sqrt
-from scipy import constants
-
-
-def doping2chargecarrierdensity(doping): #use doping in percentage by mass, not volume
- '''Returns the charge carrier density for a given doping in percentage of mass.'''
- return float(doping*3.31664067935517E+020*1E6) # doping[%] * charge carrier/cm^3 * 1E6 = n_e [m^-3]
-
-def doping2plasmaFrequency(doping,epsInfinity=1.): #use doping in percentage by mass, not volume
- '''Returns the plasma frequency for a given doping in percentage of mass.'''
- eps_0 = constants.value("vacuum electric permittivity")
- e = constants.value('elementary charge')
- m_e = constants.value('electron mass')
- c = constants.value('speed of light in vacuum')
- eps_r = epsInfinity
- n_e = doping2chargecarrierdensity(doping)
-
- w_P = np.sqrt(n_e*e*e/(eps_r*eps_0*m_e)) #Hz
- w_P = w_P /1E+12 * 33.35641 #Hz -> THz -> cm^-1
- return w_P
-
-
-def doping2surfacePlasma(doping,epsInfinity=1.):
- '''Returns the surface plasma frequency in cm-1 for a doping
- given by the argument in cm-3'''
- return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(1+epsInfinity)
-
-def loss(eps):
- '''Returns the loss function for a given eps'''
- return np.imag(eps) #HHE changed former sign
-
-def surfaceLoss(eps):
- '''Returns the surface loss function for a given eps'''
- return np.imag(-1/(1+eps)) #HHE changed former sign, according to Phd of FSc
-
-def reflectivity(eps): # IR reflectivity
- '''Returns the reflectivity for a given eps'''
- sq = np.sqrt(eps)
- a = np.abs((sq-1)/(sq+1))
- return a*a
-
-def sigma(eps,w):
- '''Returns the complex conductivity for a given eps'''
- return (eps-1)*w/1j
-
-def plotDielectrics(x,eps, title=" ", plot_show=True):
- '''This method plots a given (x,eps)-dataset as Re/Im(eps), SurfaceLoss and Reflectivity
- '''
- fig, axs = plt.subplots(3, 1, sharex=True)
- fig.suptitle(str(title), fontsize=16)
- axs[0].plot(x, np.imag(eps), label='Im( $\epsilon (\omega )$ )')
- axs[0].plot(x, np.real(eps), label='Re( $\epsilon (\omega )$ )')
- axs[0].legend()
- axs[0].set_ylabel('Dielectric Function')
-
- axs[1].plot(x, surfaceLoss(eps), lineStyle='-', label='Loss function')
- axs[1].set_ylabel('Surface Loss Function')
- axs[1].set_ylim([-0.10,None])
-
- axs[2].plot(x, reflectivity(eps))
- axs[2].set_ylabel('Reflectivity')
- axs[2].set_xlabel('Frequency')
- axs[2].set_ylim([0.0,1])
- axs[0].set_xlim(left=0, right=max(x))
-
- if plot_show==True:
- plt.show()
-
-class oscillator:
- '''This defines epsilon of an oscillator with a given (wTO, gTO, wLO, gLO)
- '''
- def __init__(self, wTO=177., gTO=20., wLO=184., gLO=20.):
- self.wTO = wTO
- self.wLO = wLO
- self.gammaTO = gTO
- self.gammaLO = gLO
-
- def eps(self,w):
- nom = self.wLO*self.wLO - w*w - 1j*w*self.gammaLO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa
- denom = self.wTO*self.wTO - w*w - 1j*w*self.gammaTO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa
- return nom/denom
-
- def __call__(self, w):
- return self.eps(w)
-
-class drude:
- '''This defines an additive Drude component to a dielectric function:
- eps_drude = drude2(omega_p, gamma_p, gamma_0)
- ... eps(omega) + eps_drude(omega)
- '''
- def __init__(self, omega_p, gamma_p, gamma_0):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
-
- print("called drude with w_P, g_P, g_0:", omega_p, gamma_p, gamma_0)
-
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p)
- denom = w*w + 1j*w*self.gamma_0 #sign changed according to Lambin notation
-
- return -(nom/denom) #HHE changed
-
- def __call__(self, w):
- return self.eps(w)
-
-
-class cole_davidson:
- '''This defines an additive Cole-Davidson component to a dielectric function:
- eps_cole_davidson = cole_davidson(omega_p, gamma_p, gamma_0)
- ... eps(omega) + eps_drude(omega)
- '''
- def __init__(self, omega_p, gamma_p, gamma_0, beta):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
- self.beta = beta
-
- print("called coledavidson with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
-
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
- #return -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)) #HHE changed
-
- #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
- #denom = (w*w - 1j*w*self.gamma_0)**self.beta #sign changed according to gervais paper, formular 7
-
- hilf = 1j * w * self.gamma_0 * (1-1j*w*(1/self.gamma_0))**(self.beta)
-
- nom = self.omega_p**2 - 1j*w*(self.gamma_p) - w*w + hilf #unfortunately also the nominator changes due to cole-davidson
- denom = hilf #sign changed according to gervais paper, formular 7 -> did not match, changed back to lambin formularism
-
-
- return -(nom/denom)
-
- def __call__(self, w):
- return self.eps(w)
-
-class cole_cole:
- '''This defines an additive Cole-Davidson component to a dielectric function:
- eps_cole_cole = cole_cole(omega_p, gamma_p, gamma_0, beta)
- ...
- '''
- def __init__(self, omega_p, gamma_p, gamma_0, beta):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
- self.beta = beta
-
- print("called colecole with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
- #denom = w*w + w*(1j*self.gamma_0)**self.beta
-
- hilf = 1j**(self.beta+1) * w**(self.beta+1) * self.gamma_0**(1-self.beta)
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) - w*w - hilf #unfortunately also the nominator changes due to cole-cole
- denom = 1j*w*(self.gamma_0) - hilf #sign changed according to lambin formularism
-
- return -(nom/denom)
-
- def __call__(self, w):
- return self.eps(w)
-
-
-class simpleOscillator(oscillator):
- '''This defines epsilon of an oscillator with a given (wPl, Q and gTO)
- '''
- def __init__(self, wPl, Q, gTO=20):
- self.wPl = wPl
- self.Q = Q
- self.gammaTO = gTO
-
- def eps(self,w):
- nom = self.wPl*self.wPl*self.Q
- denom = self.wPl*self.wPl - w*w - 1j*self.gammaTO*w #HHE changed, wrong sign of Im
- return nom/denom
-
-class simpleDrude:
- '''This defines epsilon of an oscillator with a given (wPl, gamma). gamma = gammaTO = gammaLO
- '''
- def __init__(self, wPL, gamma):
- self.wPL = wPL
- self.gamma = gamma
-
- def eps(self,w):
- return -(self.wPL*self.wPL)/(w*w-1j*w*self.gamma) #HHE changed
-
- def __call__(self, w):
- return self.eps(w)
-
-
-
-def myMain():
- # BaTiO3
- simple_data = [
- [178.6, 8.02, 3.09800E-02],
- [270.6, 24.00, 10.0800E-02],
- [522.9, 1.20, 6.97200E-02]]
-
- data = [
- [177.4 , 1.9262 , 184.003 , 9.718],
- [272.84 , 93.3091 , 470.972 , 14.31 ],
- [506.379 , 42.60 , 739.651 , 33.044],
- ]
-
- x = np.linspace(0,1000,num=1200)
- epsInfinity = 5.25
-
- simpleoscis = [simpleOscillator(TO, Q, f*TO) for (TO, Q, f) in simple_data]
- eps = epsInfinity
- for each in simpleoscis:
- eps += each(x) #sum
- plotDielectrics(x, eps, title="Eps of BaTiO calced with simpleOscillator", plot_show=False)
-
- oscis = [oscillator(wTO=a, wLO=c, gTO=b, gLO=d) for (a, b, c, d) in data]
- eps2 = epsInfinity
- for each in oscis:
- eps2 *= each(x) #kurosawa
- plotDielectrics(x, eps2, title="Eps of BaTiO calced with oscillator", plot_show=False)
-
- plt.show()
-
-if __name__ == '__main__':
- print(doping2plasmaFrequency(0.05,3.024434))
- print(doping2plasmaFrequency(0.1,4.57655596))
-
- myMain()
diff --git a/Python code/libhreels/dielectrics2_save.py b/Python code/libhreels/dielectrics2_save.py
deleted file mode 100755
index e4dabbd..0000000
--- a/Python code/libhreels/dielectrics2_save.py
+++ /dev/null
@@ -1,250 +0,0 @@
-# Testing dielectric models and their surface loss function
-import matplotlib.pyplot as plt
-import numpy as np
-from numpy import real, imag, sqrt
-from scipy import constants
-
-
-def doping2chargecarrierdensity(doping): #use doping in percentage by mass, not volume
- '''Returns the charge carrier density for a given doping in percentage of mass.'''
- return float(doping*3.31664067935517E+020*1E6) # doping[%] * charge carrier/cm^3 * 1E6 = n_e [m^-3]
-
-def doping2plasmaFrequency(doping,epsInfinity=1.): #use doping in percentage by mass, not volume
- '''Returns the plasma frequency for a given doping in percentage of mass.'''
- eps_0 = constants.value("vacuum electric permittivity")
- e = constants.value('elementary charge')
- m_e = constants.value('electron mass')
- c = constants.value('speed of light in vacuum')
- eps_r = epsInfinity
- n_e = doping2chargecarrierdensity(doping)
-
- w_P = np.sqrt(n_e*e*e/(eps_r*eps_0*m_e)) #Hz
- w_P = w_P /1E+12 * 33.35641 #Hz -> THz -> cm^-1
- return w_P
-
-
-def doping2surfacePlasma(doping,epsInfinity=1.):
- '''Returns the surface plasma frequency in cm-1 for a doping
- given by the argument in cm-3'''
- return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(1+epsInfinity)
-
-def loss(eps):
- '''Returns the loss function for a given eps'''
- return np.imag(eps) #HHE changed former sign
-
-def surfaceLoss(eps):
- '''Returns the surface loss function for a given eps'''
- return np.imag(-1/(1+eps)) #HHE changed former sign, according to Phd of FSc
-
-def reflectivity(eps): # IR reflectivity
- '''Returns the reflectivity for a given eps'''
- sq = np.sqrt(eps)
- a = np.abs((sq-1)/(sq+1))
- return a*a
-
-def sigma(eps,w):
- '''Returns the complex conductivity for a given eps'''
- return (eps-1)*w/1j
-
-def plotDielectrics(x,eps, title=" ", plot_show=True):
- '''This method plots a given (x,eps)-dataset as Re/Im(eps), SurfaceLoss and Reflectivity
- '''
- fig, axs = plt.subplots(3, 1, sharex=True)
- fig.suptitle(str(title), fontsize=16)
- axs[0].plot(x, np.imag(eps), label='Im( $\epsilon (\omega )$ )')
- axs[0].plot(x, np.real(eps), label='Re( $\epsilon (\omega )$ )')
- axs[0].legend()
- axs[0].set_ylabel('Dielectric Function')
-
- axs[1].plot(x, surfaceLoss(eps), lineStyle='-', label='Loss function')
- axs[1].set_ylabel('Surface Loss Function')
- axs[1].set_ylim([-0.10,None])
-
- axs[2].plot(x, reflectivity(eps))
- axs[2].set_ylabel('Reflectivity')
- axs[2].set_xlabel('Frequency')
- axs[2].set_ylim([0.0,1])
- axs[0].set_xlim(left=0, right=max(x))
-
- if plot_show==True:
- plt.show()
-
-class oscillator:
- '''This defines epsilon of an oscillator with a given (wTO, gTO, wLO, gLO)
- '''
- def __init__(self, wTO=177., gTO=20., wLO=184., gLO=20.):
- self.wTO = wTO
- self.wLO = wLO
- self.gammaTO = gTO
- self.gammaLO = gLO
-
- def eps(self,w):
- nom = self.wLO*self.wLO - w*w - 1j*w*self.gammaLO #HHE changed, wrong sign of Im
- denom = self.wTO*self.wTO - w*w - 1j*w*self.gammaTO #HHE changed, wrong sign of Im
- return nom/denom
-
- def __call__(self, w):
- return self.eps(w)
-
-class drude:
- '''This defines an additive Drude component to a dielectric function:
- eps_drude = drude2(omega_p, gamma_p, gamma_0)
- ... eps(omega) + eps_drude(omega)
- '''
- def __init__(self, omega_p, gamma_p, gamma_0):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
-
- print("called drude with w_P, g_P, g_0:", omega_p, gamma_p, gamma_0)
-
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_p-self.gamma_0)
- denom = w*w - 1j*w*self.gamma_0 #sign changed according to gervais paper, formular 7
-
- return -(nom/denom) #HHE changed
-
- def __call__(self, w):
- return self.eps(w)
-
-
-class cole_davidson:
- '''This defines an additive Cole-Davidson component to a dielectric function:
- eps_cole_davidson = cole_davidson(omega_p, gamma_p, gamma_0)
- ... eps(omega) + eps_drude(omega)
- '''
- def __init__(self, omega_p, gamma_p, gamma_0, beta):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
- self.beta = beta
-
- print("called coledavidson with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
-
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)))
- #return -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)) #HHE changed
-
- #nom = self.omega_p**2 * self.gamma_p**(1-self.beta)
- #denom = 1j*w*(self.gamma_p-1j*w)**self.beta
- #return -(nom/denom)
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p)**(1-self.beta) #just a try to recreate the drude with two gammas -> works only for gamma-0=gamma_p
- denom = 1j*w*(self.gamma_0-1j*w)**self.beta
- #return -(nom/denom)
-
- #nom = self.omega_p**2 + 1j*w*(self.gamma_0**(1/self.beta)+self.gamma_p**(1/self.beta))**(self.beta) #seems to be correct analytically, but does not work at all
- #denom = 1j*w*(self.gamma_0**(1/self.beta)+self.gamma_p**(1/self.beta))**self.beta
-
- nom = self.omega_p**2 + +1j*w*self.gamma_p*(1-1j*w*self.gamma_p)**(1-self.beta) #lacking second gamma again, but works ok
- denom = 1j*w*self.gamma_p*(1-1j*w*self.gamma_p)**(1-self.beta)
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
- denom = (w*w - 1j*w*self.gamma_0)**self.beta #sign changed according to gervais paper, formular 7
- return -(nom/denom)
-
- def __call__(self, w):
- return self.eps(w)
-
-class cole_cole:
- '''This defines an additive Cole-Davidson component to a dielectric function:
- eps_cole_cole = cole_cole(omega_p, gamma_p, gamma_0, beta)
- ...
- '''
- def __init__(self, omega_p, gamma_p, gamma_0, beta):
-
- self.omega_p = omega_p
- self.gamma_p = gamma_p
- self.gamma_0 = gamma_0
- self.beta = beta
-
- print("called colecole with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta)
-
- def eps(self,w):
- newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero
- w = newW
- #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P)
- #denom = w*w + w*(1j*self.gamma_0)**self.beta
-
- nom = self.omega_p**2 + 1j*w*(self.gamma_p-self.gamma_0) - w*w + w**(1+self.beta)*self.gamma_0**(1-self.beta) #unfortunately also the nominator changes due to cole-cole
- denom = w**(1+self.beta)*self.gamma_0**(1-self.beta) - 1j*w*self.gamma_0 #sign changed according to gervais paper, formular 7
-
- return -(nom/denom)
-
- def __call__(self, w):
- return self.eps(w)
-
-
-class simpleOscillator(oscillator):
- '''This defines epsilon of an oscillator with a given (wPl, Q and gTO)
- '''
- def __init__(self, wPl, Q, gTO=20):
- self.wPl = wPl
- self.Q = Q
- self.gammaTO = gTO
-
- def eps(self,w):
- nom = self.wPl*self.wPl*self.Q
- denom = self.wPl*self.wPl - w*w - 1j*self.gammaTO*w #HHE changed, wrong sign of Im
- return nom/denom
-
-class simpleDrude:
- '''This defines epsilon of an oscillator with a given (wPl, gamma). gamma = gammaTO = gammaLO
- '''
- def __init__(self, wPL, gamma):
- self.wPL = wPL
- self.gamma = gamma
-
- def eps(self,w):
- return -(self.wPL*self.wPL)/(w*w-1j*w*self.gamma) #HHE changed
-
- def __call__(self, w):
- return self.eps(w)
-
-
-
-def myMain():
- # BaTiO3
- simple_data = [
- [178.6, 8.02, 3.09800E-02],
- [270.6, 24.00, 10.0800E-02],
- [522.9, 1.20, 6.97200E-02]]
-
- data = [
- [177.4 , 1.9262 , 184.003 , 9.718],
- [272.84 , 93.3091 , 470.972 , 14.31 ],
- [506.379 , 42.60 , 739.651 , 33.044],
- ]
-
- x = np.linspace(0,1000,num=1200)
- epsInfinity = 5.25
-
- simpleoscis = [simpleOscillator(TO, Q, f*TO) for (TO, Q, f) in simple_data]
- eps = epsInfinity
- for each in simpleoscis:
- eps += each(x) #sum
- plotDielectrics(x, eps, title="Eps of BaTiO calced with simpleOscillator", plot_show=False)
-
- oscis = [oscillator(wTO=a, wLO=c, gTO=b, gLO=d) for (a, b, c, d) in data]
- eps2 = epsInfinity
- for each in oscis:
- eps2 *= each(x) #kurosawa
- plotDielectrics(x, eps2, title="Eps of BaTiO calced with oscillator", plot_show=False)
-
- plt.show()
-
-if __name__ == '__main__':
- print(doping2plasmaFrequency(0.05,3.024434))
- print(doping2plasmaFrequency(0.1,4.57655596))
-
- myMain()
diff --git a/Python code/nice_fitting_code.py b/Python code/nice_fitting_code.py
index c52ca8e..cf5886b 100755
--- a/Python code/nice_fitting_code.py
+++ b/Python code/nice_fitting_code.py
@@ -2,7 +2,7 @@
from os import truncate
import matplotlib.pyplot as plt
import numpy as np
-from libhreels import calcHREELS4 as ch4
+from libhreels import calcHREELS as ch
from libhreels.HREELS import HREELS
from copy import deepcopy
import lmfit
@@ -58,7 +58,7 @@ def lossFunction(x, **par):
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 = ch.lambin(film=extended_stack, setup=setup)
film1.gauss = par['gauss']
film1.width = par['width']
if FlagAsymSetOneMinusGauss==True:
--
GitLab