diff --git a/Examples/KTaO3.png b/Examples/KTaO3.png
index 743b6b5ef49d3800a592bad8f463a2ed791e3b30..1f6f459950d8d9efdeb04d22d45bb6c78a73d9fc 100644
Binary files a/Examples/KTaO3.png and b/Examples/KTaO3.png differ
diff --git a/Examples/KTaO3.py b/Examples/KTaO3.py
index 1928b2855f647b0a6b01090a279735780981a7be..08758eadba4940ec36d4efbd6daebe4a2bf6dcb6 100644
--- a/Examples/KTaO3.py
+++ b/Examples/KTaO3.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
-from libhreels.calcHREELS import lambin, importMaterials
+from libhreels.calcHREELS20 import lambin, importMaterials
import os
###### Define Strontium titanate film on Silver ####
diff --git a/Examples/calcHREELS1.png b/Examples/calcHREELS1.png
index f564843722d10c1a8498cbfd8d716208aab91fe4..f2ebeafe6c39b6b5c8b632a71cfd96f930b0f5d0 100644
Binary files a/Examples/calcHREELS1.png and b/Examples/calcHREELS1.png differ
diff --git a/Examples/calcHREELS1.py b/Examples/calcHREELS1.py
index 246afad550fd4438ecc6568efdcfdf55da6aadbd..46639addb4ee39d8be495ab3e3678e483e375a5a 100644
--- a/Examples/calcHREELS1.py
+++ b/Examples/calcHREELS1.py
@@ -1,24 +1,30 @@
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
-from libhreels import calcHREELS
+from libhreels import calcHREELS20
import os
###### Define Strontium titanate film on Silver ####
## (Thickness is given in Angstroem as second parm)#
-thinFilm = [[calcHREELS.importMaterials('STO'),200.],
- [calcHREELS.importMaterials('Ag'),1000.]]
-
-d = calcHREELS.lambin(thinFilm)
+thinFilm = [[{
+ "eps": 5.25,
+ "wTO": [393.7],
+ "gTO": [10.8],
+ "wLO": [584.7],
+ "gLO": [10.8],
+ "reference": "Schumann'22",
+ "name": "NiO"
+ }, 300000.]]
+d = calcHREELS20.lambin(thinFilm)
x = np.linspace(-100,1000,1100)
x,y = d.calcHREELS(x)
###### Plotting ###########################
-plt.plot(x,y, label=calcHREELS.importMaterials('STO')['reference'])
+plt.plot(x,y, label='NiO(001)')
plt.ylabel('Relative HREELS signal')
plt.xlabel('Energy Loss (cm$^{-1}$)')
-plt.title(r'SrTiO$_3$')
+plt.title(r"F.O.Schumann'22")
plt.ylim(bottom=0.)
plt.legend()
diff --git a/Examples/calcHREELS2.py b/Examples/calcHREELS2.py
index 466b1e51676719af7d5dcb600a67e5ca899923ea..871bd1ae198b8ba24e6b478a61f2a2af334f991d 100644
--- a/Examples/calcHREELS2.py
+++ b/Examples/calcHREELS2.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
-from libhreels.calcHREELS import lambin, importMaterials
+from libhreels.calcHREELS20 import lambin, importMaterials
###### Define Barium titanate film on Silver ####
## (Thickness is given in Angstroem as second parm)#
diff --git a/Examples/calcHREELS3.png b/Examples/calcHREELS3.png
index 432238c38578ff6c80451899e7f646c20df504b4..5329753edf9e8812a30768030d2a0260a8165b82 100644
Binary files a/Examples/calcHREELS3.png and b/Examples/calcHREELS3.png differ
diff --git a/Examples/calcHREELS3.py b/Examples/calcHREELS3.py
index eb6310c150e3f347e6e11044b7c0dafd2a5eb773..7d3a62093b03b924eb119a1853c1efb7734b07b0 100644
--- a/Examples/calcHREELS3.py
+++ b/Examples/calcHREELS3.py
@@ -1,7 +1,7 @@
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
-from libhreels.calcHREELS import lambin, importMaterials
+from libhreels.calcHREELS20 import lambin
import os
# Experimental setup as dictionary:
diff --git a/Examples/listDatabase.py b/Examples/listDatabase.py
index 5c5ed33a6fadc91c50012a522e685a3580fe0ba0..5a7e2b39f3e3779f8509de250138da5d966548de 100644
--- a/Examples/listDatabase.py
+++ b/Examples/listDatabase.py
@@ -1,12 +1,12 @@
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
-import libhreels.calcHREELS as hh
+import libhreels.calcHREELS20 as hh
import os, json
libDir = os.path.dirname(os.path.realpath(hh.__file__))
-file = os.path.join(libDir,'materials.json')
+file = os.path.join(libDir,'materials20.json')
with open(file) as json_file:
materials = json.load(json_file)
diff --git a/Examples/readAuger.py b/Examples/readAuger.py
index 5a7d44124d5934132226d743032692504ba26895..bd756faec416d5b69163777ce6bf58c365c525a0 100644
--- a/Examples/readAuger.py
+++ b/Examples/readAuger.py
@@ -1,8 +1,8 @@
import matplotlib.pyplot as plt
# %matplotlib qt
import libhreels.Auger as h
-# path="./Examples/data"
-path="./data"
+path="./Examples/data"
+#path="./data"
plt.rcParams['figure.figsize'] = (6,4)
plt.rcParams['figure.dpi'] = 200
plt.rcParams['font.size'] = 12
diff --git a/Examples/readScans.png b/Examples/readScans.png
index 5b91ee8036bdcaf9c3017a69708b5743aa897d5e..e28f32651535cce2de1830b322d0188911e7343e 100644
Binary files a/Examples/readScans.png and b/Examples/readScans.png differ
diff --git a/Examples/showHREELS1.png b/Examples/showHREELS1.png
index 4b8ed81c25365479e7b7baa9eafcf1d2110270cf..eb8543a376835c67a2212f27b96166f4e7c1a8f8 100644
Binary files a/Examples/showHREELS1.png and b/Examples/showHREELS1.png differ
diff --git a/Examples/showHREELS1.py b/Examples/showHREELS1.py
index 3292303575d31aee14490951cd0c4639f850acc8..6fc32b5a2f6cdbdeba907238c03eda09c2f33765 100644
--- a/Examples/showHREELS1.py
+++ b/Examples/showHREELS1.py
@@ -3,8 +3,8 @@ import os
import matplotlib.pyplot as plt
import numpy as np
-# datapath = "./Examples/data"
-datapath = "./data"
+datapath = "./Examples/data"
+#datapath = "./data"
d2 = hh.HREELS('c1a02', datapath=datapath)
d1 = hh.HREELS('c1a03', datapath=datapath)
diff --git a/Examples/showHREELS2.png b/Examples/showHREELS2.png
index 62bd86d21e5c5487f57aaca03c5ee3dad14328f8..cdfff9fa99952b0d2bff2f516ff1a85e9ab992cc 100644
Binary files a/Examples/showHREELS2.png and b/Examples/showHREELS2.png differ
diff --git a/Examples/showHREELS2.py b/Examples/showHREELS2.py
index 6bedd70595bf3b41ae2223c6f9e7684ad23e5ef8..6dd4e9291f78e036a5e9a1ee2ee56b061e489f03 100644
--- a/Examples/showHREELS2.py
+++ b/Examples/showHREELS2.py
@@ -2,8 +2,8 @@ from libhreels import HREELS
import matplotlib.pyplot as plt
import os
-d = './data'
-# d = './Examples/data'
+#d = './data'
+d = './Examples/data'
file = "c1a03"
factor = 20
diff --git a/Examples/startAuger_Browser.py b/Examples/startAuger_Browser.py
index dcf8883dd4268921b95f95da19aee40662b4c62c..bbc7b798bb7a1bab9a735099245db7829b639e4c 100644
--- a/Examples/startAuger_Browser.py
+++ b/Examples/startAuger_Browser.py
@@ -1,9 +1,6 @@
from libhreels import ViewAuger
-import sys
-import matplotlib.pyplot as plt
-import numpy as np
-# datapath = ".\\data"
-datapath = "./data"
+# datapath = "./data"
+datapath = "./Examples/data"
ViewAuger.runViewer(datapath=datapath, startWithFile='181214_A_01')
\ No newline at end of file
diff --git a/README.md b/README.md
index bad16ed40072cde768a1c9a52e3bfa89b471216b..a0c067b5d7514dff255756f948d0465e5295fab6 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,11 @@
# HREELS data handling and simulation
-Set of data analysis routines for surface vibrational spectroscopy based on a Delta05 spectrometer. Simple data reading and plotting routines are provided, but also advanced routines for in depth analysis.
+Set of data analysis routines for surface vibrational spectroscopy based on a Delta05 spectrometer with software from VSI or SPECS.
+Simple data reading and plotting routines are provided, but also advanced routines for in depth analysis.
A fast data browser as graphical user interface is included as **ViewHREELS.py**. It can also take command line arguments.
-Additional a python interface to a Fortran-based calculation of full HREEL spectra is provided in **class lambin** (in **calcHREELS.py**). Note that this part requires a local compilation (via f2py3) from the Fortran90 files. Only for the Linux-based python version 3.6, there are complied files provided. For details see below.
+Additional a python interface to a Fortran-based calculation of full HREEL spectra is provided in **class lambin** (in **calcHREELS20.py**). Note that this part requires a local compilation (via "f2py3 -c -m myEels20 myEels20.f90") from the Fortran90 files. Only for the Linux-based python version 3.8, there are complied files provided. For details see below.
Most of the routines are within the **class HREELS** (HREELS.py). A simple command line program is "showHREELS.py", which reads one data file and plots the spectrum with a second amplified trace.
@@ -50,9 +51,14 @@ The following methods are defined within the class:
plotWaterFall(...):
-# calcHREELS
+# calcHREELS20
These routines are based on the publication "Computation of the surface electron-energy-loss spectrum in specular geometry for an arbitrary plane-stratified medium" by P. Lambin, J.-P. Vigneron, and A. A. Lucas, in the Journal "Computer Physics Communications 60, 351-64(1990)".
+The code is modified to comply with Fortran90 and wraped to python function. See the example in calcHREELS1.py for NiO(001). The parameters for epsilon_infinity "eps": 5.25, as well as frequency and width for the TO phonon, "wTO": [393.7], "gTO": [10.8], and the LO phonon, "wLO": [584.7], "gLO": [10.8], need to be specified.
+
+Complex calculations for perovskite oxides are provided in the examples calcHREELS2.py and calcHREELS3.py.
+
+
diff --git a/dist/libhreels-2.0.0-py3-none-any.whl b/dist/libhreels-2.0.0-py3-none-any.whl
index 659b80607235f465bb4c841d9b1a4e2996bf3461..c4312cd3ed7bb27ffffeba9ecd3a7af8bdcd2db2 100644
Binary files a/dist/libhreels-2.0.0-py3-none-any.whl and b/dist/libhreels-2.0.0-py3-none-any.whl differ
diff --git a/dist/libhreels-2.0.0.tar.gz b/dist/libhreels-2.0.0.tar.gz
index 874feadb8c724a331e47cfd2b0198ed836bcf362..ddebd6e475865594105418d502797bae4dd01a2a 100644
Binary files a/dist/libhreels-2.0.0.tar.gz and b/dist/libhreels-2.0.0.tar.gz differ
diff --git a/libhreels/ViewAuger.py b/libhreels/ViewAuger.py
index 7f384a30ce14c0b3d5820467fd86a5e76a768483..d6757db61600678d49240ecd4a637836729f333c 100644
--- a/libhreels/ViewAuger.py
+++ b/libhreels/ViewAuger.py
@@ -11,8 +11,8 @@ from matplotlib.backends.backend_qt5agg import \
FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt5agg import \
NavigationToolbar2QT as NavigationToolbar
-#from PyQt5 import QtCore, QtGui, QtWidgets, uic
-from PyQt6 import QtCore, QtGui, QtWidgets, uic
+from PyQt5 import QtCore, QtGui, QtWidgets, uic
+#from PyQt6 import QtCore, QtGui, QtWidgets, uic
import libhreels as hh
hhPath = hh.__path__[0]
from libhreels.Auger import Auger
diff --git a/libhreels/ViewHREELS.py b/libhreels/ViewHREELS.py
index 956ac06096e08f12a3539501ee020021f5c7da77..8b9428a14a6a4e2a385912d7c587be2b1c838bc7 100644
--- a/libhreels/ViewHREELS.py
+++ b/libhreels/ViewHREELS.py
@@ -9,15 +9,15 @@ from matplotlib.backends.backend_qt5agg import \
FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt5agg import \
NavigationToolbar2QT as NavigationToolbar
-#from PyQt5 import QtCore, QtGui, QtWidgets, uic
-from PyQt6 import QtCore, QtGui, QtWidgets, uic
+from PyQt5 import QtCore, QtGui, QtWidgets, uic
+#from PyQt6 import QtCore, QtGui, QtWidgets, uic
import libhreels as hh
hhPath = os.path.dirname(__file__)
from libhreels.HREELS import HREELS, myPath
from datetime import datetime
import argparse
-import libhreels.expLogbook as lgb
-# import libhreels.NoExpLogbook as lgb
+#import libhreels.expLogbook as lgb
+import libhreels.NoExpLogbook as lgb
# fix HighRes Displays
QtWidgets.QApplication.setAttribute(QtCore.Qt.AA_EnableHighDpiScaling, True)
diff --git a/libhreels/calcHREELS20.png b/libhreels/calcHREELS20.png
new file mode 100644
index 0000000000000000000000000000000000000000..2bac0b5884373fa7917565db9fe2a7b59e5f5d7a
Binary files /dev/null and b/libhreels/calcHREELS20.png differ
diff --git a/libhreels/calcHREELS20.py b/libhreels/calcHREELS20.py
index 2c0715aa067994b4ae05973c786e334a9689f8cd..273908805d926669527625d065894cce5816dbf6 100755
--- a/libhreels/calcHREELS20.py
+++ b/libhreels/calcHREELS20.py
@@ -42,7 +42,7 @@ def importMaterials(string='', path=libDir):
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')
+ file = os.path.join(myPath(path),'materials20.json')
with open(file) as json_file:
materials = json.load(json_file)
try:
@@ -195,7 +195,6 @@ class lambin:
def myMain():
import matplotlib.pyplot as plt
import numpy as np
- from copy import deepcopy
import os
x = np.linspace(-100.,1000,400)
diff --git a/libhreels/materials.json b/libhreels/materials.json
deleted file mode 100755
index bc95115f29de261cc59ee045862b0d8d2cba7320..0000000000000000000000000000000000000000
--- a/libhreels/materials.json
+++ /dev/null
@@ -1,243 +0,0 @@
-{
- "NiO": {
- "eps": 5.25,
- "wTO": [393.7],
- "gTO": [10.8],
- "wLO": [584.7],
- "gLO": [10.8],
- "reference": "Schumann'22",
- "name": "NiO"
- },
- "BaO": {
- "eps": 3.7,
- "wTO": [145.0],
- "gTO": [46.0],
- "wLO": [422.0],
- "gLO": [56.0],
- "reference": "Goian'18",
- "name": "BaO"
- },
- "BaO_2": {
- "eps": 3.7,
- "wTO": [145.0],
- "gTO": [46.0],
- "wLO": [414.0],
- "gLO": [56.0],
- "reference": "...",
- "name": "BaO"
- },
- "Willet": {
- "eps": 4.12,
- "wTO": [188.0, 427.0, 495.72, 650.79, 708.2],
- "gTO": [0.4, 5.0, 3.8, 22.5, 55.3],
- "wLO": [276.4, 596.1, 495.5, 744.1, 702.2],
- "gLO": [3.7, 7.2, 3.8, 12.1, 66.0],
- "reference": "Willet-Gies'14",
- "name": "LaAlO_3"
- },
- "STO": {
- "eps": 5.14285,
- "wTO": [ 98.0 , 174.367 , 548.821 ],
- "gTO": [ 62.488 , 9.89216, 25.0 ],
- "wLO": [171.679 , 475.624 , 817.369 ],
- "gLO": [ 6.11124, 29.9772 , 60.0 ],
- "reference": "Crandles'99",
- "name": "SrTiO_3"
- },
- "Crandles99_2e18_80K": {
- "eps": 5.16,
- "wTO": [ 50.5 , 174.2 , 548.9 ],
- "gTO": [ 16 , 6.0 , 10.0 ],
- "wLO": [171.56 , 481.0 , 805.7 ],
- "gLO": [ 4.8 , 5.35 , 29.0 ],
- "reference": "Crandles'99 2e18 80K",
- "name": "SrTiO_3"
- },
- "Crandles99_2e18_300K": {
- "eps": 5.14,
- "wTO": [93.9831, 173.175, 543.346 ],
- "gTO": [25.8272, 8.42197, 17.5385 ],
- "wLO": [169.788, 474.634, 795.385 ],
- "gLO": [ 4.9, 7.07977, 28.3571 ],
- "reference": "Crandles'99 2e18 300K",
- "name": "SrTiO_3"
- },
- "Crandles99_2e19_80K": {
- "eps": 5.1,
- "wTO": [50.0211, 173.071, 548.723],
- "gTO": [34.8738, 5.00001, 10.0],
- "wLO": [ 170.94, 493.881, 823.411],
- "gLO": [4.69999, 15.3783, 52.6903],
- "reference": "Crandles'99 2e19 80K",
- "name": "SrTiO_3"
- },
- "Gervais": {
- "eps": 5.2,
- "wTO": [ 88.5533 , 175.999 , 543.0 ],
- "gTO": [ 27.9999 , 5.99999, 18.7496 ],
- "wLO": [171.967 , 475.999 , 796.556 ],
- "gLO": [ 3.09999, 4.62163, 25.0 ],
- "reference": "Gervais'93",
- "comment": "300 K",
- "name": "Gervais_SrTiO_3"
- },
- "Gervais_doped": {
- "eps": 5.01,
- "wTO": [ 93.0 , 175.99 , 547.99 ],
- "gTO": [ 42.874 , 8.9935, 18.999 ],
- "wLO": [172.79 , 472.53 , 795.0 ],
- "gLO": [ 4.1918, 9.9999, 26.0 ],
- "reference": "Gervais'93",
- "comment": "300 K",
- "name": "SrTiO_3"
- },
- "J0B04": {
- "eps": 5.2,
- "wTO": [60, 171.8, 531.97],
- "gTO": [80, 30, 17.706],
- "wLO": [170.9, 481, 814.92],
- "gLO": [53.6, 13.3 , 36.497],
- "reference": "Fit to J0B04",
- "name": "SrTiO_3"
- },
- "Fitted_EBE02" : {
- "eps": 5.2,
- "wTO": [94.85, 178.28, 528.921],
- "gTO": [2.00, 9.14, 18.73],
- "wLO": [176.49, 476.26, 833.7236],
- "gLO": [10.43, 18.26 , 52.63],
- "reference": "Fit by HHE",
- "name": "Fitted_EBE02_300K"
- },
- "Fitted_EBE02_210914" : {
- "eps": 3.024434,
- "wTO": [147.877698, 175.0514, 531.3065],
- "gTO": [5.02880219, 13.73305, 18.76442],
- "wLO": [168.951344, 475.621265, 833.761677],
- "gLO": [10.0004693, 19.75355, 60.70802],
- "reference": "Fit by HHE, 0.05% doped, 300K, fitted with calcHREELS3, last updated: 210916",
- "name": "Fitted_EBE02_300K"
- },
- "Fitted_K7701_210916" : {
- "eps": 4.57655596,
- "wTO": [50.0211, 174.060581, 534.354046],
- "gTO": [1.56852054, 8.17815492, 12.6232742],
- "wLO": [171.56, 482.400957, 820.342514],
- "gLO": [8.52832484, 21.0070097, 58.6519913],
- "reference": "Fit by HHE, 0.05% doped, 114K, fitted with calcHREELS3, last updated: 210917",
- "name": "Fitted_K7701_100K"
- },
-
-
- "STO_DEG1": {
- "eps": 5.14285,
- "wTO": [ 98.0 , 174.367 , 548.821, 0.0 ],
- "gTO": [ 62.488 , 9.89216, 25.0 , -100 ],
- "wLO": [171.679 , 475.624 , 817.369, 200 ],
- "gLO": [ 6.11124, 29.9772 , 60.0 , 0.0 ],
- "reference": "not working",
- "name": "SrTiO_3"
- },
- "KTO": {
- "eps": 5.12878,
- "wTO": [91.0, 201.0, 548.0, 710.0, 752.0],
- "gTO": [6.5, 6.0, 13.0, 130.0, 20.0],
- "wLO": [186.0, 423.0,701.0, 750.0, 821.0],
- "gLO": [6.0, 6.5, 130.0, 20.0, 13.0],
- "reference": "Jandl'91",
- "name": "KTaO_3"
- },
- "Massa97": {
- "eps": 5.12878,
- "wTO": [1317.9, 0.0, 181.8, 209.9, 330.5],
- "gTO": [1500.5, 1170.0, 185.1, 321.8, 391.5],
- "wLO": [1618.5, 900.3, 12.8, 829.7, 49.3],
- "gLO": [3967.7, 689.4, 5.0, 41.9, 85.9],
- "reference": "Massa'97",
- "name": "NdNiO_3"
- },
- "Massa97_77K": {
- "eps": 5.12878,
- "wTO": [ 398.0, 430.4, 470.0, 561.1, 729.9, 0.0],
- "gTO": [ 428.8, 435.8, 542.1, 567.8, 1108.5, 146.6],
- "wLO": [ 149.2, 31.5, 254.5, 107.3, 1155.6, 1705.3],
- "gLO": [ 887.6, 28.0, 290.8, 110.8, 1504.6, 2087.1],
- "reference": "Massa'97",
- "name": "NdNiO_3"
- },
- "BTO": {
- "eps": 5.1295,
- "wTO": [177.4 , 272.84 , 506.379],
- "gTO": [ 1.9262 , 93.3091 , 42.60 ],
- "wLO": [184.003 , 470.972 , 739.651],
- "gLO": [ 9.718 , 14.31 , 33.044],
- "reference": "BTO323",
- "name": "BaTiO_3"
- },
- "S2R1O4": {
- "eps": 5.5,
- "wTO": [195, 369, 412.5, 466 ],
- "gTO": [9.2, 18, 17.1, 18 ],
- "wLO": [250, 400, 450, 600],
- "gLO": [9.2, 18, 17.1, 18 ],
- "reference": "Farbe der Ruthenate",
- "name": "SRO-214"
- },
- "Ag": {
- "eps": 1.0,
- "wTO": [-31460.0],
- "gTO": [-161.0],
- "wLO": [1],
- "gLO": [1],
- "reference": "Irani'71, Ehrenreich'62, Furtak'75",
- "name": "Silver"
- },
- "Ag4": {
- "eps": 1.0,
- "wTO": [0],
- "gTO": [161.0],
- "wLO": [31460],
- "gLO": [161.0],
- "reference": "Irani'71, Ehrenreich'62, Furtak'75",
- "name": "Silver new"
- },
- "Pt": {
- "eps": 1.0,
- "wTO": [-185541.0],
- "gTO": [-16134.0],
- "wLO": [1],
- "gLO": [1],
- "reference": "Weaver'75, Seignac'72",
- "name": "Platinum"
- },
- "Pt4": {
- "eps": 1.0,
- "wTO": [0],
- "gTO": [16134.0],
- "wLO": [185541.0],
- "gLO": [16134.0],
- "reference": "Weaver'75, Seignac'72",
- "name": "Platinum new"
- },
-
- "test": {
- "eps": 5.25,
- "wTO": [493.7, -10.0],
- "gTO": [10.8, -5.0],
- "wLO": [584.7, 1],
- "gLO": [10.8, 1],
- "reference": "for playing",
- "name": "test"
- },
- "Vacuum": {
- "eps": 1.0,
- "wTO": [1000],
- "gTO": [100],
- "wLO": [1001],
- "gLO": [100],
- "reference": "...",
- "name": "vacuum"
- }
-}
-
-
diff --git a/libhreels/mod_doeels.mod b/libhreels/mod_doeels.mod
deleted file mode 100644
index 781cf9ba42d1247fa1ce9514c5fd5045af2c771a..0000000000000000000000000000000000000000
Binary files a/libhreels/mod_doeels.mod and /dev/null differ
diff --git a/libhreels/myEels2.f90 b/libhreels/myEels2.f90
deleted file mode 100644
index 5809129bb9015fcd6714a93aa8e05fe9b6fe22e1..0000000000000000000000000000000000000000
--- a/libhreels/myEels2.f90
+++ /dev/null
@@ -1,997 +0,0 @@
-module mod_doeels
- contains
-subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, &
- layers, neps, nper, name, name_size, thick, epsinf, nos, osc, osc_size,&
- contrl, mode, wn_array, debug, f_array, wn_array_size)
-
-! ******************************************************************
-! * *
-! * compute the classical eels spectrum of an arbitrary plane- *
-! * statified medium made from isotropic materials in specular *
-! * geometry using the dielectric theory of eels. *
-! *
-! * It is based on the work of Lambin'90 and modified for the use *
-! * within python *
-! * (KMS and WFW, Martin-Luther-Universität Halle-Wittenberg) * *
-! ******************************************************************
-
- implicit none
-
- integer, parameter :: nt = 5
-
- double precision, intent(in) :: e0, theta, phia, phib, wmin, wmax, dw
- double precision, intent(in) :: thick(name_size), epsinf(name_size), osc(3, osc_size)
- character*10, intent(in) :: name(name_size)
- character*10, intent(in) :: contrl, mode
- integer, intent(in) :: name_size, osc_size, neps
- integer, intent(in) :: wn_array_size
- integer, intent(in) :: layers, nper, nos(name_size)
- double precision, intent(out) :: wn_array(wn_array_size), f_array(wn_array_size)
- logical, intent(in), optional :: debug
-
- logical :: ration, user
- integer :: i, iw, lstart, nofu, nout, nw, lmax, flag
- double precision :: a, acoef, aerr, alpha, argmin, argmax, b, bcoef, beta, &
- c1, c2, ccoef, cospsi, dlimf, dx, elleps, ener, epsmac, facru, f, f0, &
- f1, fpic, pi, prefac, psia, psii, rerr, ru, sinpsi, t, &
- tanpsi, table, um, widt, wn, wpic, x, xmin, xmax, z, z1, z2
- double complex, dimension(30) :: eps
- dimension table(nt)
- logical debugFirstRun
- common /control/ debugFirstRun
-
- common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
- common / mulayr / argmin, argmax, epsmac, flag
-
- !if (present(debug)) debug = .True.
- debugFirstRun = .True.
-
- data aerr / 0.0d0 /, rerr / 1.0d-06 /, f / 0.0d0 /, f1 / 0.0d0 /
-
-! *** machine-dependent constants
-! *** epsmac + 1.0 = epsmac , cosh(argmin) = 1.0 , tanh(argmax) = 1.0
- flag = 1
-
- if (debug) write(*,*) 'doeels:'
- if (debug) write(*,*) 'thick: ', size(thick)
- if (debug) write(*,*) 'epsinf: ', size(epsinf), epsinf
- if (debug) write(*,*) 'osc: ', size(osc), size(osc, 1), size(osc, 2)
- if (debug) write(*,*) 'osc(1,1):', osc(1,1)
- if (debug) write(*,*) 'nos: ', size(nos)
- if (debug) write(*,*) 'wn_array: ', size(wn_array)
- if (debug) write(*,*) 'f_array: ', size(f_array)
-
- pi = 4.0d0 * atan(1.0d0)
- epsmac = 1.0d0
-1 epsmac = epsmac / 2.0d0
- if (1.0d0 + epsmac > 1.0d0) goto 1
- argmin = sqrt(2.0d0 * epsmac)
- argmax = 0.5d0 * log(2.0d0 / epsmac)
-
- dlimf = 0.0d0
- ration = .false.
-
-! *** read target specifications
-
- user = layers == 0
- if (user) then
-
- if (layers == 1) ration = .true.
- lstart = layers - nper + 1
- endif
-
-! *** initialize constants
-
- lmax = size(thick)
- nw = size(wn_array)
- ! if (debug) write(*,*) 'lmax: ', lmax
- ! allocate(eps(lmax))
- ! if (debug) write(*,*) 'eps: ', size(eps)
- nout = 1 + nw / 20
- ener = 8.065d+03 * e0
- psia = phia / 180.0d0 * pi
- psii = theta / 180.0d0 * pi
- cospsi = cos(psii)
- sinpsi = sin(psii)
- tanpsi = tan(psii)
- prefac = sqrt(2.555d+05 / e0)/(1.37d+02 * cospsi)
- facru = psia / cospsi * sqrt(0.2624664d0 * e0)
- elleps = (1.0d0 - phia / phib) * (1.0d0 + phia / phib)
- acoef = sinpsi**2 + elleps * cospsi**2
- bcoef = sinpsi * cospsi
- if (dlimf > 0.0d0) then
- ration = .false.
- write(*,*) ' = > electron attracted by an image charge = ', dlimf
-! *** dlimf : half the length unit imposed by the image force
- dlimf = 1.80d0 * dlimf/(e0 * cospsi**2)
- endif
- ! if (debug) write(*,*) 'ration: ', ration
- if (.not. ration) goto 35
-
-! *** set up coefficients for the rational approximation to the integral
-
- write(*,*) ' = > set up a rational approximation to the integral'
- call quanc8(fun, 0.0d0, pi / 2.0d0, aerr, rerr, alpha, c1, nofu, c2, eps, thick, layers, nper)
- alpha = (2.0d0 / pi)**2 * alpha
- c1 = 2.0d0 / pi / sqrt(1.0d0 - elleps) * sinpsi * alpha**2
- if (c1 > 0.99d0) goto 30
- c2 = 3.0d0 * alpha**2 / (1.0d0 - c1)
- c1 = c1 * c2
- xmin = wmin / (2.0d0 * ener * psia)
- xmax = wmax / (2.0d0 * ener * psia)
- if (xmin <= 0.0d0) xmin = 0.0d0
- dx = max(0.02d0, (xmax - xmin) / nt)
- z1 = 0.0d0
- z2 = 0.0d0
- do i = 1, nt
- x = xmin + i * dx
- call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
- table(i) = f
- f = f * (1.0d0 + alpha * x)**2
- if (abs(c2 * f - c1) < c2 * rerr) cycle
- z = (1.0d0 - f) / (c2 * f - c1)
- if (z <= 0.0d0) cycle
- z1 = z1 + x * z * (x**2 - z)
- z2 = z2 + (x * z)**2
- enddo
- if (z2 == 0.0d0) goto 30
- beta = z1 / z2
- z = 0.0d0
- do i = 1, nt
- x = xmin + i * dx
- z = z + (table(i) - qrat(x,alpha,beta,c1,c2))**2
- enddo
- z = sqrt(z) / nt
- if (z > 5.0d-03) goto 30
- write(*, 109) alpha, c1, c2, beta, z
- goto 35
-30 write(*, *) ' ===> cannot do it'
- ration = .false.
-
-! *** loop over the energy losses
-
-35 if (debug) write(*, 110)
- do iw = 1, nw
- f0 = f1
- f1 = f
- f = 0.0d0
- wn = wmin + (iw - 1) * dw
-! if (debug) write(*,*) 'wn: ', wn
- if (wn < 0.0d0) goto 45
- if (wn == 0.0d0) goto 40
- if (.not. user) call seteps(nos, osc_size, osc, epsinf, wn, name_size, eps)
- ! subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
-
- x = wn / (2.0d0 * ener * psia)
- if (ration) then
- f = qrat(x,alpha,beta,c1,c2) * aimag(-2.0 / (1.0 + eps(1)))
- else
- call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
- endif
- f = prefac * f / wn
-40 continue
- wn_array(iw) = wn
- f_array(iw) = f
-! *** localize a peak using a parabolic interpolation
- if (iw < 3) goto 45
- if (f1 - f0 <= aerr) goto 45
- if (f1 - f <= aerr) goto 45
- a = (f1 - f0) + (f1 - f)
- if (a <= 4.0d0 * rerr * f1) goto 45
- b = 0.5d0 * (f1 - f0 + 3.0d0 * (f1 - f))
- t = b / a
- wpic = wn - t * dw
- fpic = f + 0.5d0 * b * t
- widt = sqrt(8.0d0 * fpic / a) * dw
- if (debug) write(*, 112) wpic, fpic, widt
-45 if (mod(iw, nout) == 0) then
- if (debug) write(*, 113) 100.0 * iw / nw, wn, f
- endif
- enddo
- return
-109 format(5x, 'alpha = ', f9.4, 4x, 'c1 = ', f9.4, 4x, 'c2 = ', f9.4, 4x, &
- 'beta = ', f9.4/5x, 'accuracy = ', e9.2)
-110 format(//' run (%) wn (cm**-1) pcl(wn) (cm) |', &
- ' peak location amplitude width')
-112 format(40x, f10.2, d12.4, f10.2)
-113 format(2x, f5.1, 3x, f11.3, d14.5)
-end subroutine doeels
-! ####################################################################
-
-double precision function qrat(x, alpha, beta, c1, c2)
- double precision, intent(in) :: x, alpha, beta, c1, c2
- qrat = (1.0d0 + x * (beta + c1 * x)) / ((1.0d0 + x * (beta + c2 * x)) * (1.0d0 + alpha * x)**2)
- return
-end function qrat
-! ####################################################################
-
-subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, d, layers, nper)
-
-! estimate the integral of fun(x) from a to b
-! to a user provided tolerance.
-! an automatic adaptive routine based on
-! the 8-panel newton-cotes rule (g. forsythe et al, 1977, p. 92)
-!
-! input ..
-!
-! fun the name of the integrand function subprogram fun(x).
-! a the lower limit of integration.
-! b the upper limit of integration.(b may be less than a.)
-! relerr a relative error tolerance. (should be non-negative)
-! abserr an absolute error tolerance. (should be non-negative)
-!
-! output ..
-!
-! result an approximation to the integral hopefully satisfying the
-! least stringent of the two error tolerances.
-! errest an estimate of the magnitude of the actual error.
-! nofun the number of function values used in calculation of result.
-! flag a reliability indicator. if flag is zero, then result
-! probably satisfies the error tolerance. if flag is
-! xxx.yyy , then xxx = the number of intervals which have
-! not converged and 0.yyy = the fraction of the interval
-! left to do when the limit on nofun was approached.
-
- implicit none
- external fun
- double precision :: fun
- double precision, intent(in) :: a
- double precision, intent(in) :: b
- double precision, intent(in out) :: abserr
- double precision, intent(in) :: relerr
- double precision, intent(out) :: result
- double precision, intent(out) :: errest
- integer, intent(out) :: nofun
- double precision, intent(out) :: flag
-
- double precision, intent(in) :: d(:)
- double complex, intent(in) :: eps(:)
- integer, intent(in) :: layers, nper
-
- double precision :: w0, w1, w2, w3, w4, area, x0, f0, stone, step, cor11, temp
- double precision :: qprev, qnow, qdiff, qleft, esterr, tolerr
- double precision :: qright(31), f(16), x(16), fsave(8, 30), xsave(8, 30)
- double precision :: dabs, dmax1
-
- integer :: levmin, levmax, levout, nomax, nofin, lev, nim, i, j
-
-! *** stage 1 *** general initialization
-! set constants.
-
-! write (*,*) 'quanc8:'
-! write (*,*) 'd: ', size(d)
-! write (*,*) 'eps: ', size(eps)
-
- levmin = 1
- levmax = 30
- levout = 6
- nomax = 5000
- nofin = nomax - 8 * (levmax - levout + 2**(levout + 1))
-
-! trouble when nofun reaches nofin
-
- w0 = 3956.0d0 / 14175.0d0
- w1 = 23552.0d0 / 14175.0d0
- w2 = -3712.0d0 / 14175.0d0
- w3 = 41984.0d0 / 14175.0d0
- w4 = -18160.0d0 / 14175.0d0
-
-! initialize running sums to zero.
-
- flag = 0.0d0
- result = 0.0d0
- cor11 = 0.0d0
- errest = 0.0d0
- area = 0.0d0
- nofun = 0
- if (a == b) return
-
-! *** stage 2 *** initialization for first interval
-
- lev = 0
- nim = 1
- x0 = a
- x(16) = b
- qprev = 0.0d0
- f0 = fun(x0, eps, d, layers, nper, size(eps))
- stone = (b - a) / 16.0d0
- x(8) = (x0 + x(16)) / 2.0d0
- x(4) = (x0 + x(8)) / 2.0d0
- x(12) = (x(8) + x(16)) / 2.0d0
- x(2) = (x0 + x(4)) / 2.0d0
- x(6) = (x(4) + x(8)) / 2.0d0
- x(10) = (x(8) + x(12)) / 2.0d0
- x(14) = (x(12) + x(16)) / 2.0d0
- do j = 2, 16, 2
- f(j) = fun(x(j), eps, d, layers, nper, size(eps))
- enddo
- nofun = 9
-
- do
-
-! *** stage 3 *** central calculation
-! requires qprev, x0, x2, x4, ..., x16, f0, f2, f4, ..., f16.
-! calculates x1, x3, ...x15, f1, f3, ...f15, qleft, qright, qnow, qdiff, area.
-
- x(1) = (x0 + x(2)) / 2.0d0
- f(1) = fun(x(1), eps, d, layers, nper, size(eps))
- do j = 3, 15, 2
- x(j) = (x(j - 1) + x(j + 1)) / 2.0d0
- f(j) = fun(x(j), eps, d, layers, nper, size(eps))
- enddo
- nofun = nofun + 8
- step = (x(16) - x0) / 16.0d0
- qleft = (w0 * (f0 + f(8)) + w1 * (f(1) + f(7)) + w2 * (f(2) + f(6)) &
- + w3 * (f(3) + f(5)) + w4 * f(4)) * step
- qright(lev + 1) = (w0 * (f(8) + f(16)) + w1 * (f(9) + f(15)) + w2 * (f(10) + f(14)) &
- + w3 * (f(11) + f(13)) + w4 * f(12)) * step
- qnow = qleft + qright(lev + 1)
- qdiff = qnow - qprev
- area = area + qdiff
-
-! *** stage 4 *** interval convergence test
-
- esterr = dabs(qdiff) / 1023.0d0
- tolerr = dmax1(abserr, relerr * dabs(area)) * (step / stone)
-
- if (lev >= levmin) then
- if (lev >= levmax) then
-! current level is levmax.
- flag = flag + 1.0d0
- else
- if (nofun > nofin) then
-! *** stage 6 *** trouble section
-! number of function values is about to exceed limit.
- nofin = 2 * nofin
- levmax = levout
- flag = flag + (b - x0) / (b - a)
- else
- if (esterr > tolerr) then
-! *** stage 5 *** no convergence
-! locate next interval.
- nim = 2 * nim
- lev = lev + 1
-! store right hand elements for future use.
- do i = 1, 8
- fsave(i, lev) = f(i + 8)
- xsave(i, lev) = x(i + 8)
- enddo
-! assemble left hand elements for immediate use.
- qprev = qleft
- do i = 1, 8
- f(18 - 2 * i) = f(9 - i)
- x(18 - 2 * i) = x(9 - i)
- enddo
- cycle
- endif
- endif
- endif
-
-! *** stage 7 *** interval converged
-! add contributions into running sums.
- result = result + qnow
- errest = errest + esterr
- cor11 = cor11 + qdiff / 1023.0d0
-! locate next interval.
- do while (nim /= 2 * (nim / 2))
- nim = nim / 2
- lev = lev - 1
- enddo
- nim = nim + 1
-
- if (lev <= 0) exit
-
-! assemble elements required for the next interval.
- qprev = qright(lev)
- x0 = x(16)
- f0 = f(16)
- do i = 1, 8
- f(2*i) = fsave(i, lev)
- x(2*i) = xsave(i, lev)
- enddo
- cycle
- else
-! *** stage 5 *** no convergence
-! locate next interval.
- nim = 2 * nim
- lev = lev + 1
-! store right hand elements for future use.
- do i = 1, 8
- fsave(i, lev) = f(i + 8)
- xsave(i, lev) = x(i + 8)
- enddo
-! assemble left hand elements for immediate use.
- qprev = qleft
- do i = 1, 8
- f(18 - 2 * i) = f(9 - i)
- x(18 - 2 * i) = x(9 - i)
- enddo
- endif
-
- enddo
-
- ! *** stage 8 *** finalize and return
- result = result + cor11
-
-! make sure errest not less than roundoff level.
- if (errest /= 0.0d0) then
- temp = dabs(result) + errest
- do while (temp == dabs(result))
- errest = 2 * errest
- temp = dabs(result) + errest
- enddo
- endif
- return
-end subroutine quanc8
-! ############################################################################
-
-double precision function fun(phi)
-
-! ******************************************************************
-! * *
-! * integrand of the expression of the 1st order term in the *
-! * expansion of the eels integral for a homogeneous target. *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in out) :: phi
-
- logical :: user, ration
- double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps, ru
- double precision :: sinphi, sinpsi, tanpsi, um, wn
-
- common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
-
- sinphi = sin(phi)
- fun = sqrt((1.0d0 - elleps + elleps * sinphi**2) * &
- (1.0d0 - sinpsi * sinphi) * &
- (1.0d0 + sinpsi * sinphi))
- return
-end function fun
-! ##########################################################################
-
-subroutine queels(x, f, aerr, rerr, facru, eps, d, layers, nper)
-
-! ******************************************************************
-! * *
-! * perform q-space integration for computing the eels spectrum of *
-! * a isotropic target using polar coordinates. *
-! * *
-! * x is the dimensionless energy loss hbar*omega/(2*e0*phia) *
-! * aerr and rerr are the desired absolute and relative accuracies *
-! * facru*x is the units of wavevectors omega/v_perpendicular *
-! * f is the q-integral multiplied by (2/pi)**2 *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: x
- double precision, intent(out) :: f
- double precision, intent(in out) :: aerr
- double precision, intent(in out) :: rerr
- double precision, intent(in) :: facru
-
- double precision, intent(in) :: d(:)
- double complex, intent(in) :: eps(:)
- integer, intent(in) :: layers, nper
-
- logical :: ration, user
- double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps
- double precision :: error, flag, ru, sinpsi
- double precision :: u1, u2, um, ut, tanpsi, wn, y
- integer :: ie, nofu
- dimension error(3), flag(3)
-
- common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
-
-! write (*,*) 'queels:'
-! write (*,*) 'eps: ', size(eps)
-! write (*,*) 'd: ', size(d)
-
- f = 0.0d0
- if (x <= 0.0d0) then
- return
- endif
- ru = facru * x
- ccoef = cospsi**2 / x
- ut = ccoef - bcoef
- u1 = abs(ut)
- u2 = ccoef + bcoef
- if (ut > 0.0d0) then
- call quanc8(fint1, 0.0d0, u1, aerr, rerr, y, error(1), nofu, flag(1), eps, d, layers, nper)
- f = y
- else
- flag(1) = 0.0d0
- endif
- if (u2 > u1) then
- call quanc8(fint2, u1, u2, aerr, rerr, y, error(2), nofu, flag(2), eps, d, layers, nper)
- f = f + y
- else
- flag(2) = 0.0d0
- endif
- if (abs(acoef) > x * (1.0d0 - elleps) * bcoef) then
- um = sqrt(ccoef / x / (1.0d0 - elleps) + bcoef**2 / acoef)
- if (um > u2) then
- call quanc8(fint3, u2, um, aerr, rerr, y, error(3), nofu, flag(3), eps, d, layers, nper)
- f = f + y
- endif
- if (um < u1) then
- call quanc8(fint3, um, u1, aerr, rerr, y, error(3), nofu, flag(3), eps, d, layers, nper)
- f = f - y
- endif
- else
- flag(3) = 0.0d0
- endif
- do ie = 1, 3
- if (flag(ie) == 0.0d0) cycle
- write(*,*) ' +++ flag(', ie, ') =', flag(ie), ', error =', error(ie), ' +++'
- if (flag(ie) - aint(flag(ie)) > 0.5d-02) then
- stop '*** execution aborted ***'
- endif
- enddo
- f = (2.0d0 / 3.141592653589793238d0)**2 * f
- return
-end subroutine queels
-! ###########################################################################
-
-double precision function fint1(u, eps, d, layers, nper, eps_size)
-
-! ******************************************************************
-! * *
-! * integration over the azimutal angle from 0.0 to pi *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
- double complex, intent(in) :: eps(eps_size)
- integer, intent(in) :: layers, nper, eps_size
-
- logical :: ration, user
- double precision :: acoef, bcoef, ccoef, cospsi, den, dif, dlimf, e, elleps
- double precision :: pi, rom, rop, ru, sinpsi, sum, um, t
- double precision :: tanpsi, wn
-
- common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
-
- data pi / 3.141592653589793238d0 /
-
-! write (*,*) 'fint1:'
-! write (*,*) 'd: ', size(d)
-! write (*,*) 'eps: ', size(eps)
-
- if (u == 0.0d0) then
- fint1 = 0.0d0
- return
- endif
- e = tanpsi * u
- rom = (1.0d0 - e)**2 + u**2
- rop = (1.0d0 + e)**2 + u**2
- sum = rop + rom
- rom = sqrt(rom)
- rop = sqrt(rop)
- dif = rop - rom
- den = sqrt((2.0d0 - dif) * (2.0d0 + dif)) * rop * rom
- fint1 = pi * u**2 * (4.0d0 * sum - dif**2 * (sum - rop * rom)) / den**3
- if (ration) then
- return
- endif
- if (user) then
- fint1 = fint1 * usurlo(ru * u, wn)
- else
- fint1 = fint1 * surlos(ru * u, eps, d, layers, nper)
- if (dlimf > 0.0d0) then
- t = ru * u * dlimf
- fint1 = fint1 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
- endif
- endif
- return
-end function fint1
-! ###################################################
-
-double precision function fint2(u, eps, d, layers, nper, eps_size)
-
-! ******************************************************************
-! * *
-! * integration over the azimutal angle from 0.0 to phi < pi *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
- double complex, intent(in) :: eps(eps_size)
- integer, intent(in) :: layers, nper, eps_size
-
- logical :: ration, user
- double precision :: a, arg, b, b2, c, ccoef, cospsi, dlimf, elleps, phi
- double precision :: ru, sinpsi, um, t, tanpsi, wn, x
-
- common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
-
-! write (*,*) 'fint2:'
-! write (*,*) 'd: ', size(d)
-! write (*,*) 'eps: ', size(eps)
-
- if (u == 0.0d0) then
- fint2 = 0.0d0
- return
- endif
- b2 = b**2
- c = (1.0d0 - elleps) * (cospsi * u)**2 + (b - ccoef) * (b + ccoef)
- if (abs(a * c) > 1.0d-03 * b2) then
- x = (b - sqrt(b2 - a * c)) / a
- else
- x = a * c / b2
- x = 0.5d0 * c * (1.d0 + 0.25d0 * x * (1.d0 + 0.5d0 * x * (1.d0 + 0.625d0 * x))) / b
- endif
- arg = x / u
- if (abs(arg) > 1.0d0) then
- arg = sign(1.0d0, arg)
- endif
- phi = acos(arg)
- fint2 = phint(phi, tanpsi, u)
- if (ration) then
- return
- endif
- if (user) then
- fint2 = fint2 * usurlo(ru * u, wn)
- else
- fint2 = fint2 * surlos(ru * u, eps, d, layers, nper)
- if (dlimf > 0.0d0) then
- t = ru * u * dlimf
- fint2 = fint2 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
- endif
- endif
- return
-end function fint2
-! ########################################################################
-
-double precision function fint3(u, eps, d, layers, nper, eps_size)
-
-! ******************************************************************
-! * *
-! * integration over the azimutal angle from phi1 > 0 to phi2 < pi *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: u
-
- double precision, intent(in) :: d(eps_size)
- double complex, intent(in) :: eps(eps_size)
- integer, intent(in) :: layers, nper, eps_size
-
- logical :: ration, user
- double precision :: a, arg, b, ccoef, cospsi, dlimf, elleps, phi1, phi2
- double precision :: sinpsi, rac, ru, um, t, tanpsi, wn
-
- common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi, &
- ru, um, dlimf, wn, user, ration
-
-
-! write (*,*) 'fint3:'
-! write (*,*) 'd: ', size(d)
-! write (*,*) 'eps: ', size(eps)
-
- if (u == 0.0d0) then
- fint3 = 0.0d0
- return
- endif
- rac = sign(1.0d0, a) * cospsi * sqrt((1.0d0 - elleps) * a * (um - u) * (um + u))
- arg = (b - rac) / (u * a)
- if (abs(arg) > 1.0d0) arg = sign(1.0d0, arg)
- phi2 = acos(arg)
- fint3 = phint(phi2, tanpsi, u)
- arg = (b + rac) / (u * a)
- if (abs(arg) > 1.0d0) arg = sign(1.0d0, arg)
- phi1 = acos(arg)
- fint3 = fint3 - phint(phi1, tanpsi, u)
- if (ration) return
- if (user) then
- fint3 = fint3 * usurlo(ru * u, wn)
- else
- fint3 = fint3 * surlos(ru * u, eps, d, layers, nper)
- if (dlimf > 0.0d0) then
- t = ru * u * dlimf
- fint3 = fint3 * (1.d0 + t * log(t / (t + 0.26d0)))**2 / (1.d0 + 1.40d0 * t)
- endif
- endif
- return
-end function fint3
-! ##########################################################################
-
-double precision function usurlo(dq, wn)
-
-! ******************************************************************
-! * *
-! * user-supplied dielectric surface loss function aimag(g(dq, wn)) *
-! * input arguments : *
-! * dq : modulus of the two-dimensional surface wave vector *
-! * (angstroem**-1) *
-! * wn : frequency (cm**-1) *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: dq
- double precision, intent(in) :: wn
-
- if ((wn .GT. 0).AND.(dq.GT.0)) then
- write(*,*) 'hello, here is the user loss function'
- endif
- usurlo = 1.0d0
- return
-end function usurlo
-! #########################################################################
-
-double precision function surlos(dk, eps, d, layers, nper)
-
-! ******************************************************************
-! * *
-! * eels surface loss function for an arbitrary multilayered target*
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in) :: dk
- double complex, intent(in) :: eps(30)
- double precision, intent(in) :: d(:)
- integer, intent(in) :: layers, nper
-
- integer :: lmax, flag
- logical :: static, zero
- integer :: lstart, n
- double precision, allocatable :: arg(:)
- double precision :: argmin, argmax, cn, cnm1, epsmac, sn, snm1, t
- double complex :: a, b, csi, pnm2, pnm1, pn, pp, qnm2, qnm1, qn, qp, z
-
- common / mulayr / argmin, argmax, epsmac, flag
-
- zero(z) = (real(z) == 0.0) .and. (aimag(z) == 0.0)
-
-! write (*,*) 'surlos:'
-! write (*,*) 'd: ', size(d)
-! write (*,*) 'eps: ', size(eps)
-
- lmax = size(eps)
- allocate (arg(lmax))
- lstart = layers - nper + 1
- static = .true.
- n = 1
-1 arg(n) = dk * d(n)
- if (arg(n) > argmax .or. zero(eps(n))) goto 10
- static = .not. (n >= lstart .and. arg(n) > argmin)
- n = n + 1
- if (n <= layers) goto 1
-
-! *** periodic continued fraction, period = nper
-
- if (nper > 1) goto 2
- csi = eps(layers)
- goto 9
-2 if (static) goto 5
- cn = cosh(arg(lstart))
- sn = sinh(arg(lstart))
- pnm1 = 1.0
- pn = cn
- pp = eps(lstart) * sn
- qnm1 = 0.0
- qn = sn / eps(lstart)
- qp = pn
- do n = lstart + 1, layers
- cnm1 = cn
- snm1 = sn
- cn = cosh(arg(n))
- sn = sinh(arg(n))
- a = eps(n) * sn
- pp = cn * pp + a * pn
- qp = cn * qp + a * qn
- b = (eps(n - 1) / eps(n)) * (sn / snm1)
- a = cnm1 * b + cn
- pnm2 = pnm1
- pnm1 = pn
- qnm2 = qnm1
- qnm1 = qn
- pn = a * pnm1 - b * pnm2
- qn = a * qnm1 - b * qnm2
- enddo
- if (zero(qn)) goto 4
- a = 0.5 * (pn - qp) / qn
- b = sqrt(a**2 + pp / qn)
- pn = a - pn / qn
- if (abs(pn + b) > abs(pn - b)) then
- b = -b
- endif
- csi = a + b
- goto 9
-4 a = qp - pn
- if (zero(a)) goto 12
- csi = pp / a
- goto 9
-
-! *** small-dk limit of the periodic tail
-
-5 pn = 0.0
- qn = 0.0
- do n = lstart, layers
- pn = pn + d(n) * eps(n)
- qn = qn + d(n) / eps(n)
- enddo
- if (zero(qn)) goto 12
- csi = sqrt(pn / qn)
- if (aimag(csi) > 0.0) goto 9
- if ((aimag(csi) < 0.0) .or. (real(qn) < 0.0)) then
- csi = -csi
- endif
-9 n = lstart
- goto 11
-10 csi = eps(n)
-
-! *** backward algorithm
-
-11 n = n - 1
- if (n <= 0) goto 15
- if (arg(n) == 0.0d0) goto 11
- t = tanh(arg(n))
- b = eps(n) + csi * t
- if (zero(b)) goto 13
- csi = eps(n) * (csi + t * eps(n)) / b
- goto 11
-12 n = lstart
-13 if (n <= 1) goto 14
- n = n - 1
- csi = eps(n) / tanh(arg(n))
- goto 11
-14 surlos = 0.0d0
- return
-15 a = csi + 1.0
- if (zero(a)) then
- surlos = 2.0d0 / epsmac
- else
- surlos = aimag(-2.0 / a)
- endif
- return
-end function surlos
-! #########################################################################
-
-double precision function phint(phi, a, u)
-
-! ******************************************************************
-! * *
-! * evaluate the integral from zero to phi of *
-! * *
-! * u 2 *
-! * ( ----------------------------- ) dphi *
-! * 2 2 *
-! * (1 - a * u * cos(phi)) + u *
-! * *
-! * for 0 <= phi <= pi , u >= 0 and a >= 0 *
-! * *
-! ******************************************************************
-
- implicit none
-
- double precision, intent(in out) :: phi
- double precision, intent(in) :: a
- double precision, intent(in) :: u
-
- double precision :: ai, ar, bi, br, c, cpr, d, e, esr, pi, qr, ri, rm, root
- double precision :: rp, rr, s, spr, tm, tp, u2, x, zeta, zetai, zetar, zr
-
- pi = 3.141592653589793238d0
- c = cos(phi)
- s = sin(phi)
- u2 = u**2
- e = a*u
- if (u < 1.0d0 .and. e < 1.0d-02 * (1.0d0 + u2)) then
- zr = 1.0d0 + u2
- esr = e / zr
- phint = u2 / zr**2 * ((( (4.0d0 / 3.0d0) * (2.0d0 + c**2) * s * (5.0d0 - 3.0d0 * u2) * &
- esr + (phi + c * s) * (5.0d0 - u2)) * esr + 4.0d0 * s) * esr + phi)
- else
- rm = sqrt((1.0d0 - e)**2 + u2)
- tm = 0.5d0 * atan2(u, 1.0d0 - e)
- rp = sqrt((1.0d0 + e)**2 + u2)
- tp = 0.5d0 * atan2(u, 1.0d0 + e)
- root = sqrt(rm * rp)
- cpr = cos(tm + tp)
- spr = sin(tm + tp)
- if (c >= 0.0d0) then
- x = s / (1.0d0 + c)
- elseif (abs(s) > 1.0d-07) then
- x = (1.0d0 - c) / s
- endif
- if ((c >= 0.0d0) .or. (abs(s) > 1.0d-07)) then
- zeta = sqrt(rm / rp)
- zetar = -zeta * sin(tm - tp)
- zetai = zeta * cos(tm - tp)
- br = 0.5d0 * log(((zetar + x)**2 + zetai**2) / ((zetar - x)**2 + zetai**2))
- bi = atan2(zetai, zetar + x) - atan2(zetai, zetar - x)
- rr = -(br * spr - bi * cpr) / root
- ri = -(bi * spr + br * cpr) / root
- d = e * s / ((1.0d0 - e * c)**2 + u2)
- ar = d * (1.0d0 - e * c) - rr + u * ri
- ai = -d * u - ri - u * rr
- else
- rr = -pi / root * cpr
- ri = pi / root * spr
- ar = -rr + u * ri
- ai = -ri - u * rr
- endif
- qr = (ar * (cpr - spr) * (cpr + spr) + 2.0d0 * ai * cpr * spr) / (rm * rp)
- phint = 0.5d0 * (ri / u - qr)
- endif
- return
-end function phint
-! ################################################################
-
-subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
-
- ! ******************************************************************
- ! * set up long-wavelength dielectric functions of the layers for *
- ! * the present frequency wn (in cm**-1) *
- ! ******************************************************************
-
- ! implicit none
- integer, intent(in) :: nLayer
- integer, dimension(nLayer),intent(in) :: nos
- integer, intent(in) :: osc_size
- double precision, dimension(3,osc_size),intent(in) :: osc
- !f2py depend(osc_size) osc
- double precision, dimension(nLayer),intent(in) :: epsinf
- double precision, intent(in) :: wn
- double complex, dimension(nLayer), intent(out) :: eps
- !f2py depend(nLayer) nos, epsinf, eps
-
- double complex :: nomi, deno, addeps
- double precision :: wn2, b
- integer :: flag
- logical debugFirstRun
- common /control/ debugFirstRun
-
- j = 0
- do l = 1, nLayer ! loop over different thin film layers
- m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list
- nomi = dcmplx(1.0d0, 0.0d0)
- deno = dcmplx(1.0d0, 0.0d0)
- addeps = dcmplx(0.0d0, 0.0d0)
- wn2 = wn**2
- do k = 1, m ! loop over all TO modes
- j = j + 1
- if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosa' form
- b = wn/osc(1, j+m)
- nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) )
- deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) )
- else ! Negative TO mode means: treat as term added to epsilon
- if (osc(3,j) > 0) then ! it is a Lorentz oscillator
- addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j))
- else ! it is a Drude term
- addeps = addeps - osc(1,j)**2/dcmplx(wn2, -1*wn*osc(3,j))
- end if
- end if
- enddo
- j = j+m ! we have already looped over the LO modes, therefore increase the index
- eps(l) = epsinf(l) * nomi / deno + addeps
- enddo
- debugFirstRun = .false.
- return
-end subroutine seteps
-
-end module mod_doeels
\ No newline at end of file
diff --git a/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so b/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so
index 59bf896c1432e147d62c4476225db64664fc9923..dcb69b0d97b9e04ad9db972c13d186e86da59d08 100755
Binary files a/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so and b/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so differ