added preliminary version of matrix analysis of NSE (f3625e04) · Commits · Zolnierczuk, Piotr / pysen

misc/fix_nexus_value.py

+40 −15

Original line number	Diff line number	Diff line
		@@ -7,14 +7,34 @@ import argparse
		#
		#
		import h5py
		from h5py import string_dtype

		def copy_file(filename, newfile):
		"copy file, do not overwrite"
		with open(filename, "rb") as src, open(newfile, "xb") as dst:
		shutil.copyfileobj(src, dst)

		def replace_value_base(_df, _key, newvalue):
		"trivial replacement"
		return newvalue

		def fix_value(filename, key, newvalue):

		def replace_value_full(df, key, newvalue):
		"""'full replacement'
		recreate string key with new length
		"""
		try:
		# let's hope it is a string
		len(df[key][0]) # old length (test if throws an exception)
		new_len = len(newvalue) # new length
		del df[key]
		df.create_dataset(key, shape=1, dtype=string_dtype(length=new_len))
		except TypeError:
		# scalar value
		pass
		return newvalue

		def fix_value(filename, newvalues):
		"fix value"
		try:
		#1. make new file, do not overwrite if exists
		@@ -23,12 +43,16 @@ def fix_value(filename, key, newvalue):
		copy_file(filename, newfile)
		#2. open this new file in a write mode
		hf = h5py.File(newfile, mode='r+')
		for key, newvalue in newvalues.items():
		#3. replace the value
		oldvalue = hf[key][0]
		hf[key][0] = newvalue
		hf[key][0] = replace_value_base(hf, key, newvalue)
		#hf[key][0] = replace_value_full(hf, key, newvalue)
		print(f"{newfile}: {key}={newvalue} (was {oldvalue})")
		#4. done
		hf.close()
		print(f"{newfile}: replaced {key}={oldvalue} with {newvalue}")
		except KeyError:
		print(f"{newfile}: key {key} does not exist")
		except FileExistsError:
		print(f"{newfile}: already exists and will not be overwritten.")

		@@ -36,19 +60,20 @@ def fix_value(filename, key, newvalue):
		def main():
		"the main"
		parser = argparse.ArgumentParser(description='echo plot')
		parser.set_defaults(key='/entry/title', new_value=None)
		parser.add_argument('file', metavar='filename', help='file to process')
		parser.add_argument('--key', '-k', dest='key',
		help='set the key to fix (default=%(default)s)')
		parser.add_argument('--newvalue', '-n', dest='new_value',
		help='set the new value')
		parser.set_defaults(newvalues=None)
		parser.add_argument('file', metavar='filename', nargs='+', help='file to process')
		parser.add_argument('--key', '-k', dest='newvalues', nargs='+')

		args = parser.parse_args()
		if not args.key or args.new_value is None:
		print('need a key and new value!!!')
		return -1
		#print(vars(args))

		fix_value(args.file, key=args.key, newvalue=args.new_value)
		newvalues = {}
		for line in args.newvalues:
		key,val = line.split(':')
		newvalues[key] = val
		print(newvalues)
		for filename in args.file:
		fix_value(filename, newvalues=newvalues)
		return 0

		if __name__ == "__main__":

pyproject.toml

+1 −1

Original line number	Diff line number	Diff line
		@@ -20,7 +20,7 @@ classifiers = [
		"Topic :: Scientific/Engineering",
		"Topic :: Scientific/Engineering :: Visualization",
		]
		dependencies = ['numpy', 'scipy', 'matplotlib', 'pandas', 'h5py' , 'qtpy', 'qtconsole']
		dependencies = ['numpy', 'scipy', 'matplotlib', 'pandas', 'h5py' , 'qtpy', 'qtconsole', 'sympy']

		[project.scripts]
		nseplot = "pysen.ui.nseplot:main"

pysen/echo/matrix_analysis.py

0 → 100644

