Skip to content
Snippets Groups Projects
Commit ff95236c authored by Tobias WEBER's avatar Tobias WEBER
Browse files

updated covariance tool

parent 5701d02e
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
tascalc/cov.py
+ 353
− 140
View file @ ff95236c
......@@ -3,85 +3,158 @@
# calculate covariance from neutron events
# @author Tobias Weber <tweber@ill.fr>
# @date 30-mar-2019
# @license see 'LICENSE' file
# @license GPLv3, see 'LICENSE' file
#
# @desc For a good explanation of the covariance matrix method, see (Arens 2015), pp. 795 and 1372.
# @desc reimplements the functionality of https://github.com/McStasMcXtrace/McCode/blob/master/tools/Legacy-Perl/mcresplot.pl
#
import os
import tas
try:
import numpy as np
import numpy.linalg as la
except ImportError:
print("Numpy could not be imported!")
exit(-1)
import numpy as np
import numpy.linalg as la
try:
np.set_printoptions(
precision = 4,
floatmode = "fixed")
except TypeError:
print("Warning: Numpy print options could not be set.")
verbose = True # console outputs
plot_results = True # show plot window
plot_neutrons = True # also plot neutron events
centre_on_Q = False # centre plots on Q or zero?
ellipse_points = 128 # number of points to draw ellipses
options = {
"verbose" : True, # console outputs
"plot_results" : True, # show plot window
"plot_neutrons" : True, # also plot neutron events
"centre_on_Q" : False, # centre plots on Q or zero?
"ellipse_points" : 128, # number of points to draw ellipses
"symsize" : 32.,
"dpi" : 600,
"use_tex" : False,
# column indices in mc file
ki_start_idx = 0 # start index of ki 3-vector
kf_start_idx = 3 # start index of kf 3-vector
wi_idx = 9 # start index of ki weight factor
wf_idx = 10 # start index of kf weight factor
# column indices in ki,kf files
"ki_start_idx" : 0, # start index of ki 3-vector
"kf_start_idx" : 3, # start index of kf 3-vector
"wi_idx" : 9, # start index of ki weight factor
"wf_idx" : 10, # start index of kf weight factor
# column indices in Q,E files
"Q_start_idx" : 0,
"E_idx" : 3,
"w_idx" : 4,
"filter_eps" : 1e-4,
}
# constants
sig2hwhm = np.sqrt(2. * np.log(2.))
sig2fwhm = 2.*sig2hwhm
#import scipy.constants as co
#hbar_in_meVs = co.Planck/co.elementary_charge*1000./2./np.pi
#E_to_k2 = 2.*co.neutron_mass/hbar_in_meVs**2. / co.elementary_charge*1000. * 1e-20
E_to_k2 = 0.482596406464 # E -> k^2, calculated with scipy, using the code above
k2_to_E = 1./E_to_k2 # k^2 -> E
#
# normalises events and filters out too low probabilities
#
def filter_events(Q, E, w):
if w.size == 0:
raise ValueError("No neutron events available.")
# normalise intensity/probability
maxw = np.max(w)
if np.abs(maxw) < np.finfo(w[0].__class__).eps:
raise ValueError("Neutron probability factors are zero.")
w /= maxw
# filter out too low probabilities
beloweps = lambda d: np.abs(d) <= options["filter_eps"]
nonzero_idx = [i for i in range(len(w)) if not beloweps(w[i])]
Q = Q[nonzero_idx]
E = E[nonzero_idx]
w = w[nonzero_idx]
if w.size == 0:
raise ValueError("No neutron events left after filtering.")
return [Q, E, w]
#
# loads a list of neutron events in the [ ki_vec, kf_vec, pos_vec, wi, wf ] format
#
def load_events(filename):
def load_events_kikf(filename):
dat = np.loadtxt(filename)
ki = dat[:, ki_start_idx:ki_start_idx+3]
kf = dat[:, kf_start_idx:kf_start_idx+3]
wi = dat[:, wi_idx]
wf = dat[:, wf_idx]
ki = dat[:, options["ki_start_idx"]:options["ki_start_idx"]+3]
kf = dat[:, options["kf_start_idx"]:options["kf_start_idx"]+3]
wi = dat[:, options["wi_idx"]]
wf = dat[:, options["wf_idx"]]
w = wi * wf
Q = ki - kf
E = k2_to_E * (np.multiply(ki, ki).sum(1) - np.multiply(kf, kf).sum(1))
E = tas.k2_to_E * (np.multiply(ki, ki).sum(1) - np.multiply(kf, kf).sum(1))
return [Q, E, w]
return filter_events(Q, E, w)
#
# loads a list of neutron events in the [ h, k, l, E, w ] format
#
def load_events_QE(filename):
dat = np.loadtxt(filename)
Q = dat[:, options["Q_start_idx"]:options["Q_start_idx"]+3]
E = dat[:, options["E_idx"]]
w = dat[:, options["w_idx"]]
return filter_events(Q, E, w)
#
# calculates the covariance matrix of the (Q, E) 4-vectors
#
def calc_covar(Q, E, w):
def calc_covar(Q, E, w, Qpara, Qperp):
# make a [Q, E] 4-vector
Q4 = np.insert(Q, 3, E, axis=1)
# calculate the mean Q 4-vector
Qmean = [ np.average(Q4[:,i], weights = w) for i in range(4) ]
if verbose:
if options["verbose"]:
print("Mean (Q, E) vector in lab system:\n%s\n" % Qmean)
# get the weighted covariance matrix
Qcov = np.cov(Q4, rowvar = False, aweights = w, ddof = 0)
if verbose:
Qcov = np.cov(Q4, rowvar = False, aweights = np.abs(w), ddof = 0)
if options["verbose"]:
print("Covariance matrix in lab system:\n%s\n" % Qcov)
# the resolution is the inverse of the covariance
Qres = la.inv(Qcov)
if verbose:
if options["verbose"]:
print("Resolution matrix in lab system:\n%s\n" % Qres)
# create a matrix to transform into the coordinate system with Q along x
Qnorm = Qmean[0:3] / la.norm(Qmean[0:3])
Qup = np.array([0, 1, 0])
Qside = np.cross(Qup, Qnorm)
# choose given coordinate system
if len(Qpara) == 3 and len(Qperp) == 3:
Qnorm = Qpara / la.norm(Qpara)
Qside = Qperp / la.norm(Qperp)
Qup = np.cross(Qnorm, Qside)
else:
Qnorm = Qmean[0:3] / la.norm(Qmean[0:3])
Qup = np.array([0, 1, 0])
Qside = np.cross(Qup, Qnorm)
if options["verbose"]:
print("Qpara = %s\nQperp = %s\nQup = %s\n" % (Qnorm, Qside, Qup))
# trafo matrix
T = np.transpose(np.array([
......@@ -90,29 +163,32 @@ def calc_covar(Q, E, w):
np.insert(Qup, 3, 0),
[0, 0, 0, 1] ]))
if verbose:
if options["verbose"]:
print("Transformation into (Qpara, Qperp, Qup, E) system:\n%s\n" % T)
# transform mean Q vector
Qmean_Q = np.dot(np.transpose(T), Qmean)
if verbose:
if options["verbose"]:
print("Mean (Q, E) vector in (Qpara, Qperp, Qup, E) system:\n%s\n" % Qmean_Q)
# transform the covariance matrix
Qcov_Q = np.dot(np.transpose(T), np.dot(Qcov, T))
if verbose:
if options["verbose"]:
print("Covariance matrix in (Qpara, Qperp, Qup, E) system:\n%s\n" % Qcov_Q)
# the resolution is the inverse of the covariance
Qres_Q = la.inv(Qcov_Q)
if verbose:
if options["verbose"]:
print("Resolution matrix in (Qpara, Qperp, Qup, E) system:\n%s\n" % Qres_Q)
#[ evals, evecs ] = la.eig(Qcov_Q)
#print("Ellipsoid fwhm radii:\n%s\n" % (np.sqrt(np.abs(evals)) * sig2fwhm))
# transform all neutron events
Q4_Q = np.array([])
if plot_neutrons:
if options["plot_neutrons"]:
Q4_Q = np.dot(Q4, T)
if not centre_on_Q:
if not options["centre_on_Q"]:
Q4_Q -= Qmean_Q
......@@ -121,11 +197,13 @@ def calc_covar(Q, E, w):
#
# calculates the characteristics of a given ellipse
# calculates the characteristics of a given ellipse by principal axis trafo
#
def descr_ellipse(quadric):
[ evals, evecs ] = la.eig(quadric)
fwhms = 1./np.sqrt(evals) * sig2fwhm
#print("Evals: %s" % evals)
fwhms = 1./np.sqrt(np.abs(evals)) * sig2fwhm
angles = np.array([])
if len(quadric) == 2:
......@@ -135,114 +213,171 @@ def descr_ellipse(quadric):
#
# projects along one axis of the quadric
#
def proj_quad(_E, idx):
E = np.delete(np.delete(_E, idx, axis=0), idx, axis=1)
if np.abs(_E[idx, idx]) < 1e-8:
return E
v = (_E[idx,:] + _E[:,idx]) * 0.5
vv = np.outer(v, v) / _E[idx, idx]
vv = np.delete(np.delete(vv, idx, axis=0), idx, axis=1)
return E - vv
#
# describes the ellipsoid by a principal axis trafo and by 2d cuts
#
def calc_ellipses(Qres_Q):
if verbose:
print("4d resolution ellipsoid diagonal elements fwhm (coherent-elastic scattering) lengths:\n%s\n" \
% (1./np.sqrt(np.diag(Qres_Q)) * sig2fwhm))
# 4d ellipsoid
[fwhms, angles, rot] = descr_ellipse(Qres_Q)
if verbose:
print("4d resolution ellipsoid principal axes fwhm lengths:\n%s\n" % fwhms)
axes_to_delete = [ [2, 1], [2, 0], [1,0], [3,2] ]
slice_first = [ True, True, True, False ]
results = []
for ellidx in range(len(axes_to_delete)):
# sliced 2d ellipse
Qres = np.delete(np.delete(Qres_Q, axes_to_delete[ellidx][0], axis=0), axes_to_delete[ellidx][0], axis=1)
Qres = np.delete(np.delete(Qres, axes_to_delete[ellidx][1], axis=0), axes_to_delete[ellidx][1], axis=1)
[fwhms, angles, rot] = descr_ellipse(Qres)
# projected 2d ellipse
if slice_first[ellidx]:
Qres_proj = np.delete(np.delete(Qres_Q, axes_to_delete[ellidx][0], axis=0), axes_to_delete[ellidx][0], axis=1)
Qres_proj = proj_quad(Qres_proj, axes_to_delete[ellidx][1])
else:
Qres_proj = proj_quad(Qres_Q, axes_to_delete[ellidx][0])
Qres_proj = np.delete(np.delete(Qres_proj, axes_to_delete[ellidx][1], axis=0), axes_to_delete[ellidx][1], axis=1)
[fwhms_proj, angles_proj, rot_proj] = descr_ellipse(Qres_proj)
results.append({ "fwhms" : fwhms, "angles" : angles, "rot" : rot,
"fwhms_proj" : fwhms_proj, "angles_proj" : angles_proj, "rot_proj" : rot_proj })
if options["verbose"]:
print("4d resolution ellipsoid diagonal elements fwhm (coherent-elastic scattering) lengths:\n%s\n" \
% (1./np.sqrt(np.abs(np.diag(Qres_Q))) * sig2fwhm))
print("4d resolution ellipsoid principal axes fwhm: %s" % fwhms)
# 2d sliced ellipses
Qres_QxE = np.delete(np.delete(Qres_Q, 2, axis=0), 2, axis=1)
Qres_QxE = np.delete(np.delete(Qres_QxE, 1, axis=0), 1, axis=1)
[fwhms_QxE, angles_QxE, rot_QxE] = descr_ellipse(Qres_QxE)
if verbose:
print("2d Qx/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QxE, angles_QxE[0]))
Qres_proj = proj_quad(proj_quad(proj_quad(Qres_Q, 2), 1), 0)
print("Incoherent-elastic fwhm: %.4f meV\n" % (1./np.sqrt(np.abs(Qres_proj[0,0])) * sig2fwhm))
Qres_QyE = np.delete(np.delete(Qres_Q, 2, axis=0), 2, axis=1)
Qres_QyE = np.delete(np.delete(Qres_QyE, 0, axis=0), 0, axis=1)
[fwhms_QyE, angles_QyE, rot_QyE] = descr_ellipse(Qres_QyE)
if verbose:
print("2d Qy/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QyE, angles_QyE[0]))
print("Qx,E sliced ellipse fwhm and slope angle: %s, %.4f" % (results[0]["fwhms"], results[0]["angles"][0]))
print("Qy,E sliced ellipse fwhm and slope angle: %s, %.4f" % (results[1]["fwhms"], results[1]["angles"][0]))
print("Qz,E sliced ellipse fwhm and slope angle: %s, %.4f" % (results[2]["fwhms"], results[2]["angles"][0]))
print("Qx,Qy sliced ellipse fwhm and slope angle: %s, %.4f" % (results[3]["fwhms"], results[3]["angles"][0]))
print()
print("Qx,E projected ellipse fwhm and slope angle: %s, %.4f" % (results[0]["fwhms_proj"], results[0]["angles_proj"][0]))
print("Qy,E projected ellipse fwhm and slope angle: %s, %.4f" % (results[1]["fwhms_proj"], results[1]["angles_proj"][0]))
print("Qz,E projected ellipse fwhm and slope angle: %s, %.4f" % (results[2]["fwhms_proj"], results[2]["angles_proj"][0]))
print("Qx,Qy projected ellipse fwhm and slope angle: %s, %.4f" % (results[3]["fwhms_proj"], results[3]["angles_proj"][0]))
Qres_QzE = np.delete(np.delete(Qres_Q, 1, axis=0), 1, axis=1)
Qres_QzE = np.delete(np.delete(Qres_QzE, 0, axis=0), 0, axis=1)
[fwhms_QzE, angles_QzE, rot_QzE] = descr_ellipse(Qres_QzE)
if verbose:
print("2d Qz/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QzE, angles_QzE[0]))
Qres_QxQy = np.delete(np.delete(Qres_Q, 3, axis=0), 3, axis=1)
Qres_QxQy = np.delete(np.delete(Qres_QxQy, 2, axis=0), 2, axis=1)
[fwhms_QxQy, angles_QxQy, rot_QxQy] = descr_ellipse(Qres_QxQy)
if verbose:
print("2d Qx/Qy ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QxQy, angles_QxQy[0]))
return results
return [fwhms_QxE, rot_QxE, fwhms_QyE, rot_QyE, fwhms_QzE, rot_QzE, fwhms_QxQy, rot_QxQy]
#
# shows the 2d ellipses
#
def plot_ellipses(file, Q4, Qmean, fwhms_QxE, rot_QxE, fwhms_QyE, rot_QyE, fwhms_QzE, rot_QzE, fwhms_QxQy, rot_QxQy):
import matplotlib.pyplot as plot
def plot_ellipses(file, Q4, w, Qmean, ellis):
try:
import mpl_toolkits.mplot3d as mplot3d
import matplotlib
import matplotlib.pyplot as plot
except ImportError:
print("Matplotlib could not be imported!")
exit(-1)
matplotlib.rc("text", usetex=options["use_tex"])
thesymsize = options["symsize"] * w
themarker = "."
ellfkt = lambda rad, rot, phi, Qmean2d : \
np.dot(rot, np.array([ rad[0]*np.cos(phi), rad[1]*np.sin(phi) ])) + Qmean2d
# centre plots on zero or mean Q vector ?
QxE = np.array([[0], [0]])
QyE = np.array([[0], [0]])
QzE = np.array([[0], [0]])
QxQy = np.array([[0], [0]])
# 2d plots
fig = plot.figure()
if centre_on_Q:
QxE = np.array([[Qmean[0]], [Qmean[3]]])
QyE = np.array([[Qmean[1]], [Qmean[3]]])
QzE = np.array([[Qmean[2]], [Qmean[3]]])
QxQy = np.array([[Qmean[0]], [Qmean[1]]])
num_ellis = len(ellis)
coord_axes = [[0,3], [1,3], [2,3], [0,1]]
coord_names = ["Qpara (1/A)", "Qperp (1/A)", "Qup (1/A)", "E (meV)"]
if options["use_tex"]:
coord_names[0] = "$Q_{\parallel}$ (\AA$^{-1}$)"
coord_names[1] = "$Q_{\perp}$ (\AA$^{-1}$)"
coord_names[2] = "$Q_{up}$ (\AA$^{-1}$)"
phi = np.linspace(0, 2.*np.pi, ellipse_points)
ell_QxE = ellfkt(fwhms_QxE, rot_QxE, phi, QxE)
ell_QyE = ellfkt(fwhms_QyE, rot_QyE, phi, QyE)
ell_QzE = ellfkt(fwhms_QzE, rot_QzE, phi, QzE)
ell_QxQy = ellfkt(fwhms_QxQy, rot_QxQy, phi, QxQy)
fig = plot.figure()
subplot_QxE = fig.add_subplot(221)
subplot_QxE.set_xlabel("Qpara (1/A)")
subplot_QxE.set_ylabel("E (meV)")
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot_QxE.scatter(Q4[:, 0], Q4[:, 3], s=0.25)
subplot_QxE.plot(ell_QxE[0], ell_QxE[1], c="black")
ellplots = []
for ellidx in range(num_ellis):
# centre plots on zero or mean Q vector ?
QxE = np.array([[0], [0]])
subplot_QyE = fig.add_subplot(222)
subplot_QyE.set_xlabel("Qperp (1/A)")
subplot_QyE.set_ylabel("E (meV)")
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot_QyE.scatter(Q4[:, 1], Q4[:, 3], s=0.25)
subplot_QyE.plot(ell_QyE[0], ell_QyE[1], c="black")
if options["centre_on_Q"]:
QxE = np.array([[Qmean[coord_axes[ellidx][0]]], [Qmean[coord_axes[ellidx][0]]]])
subplot_QzE = fig.add_subplot(223)
subplot_QzE.set_xlabel("Qup (1/A)")
subplot_QzE.set_ylabel("E (meV)")
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot_QzE.scatter(Q4[:, 2], Q4[:, 3], s=0.25)
subplot_QzE.plot(ell_QzE[0], ell_QzE[1], c="black")
subplot_QxQy = fig.add_subplot(224)
subplot_QxQy.set_xlabel("Qpara (1/A)")
subplot_QxQy.set_ylabel("Qperp (meV)")
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot_QxQy.scatter(Q4[:, 0], Q4[:, 1], s=0.25)
subplot_QxQy.plot(ell_QxQy[0], ell_QxQy[1], c="black")
phi = np.linspace(0, 2.*np.pi, options["ellipse_points"])
ell_QxE = ellfkt(ellis[ellidx]["fwhms"]*0.5, ellis[ellidx]["rot"], phi, QxE)
ell_QxE_proj = ellfkt(ellis[ellidx]["fwhms_proj"]*0.5, ellis[ellidx]["rot_proj"], phi, QxE)
ellplots.append({"sliced":ell_QxE, "proj":ell_QxE_proj})
subplot_QxE = fig.add_subplot(221 + ellidx)
subplot_QxE.set_xlabel(coord_names[coord_axes[ellidx][0]])
subplot_QxE.set_ylabel(coord_names[coord_axes[ellidx][1]])
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot_QxE.scatter(Q4[:, coord_axes[ellidx][0]], Q4[:, coord_axes[ellidx][1]], marker=themarker, s=thesymsize)
subplot_QxE.plot(ell_QxE[0], ell_QxE[1], c="black", linestyle="dashed")
subplot_QxE.plot(ell_QxE_proj[0], ell_QxE_proj[1], c="black", linestyle="solid")
plot.tight_layout()
# 3d plot
fig3d = plot.figure()
subplot3d = fig3d.add_subplot(111, projection="3d")
subplot3d.set_xlabel(coord_names[0])
subplot3d.set_ylabel(coord_names[1])
subplot3d.set_zlabel(coord_names[3])
if len(Q4.shape)==2 and len(Q4)>0 and len(Q4[0])==4:
subplot3d.scatter(Q4[:,0], Q4[:,1], Q4[:,3], marker=themarker, s=thesymsize)
# xE
subplot3d.plot(ellplots[0]["sliced"][0], ellplots[0]["sliced"][1], zs=0., zdir="y", c="black", linestyle="dashed")
subplot3d.plot(ellplots[0]["proj"][0], ellplots[0]["proj"][1], zs=0., zdir="y", c="black", linestyle="solid")
# yE
subplot3d.plot(ellplots[1]["sliced"][0], ellplots[1]["sliced"][1], zs=0., zdir="x", c="black", linestyle="dashed")
subplot3d.plot(ellplots[1]["proj"][0], ellplots[1]["proj"][1], zs=0., zdir="x", c="black", linestyle="solid")
# xy
subplot3d.plot(ellplots[3]["sliced"][0], ellplots[3]["sliced"][1], zs=0., zdir="z", c="black", linestyle="dashed")
subplot3d.plot(ellplots[3]["proj"][0], ellplots[3]["proj"][1], zs=0., zdir="z", c="black", linestyle="solid")
if file != "":
if verbose:
print("Saving plot to \"%s\"." % file)
plot.savefig(file)
splitext = os.path.splitext(file)
file3d = splitext[0] + "_3d" + splitext[1]
if options["verbose"]:
print("Saving 2d plot to \"%s\"." % file)
print("Saving 3d plot to \"%s\"." % file3d)
fig.savefig(file, dpi=options["dpi"])
fig3d.savefig(file3d, dpi=options["dpi"])
if plot_results:
if options["plot_results"]:
plot.show()
......@@ -261,45 +396,123 @@ def check_versions():
#
# main
# entry point
#
if __name__ == "__main__":
def run_cov():
print("This is a covariance matrix calculator using neutron events,\n\twritten by T. Weber <tweber@ill.fr>, 30 March 2019.\n")
check_versions()
import argparse as arg
try:
import argparse as arg
except ImportError:
print("Argparse could not be imported!")
exit(-1)
args = arg.ArgumentParser(description="Calculates the covariance matrix of neutron scattering events.")
args.add_argument("file", type=str, help="input file")
args.add_argument("-s", "--save", default="", type=str, nargs="?", help="save plot to file")
args.add_argument("--ellipse", default=ellipse_points, type=int, nargs="?", help="number of points to draw ellipses")
args.add_argument("--ki", default=ki_start_idx, type=int, nargs="?", help="index of ki vector's first column in file")
args.add_argument("--kf", default=kf_start_idx, type=int, nargs="?", help="index of kf vector's first column in file")
args.add_argument("--wi", default=wi_idx, type=int, nargs="?", help="index of ki weight factor column in file")
args.add_argument("--wf", default=wf_idx, type=int, nargs="?", help="index of kf weight factor column in file")
args.add_argument("--ellipse", default=options["ellipse_points"], type=int, nargs="?", help="number of points to draw ellipses")
args.add_argument("--ki", default=options["ki_start_idx"], type=int, nargs="?", help="index of ki vector's first column in kikf file")
args.add_argument("--kf", default=options["kf_start_idx"], type=int, nargs="?", help="index of kf vector's first column in kikf file")
args.add_argument("--wi", default=options["wi_idx"], type=int, nargs="?", help="index of ki weight factor column in kikf file")
args.add_argument("--wf", default=options["wf_idx"], type=int, nargs="?", help="index of kf weight factor column in kikf file")
args.add_argument("--w", default=options["w_idx"], type=int, nargs="?", help="index of neutron weight factor column in QE file")
args.add_argument("--Q", default=options["Q_start_idx"], type=int, nargs="?", help="index of Q vector's first column in QE file")
args.add_argument("--E", default=options["E_idx"], type=int, nargs="?", help="index of E column in QE file")
args.add_argument("--QEfile", action="store_true", help="use the QE file type")
args.add_argument("--centreonQ", action="store_true", help="centre plots on mean Q")
args.add_argument("--noverbose", action="store_true", help="don't show console logs")
args.add_argument("--noplot", action="store_true", help="don't show plot window")
args.add_argument("--noneutrons", action="store_true", help="don't show mc neutrons in plot")
args.add_argument("--noplot", action="store_true", help="don't show any plot windows")
args.add_argument("--noneutrons", action="store_true", help="don't show neutron events in plots")
args.add_argument("--symsize", default=options["symsize"], type=float, nargs="?", help="size of the symbols in plots")
args.add_argument("--ax", default=None, type=float, nargs="?", help="x component of first orientation vector")
args.add_argument("--ay", default=None, type=float, nargs="?", help="y component of first orientation vector")
args.add_argument("--az", default=None, type=float, nargs="?", help="z component of first orientation vector")
args.add_argument("--bx", default=None, type=float, nargs="?", help="x component of second orientation vector")
args.add_argument("--by", default=None, type=float, nargs="?", help="y component of second orientation vector")
args.add_argument("--bz", default=None, type=float, nargs="?", help="z component of second orientation vector")
args.add_argument("--a", "--as", "-a", default=None, type=float, nargs="?", help="lattice constant a (only needed in case data is in rlu)")
args.add_argument("--b", "--bs", "-b", default=None, type=float, nargs="?", help="lattice constant b (only needed in case data is in rlu)")
args.add_argument("--c", "--cs", "-c", default=None, type=float, nargs="?", help="lattice constant c (only needed in case data is in rlu)")
args.add_argument("--aa", "--alpha", default=90., type=float, nargs="?", help="lattice angle alpha (only needed in case data is in rlu)")
args.add_argument("--bb", "--beta", default=90., type=float, nargs="?", help="lattice angle beta (only needed in case data is in rlu)")
args.add_argument("--cc", "--gamma", default=90., type=float, nargs="?", help="lattice angle gamma (only needed in case data is in rlu)")
args.add_argument("--filtereps", default=options["filter_eps"], type=float, nargs="?", help="epsilon probability below which neutron events are filtered out")
args.add_argument("--dpi", default=options["dpi"], type=int, nargs="?", help="DPI of output plot file")
args.add_argument("--tex", action="store_true", help="use tex in plots")
argv = args.parse_args()
verbose = (argv.noverbose==False)
plot_results = (argv.noplot==False)
plot_neutrons = (argv.noneutrons==False)
centre_on_Q = argv.centreonQ
options["verbose"] = (argv.noverbose==False)
options["plot_results"] = (argv.noplot==False)
options["plot_neutrons"] = (argv.noneutrons==False)
options["centre_on_Q"] = argv.centreonQ
options["dpi"] = argv.dpi
options["use_tex"] = argv.tex
B = []
if argv.a!=None and argv.b!=None and argv.c!=None and argv.aa!=None and argv.bb!=None and argv.cc!=None:
lattice = np.array([ argv.a, argv.b, argv.c, ])
angles = np.array([ argv.aa, argv.bb, argv.cc ]) / 180.*np.pi
B = tas.get_B(lattice, angles)
if options["verbose"]:
print("Crystal B matrix:\n%s\n" % B)
infile = argv.file
outfile = argv.save
ellipse_points = argv.ellipse
ki_start_idx = argv.ki
kf_start_idx = argv.kf
wi_idx = argv.wi
wf_idx = argv.wf
options["ellipse_points"] = argv.ellipse
options["ki_start_idx"] = argv.ki
options["kf_start_idx"] = argv.kf
options["wi_idx"] = argv.wi
options["wf_idx"] = argv.wf
options["filter_eps"] = argv.filtereps
options["symsize"] = argv.symsize
avec = [ argv.az, argv.ay, argv.az ]
bvec = [ argv.bx, argv.by, argv.bz ]
try:
# input file is in h k l E w format?
if argv.QEfile:
[Q, E, w] = load_events_QE(infile)
# convert rlu to 1/A
if len(B) != 0:
Q = np.dot(Q, np.transpose(B))
# input file is in the kix kiy kiz kfx kfy kfz wi wf format?
else:
[Q, E, w] = load_events_kikf(infile)
except OSError:
print("Could not load input file %s." % infile)
exit(-1)
except NameError:
print("Error processing input file %s." % infile)
exit(-1)
if avec[0]!=None and avec[1]!=None and avec[2]!=None and bvec[0]!=None and bvec[1]!=None and bvec[2]!=None:
Qpara = np.array(avec)
Qperp = np.array(bvec)
# convert rlu to 1/A
if len(B) != 0:
Qpara = np.dot(B, Qpara)
Qperp = np.dot(B, Qperp)
else:
Qpara = np.array([])
Qperp = np.array([])
[Qres, Q4, Qmean] = calc_covar(Q, E, w, Qpara, Qperp)
calcedellis = calc_ellipses(Qres)
print("")
[Q, E, w] = load_events(infile)
if options["plot_results"] or outfile!="":
plot_ellipses(outfile, Q4, w, Qmean, calcedellis)
[Qres, Q4, Qmean] = calc_covar(Q, E, w)
[fwhms_QxE, rot_QxE, fwhms_QyE, rot_QyE, fwhms_QzE, rot_QzE, fwhms_QxQy, rot_QxQy] = calc_ellipses(Qres)
if plot_results or outfile!="":
plot_ellipses(outfile, Q4, Qmean, fwhms_QxE, rot_QxE, fwhms_QyE, rot_QyE, fwhms_QzE, rot_QzE, fwhms_QxQy, rot_QxQy)
#
# main
#
if __name__ == "__main__":
run_cov()
tascalc/qtas.py deleted 100644 → 0
+ 0
− 517
View file @ 5701d02e
#
# tascalc gui
# @author Tobias Weber <tweber@ill.fr>
# @date 24-oct-18
# @license see 'LICENSE' file
#
# __ident__ : TAS Calculator
# __descr__ : Calculates triple-axis angles.
#
# -----------------------------------------------------------------------------
# dependencies
import sys
import tascalc as tas
import numpy as np
import numpy.linalg as la
# try to import qt5...
try:
import PyQt5 as qt
import PyQt5.QtCore as qtc
import PyQt5.QtWidgets as qtw
qt_ver = 5
except ImportError:
# ...and if not possible try to import qt4 instead
try:
import PyQt4 as qt
import PyQt4.QtCore as qtc
import PyQt4.QtGui as qtw
qt_ver = 4
except ImportError:
print("Error: No suitable version of Qt was found!")
exit(-1)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# main application
np.set_printoptions(suppress=True, precision=4)
app = qtw.QApplication(sys.argv)
app.setApplicationName("qtas")
#app.setStyle("Fusion")
sett = qtc.QSettings("tobis_stuff", "in20tool")
if sett.contains("mainwnd/theme"):
app.setStyle(sett.value("mainwnd/theme"))
tabs = qtw.QTabWidget()
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# variables
B = np.array([[1,0,0], [0,1,0], [0,0,1]])
orient_rlu = np.array([1,0,0])
#orient2_rlu = np.array([0,1,0])
orient_up_rlu = np.array([0,0,1])
g_eps = 1e-4
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# helpers
def getfloat(str):
try:
return float(str)
except ValueError:
return 0.
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# crystal tab
xtalpanel = qtw.QWidget()
xtallayout = qtw.QGridLayout(xtalpanel)
scpanel = xtalpanel
sclayout = xtallayout
editA = qtw.QLineEdit(xtalpanel)
editB = qtw.QLineEdit(xtalpanel)
editC = qtw.QLineEdit(xtalpanel)
editAlpha = qtw.QLineEdit(xtalpanel)
editBeta = qtw.QLineEdit(xtalpanel)
editGamma = qtw.QLineEdit(xtalpanel)
separatorXtal = qtw.QFrame(xtalpanel)
separatorXtal.setFrameStyle(qtw.QFrame.HLine)
editAx = qtw.QLineEdit(scpanel)
editAy = qtw.QLineEdit(scpanel)
editAz = qtw.QLineEdit(scpanel)
editBx = qtw.QLineEdit(scpanel)
editBy = qtw.QLineEdit(scpanel)
editBz = qtw.QLineEdit(scpanel)
editBMat = qtw.QPlainTextEdit(xtalpanel)
editBMat.setReadOnly(True)
def xtalChanged():
global B, orient_rlu, orient_up_rlu
lattice = np.array([getfloat(editA.text()), getfloat(editB.text()), getfloat(editC.text())])
angles = np.array([getfloat(editAlpha.text()), getfloat(editBeta.text()), getfloat(editGamma.text())])
orient_rlu = np.array([getfloat(editAx.text()), getfloat(editAy.text()), getfloat(editAz.text())])
orient2_rlu = np.array([getfloat(editBx.text()), getfloat(editBy.text()), getfloat(editBz.text())])
try:
B = tas.get_B(lattice, angles/180.*np.pi)
invB = la.inv(B)
metric = tas.get_metric(B)
ang = tas.angle(orient_rlu, orient2_rlu, metric)
orient_up_rlu = tas.cross(orient_rlu, orient2_rlu, B)
orient_up_rlu_norm = orient_up_rlu / la.norm(orient_up_rlu)
UB = tas.get_UB(B, orient_rlu, orient2_rlu, orient_up_rlu)
invUB = la.inv(UB)
editBMat.setPlainText("Scattering plane normal: %s rlu.\n" % str(orient_up_rlu_norm) \
+"Angle between orientation vectors 1 and 2: %.4g deg.\n" % (ang/np.pi*180.) \
+"\nB =\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n" \
% (B[0,0],B[0,1],B[0,2], B[1,0],B[1,1],B[1,2], B[2,0],B[2,1],B[2,2]) \
+"\nB^(-1) =\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n" \
% (invB[0,0],invB[0,1],invB[0,2], invB[1,0],invB[1,1],invB[1,2], invB[2,0],invB[2,1],invB[2,2]) \
+"\nUB =\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n" \
% (UB[0,0],UB[0,1],UB[0,2], UB[1,0],UB[1,1],UB[1,2], UB[2,0],UB[2,1],UB[2,2]) \
+"\n(UB)^(-1) =\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n%10.4f %10.4f %10.4f\n" \
% (invUB[0,0],invUB[0,1],invUB[0,2], invUB[1,0],invUB[1,1],invUB[1,2], invUB[2,0],invUB[2,1],invUB[2,2]) \
)
except (ArithmeticError, la.LinAlgError) as err:
editBMat.setPlainText("invalid")
QChanged()
def planeChanged():
xtalChanged()
#QChanged()
editA.textEdited.connect(xtalChanged)
editB.textEdited.connect(xtalChanged)
editC.textEdited.connect(xtalChanged)
editAlpha.textEdited.connect(xtalChanged)
editBeta.textEdited.connect(xtalChanged)
editGamma.textEdited.connect(xtalChanged)
editAx.textEdited.connect(planeChanged)
editAy.textEdited.connect(planeChanged)
editAz.textEdited.connect(planeChanged)
editBx.textEdited.connect(planeChanged)
editBy.textEdited.connect(planeChanged)
editBz.textEdited.connect(planeChanged)
editA.setText("%.6g" % sett.value("qtas/a", 5., type=float))
editB.setText("%.6g" % sett.value("qtas/b", 5., type=float))
editC.setText("%.6g" % sett.value("qtas/c", 5., type=float))
editAlpha.setText("%.6g" % sett.value("qtas/alpha", 90., type=float))
editBeta.setText("%.6g" % sett.value("qtas/beta", 90., type=float))
editGamma.setText("%.6g" % sett.value("qtas/gamma", 90., type=float))
editAx.setText("%.6g" % sett.value("qtas/ax", 1., type=float))
editAy.setText("%.6g" % sett.value("qtas/ay", 0., type=float))
editAz.setText("%.6g" % sett.value("qtas/az", 0., type=float))
editBx.setText("%.6g" % sett.value("qtas/bx", 0., type=float))
editBy.setText("%.6g" % sett.value("qtas/by", 1., type=float))
editBz.setText("%.6g" % sett.value("qtas/bz", 0., type=float))
xtallayout.addWidget(qtw.QLabel(u"a (\u212b):", xtalpanel), 0,0, 1,1)
xtallayout.addWidget(editA, 0,1, 1,3)
xtallayout.addWidget(qtw.QLabel(u"b (\u212b):", xtalpanel), 1,0, 1,1)
xtallayout.addWidget(editB, 1,1, 1,3)
xtallayout.addWidget(qtw.QLabel(u"c (\u212b):", xtalpanel), 2,0, 1,1)
xtallayout.addWidget(editC, 2,1, 1,3)
xtallayout.addWidget(qtw.QLabel(u"\u03b1 (deg):", xtalpanel), 3,0, 1,1)
xtallayout.addWidget(editAlpha, 3,1, 1,3)
xtallayout.addWidget(qtw.QLabel(u"\u03b2 (deg):", xtalpanel), 4,0, 1,1)
xtallayout.addWidget(editBeta, 4,1, 1,3)
xtallayout.addWidget(qtw.QLabel(u"\u03b3 (deg):", xtalpanel), 5,0, 1,1)
xtallayout.addWidget(editGamma, 5,1, 1,3)
xtallayout.addWidget(separatorXtal, 6,0, 1,4)
sclayout.addWidget(qtw.QLabel("Orient. 1 (rlu):", scpanel), 7,0, 1,1)
sclayout.addWidget(editAx, 7,1, 1,1)
sclayout.addWidget(editAy, 7,2, 1,1)
sclayout.addWidget(editAz, 7,3, 1,1)
sclayout.addWidget(qtw.QLabel("Orient. 2 (rlu):", scpanel), 8,0, 1,1)
sclayout.addWidget(editBx, 8,1, 1,1)
sclayout.addWidget(editBy, 8,2, 1,1)
sclayout.addWidget(editBz, 8,3, 1,1)
xtallayout.addWidget(editBMat, 9,0, 2,4)
tabs.addTab(xtalpanel, "Crystal")
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# tas tab
taspanel = qtw.QWidget()
taslayout = qtw.QGridLayout(taspanel)
Qpanel = taspanel
Qlayout = taslayout
editA1 = qtw.QLineEdit(taspanel)
editA2 = qtw.QLineEdit(taspanel)
editA3 = qtw.QLineEdit(taspanel)
editA4 = qtw.QLineEdit(taspanel)
editA5 = qtw.QLineEdit(taspanel)
editA6 = qtw.QLineEdit(taspanel)
checkA4Sense = qtw.QCheckBox(taspanel)
editDm = qtw.QLineEdit(taspanel)
editDa = qtw.QLineEdit(taspanel)
edith = qtw.QLineEdit(Qpanel)
editk = qtw.QLineEdit(Qpanel)
editl = qtw.QLineEdit(Qpanel)
editE = qtw.QLineEdit(Qpanel)
editKi = qtw.QLineEdit(Qpanel)
editKf = qtw.QLineEdit(Qpanel)
editQAbs = qtw.QLineEdit(Qpanel)
editQAbs.setReadOnly(True)
tasstatus = qtw.QLabel(taspanel)
separatorTas = qtw.QFrame(Qpanel)
separatorTas.setFrameStyle(qtw.QFrame.HLine)
separatorTas2 = qtw.QFrame(Qpanel)
separatorTas2.setFrameStyle(qtw.QFrame.HLine)
separatorTas3 = qtw.QFrame(Qpanel)
separatorTas3.setFrameStyle(qtw.QFrame.HLine)
def TASChanged():
global orient_rlu, orient_up_rlu
a1 = getfloat(editA1.text()) / 180. * np.pi
a2 = a1 * 2.
a3 = getfloat(editA3.text()) / 180. * np.pi
a4 = getfloat(editA4.text()) / 180. * np.pi
a5 = getfloat(editA5.text()) / 180. * np.pi
a6 = a5 * 2.
dmono = getfloat(editDm.text())
dana = getfloat(editDa.text())
sense_sample = 1.
if checkA4Sense.isChecked() == False:
sense_sample = -1.
editA2.setText("%.6g" % (a2 / np.pi * 180.))
editA6.setText("%.6g" % (a6 / np.pi * 180.))
try:
ki = tas.get_monok(a1, dmono)
kf = tas.get_monok(a5, dana)
E = tas.get_E(ki, kf)
Qlen = tas.get_Q(ki, kf, a4)
Qvec = tas.get_hkl(ki, kf, a3, Qlen, orient_rlu, orient_up_rlu, B, sense_sample)
edith.setText("%.6g" % Qvec[0])
editk.setText("%.6g" % Qvec[1])
editl.setText("%.6g" % Qvec[2])
editQAbs.setText("%.6g" % Qlen)
editKi.setText("%.6g" % ki)
editKf.setText("%.6g" % kf)
editE.setText("%.6g" % E)
except (ArithmeticError, la.LinAlgError) as err:
edith.setText("invalid")
editk.setText("invalid")
editl.setText("invalid")
editKi.setText("invalid")
editKf.setText("invalid")
editE.setText("invalid")
def A2Changed():
a2 = getfloat(editA2.text()) / 180. * np.pi
editA1.setText("%.6g" % (0.5*a2 / np.pi * 180.))
TASChanged()
def A6Changed():
a6 = getfloat(editA6.text()) / 180. * np.pi
editA5.setText("%.6g" % (0.5*a6 / np.pi * 180.))
TASChanged()
def DChanged():
QChanged()
def QChanged():
global orient_rlu, orient_up_rlu, g_eps
Q_rlu = np.array([getfloat(edith.text()), getfloat(editk.text()), getfloat(editl.text())])
ki = getfloat(editKi.text())
kf = getfloat(editKf.text())
try:
[a1, a2] = tas.get_a1a2(ki, getfloat(editDm.text()))
editA1.setText("%.6g" % (a1 / np.pi * 180.))
editA2.setText("%.6g" % (a2 / np.pi * 180.))
except (ArithmeticError, la.LinAlgError) as err:
editA1.setText("invalid")
editA2.setText("invalid")
try:
[a5, a6] = tas.get_a1a2(kf, getfloat(editDa.text()))
editA5.setText("%.6g" % (a5 / np.pi * 180.))
editA6.setText("%.6g" % (a6 / np.pi * 180.))
except (ArithmeticError, la.LinAlgError) as err:
editA5.setText("invalid")
editA6.setText("invalid")
try:
sense_sample = 1.
if checkA4Sense.isChecked() == False:
sense_sample = -1.
[a3, a4, dist_Q_plane] = tas.get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, sense_sample)
Qlen = tas.get_Q(ki, kf, a4)
Q_in_plane = np.abs(dist_Q_plane) < g_eps
editA3.setText("%.6g" % (a3 / np.pi * 180.))
editA4.setText("%.6g" % (a4 / np.pi * 180.))
editQAbs.setText("%.6g" % Qlen)
if Q_in_plane:
tasstatus.setText("")
else:
tasstatus.setText(u"WARNING: Q is out of the plane by %.4g / \u212b!" % dist_Q_plane)
except (ArithmeticError, la.LinAlgError) as err:
editA3.setText("invalid")
editA4.setText("invalid")
def KiKfChanged():
ki = getfloat(editKi.text())
kf = getfloat(editKf.text())
try:
E = tas.get_E(ki, kf)
editE.setText("%.6g" % E)
QChanged()
except (ArithmeticError, la.LinAlgError) as err:
editE.setText("invalid")
def EChanged():
E = getfloat(editE.text())
kf = getfloat(editKf.text())
try:
ki = tas.get_ki(kf, E)
editKi.setText("%.6g" % ki)
QChanged()
except (ArithmeticError, la.LinAlgError) as err:
editKi.setText("invalid")
editA1.textEdited.connect(TASChanged)
editA3.textEdited.connect(TASChanged)
editA4.textEdited.connect(TASChanged)
editA5.textEdited.connect(TASChanged)
editA2.textEdited.connect(A2Changed)
editA6.textEdited.connect(A6Changed)
editDm.textEdited.connect(DChanged)
editDa.textEdited.connect(DChanged)
edith.textEdited.connect(QChanged)
editk.textEdited.connect(QChanged)
editl.textEdited.connect(QChanged)
editKi.textEdited.connect(KiKfChanged)
editKf.textEdited.connect(KiKfChanged)
editE.textEdited.connect(EChanged)
editDm.setText("%.6g" % sett.value("qtas/dm", 3.355, type=float))
editDa.setText("%.6g" % sett.value("qtas/da", 3.355, type=float))
edith.setText("%.6g" % sett.value("qtas/h", 1., type=float))
editk.setText("%.6g" % sett.value("qtas/k", 0., type=float))
editl.setText("%.6g" % sett.value("qtas/l", 0., type=float))
#editE.setText("%.6g" % sett.value("qtas/E", 0., type=float))
editKi.setText("%.6g" % sett.value("qtas/ki", 2.662, type=float))
editKf.setText("%.6g" % sett.value("qtas/kf", 2.662, type=float))
checkA4Sense.setText("a4 sense is counter-clockwise")
checkA4Sense.setChecked(sett.value("qtas/a4_sense", 1, type=bool))
checkA4Sense.stateChanged.connect(QChanged)
Qlayout.addWidget(qtw.QLabel("h (rlu):", Qpanel), 0,0, 1,1)
Qlayout.addWidget(edith, 0,1, 1,2)
Qlayout.addWidget(qtw.QLabel("k (rlu):", Qpanel), 1,0, 1,1)
Qlayout.addWidget(editk, 1,1, 1,2)
Qlayout.addWidget(qtw.QLabel("l (rlu):", Qpanel), 2,0, 1,1)
Qlayout.addWidget(editl, 2,1, 1,2)
Qlayout.addWidget(qtw.QLabel("E (meV):", Qpanel), 3,0, 1,1)
Qlayout.addWidget(editE, 3,1, 1,2)
Qlayout.addWidget(qtw.QLabel(u"ki, kf (1/\u212b):", Qpanel), 4,0, 1,1)
Qlayout.addWidget(editKi, 4,1, 1,1)
Qlayout.addWidget(editKf, 4,2, 1,1)
Qlayout.addWidget(qtw.QLabel(u"|Q| (1/\u212b):", Qpanel), 5,0, 1,1)
Qlayout.addWidget(editQAbs, 5,1, 1,2)
Qlayout.addWidget(separatorTas, 6,0,1,3)
taslayout.addWidget(qtw.QLabel("a1, a2 (deg):", taspanel), 7,0, 1,1)
taslayout.addWidget(editA1, 7,1, 1,1)
taslayout.addWidget(editA2, 7,2, 1,1)
taslayout.addWidget(qtw.QLabel("a3, a4 (deg):", taspanel), 8,0, 1,1)
taslayout.addWidget(editA3, 8,1, 1,1)
taslayout.addWidget(editA4, 8,2, 1,1)
taslayout.addWidget(qtw.QLabel("a5, a6 (deg):", taspanel), 9,0, 1,1)
taslayout.addWidget(editA5, 9,1, 1,1)
taslayout.addWidget(editA6, 9,2, 1,1)
taslayout.addWidget(separatorTas2, 10,0, 1,3)
taslayout.addWidget(qtw.QLabel("Sense:", taspanel), 11,0, 1,1)
taslayout.addWidget(checkA4Sense, 11,1, 1,2)
taslayout.addWidget(separatorTas3, 12,0, 1,3)
taslayout.addWidget(qtw.QLabel("Mono., Ana. d (A):", taspanel), 13,0, 1,1)
taslayout.addWidget(editDm, 13,1, 1,1)
taslayout.addWidget(editDa, 13,2, 1,1)
taslayout.addItem(qtw.QSpacerItem(16,16, qtw.QSizePolicy.Minimum, qtw.QSizePolicy.Expanding), 14,0, 1,3)
taslayout.addWidget(tasstatus, 15,0, 1,3)
tabs.addTab(taspanel, "TAS")
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# info/settings tab
infopanel = qtw.QWidget()
infolayout = qtw.QGridLayout(infopanel)
comboA3 = qtw.QComboBox(infopanel)
comboA3.addItems(["Takin", "NOMAD", "SICS", "NICOS"])
a3_offs = [np.pi/2., np.pi, 0., 0.]
comboA3.setCurrentIndex(sett.value("qtas/a3_conv", 1, type=int))
def comboA3ConvChanged():
idx = comboA3.currentIndex()
tas.set_a3_offs(a3_offs[idx])
QChanged()
comboA3.currentIndexChanged.connect(comboA3ConvChanged)
separatorInfo = qtw.QFrame(infopanel)
separatorInfo.setFrameStyle(qtw.QFrame.HLine)
infolayout.addWidget(qtw.QLabel("TAS Calculator.", infopanel), 0,0, 1,2)
infolayout.addWidget(qtw.QLabel("Written by Tobias Weber <tweber@ill.fr>.", infopanel), 1,0, 1,2)
infolayout.addWidget(qtw.QLabel("Date: October 24, 2018.", infopanel), 2,0, 1,2)
infolayout.addWidget(separatorInfo, 3,0, 1,2)
infolayout.addWidget(qtw.QLabel("Interpreter Version: " + sys.version + ".", infopanel), 4,0, 1,2)
infolayout.addWidget(qtw.QLabel("Numpy Version: " + np.__version__ + ".", infopanel), 5,0, 1,2)
infolayout.addWidget(qtw.QLabel("Qt Version: " + qtc.QT_VERSION_STR + ".", infopanel), 6,0, 1,2)
infolayout.addItem(qtw.QSpacerItem(16,16, qtw.QSizePolicy.Minimum, qtw.QSizePolicy.Expanding), 7,0, 1,2)
infolayout.addWidget(qtw.QLabel("A3 Convention:", infopanel), 8,0, 1,1)
infolayout.addWidget(comboA3, 8,1, 1,1)
tabs.addTab(infopanel, "Infos")
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# main dialog window
dlg = qtw.QDialog()
dlg.setWindowTitle("TAS Calculator")
mainlayout = qtw.QGridLayout(dlg)
mainlayout.addWidget(tabs)
if sett.contains("qtas/geo"):
geo = sett.value("qtas/geo")
if qt_ver == 4:
geo = geo.toByteArray()
dlg.restoreGeometry(geo)
xtalChanged()
KiKfChanged()
comboA3ConvChanged()
#QChanged()
dlg.show()
app.exec_()
# save settings
sett.setValue("qtas/a", getfloat(editA.text()))
sett.setValue("qtas/b", getfloat(editB.text()))
sett.setValue("qtas/c", getfloat(editC.text()))
sett.setValue("qtas/alpha", getfloat(editAlpha.text()))
sett.setValue("qtas/beta", getfloat(editBeta.text()))
sett.setValue("qtas/gamma", getfloat(editGamma.text()))
sett.setValue("qtas/ax", getfloat(editAx.text()))
sett.setValue("qtas/ay", getfloat(editAy.text()))
sett.setValue("qtas/az", getfloat(editAz.text()))
sett.setValue("qtas/bx", getfloat(editBx.text()))
sett.setValue("qtas/by", getfloat(editBy.text()))
sett.setValue("qtas/bz", getfloat(editBz.text()))
sett.setValue("qtas/dm", getfloat(editDm.text()))
sett.setValue("qtas/da", getfloat(editDa.text()))
sett.setValue("qtas/h", getfloat(edith.text()))
sett.setValue("qtas/k", getfloat(editk.text()))
sett.setValue("qtas/l", getfloat(editl.text()))
#sett.setValue("qtas/E", getfloat(editE.text()))
sett.setValue("qtas/ki", getfloat(editKi.text()))
sett.setValue("qtas/kf", getfloat(editKf.text()))
sett.setValue("qtas/a3_conv", comboA3.currentIndex())
sett.setValue("qtas/a4_sense", checkA4Sense.isChecked())
sett.setValue("qtas/geo", dlg.saveGeometry())
# -----------------------------------------------------------------------------
tascalc/tas.py 0 → 100644
+ 1083
− 0
View file @ ff95236c
This diff is collapsed.
tascalc/tascalc.py deleted 100644 → 0
+ 0
− 225
View file @ 5701d02e
#
# calculates TAS angles from rlu (part of in20tools)
# @author Tobias Weber <tweber@ill.fr>
# @date 1-aug-18
# @license see 'LICENSE' file
#
import numpy as np
import numpy.linalg as la
use_scipy = False
# -----------------------------------------------------------------------------
# choose an a3 convention
#a3_offs = 0. # for sics: a3 is angle between ki and orient1
#a3_offs = np.pi/2. # for takin: Q along orient1 => a3:=a4/2
a3_offs = np.pi # for nomad: a3 is angle between ki and orient1 + 180 deg
def set_a3_offs(offs):
global a3_offs
a3_offs = offs
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# rotate a vector around an axis using Rodrigues' formula
# see: https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
def rotate(_axis, vec, phi):
axis = _axis / la.norm(_axis)
s = np.sin(phi)
c = np.cos(phi)
return c*vec + (1.-c)*np.dot(vec, axis)*axis + s*np.cross(axis, vec)
# get metric from crystal B matrix
# basis vectors are in the columns of B, i.e. the second index
def get_metric(B):
#return np.einsum("ij,ik -> jk", B, B)
return np.dot(np.transpose(B), B)
# cross product in fractional coordinates
def cross(a, b, B):
# levi-civita in fractional coordinates
def levi(i,j,k, B):
M = np.array([B[:,i], B[:,j], B[:,k]])
return la.det(M)
metric_inv = la.inv(get_metric(B))
eps = [[[ levi(i,j,k, B) for k in range(0,3) ] for j in range(0,3) ] for i in range(0,3) ]
return np.einsum("ijk,j,k,li -> l", eps, a, b, metric_inv)
# dot product in fractional coordinates
def dot(a, b, metric):
return np.dot(a, np.dot(metric, b))
# angle between peaks in fractional coordinates
def angle(a, b, metric):
len_a = np.sqrt(dot(a, a, metric))
len_b = np.sqrt(dot(b, b, metric))
c = dot(a, b, metric) / (len_a * len_b)
return np.arccos(c)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
if use_scipy:
import scipy as sp
import scipy.constants as co
hbar_in_meVs = co.Planck/co.elementary_charge*1000./2./np.pi
E_to_k2 = 2.*co.neutron_mass/hbar_in_meVs**2. / co.elementary_charge*1000. * 1e-20
else:
E_to_k2 = 0.482596406464 # calculated with scipy, using the formula above
#print(1./E_to_k2)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# mono (or ana) k -> A1 & A2 angles (or A5 & A6)
def get_a1a2(k, d):
s = np.pi/(d*k)
a1 = np.arcsin(s)
return [a1, 2.*a1]
# a1 angle (or a5) -> mono (or ana) k
def get_monok(theta, d):
s = np.sin(theta)
k = np.pi/(d*s)
return k
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# scattering angle a4
def get_a4(ki, kf, Q):
c = (ki**2. + kf**2. - Q**2.) / (2.*ki*kf)
return np.arccos(c)
# get |Q| from ki, kf and a4
def get_Q(ki, kf, a4):
c = np.cos(a4)
return np.sqrt(ki**2. + kf**2. - c*(2.*ki*kf))
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# angle enclosed by ki and Q
def get_psi(ki, kf, Q, sense=1.):
c = (ki**2. + Q**2. - kf**2.) / (2.*ki*Q)
return sense*np.arccos(c)
# crystallographic A matrix converting fractional to lab coordinates
# see: https://de.wikipedia.org/wiki/Fraktionelle_Koordinaten
def get_A(lattice, angles):
cs = np.cos(angles)
s2 = np.sin(angles[2])
a = lattice[0] * np.array([1, 0, 0])
b = lattice[1] * np.array([cs[2], s2, 0])
c = lattice[2] * np.array([cs[1], \
(cs[0]-cs[1]*cs[2]) / s2, \
(np.sqrt(1. - np.dot(cs,cs) + 2.*cs[0]*cs[1]*cs[2])) / s2])
# testing equality with own derivation
#print((np.sqrt(1. - np.dot(cs,cs) + 2.*cs[0]*cs[1]*cs[2])) / s2)
#print(np.sqrt(1. - cs[1]*cs[1] - ((cs[0] - cs[2]*cs[1])/s2)**2.))
# the real-space basis vectors form the columns of the A matrix
return np.transpose(np.array([a, b, c]))
# crystallographic B matrix converting rlu to 1/A
# the reciprocal-space basis vectors form the columns of the B matrix
def get_B(lattice, angles):
A = get_A(lattice, angles)
B = 2.*np.pi * np.transpose(la.inv(A))
return B
# UB orientation matrix
def get_UB(B, orient1_rlu, orient2_rlu, orientup_rlu):
orient1_invA = np.dot(B, orient1_rlu)
orient2_invA = np.dot(B, orient2_rlu)
orientup_invA = np.dot(B, orientup_rlu)
orient1_invA = orient1_invA / la.norm(orient1_invA)
orient2_invA = orient2_invA / la.norm(orient2_invA)
orientup_invA = orientup_invA / la.norm(orientup_invA)
U_invA = np.array([orient1_invA, orient2_invA, orientup_invA])
UB = np.dot(U_invA, B)
return UB
# a3 & a4 angles
def get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, sense_sample=1.):
metric = get_metric(B)
# angle xi between Q and orientation reflex
xi = angle(Q_rlu, orient_rlu, metric)
# sign of xi
if dot(cross(orient_rlu, Q_rlu, B), orient_up_rlu, metric) < 0.:
xi = -xi
# length of Q
Qlen = np.sqrt(dot(Q_rlu, Q_rlu, metric))
# distance to plane
up_len = np.sqrt(dot(orient_up_rlu, orient_up_rlu, metric))
dist_Q_plane = dot(Q_rlu, orient_up_rlu, metric) / up_len
# angle psi enclosed by ki and Q
psi = get_psi(ki, kf, Qlen, sense_sample)
a3 = - psi - xi + a3_offs
a4 = get_a4(ki, kf, Qlen)
#print("xi = " + str(xi/np.pi*180.) + ", psi = " + str(psi/np.pi*180.))
return [a3, a4, dist_Q_plane]
def get_hkl(ki, kf, a3, Qlen, orient_rlu, orient_up_rlu, B, sense_sample=1.):
B_inv = la.inv(B)
# angle enclosed by ki and Q
psi = get_psi(ki, kf, Qlen, sense_sample)
# angle between Q and orientation reflex
xi = - a3 + a3_offs - psi
Q_lab = rotate(np.dot(B, orient_up_rlu), np.dot(B, orient_rlu*Qlen), xi)
Q_lab *= Qlen / la.norm(Q_lab)
Q_rlu = np.dot(B_inv, Q_lab)
return Q_rlu
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# get ki from kf and energy transfer
def get_ki(kf, E):
return np.sqrt(kf**2. + E_to_k2*E)
# get kf from ki and energy transfer
def get_kf(ki, E):
return np.sqrt(ki**2. - E_to_k2*E)
# get energy transfer from ki and kf
def get_E(ki, kf):
return (ki**2. - kf**2.) / E_to_k2
# -----------------------------------------------------------------------------
tascalc/tst.py
+ 1
− 1
View file @ ff95236c
......@@ -5,7 +5,7 @@
# @license see 'LICENSE' file
#
import tascalc as tas
import tas
import numpy as np
import numpy.linalg as la
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment