#/*##########################################################################
#
# The PyMca X-Ray Fluorescence Toolkit
#
# Copyright (c) 2004-2014 European Synchrotron Radiation Facility
#
# This file is part of the PyMca X-ray Fluorescence Toolkit developed at
# the ESRF by the Software group.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
#############################################################################*/
__author__ = "V.A. Sole - ESRF Data Analysis"
__contact__ = "sole@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
__doc__ = """
Methods to convert single point or complete images to reciprocal space.
It is fully vectorized and therefore very fast for converting complete
images.
"""
import numpy
cos = numpy.cos
sin = numpy.sin
[docs]class SixCircle(object):
def __init__(self):
self._energy = None
self._lambda = None
self._K = 1.0
self._ub = None
self.setLambda(1.0)
self.setUB([1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0])
[docs] def setUB(self, ublist):
"""
:param ublist: the ub matrix element values
:type ublist: list, tuple or array to convert to a 3x3 matrix
"""
self._ub = numpy.array(ublist, copy=True, dtype=numpy.float)
self._ub.shape = 3, 3
[docs] def getUB(self):
"""
:return: the ub matrix element values
:rtype: list(float)
"""
a = self._ub * 1
a.shape = -1
return a.tolist()
[docs] def setEnergy(self, energy):
"""
:param energy: the energy to set in KeV
:type energy: float
"""
self._lambda = 12.39842 / energy
self._energy = energy
self.update()
[docs] def getEnergy(self):
"""
:return: the energy in KeV
:rtype: float
"""
return self._energy
[docs] def setLambda(self, value):
"""
:param value: the wavelength to set in Angstroms
:type value: float
"""
self._lambda = value
self._energy = 12.39842 / value
self.update()
[docs] def getLambda(self):
"""
:return: the wavelength in Angstroms
:rtype: float
"""
return self._lambda
[docs] def update(self):
"""
compute K from the wavelength value
"""
self._K = (2 * numpy.pi) / self._lambda
[docs] def getPhiMatrix(self, phi):
"""
:param phi: the phi angle in degree
:type phi: float
:return: the rotation matrix of the phi axis for a given angle
:rtype: numpy.ndarray
"""
angle = numpy.radians(phi)
cphi = cos(angle)
sphi = sin(angle)
return numpy.array([[cphi, sphi, 0.0],
[-sphi, cphi, 0.0],
[0.0, 0.0, 1.0]], numpy.float)
[docs] def getChiMatrix(self, chi):
"""
:param chi: the chi angle in degree
:type chi: float
:return: the rotation matrix of the chi
:rtype: numpy.ndarray
"""
angle = numpy.radians(chi)
cchi = cos(angle)
schi = sin(angle)
return numpy.array([[cchi, 0.0, schi],
[0.0, 1.0, 0.0],
[-schi, 0.0, cchi]], numpy.float)
[docs] def getThetaMatrix(self, th):
"""
:param th: the theta angle in Degree
:type th: float
:return: the rotation matrix of the theta axis
:rtype: numpy.ndarray
"""
angle = numpy.radians(th)
cth = cos(angle)
sth = sin(angle)
return numpy.array([[cth, sth, 0],
[-sth, cth, 0],
[0, 0, 1]], numpy.float)
[docs] def getDeltaMatrix(self, delta):
"""
:param delta: the delta angle in Degree
:type delta: float
:return: the rotation matrix of the delta axis
:rtype: numpy.ndarray
"""
angle = numpy.radians(delta)
cdel = cos(angle)
sdel = sin(angle)
return numpy.array([[cdel, sdel, 0],
[-sdel, cdel, 0],
[0, 0, 1]], numpy.float)
[docs] def getGammaMatrix(self, gamma):
"""
:param gamma: the gamma angle in Degree
:type gamma: float
:return: the rotation matrix of the gamma axis
:rtype: numpy.ndarray
"""
angle = numpy.radians(gamma)
cgam = cos(angle)
sgam = sin(angle)
return numpy.array([[1.0, 0.0, 0.0],
[0.0, cgam, -sgam],
[0.0, sgam, cgam]], numpy.float)
[docs] def getMuMatrix(self, mu):
"""
:param mu: the mu angle in degree
:type mu: float
:return: the rotation matrix of the mu axis
:rtype: numpy.ndarray
"""
angle = numpy.radians(mu)
cmu = cos(angle)
smu = sin(angle)
return numpy.array([[1.0, 0.0, 0.0],
[0.0, cmu, -smu],
[0.0, smu, cmu]], numpy.float)
def _getDeltaDotGammaMatrix(self, delta, gamma, gamma_first=False):
"""
:param delta: the delta angles in Degrees
:type delta: numpy.ndarray (1D)
:param gamma: the gamma values in Degrees
:type gamma: numpy.ndarray (1D)
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: all the rotation matrix of all the delta, gamma combinations
:rtype: numpy.ndarray (3x3, len(delta) * len(gamma))
"""
delr = numpy.radians(delta)
gamr = numpy.radians(gamma)
if gamma_first:
cgam, cdel = numpy.meshgrid(numpy.cos(gamr), numpy.cos(delr))
sgam, sdel = numpy.meshgrid(numpy.sin(gamr), numpy.sin(delr))
else:
#this is to give the same result as Didier and not the transpose
cdel, cgam = numpy.meshgrid(numpy.cos(delr), numpy.cos(gamr))
sdel, sgam = numpy.meshgrid(numpy.sin(delr), numpy.sin(gamr))
deltaDotGamma = numpy.zeros((3, 3, len(delta), len(gamma)),
numpy.float)
# 1st row of dot(deltamatrix, gammaMatrix)
deltaDotGamma[0, 0, :] = cdel
deltaDotGamma[0, 1, :] = (sdel * cgam)[:]
deltaDotGamma[0, 2, :] = -sdel * sgam
# 2nd row of dot(deltaMatrix, gammaMatrix)
deltaDotGamma[1, 0, :] = -sdel
deltaDotGamma[1, 1, :] = cdel * cgam
deltaDotGamma[1, 2, :] = -cdel * sgam
# 3rd row of dot(deltaMatrix, gammaMatrix)
deltaDotGamma[2, 0, :] = 0.0
deltaDotGamma[2, 1, :] = sgam
deltaDotGamma[2, 2, :] = cgam
deltaDotGamma.shape = 3, 3, len(delta) * len(gamma)
return deltaDotGamma
[docs] def getQMu(self, phi=0., chi=0., theta=0., mu=0.,
delta=0., gamma=0., gamma_first=False):
"""
:param phi: angle in Degrees
:type phi: float
:param chi: angle in Degrees
:type chi: float
:param theta: angle in Degrees
:type theta: float
:param mu: angle in Degrees
:type mu: float
:param delta: angle in Degrees
:type delta: float or numpy.ndarray
:param gamma: angle in Degrees
:type gamma: float or numpy.ndarray
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: Q coordinates for all the given delta, gamma values
:rtype: numpy.ndarray (len(delta), len(gamma), 3)
"""
PHIi = self.getPhiMatrix(phi).T
CHIi = self.getChiMatrix(chi).T
THi = self.getThetaMatrix(theta).T
MUi = self.getMuMatrix(mu).T
tmpArray = numpy.dot(PHIi, numpy.dot(CHIi, numpy.dot(THi, MUi)))
Q = self.getQLab(mu=mu, delta=delta, gamma=gamma, gamma_first=gamma_first)
Q.shape = 3, -1
Q = numpy.transpose(numpy.dot(tmpArray, Q))
if type(delta) in [type(1.0), type(1)]:
lendelta = 1
else:
lendelta = len(delta)
if type(gamma) in [type(1.0), type(1)]:
lengamma = 1
else:
lengamma = len(gamma)
Q.shape = lengamma, lendelta, 3
return Q
[docs] def getQSurface(self, phi=0., chi=0., theta=0., mu=0.,
delta=0., gamma=0., gamma_first=False):
"""
:param phi: angle in Degrees
:type phi: float
:param chi: angle in Degrees
:type chi: float
:param theta: angle in Degrees
:type theta: float
:param mu: angle in Degrees
:type mu: float
:param delta: angle in Degrees
:type delta: float or numpy.ndarray
:param gamma: angle in Degrees
:type gamma: float or numpy.ndarray
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: Q values for all the given delta, gamma values
This is only true if the diffractometer has been properly aligned.
"""
PHIi = self.getPhiMatrix(phi).T
CHIi = self.getChiMatrix(chi).T
THi = self.getThetaMatrix(theta).T
MUi = self.getMuMatrix(mu).T
tmpArray = numpy.dot(PHIi, numpy.dot(CHIi, numpy.dot(THi, MUi)))
Q = self.getQLab(mu=mu, delta=delta, gamma=gamma, gamma_first=gamma_first)
Q.shape = 3, -1
return (numpy.dot(tmpArray, Q))
[docs] def getQLab(self, mu=0.0, delta=0.0, gamma=0.0, gamma_first=False):
"""
:param mu: angle in Degrees
:type mu: float
:param delta: angle in Degrees
:type delta: float or numpy.ndarray
:param gamma: angle in Degrees
:type gamma: float or numpy.ndarray
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: the Q coordinates in the Lab system
:rtype: numpy.ndarray ()
Q = Kf - Ki = (2 * pi / lambda) * (MU DELTA GAMMA - I) * (0, 1, 0)
This gives (transforming angles to radians):
(2*pi/lambda) * ( sin(delta) cos(gamma),
cos(mu) cos(delta) cos(gamma) - sin(mu) sin(gamma) - 1,
sin(mu) cos(delta) cos(gamma) + cos(mu) sin(gamma))
or, in terms of DG = numpy.dot(DELTA, GAMMA):
(2*pi/lambda) * ( DG[0,1],
cos(mu)* DG[1,1] - sin(mu) * DG[2,1] - 1
sin(mu)* DG[1,1] + cos(mu) * DG[2,1])
"""
alpha = numpy.radians(mu)
cmu = cos(alpha)
smu = sin(alpha)
alpha = numpy.radians(delta)
cdel = cos(alpha)
sdel = sin(alpha)
alpha = numpy.radians(gamma)
cgam = cos(alpha)
sgam = sin(alpha)
if isinstance(delta, numpy.ndarray) or \
isinstance(gamma, numpy.ndarray):
if gamma_first:
cgam, cdel = numpy.meshgrid(cgam, cdel)
sgam, sdel = numpy.meshgrid(sgam, sdel)
else:
# this is to give the same result as Didier and not the transpose
cdel, cgam = numpy.meshgrid(cdel, cgam)
sdel, sgam = numpy.meshgrid(sdel, sgam)
Q = numpy.zeros((3, sdel.shape[0], sdel.shape[1]), numpy.float)
Q[0, :, :] = sdel * cgam
Q[1, :, :] = cmu * cdel * cgam - smu * sgam - 1
Q[2, :, :] = smu * cdel * cgam + cmu * sgam
else:
Q = numpy.zeros((3, 1), numpy.float)
Q[0, 0] = sdel * cgam
Q[1, 0] = cmu * cdel * cgam - smu * sgam - 1
Q[2, 0] = smu * cdel * cgam + cmu * sgam
return Q * self._K
[docs] def getHKL(self, phi=0., chi=0., theta=0., mu=0.,
delta=0., gamma=0., gamma_first=False):
"""
:param phi: angle in Degrees
:type phi: float
:param chi: angle in Degrees
:type chi: float
:param theta: angle in Degrees
:type theta: float
:param mu: angle in Degrees
:type mu: float
:param delta: angle in Degrees
:type delta: float or numpy.ndarray
:param gamma: angle in Degrees
:type gamma: float or numpy.ndarray
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: HKL values for all the given delta, gamma values
"""
PHIi = self.getPhiMatrix(phi).T
CHIi = self.getChiMatrix(chi).T
THi = self.getThetaMatrix(theta).T
MUi = self.getMuMatrix(mu).T
UBi = numpy.linalg.inv(self._ub)
tmpArray = numpy.dot(UBi,
numpy.dot(PHIi,
numpy.dot(CHIi,
numpy.dot(THi, MUi))))
Q = self.getQLab(mu=mu, delta=delta, gamma=gamma, gamma_first=gamma_first)
Q.shape = 3, -1
return (numpy.dot(tmpArray, Q))
[docs]def getHKL(wavelength, ub, phi=0., chi=0., theta=0., mu=0.,
delta=0., gamma=0., gamma_first=False):
"""
A convenience function that takes the whole input in one go.
:param wavelength: the wavelength in Angstroms
:type wavelength: float
:param ub: the ub matrix element values
:type ub: list(float)
:param phi: angle in Degrees
:type phi: float
:param chi: angle in Degrees
:type chi: float
:param theta: angle in Degrees
:type theta: float
:param mu: angle in Degrees
:type mu: float
:param delta: angle in Degrees
:type delta: float or numpy.ndarray
:param gamma: angle in Degrees
:type gamma: float or numpy.ndarray
:param gamma_first: if delta and gamma are arrays, which one variates first.
:type gamma_first: boolean
:return: HKL values for all the given delta, gamma values
"""
a = SixCircle()
a.setLambda(wavelength)
a.setUB(ub)
return a.getHKL(delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma, gamma_first=gamma_first)
[docs]def main():
wavelength = 0.363504
UB = [1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0]
UB[0] = -4.080
UB[1] = 0.000
UB[2] = 0.000
UB[3] = 0.000
UB[4] = 4.080
UB[5] = 0.000
UB[6] = 0.000
UB[7] = 0.000
UB[8] = -4.080
d = SixCircle()
d.setLambda(wavelength)
d.setUB(UB)
print("H = 0 K = 0 L = 1")
delta, theta, chi, phi, mu, gamma = 13.5558, 6.77779, -90, 0.0, 0.0, 0.0
print(d.getHKL(delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma))
print("H = 0 K = 1 L = 0")
delta, theta, chi, phi, mu, gamma = 13.5558, 96.77779, -90, 0.0, 0.0, 0.0
print(d.getHKL(delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma))
print("H = 1 K = 1 L = 1")
delta, theta, chi, phi, mu, gamma = 23.5910, 47.0595, -135., 0.0, 0.0, 0.0
print(d.getHKL(delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma))
print("H = 2 K = -1 L = 0")
delta, theta, chi, phi, mu, gamma = 30.6035, -11.2635, 180.0, 0.0, 0.0, 0.0
print(d.getHKL(delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma))
print("H = 2 K = -1 L = 0")
print(getHKL(wavelength, UB, delta=delta, theta=theta, chi=chi, phi=phi,
mu=mu, gamma=gamma))
if __name__ == "__main__":
main()