+156 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python
		# coding: utf-8
		r"""
		Matrix Analysis of Neutron Spin Echo

		References:
		[1] Hayter, J.B. Matrix analysis of neutron spin-echo. Z Physik B 31, 117–125 (1978)
		[2] Ohl


		Beam Coordinate System (BCS) Sample Coordinate System (SCS)
		Hayter in [1] Monkenbusch in [2]

		^ x-axis z-axis ^ ^
		\| (up) (up) \| / y-axis (q)
		\| \| /
		\| \| /
		\| \|/
		+--------- > z-axis (beam) +------- > x-axis (beam)
		/
		/ y-axis
		v ('anti-q')

		"""
		__all__ = ['Rx', 'Ry', 'Rz', 'B2S', 'M_BS', 'T_op', 'Fx', 'Fy', 'Cx', 'Lz' ]


		from sympy import sin, cos, pi, Matrix, trigsimp

		#pylint: disable=invalid-name

		# Rotation Matrices
		# Vector rotation maxtrices about x, y and z-axis ("active" rotation)
		# To rotate system of coordinates ("passive" rotation) simply replace $\alpha$ by $-\alpha$.

		def Rx(alpha):
		"rotation about x-axis"
		return Matrix([[ 1, 0 , 0],
		[ 0, cos(alpha), -sin(alpha)],
		[ 0, sin(alpha), cos(alpha)]])

		def Ry(beta):
		"rotation about y-axis"
		return Matrix([[ cos(beta), 0 , sin(beta)],
		[ 0, 1 , 0],
		[ -sin(beta), 0 , cos(beta)]])

		def Rz(gamma):
		"rotation about z-axis"
		return Matrix([[ cos(gamma), -sin(gamma), 0],
		[ sin(gamma), cos(gamma), 0],
		[ 0, 0 , 1]])



		M_BS= Matrix([[0, 0, 1],
		[0, -1, 0],
		[1, 0, 0]])
		def B2S(Op):
		"convert between BCS and SCS"
		return M_BSOpM_BS

		def T_op(theta, psi):
		r"""
		$\vec{B}$ at an angle $\theta$ to x-axis (vertical direction)

		1. Rotate system of coordinates about y axis to make z parallel to $\vec{B}$.
		2. Rotate polarization vector about the new z-axis ($\vec{B}$), angle is $\psi$.
		3. Rotate system of coordinates about y axis back.

		$\mathbf{T}=R_y(\varphi) R_z(\psi) R^_y(-\varphi)$

		where $\varphi=\pi/2-\theta$

		We get spin transfer operator, Equation (12) (HCS)."""
		op = Ry(+(pi/2-theta))Rz(psi)Ry(-(pi/2-theta))
		return trigsimp(op)



		def Fx(psi, positive=True):
		r"""The $\pi$ Spin Turn Operator
		Equation (16) in [1]

		$\mathbf{\hat F}_{\pm x} =
		\begin{pmatrix}
		1 & 0 & 0 \\
		0 & \cos\psi & \mp\sin\psi \\
		0 & \pm\sin\psi & \cos\psi
		\end{pmatrix}$"""
		if positive:
		return T_op(0, psi)
		return T_op(pi, psi)


		def Fy(psi, positive=True):
		r"""The $\pi$ Spin Turn Operator
		spin turn about y axis
		Equation (17) in [1] (corrected)

		$\mathbf{\hat F}_{\pm y} =
		\begin{pmatrix}
		\cos\psi & 0 & \pm\sin\psi \\
		0 & 1 & 0 \\
		\mp\sin\psi & 0 & \cos\psi
		\end{pmatrix}$"""
		if positive:
		return Ry(psi)
		return Ry(-psi)

		def Cx(psi, positive=True):
		r"""The $\pi/2$ Spin Turn Operator

		Equation (28) in [1]

		$\mathbf{\hat C}_{\pm x} = \frac{1}{2}
		\begin{pmatrix}
		1+\cos\psi & -\sqrt{2}\sin\psi & \pm(1-\cos\psi) \\
		\sqrt{2}\sin\psi & 2\cos\psi & \mp\sqrt{2}\sin\psi \\
		\pm(1-\cos\psi) & \pm\sqrt{2}\sin\psi & 1+\cos\psi
		\end{pmatrix}$"""
		if positive:
		return T_op(+pi/4, psi)
		#return T_op(-pi/4, -psi)
		return T_op(3*pi/4, psi)

		def CxCx(psi, positive=True):
		r"""

		Equation (33) in [1]

		$\mathbf{\hat C}_{\pm x} mathbf{\hat C}_{\mp x} = \frac{1}{2}
		\begin{pmatrix}
		2\cos\psi - \sin^2\psi & -\sqrt{2}\sin\psi (1+\cos\psi) & \pm\sin^2\psi \\
		\sqrt{2}\sin\psi (1+\cos\psi) & 2\cos^2\psi & \sqrt{2}\sin\psi (1-\cos\psi) \\
		\pm\sin^2\psi & \sqrt{2}\sin\psi (1-\cos\psi) & 2\cos\psi + \sin^2\psi
		\end{pmatrix}$"""
		if positive:
		return trigsimp(Cx(psi, positive=True)*Cx(psi, positive=False))
		return trigsimp(Cx(psi, positive=False)*Cx(psi, positive=True))

		def Lz(omega, positive=True):
		r"""Precession About Beam Axis

		Equation (34) in [1]

		$\mathbf{\hat L}_{\pm z} =
		\begin{pmatrix}
		\cos\Omega & \mp\sin\Omega & 0 \\
		\pm\sin\Omega & \cos\Omega & 0 \\
		0 & 0 & 1
		\end{pmatrix}$"""

		if positive:
		return Rz(+omega)
		return Rz(-omega)

