-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMS_FeSe.py
133 lines (126 loc) · 5.33 KB
/
MS_FeSe.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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
import sys
sys.path.append('/dls_sw/i16/software/python/DMS_Code')
sys.path.append('/dls_sw/apps/scisoftpy/2.7')
#===============================================================================
# Gareth Nisbet Diamond Light Source - 25 June 2013
#===============================================================================
import os
import sys
import subprocess
from threading import Thread
import sys
from mpl_toolkits.mplot3d import Axes3D
dir = os.path.dirname(__file__)
sys.path.append(os.path.dirname(__file__))
sys.path.append(os.path.dirname(__file__)+'calcms')
import matplotlib.pyplot as plt
from numpy import linalg as LA
from matplotlib.pyplot import step, legend, xlim, ylim
from calcms import ts
import numpy as np
import scisoftpy.io as do
from scipy import ndimage
cif_file='FeSe_Cmma.cif'
hkllabels=0
hkl=[[3,3,0]]
hklint=hkl
azir=[1,0,0]
PV=[1,0]
energyrange=[7,7.2]
xmin, xmax =[-90,90]
numsteps=50
#===============================================================================
# DMS Calculation
#===============================================================================
mslist=[[np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN]] #full
################## Generate Reflist from Cif ###################################
SF, reflist, lattice , structure= ts.loadcif(cif_file,energyrange[-1])
refindex=~np.isnan(ts.vfind(reflist,np.round(hkl)-reflist).vindex())
SF=SF[refindex]
reflist=reflist[refindex]
SF=SF[ts.vfind(reflist,np.round(hkl)-reflist).vindex()]*SF
reflist=np.squeeze(reflist[np.where(SF!=0),:])
SF=np.squeeze(SF[np.where(SF!=0),:])
SF2=SF[ts.vfind(reflist,np.round(hkl)-reflist).vindex()]
loopnum=1
full, pv1, pv2, sfonly, pv1xsf1 = 0,0,0,1,0
#------------------------------------------------------------------------------
if pv1 + pv2 + sfonly + full >= 2:
print('Choose only one intensity option')
sys.exit()
elif pv1 + pv2 + sfonly + full + pv1xsf1 == 0:
print('Geometry Only')
mslist=[[np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN,np.NAN]]
for enval in np.linspace(energyrange[0],energyrange[1],numsteps):
print(str(loopnum) + ' of ' +str(numsteps))
#===========================================================================
# SF0*Gauss*SF1*SF2*PV2
#===========================================================================
if full:
calcstr='SF1*SF2*PV2'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir,SF,SF2)#[:,[3,4,5]]
polfull=ms.polfull(PV)
mslist=np.concatenate((mslist,ms.polfull(PV)),0)
#===========================================================================
# PV1 only
#===========================================================================
elif pv1:
calcstr='PV1'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir,SF,SF2)
mslist=np.concatenate((mslist,ms.pol1only(PV)),0)
#===========================================================================
# PV2 only
#===========================================================================
elif pv2:
calcstr='PV2'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir,SF,SF2)
mslist=np.concatenate((mslist,ms.pol2only(PV)),0)
#===========================================================================
# SF only
#===========================================================================
elif sfonly:
calcstr='SF1*SF2'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir,SF,SF2)
mslist=np.concatenate((mslist,ms.sfonly()),0)
#===========================================================================
# SF only
#===========================================================================
elif pv1xsf1:
calcstr='SF1*PV1'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir,SF)
mslist=np.concatenate((mslist,ms.pv1xsf1(PV)),0)
#===========================================================================
# Geometry only - no structure factors
#===========================================================================
else:
calcstr='Geometry Only'
ms=ts.calcms(lattice,hkl,hklint,reflist,enval,azir)
mslist=np.concatenate((mslist,ms.geometry()),0)
loopnum=loopnum+1
keepindex=np.where([~np.isnan(mslist).any(1)])[1]
mslist=mslist[keepindex,:]
#===============================================================================
# MS Plotting
#===============================================================================
fig = plt.figure(figsize=(6, 3),dpi=130)
ax = fig.add_subplot(111,axisbg=[150/255.0,243/255.0,86/255.0])
mslist=np.array(mslist)
if pv1 + pv2 + sfonly + full + pv1xsf1 != 0:
plt.scatter(mslist[:,3], mslist[:,7],c=mslist[:,-1],s=2,cmap=plt.cm.gray_r,lw = 0)
plt.scatter(mslist[:,4], mslist[:,7],c=mslist[:,-1],s=2,cmap=plt.cm.gray_r,lw = 0)
plt.colorbar()
else:
plt.scatter(mslist[:,3], mslist[:,-1],s=2,lw = 0)
plt.scatter(mslist[:,4], mslist[:,-1],s=2,lw = 0)
xlim(xmin,xmax)
ylim(energyrange[0],energyrange[1])
ax.set_xlim(xmin, xmax)
plt.title('CoCO3, hkl = ' + str(hkl)+' Incident polarization vector ' + str(str(PV)))
#plt.xlabel('Psi (deg)')
plt.xlabel(r'$\psi$ (deg)',fontsize=10)
plt.ylabel('Energy (keV)')
fig.subplots_adjust(bottom=0.2)
# plt.savefig('FeSe.svg',format='svg')
#plt.savefig('FeSe.pdf',format='pdf')
plt.show()
# ts.printfig('ps')