-
Notifications
You must be signed in to change notification settings - Fork 0
/
DLS_Energy_dispersive_curve_loader.py
109 lines (87 loc) · 3.68 KB
/
DLS_Energy_dispersive_curve_loader.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#----------------------BEGINNING OF IMPORTS--------------------------
import numpy as np
import h5py
import matplotlib.pyplot as plt
from scipy import constants
import scipy
from scipy.interpolate import interp1d
#--------------------------END OF IMPORTS-----------------------------
#Defining physical constants
pi = scipy.pi
Me = constants.value(u'electron mass')
hbar = constants.value(u'Planck constant over 2 pi')
e = constants.value(u'elementary charge')
#Set fermi energy (can be done by seeing where the energy dispersive curve reaches a minimum)
fermiE = 0
#Set filename (must be in same directory)
filename = "i05-1-14626.nxs"
f = h5py.File(filename, "r")
dsIntensity = f["/entry1/analyser/data"][0]
dsEnergyValues = f["/entry1/analyser/energies"]
dsAngles = f["/entry1/analyser/angles"]
#Fermi energy is repeated across array so that it can be subtracted from the energy values dataset
dsFermiEnergy = np.repeat(fermiE, len(dsEnergyValues))
dsEnergyValuesShifted = dsEnergyValues - dsFermiEnergy
#Mean energy is calculated so an estimate of the momentum can be calculated
meanE = np.mean(dsEnergyValues)
dsk = (2*Me*e*meanE)**(0.5) * 0.0173648 * dsAngles[:] * ( 3.925 * (10 ** -10) /pi) /hbar
#Intensity data is integrated across all angles/momentum
dsDOSfull = np.mean(dsIntensity, axis=0)
#Intensity data is normalised so maximum is 1
dsDOSfull = dsDOSfull / np.max(dsDOSfull)
#The data is selected to be plotted
xs = dsEnergyValuesShifted
ys = dsDOSfull
#Energy dispersive curve is plotted across full angular range
fig, (ax0) = plt.subplots(nrows=1,figsize=(8,5))
ax0.set_title("Energy Dispersive Curve integrated along the entirety of the Brillouin Zone cut")
plt.plot(xs, ys, color="xkcd:turquoise")
plt.plot((0,0),(0, 1.1), linestyle = 'dashed', color="black")
plt.ylabel("Normalised Intensity")
plt.xlabel("$E-E_{F}$ ($eV$)")
#Following lines can be included to improve aesthetics of plot
#plt.xlim((-0.8,0.2))
#plt.ylim((0,1.1))
#Following lines can be included to improve aesthetics of x ticks
#xticks = np.arange(-0.8,0.21,0.1)
#plt.xticks(xticks)
#plt.minorticks_on()
plt.tight_layout()
plt.show()
#Figure is created to be used for different momentum cuts
fig, (ax0) = plt.subplots(nrows=1,figsize=(8,5))
#Momentum range is selected by the lowest value, highest value and the step
kmin = 0
kmax = 2
kstep = 0.4
kRange = np.arange(kmin,kmax,kstep)
#The momentum range is iterated over, each time calculating the energy dispersive curve for that range
#These curves are then plotted on top of each other on a single figure
for i in kRange:
kMid = np.around(i,1)
kHalfRange = kstep/2
kMin = np.around(kMid - kHalfRange, 1)
kMax = np.around(kMid + kHalfRange, 1)
kMinIndex = (np.abs(dsk - kMin)).argmin()
kMaxIndex = (np.abs(dsk - kMax)).argmin()
dsIntensitySlice = dsIntensity[kMinIndex:kMaxIndex,]
dsDOSslice = np.mean(dsIntensitySlice, axis=0)
xs = dsEnergyValuesShifted
ys = dsDOSslice
plt.plot(xs, ys, label=(str(kMin) + " to " + str(kMax)))
#Aesthetics are selected for plot of different momentum cuts
plt.title("Energy Dispersive Curves integrated along different regions of the Brillouin Zone cut")
plt.ylabel("Normalised Intensity")
plt.plot((0,0),(0, 1.1), linestyle = 'dashed', color="black")
plt.xlabel("$E-E_{F}$ ($eV$)")
plt.legend(title="Range of k ($\pi / a$):")
#Following lines can be included to improve aesthetics of plot
#plt.xlim((-0.8,0.2))
#plt.ylim((0,1.1))
#Following lines can be included to improve aesthetics of x ticks
#xticks = np.arange(-0.8,0.21,0.1)
#plt.xticks(xticks)
#plt.minorticks_on()
fig.tight_layout()
plt.tight_layout()
plt.show()