pysen/revision.py

+2 −2

Original line number	Diff line number	Diff line
		@@ -2,8 +2,8 @@
		PySEN revision module
		"""
		import sys
		__version__ = "2.1.0.dev2"
		__date__ = "Nov 25, 2025"
		__version__ = "2.1.0.dev3"
		__date__ = "Dec 4, 2025"

		def version(full=False):
		"get pysen version number"

test/test_echo_matrix_analysis.py

0 → 100644

+186 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python
		# coding: utf-8
		"""
		Python unit test cases for pysen.echo.matrix_analysis module
		"""

		import unittest

		from sympy import symbols, sqrt, sin, cos, pi, Matrix, trigsimp
		from sympy import init_printing, pprint

		from pysen.echo.matrix_analysis import Fx, Fy, Cx, Lz, T_op, B2S


		#pylint: disable=invalid-name


		class EchoMatrixAnalysisTestCase(unittest.TestCase):
		"unit test case class"

		def setUp(self):
		"setup method"
		init_printing()
		self.verbose = False


		def tearDown(self):
		"tear down method"



		def test_t_operators(self):
		"test T_op operator"
		theta, psi = symbols('theta psi')
		#
		texp = Matrix([[ cos(theta)2 + sin(theta)2cos(psi), -sin(theta)sin(psi),
		sin(theta)cos(theta)(1-cos(psi)) ],
		[ sin(theta)*sin(psi) , cos(psi) ,
		-cos(theta)*sin(psi) ],
		[ sin(theta)cos(theta)(1-cos(psi)) , cos(theta)*sin(psi),
		sin(theta)2 + cos(theta)2*cos(psi)]])
		#
		tact = T_op(theta,psi)
		if self.verbose:
		print("\ntesting T operator (BCS)")
		print("T_op(exp)")
		pprint(texp)
		print("T_op(act)")
		pprint(tact)
		self.assertEqual(texp, tact, 'T_op(act)==T_op(exp)')

		if self.verbose:
		print("\ntesting T operator (SCS)")
		print("T_op(exp)")
		pprint(B2S(texp))
		print("T_op(act)")
		pprint(B2S(tact))
		self.assertEqual(B2S(texp), B2S(tact), 'T_op(act)==T_op(exp) [B2S]')

		def test_f_operators(self):
		"test F operator"
		psi = symbols('psi')
		#
		tpos = Matrix([[ 1, 0, 0],
		[ 0 , cos(psi), -sin(psi)],
		[ 0 ,+sin(psi), cos(psi)]])

		tneg = Matrix([[ 1, 0, 0],
		[ 0 , cos(psi), +sin(psi)],
		[ 0 ,-sin(psi), cos(psi)]])
		tact = Fx(psi, positive=True)
		#
		if self.verbose:
		print("\ntesting F_+x operator")
		print("T_op(exp)")
		pprint(tpos)
		print("T_op(act)")
		pprint(tact)
		self.assertEqual(tpos, tact,'F_+x(act)==F_+x(exp)')

		tact = Fx(psi, positive=False)
		if self.verbose:
		print("\ntesting F_-x operator")
		print("T_op(exp)")
		pprint(tneg)
		print("T_op(act)")
		pprint(tact)
		self.assertEqual(tneg, tact,'F_-x(act)==F_-x(exp)')

		def test_c_operators(self):
		"test C operator"
		psi = symbols('psi')
		#
		tpos = Matrix([[ (1+cos(psi)) , -sqrt(2)*sin(psi), +(1-cos(psi)) ],
		[ sqrt(2)sin(psi), 2cos(psi), -sqrt(2)*sin(psi)],
		[ (1-cos(psi)) , +sqrt(2)*sin(psi), (1+cos(psi)) ]])/2

		tneg = Matrix([[ (1+cos(psi)) , -sqrt(2)*sin(psi), -(1-cos(psi)) ],
		[ sqrt(2)sin(psi), 2cos(psi), +sqrt(2)*sin(psi)],
		[-(1-cos(psi)) , -sqrt(2)*sin(psi), (1+cos(psi)) ]])/2

		tact = Cx(psi, positive=True)
		#
		if self.verbose:
		print("\ntesting C_+x operator")
		print("T_op(exp)")
		pprint(tpos*2)
		print("T_op(act)")
		pprint(tact*2)
		#self.assertEqual(tpos, tact,'C_+x(act)==C_+x(exp)')

		tact = Cx(psi, positive=False)
		#
		if self.verbose:
		print("\ntesting C_-x operator")
		print("T_op(exp)")
		pprint(tneg*2)
		print("T_op(act)")
		pprint(tact*2)
		#self.assertEqual(tneg, tact,'C_-x(act)==C_-x(exp)')
		pprint('C-DIRECT')
		Cp = T_op( +pi/4, psi)
		Cn = T_op(3*pi/4, psi)

		self.assertEqual(tpos, Cp,'C+')
		self.assertEqual(tneg, Cn,'C-')

		TpnH = Matrix([[ 2cos(psi)-sin(psi)2 , -sqrt(2)sin(psi)(1+cos(psi)), +sin(psi)*2],
		[ sqrt(2)sin(psi)(1+cos(psi)), 2cos(psi)2 , -sqrt(2)sin(psi)*(1-cos(psi))],
		[ +sin(psi)*2 , +sqrt(2)sin(psi)(1-cos(psi)), 2cos(psi) + sin(psi)**2]])/2
		Tpn = Matrix([[ 2cos(psi)-sin(psi)2 , -sqrt(2)sin(psi)(1+cos(psi)), -sin(psi)*2],
		[ sqrt(2)sin(psi)(1+cos(psi)), 2cos(psi)2 , -sqrt(2)sin(psi)*(1-cos(psi))],
		[ +sin(psi)*2 , -sqrt(2)sin(psi)(1-cos(psi)), 2cos(psi) + sin(psi)**2]])/2
		TnpH = Matrix([[ 2cos(psi)-sin(psi)2 , -sqrt(2)sin(psi)(1+cos(psi)), -sin(psi)*2],
		[ sqrt(2)sin(psi)(1+cos(psi)), 2cos(psi)2 , -sqrt(2)sin(psi)*(1-cos(psi))],
		[ -sin(psi)*2 , +sqrt(2)sin(psi)(1-cos(psi)), 2cos(psi) + sin(psi)**2]])/2
		Tnp = Matrix([[ 2cos(psi)-sin(psi)2 , -sqrt(2)sin(psi)(1+cos(psi)), +sin(psi)*2],
		[ sqrt(2)sin(psi)(1+cos(psi)), 2cos(psi)2 , +sqrt(2)sin(psi)*(1-cos(psi))],
		[ -sin(psi)*2 , +sqrt(2)sin(psi)(1-cos(psi)), 2cos(psi) + sin(psi)**2]])/2

		Cpn = Cp*Cn
		Cnp = Cn*Cp
		print('Cpn')
		#pprint(trigsimp(Tpn*2))
		#pprint(trigsimp(Cpn*2))
		pprint(trigsimp(Cpn-Tpn))

		print('Cnp')
		#pprint(trigsimp(Tnp*2))
		#pprint(trigsimp(Cnp*2))
		pprint(trigsimp(Cnp-Tnp))

		#pprint('Terr1')
		#pprint(trigsimp(Tpn-TpnH))
		#pprint('Terr2')
		#pprint(trigsimp(Tnp-TnpH))

		pprint('Tsum')
		pprint(trigsimp(Tpn+Tnp))
		pprint('Tdif')
		pprint(trigsimp(Tpn-Tnp))



		def test_se_operators(self):
		"test function"
		Omega = symbols('Omega')
		mx = Matrix([[1, 0, 0], [0, -1, 0], [0, 0, -1]])
		my = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])


		Tx_SE = trigsimp(Cx(pi, False)Lz(Omega)Fx(pi)Lz(Omega)Cx(pi, True))
		Ty_SE = trigsimp(Cx(pi, False)Lz(Omega)Fy(pi)Lz(Omega)Cx(pi, True))

		if self.verbose:
		print("\ntesting T_SE operators")
		print("Tx_SE(exp)")
		pprint(Tx_SE)
		print("Ty_SE(exp)")
		pprint(Ty_SE)
		self.assertEqual(Tx_SE, mx, "Tx_SE")
		self.assertEqual(Ty_SE, my, "Ty_SE")


		if __name__ == "__main__":
		unittest.main(exit=False, verbosity=2)