Commit f25b12b990846d18ab4b1a38107f75c2769e837c - xrayutilities

Imported Upstream version 1.0.4 Eugen Wintersberger 10 years ago

29 changed file(s) with 1204 addition(s) and 378 deletion(s). Raw diff Collapse all Expand all

+11

-0

CHANGES.txt less more

	0	v1.0.4, 2013-12-22 "the final Linz edition"
	1
	2	* bugfix in GID experimental class: qconv argument so far not used correctly!
	3	* enhanced fitting routines
	4	* spec files can now be parsed while gzipped
	5	* CBF file parser for detector images
	6	* compile fixes for Microsoft compilers
	7	* compile fixes on Mac
	8	* new predefined materials (NaCl, FeO, Fe3O4, Co3O4,...)
	9	* some minor bug fixes
	10
0	11	v1.0.3, 2013-10-25
1	12
2	13	* build_doc target for setup.py to help build docs (thanks to F. Picca)

-1

PKG-INFO less more

0	0	Metadata-Version: 1.1
1	1	Name: xrayutilities
2		Version: 1.0.3
	2	Version: 1.0.4
3	3	Summary: package for x-ray diffraction data evaluation
4	4	Home-page: http://xrayutilities.sourceforge.net
5	5	Author: Dominik Kriegner

-1

doc/source/conf.py less more

49	49	# The short X.Y version.
50	50	version = '1.0'
51	51	# The full version, including alpha/beta/rc tags.
52		release = '1.0.3'
	52	release = '1.0.4'
53	53
54	54	# The language for content autogenerated by Sphinx. Refer to documentation
55	55	# for a list of supported languages.

+24

-0

doc/source/xrayutilities.io.rst less more

8	8	:undoc-members:
9	9	:show-inheritance:
10	10
	11	:mod:`cbf` Module
	12	-----------------
	13
	14	.. automodule:: xrayutilities.io.cbf
	15	:members:
	16	:undoc-members:
	17	:show-inheritance:
	18
	19	:mod:`desy_tty08` Module
	20	------------------------
	21
	22	.. automodule:: xrayutilities.io.desy_tty08
	23	:members:
	24	:undoc-members:
	25	:show-inheritance:
	26
11	27	:mod:`edf` Module
12	28	-----------------
13	29
14	30	.. automodule:: xrayutilities.io.edf
	31	:members:
	32	:undoc-members:
	33	:show-inheritance:
	34
	35	:mod:`helper` Module
	36	--------------------
	37
	38	.. automodule:: xrayutilities.io.helper
15	39	:members:
16	40	:undoc-members:
17	41	:show-inheritance:

-2

setup.py less more

84	84	os.path.join('xrayutilities','src','gridder2d.c'),
85	85	os.path.join('xrayutilities','src','block_average.c'),
86	86	os.path.join('xrayutilities','src','qconversion.c'),
87		os.path.join('xrayutilities','src','gridder3d.c')],
	87	os.path.join('xrayutilities','src','gridder3d.c'),
	88	os.path.join('xrayutilities','src','file_io.c')],
88	89	define_macros = user_macros)
89	90
90	91	try:

113	114	pass
114	115
115	116	setup(name="xrayutilities",
116		version="1.0.3",
	117	version="1.0.4",
117	118	author="Eugen Wintersberger, Dominik Kriegner",
118	119	description="package for x-ray diffraction data evaluation",
119	120	classifiers=["Topic :: Scientific/Engineering :: Physics",

-1

xrayutilities/__init__.py less more

38	38	from .experiment import HXRD
39	39	from .experiment import NonCOP
40	40	from .experiment import GID
41		from .experiment import GID_ID10B
42	41	from .experiment import GISAXS
43	42	from .experiment import Powder
44	43	from .experiment import QConversion

+65

-38

xrayutilities/analysis/sample_align.py less more

352	352	## detector parameter calculation from scan with
353	353	## area detector (determine maximum by center of mass)
354	354	######################################################
355		def area_detector_calib(angle1,angle2,ccdimages,detaxis,r_i,plot=True,cut_off = 0.7,start = (0,0,0,0), fix = (False,False,False,False), fig=None,wl=None):
	355	def area_detector_calib(angle1,angle2,ccdimages,detaxis,r_i,plot=True,cut_off = 0.7,start = (0,0,0,0), fix = (False,False,False,False), fig=None,wl=None,plotlog=False,debug=False):
356	356	"""
357	357	function to calibrate the detector parameters of an area detector
358	358	it determines the detector tilt possible rotations and offsets in the

382	382	default: None (creates own figure)
383	383	wl ...... wavelength of the experiment in Angstrom (default: config.WAVELENGTH)
384	384	value does not matter here and does only affect the scaling of the error
385		"""
386
387		debug=False
388		plotlog=False
	385	plotlog . flag to specify if the created error plot should be on log-scale
	386	debug ... flag to specify that you want to see verbose output and saving of images to show if the CEN determination works
	387
	388	"""
	389
389	390	if plot:
390	391	try: plt.__name__
391	392	except NameError:

422	423	img = ccdimages[i]
423	424	if numpy.sum(img) > cut_off*avg:
424	425	[cen1,cen2] = center_of_mass(img)
	426	if debug:
	427	plt.figure("_ccd")
	428	plt.imshow(utilities.maplog(img),origin='low')
	429	plt.plot(cen2,cen1,"wo",mfc='none')
	430	plt.axis([cen2-25,cen2+25,cen1-25,cen1+25])
	431	plt.savefig("xu_calib_ccd_img%d.png"%i)
	432	plt.close("_ccd")
	433
425	434	n1 = numpy.append(n1,cen1)
426	435	n2 = numpy.append(n2,cen2)
427	436	ang1 = numpy.append(ang1,angle1[i])

438	447
439	448	epslist = []
440	449	paramlist = []
441		epsmin = 1.
	450	epsmin = numpy.inf
442	451	fitmin = None
443	452
444	453	print("tiltaz tilt detrot offset: error (relative) (fittime)")

620	629
621	630	return detdir1,detdir2
622	631
623		def _area_detector_calib_fit(ang1,ang2,n1,n2, detaxis, r_i, detdir1, detdir2, start = (0,0,0,0), fix = (False,False,False,False),full_output=False,wl=1.):
	632	def _area_detector_calib_fit(ang1,ang2,n1,n2, detaxis, r_i, detdir1, detdir2, start = (0,0,0,0), fix = (False,False,False,False),full_output=False,wl=1.,debug=False):
624	633	"""
625	634	INTERNAL FUNCTION
626	635	function to calibrate the detector parameters of an area detector

648	657	full_output flag to tell if to return fit object with final parameters and detector directions
649	658	wl ...... wavelength of the experiment in Angstrom (default: 1)
650	659	value does not matter here and does only affect the scaling of the error
	660	debug ... flag to tell if you want to see debug output of the script (switch this to true only if you can handle it :))
651	661
652	662	returns
653	663	-------

656	666	if full_output:
657	667	eps,param,fit
658	668	"""
659
660		debug=False
661	669
662	670	def areapixel(params,detectorDir1,detectorDir2,r_i,detectorAxis,args,*kwargs):
663	671	"""

965	973	## detector parameter calculation from scan with
966	974	## area detector (determine maximum by center of mass)
967	975	######################################################
968		def area_detector_calib_hkl(sampleang,angle1,angle2,ccdimages,hkls,experiment,material,detaxis,r_i,plot=True,cut_off = 0.1,start = (0,0,0,0,0,0,'config'), fix = (False,False,False,False,False,False,False), fig=None):
	976	def area_detector_calib_hkl(sampleang,angle1,angle2,ccdimages,hkls,experiment,material,detaxis,r_i,plot=True,cut_off = 0.1,start = (0,0,0,0,0,0,'config'), fix = (False,False,False,False,False,False,False), fig=None, plotlog=False, debug=False):
969	977	"""
970	978	function to calibrate the detector parameters of an area detector
971	979	it determines the detector tilt possible rotations and offsets in the

983	991	angle2 ..... inner detector arm angle
984	992	ccdimages .. images of the ccd taken at the angles given above
985	993	hkls ....... array/list of hkl values for every image
	994	experiment . Experiment class object needed to get the UB matrix for the hkl peak treatment
986	995	material ... material used as reference crystal
987	996	detaxis .... detector arm rotation axis
988	997	default: ['z+','y-']

1001	1010	fix ..... fix parameters of start (default: (False,False,False,False,False,False,False))
1002	1011	fig ..... matplotlib figure used for plotting the error
1003	1012	default: None (creates own figure)
1004		"""
1005
1006		debug=False
1007		plotlog=False
	1013	plotlog . flag to specify if the created error plot should be on log-scale
	1014	debug ... flag to tell if you want to see debug output of the script (switch this to true only if you can handle it :))
	1015	"""
	1016
1008	1017	if plot:
1009	1018	try: plt.__name__
1010	1019	except NameError:
1011		print("XU.analyis.area_detector_calib: Warning: plot functionality not available")
	1020	print("XU.analyis.area_detector_calib_hkl: Warning: plot functionality not available")
1012	1021	plot = False
1013	1022
1014	1023	if start[-1]=='config':

1023	1032
1024	1033	# determine center of mass position from detector images
1025	1034	# also use only images with an intensity larger than 70% of the average intensity
	1035	# the image selection is only performed for images in the primary beam
1026	1036	n1 = numpy.zeros(0,dtype=numpy.double)
1027	1037	n2 = n1
1028	1038	ang1 = n1

1031	1041	usedhkls = []
1032	1042
1033	1043	avg = 0
	1044	imgpbcnt = 0
1034	1045	for i in range(Npoints):
1035		avg += numpy.sum(ccdimages[i])
1036		avg /= float(Npoints)
	1046	if (numpy.all(hkls[i] == (0,0,0))):
	1047	avg += numpy.sum(ccdimages[i])
	1048	imgpbcnt+=1
	1049
	1050	if imgpbcnt > 0:
	1051	avg /= float(imgpbcnt)
	1052	else:
	1053	avg = 0
1037	1054	(N1,N2) = ccdimages[0].shape
1038	1055
1039	1056	if debug:

1041	1058
1042	1059	for i in range(Npoints):
1043	1060	img = ccdimages[i]
1044		if numpy.sum(img) > cut_off*avg:
	1061	if ((numpy.sum(img) > cut_off*avg) or (numpy.all(hkls[i] != (0,0,0)))):
1045	1062	[cen1,cen2] = center_of_mass(img)
1046		# if True:
1047		# plt.figure("_ccd")
1048		# plt.imshow(utilities.maplog(img),origin='low')
1049		# plt.plot(cen2,cen1,"wo",mfc='none')
1050		# plt.axis([cen2-25,cen2+25,cen1-25,cen1+25])
1051		# plt.savefig("_ccd/img%d.png"%i)
1052		# plt.close("_ccd")
	1063	if debug:
	1064	plt.figure("_ccd")
	1065	plt.imshow(utilities.maplog(img),origin='low')
	1066	plt.plot(cen2,cen1,"wo",mfc='none')
	1067	plt.axis([cen2-25,cen2+25,cen1-25,cen1+25])
	1068	plt.savefig("xu_calib_hkl_ccd_img%d.png"%i)
	1069	plt.close("_ccd")
1053	1070	n1 = numpy.append(n1,cen1)
1054	1071	n2 = numpy.append(n2,cen2)
1055	1072	ang1 = numpy.append(ang1,angle1[i])

1080	1097
1081	1098	epslist = []
1082	1099	paramlist = []
1083		epsmin = 1.
	1100	epsmin = numpy.inf
1084	1101	fitmin = None
1085	1102
1086	1103	print("tiltaz tilt detrot offset sampletilt+azimuth wavelength: error (relative) (fittime)")

1178	1195	return (cch1,cch2,pwidth1,pwidth2,tiltazimuth,tilt,detrot,outerangle_offset),eps
1179	1196
1180	1197
1181		def _area_detector_calib_fit2(sang,ang1,ang2,n1,n2, hkls, experiment, material, detaxis, r_i, detdir1, detdir2, start = (0,0,0,0,0,0,1.0), fix = (False,False,False,False,False,False,False),full_output=False):
	1198	def _area_detector_calib_fit2(sang,ang1,ang2,n1,n2, hkls, experiment, material, detaxis, r_i, detdir1, detdir2, start = (0,0,0,0,0,0,1.0), fix = (False,False,False,False,False,False,False),full_output=False,debug=False):
1182	1199	"""
1183	1200	INTERNAL FUNCTION
1184	1201	function to calibrate the detector parameters of an area detector

1192	1209	angle2 ..... inner detector arm angle
1193	1210	n1,n2 ...... pixel number at which the beam was observed
1194	1211	hkls ....... Miller indices of the reflection were the images were taken (use (0,0,0)) for primary beam
	1212	experiment . Experiment class object needed to get the UB matrix needed for the hkl peak treatment
1195	1213	material ... material used as reference crystal
1196	1214	detaxis .... detector arm rotation axis
1197	1215	default: ['z+','y-']

1207	1225	By default: (0,0,0,0,0,0,1.0)
1208	1226	fix ..... fix parameters of start
1209	1227	full_output flag to tell if to return fit object with final parameters and detector directions
1210
	1228	debug ... flag to tell if you want to see debug output of the script (switch this to true only if you can handle it :))
	1229
1211	1230	returns
1212	1231	-------
1213	1232	eps final epsilon of the fit

1215	1234	if full_output:
1216	1235	eps,param,fit
1217	1236	"""
1218
1219		debug=False
1220	1237
1221	1238	def areapixel2(params,detectorDir1,detectorDir2,r_i,detectorAxis,args,*kwargs):
1222	1239	"""

1739	1756	# correct substrate Bragg peak position in
1740	1757	# reciprocal space maps
1741	1758	#################################################
1742		def fit_bragg_peak(om,tt,psd,omalign,ttalign,exphxrd,frange=(0.03,0.03),plot=True):
	1759	def fit_bragg_peak(om,tt,psd,omalign,ttalign,exphxrd,frange=(0.03,0.03),peaktype='Gauss',plot=True):
1743	1760	"""
1744	1761	helper function to determine the Bragg peak position in a reciprocal
1745	1762	space map used to obtain the position needed for correction of the data.
1746	1763	the determination is done by fitting a two dimensional Gaussian
1747		(xrayutilities.math.Gauss2d)
	1764	(xrayutilities.math.Gauss2d) or Lorentzian
	1765	(xrayutilities.math.Lorentz2d)
1748	1766
1749	1767	PLEASE ALWAYS CHECK THE RESULT CAREFULLY!
1750	1768

1761	1779	reciprocal space.
1762	1780	frange: data range used for the fit in both directions
1763	1781	(see above for details default:(0.03,0.03) unit: \AA^{-1})
1764		plot: if True (default) function will plot the result of the fit in comparison
1765		with the measurement.
	1782	peaktype: can be 'Gauss' or 'Lorentz' to fit either of the two peak
	1783	shapes
	1784	plot: if True (default) function will plot the result of the fit in
	1785	comparison with the measurement.
1766	1786
1767	1787	Returns
1768	1788	-------
1769	1789	omfit,ttfit,params,covariance: fitted angular values, and the fit
1770		parameters (of the Gaussian) as well as their errors
1771		"""
	1790	parameters (of the Gaussian/Lorentzian) as well as their errors
	1791	"""
	1792	if peaktype=='Gauss':
	1793	func = math.Gauss2d
	1794	elif peaktype=='Lorentz':
	1795	func = math.Lorentz2d
	1796	else:
	1797	raise InputError("peaktype must be either 'Gauss' or 'Lorentz'")
	1798
1772	1799	if om.size != psd.size:
1773	1800	[qx,qy,qz] = exphxrd.Ang2Q.linear(om,tt)
1774	1801	else:
1775	1802	[qx,qy,qz] = exphxrd.Ang2Q(om,tt)
1776	1803	[qxsub,qysub,qzsub] = exphxrd.Ang2Q(omalign,ttalign)
1777	1804	params = [qysub,qzsub,0.001,0.001,psd.max(),0,0.]
1778		params,covariance = math.fit_peak2d(qy.flatten(),qz.flatten(),psd.flatten(),params,[qysub-frange[0],qysub+frange[0],qzsub-frange[1],qzsub+frange[1]],math.Gauss2d,maxfev=10000)
	1805	params,covariance = math.fit_peak2d(qy.flatten(),qz.flatten(),psd.flatten(),params,[qysub-frange[0],qysub+frange[0],qzsub-frange[1],qzsub+frange[1]],func,maxfev=10000)
1779	1806	# correct params
1780	1807	params[6] = params[6]%(numpy.pi)
1781	1808	if params[5]<0 : params[5] = 0

1794	1821	INT = utilities.maplog(gridder.gdata.transpose(),4,0)
1795	1822	QXm = gridder.xmatrix
1796	1823	QZm = gridder.ymatrix
1797		cl = plt.contour(gridder.xaxis,gridder.yaxis,utilities.maplog(math.Gauss2d(QXm,QZm,*params),4,0).transpose(),8,colors='k',linestyles='solid')
	1824	cl = plt.contour(gridder.xaxis,gridder.yaxis,utilities.maplog(func(QXm,QZm,*params),4,0).transpose(),8,colors='k',linestyles='solid')
1798	1825	cf = plt.contourf(gridder.xaxis, gridder.yaxis, INT,35)
1799	1826	cf.collections[0].set_label('data')
1800	1827	cl.collections[0].set_label('fit')

-100

xrayutilities/experiment.py less more

1526	1526	Experiment.__init__(self,idir,ndir,**keyargs)
1527	1527
1528	1528	# initialize Ang2Q conversion
1529		self._A2QConversion = QConversion(['x+','y+','z-'],'x+',[0,1,0],wl=self._wl) # 3S+1D goniometer (as in the MRD, omega,chi,phi,theta)
1530		self.Ang2Q = self._A2QConversion
	1529	if "qconv" not in keyargs:
	1530	self._A2QConversion = QConversion(['x+','y+','z-'],'x+',[0,1,0],wl=self._wl) # 3S+1D goniometer (standard four-circle goniometer, omega,chi,phi,theta)
	1531	self.Ang2Q = self._A2QConversion
1531	1532
1532	1533	def Ang2Q(self,om,chi,phi,tt,**kwargs):
1533	1534	"""

1686	1687	Experiment.__init__(self,idir,ndir,**keyargs)
1687	1688
1688	1689	# initialize Ang2Q conversion
1689		self._A2QConversion = QConversion(['z-','x+'],['x+','z-'],[0,1,0],wl=self._wl) # 2S+2D goniometer
1690		self.Ang2Q = self._A2QConversion
	1690	if "qconv" not in keyargs:
	1691	self._A2QConversion = QConversion(['z-','x+'],['x+','z-'],[0,1,0],wl=self._wl) # 2S+2D goniometer
	1692	self.Ang2Q = self._A2QConversion
1691	1693
1692	1694	def _set_transform(self,v1,v2,v3):
1693	1695	"""

1805	1807	# dummy function to have some documentation string available
1806	1808	# the real function is generated dynamically in the __init__ routine
1807	1809	pass
1808
1809		class GID_ID10B(GID):
1810		"""
1811		class describing grazing incidence x-ray diffraction experiments
1812		the class helps with calculating the angles of Bragg reflections
1813		as well as it helps with analyzing measured data
1814
1815		the class describes a four circle (theta,omega,delta,gamma)
1816		goniometer to help with GID experiments at ID10B / ESRF.
1817		3D data can be treated with the use of linear and area detectors.
1818		see help self.Ang2Q
1819		"""
1820		def __init__(self,idir,ndir,**keyargs):
1821		"""
1822		initialization routine for the GID Experiment class
1823
1824		idir defines the inplane reference direction (idir points into the PB
1825		direction at zero angles)
1826
1827		Parameters
1828		----------
1829		same as for the Experiment base class
1830
1831		"""
1832		Experiment.__init__(self,idir,ndir,**keyargs)
1833
1834		# initialize Ang2Q conversion
1835		self._A2QConversion = QConversion(['x+','z-'],['x+','z-'],[0,1,0],wl=self._wl) # 2S+2D goniometer
1836		self.Ang2Q = self._A2QConversion
1837
1838		def Ang2Q(self,th,om,delta,gamma,**kwargs):
1839		"""
1840		angular to momentum space conversion for a point detector. Also see
1841		help GID_ID10B.Ang2Q for procedures which treat line and area detectors
1842
1843		Parameters
1844		----------
1845		th,om,delta,gamma: sample and detector angles as numpy array, lists or Scalars
1846		must be given. all arguments must have the same shape or
1847		length
1848
1849		**kwargs: optional keyword arguments
1850		delta: giving delta angles to correct the given ones for misalignment
1851		delta must be an numpy array or list of length 4.
1852		used angles are than th,om,delta,gamma - delta
1853		UB: matrix for conversion from (hkl) coordinates to Q of sample
1854		used to determine not Q but (hkl)
1855		(default: identity matrix)
1856		wl: x-ray wavelength in angstroem (default: self._wl)
1857		deg: flag to tell if angles are passed as degree (default: True)
1858
1859		Returns
1860		-------
1861		reciprocal space positions as numpy.ndarray with shape ( * , 3 )
1862		where * corresponds to the number of points given in the input
1863		"""
1864		# dummy function to have some documentation string available
1865		# the real function is generated dynamically in the __init__ routine
1866		pass
1867
1868		def Q2Ang(self,Q,trans=True,deg=True,**kwargs):
1869		"""
1870		calculate the GID angles needed in the experiment
1871		the inplane reference direction defines the direction were
1872		the reference direction is parallel to the primary beam
1873		(i.e. lattice planes perpendicular to the beam)
1874
1875		Parameters
1876		----------
1877		Q: a list or numpy array of shape (3) with
1878		q-space vector components
1879
1880		optional keyword arguments:
1881		trans: True/False apply coordinate transformation on Q
1882		deg: True/Flase (default True) determines if the
1883		angles are returned in radians or degrees
1884
1885		Returns
1886		-------
1887		a numpy array of shape (4) with the four GID scattering angles which are
1888		(theta,omega,delta,gamma)
1889
1890		theta: incidence angle to surface (at the moment always 0)
1891		omega: sample rotation with respect to the inplane reference direction
1892		delta: exit angle from surface (at the moment always 0)
1893		gamma: scattering angle
1894		"""
1895
1896		for k in kwargs.keys():
1897		if k not in ['trans','deg']:
1898		raise Exception("unknown keyword argument given: allowed are 'trans': coordinate transformation flag, 'deg': degree-flag")
1899
1900		[ai,azi,tt,beta] = GID.Q2Ang(self, Q, trans,deg, **kwargs)
1901
1902		return [0,azi,0,tt]
1903
1904	1810
1905	1811	class GISAXS(Experiment):
1906	1812	"""

-0

xrayutilities/io/__init__.py less more

14	14	#
15	15	# Copyright (C) 2009-2010 Eugen Wintersberger <eugen.wintersberger@desy.de>
16	16	# Copyright (C) 2009-2013 Dominik Kriegner <dominik.kriegner@gmail.com>
	17
	18	from .helper import xu_open
17	19
18	20	from .radicon import rad2hdf5
19	21	from .radicon import hst2hdf5

37	39	from .spec import geth5_scan as geth5_map
38	40
39	41	from .edf import EDFFile
	42	from .edf import EDFDirectory
	43	from .cbf import CBFFile
	44	from .cbf import CBFDirectory
40	45
41	46	from .spectra import Spectra
42	47

45	50
46	51	# parser for the alignment log file of the rotating anode
47	52	from .rotanode_alignment import RA_Alignment
	53
	54	from .desy_tty08 import tty08File
	55	from .desy_tty08 import gettty08_scan

+220

-0

xrayutilities/io/cbf.py less more

	0	# This file is part of xrayutilities.
	1	#
	2	# xrayutilities is free software; you can redistribute it and/or modify
	3	# it under the terms of the GNU General Public License as published by
	4	# the Free Software Foundation; either version 2 of the License, or
	5	# (at your option) any later version.
	6	#
	7	# This program is distributed in the hope that it will be useful,
	8	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	9	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	10	# GNU General Public License for more details.
	11	#
	12	# You should have received a copy of the GNU General Public License
	13	# along with this program; if not, see <http://www.gnu.org/licenses/>.
	14	#
	15	# Copyright (C) 2013 Dominik Kriegner <dominik.kriegner@gmail.com>
	16
	17	#module for handling files stored in the CBF data format
	18
	19	import numpy
	20	import os
	21	import os.path
	22	import gzip
	23	import re
	24
	25	from .. import cxrayutilities
	26	from .. import config
	27
	28	cbf_name_start_num=re.compile(r"^\d")
	29
	30	class CBFFile(object):
	31	def __init__(self,fname,nxkey="X-Binary-Size-Fastest-Dimension",nykey="X-Binary-Size-Second-Dimension",dtkey="DataType",path=None):
	32	"""
	33	CBF detector image parser
	34
	35	required arguments:
	36	fname ................ name of the CBF file of type .cbf or .cbf.gz
	37
	38	keyword arguments:
	39	nxkey ................ name of the header key that holds the number of points in x-direction
	40	nykey ................ name of the header key that holds the number of points in y-direction
	41	dtkey ................ name of the header key that holds the datatype for the binary data
	42	path ................. path to the CBF file
	43	"""
	44
	45	self.filename = fname
	46	if path:
	47	self.full_filename = os.path.join(path,fname)
	48	else:
	49	self.full_filename = self.filename
	50
	51	#evaluate keyword arguments
	52	self.nxkey = nxkey
	53	self.nykey = nykey
	54	self.dtkey = dtkey
	55
	56	#create attributes for holding data
	57	self.data = None
	58	self._open()
	59	self.ReadData()
	60	self.fid.close()
	61
	62	def _open(self):
	63	"""
	64	open data file for reading
	65	"""
	66	try:
	67	if os.path.splitext(self.full_filename)[-1] == '.gz':
	68	self.fid = gzip.open(self.full_filename,"rb")
	69	else :
	70	self.fid = open(self.full_filename,"rb")
	71	except:
	72	raise IOError("cannot open file %s" %(self.full_filename))
	73
	74	def ReadData(self):
	75	"""
	76	Read the CCD data into the .data object
	77	this function is called by the initialization
	78	"""
	79	wasclosed = False
	80	if self.fid.closed:
	81	wasclosed = True
	82	self._open()
	83
	84	tmp = numpy.fromfile(file=self.fid,dtype="u1").tostring()
	85
	86	# read header information
	87	pos = tmp.index(self.nxkey+':')+len(self.nxkey+':')
	88	self.xdim = int(tmp[pos:pos+6].strip())
	89	pos = tmp.index(self.nykey+':')+len(self.nykey+':')
	90	self.ydim = int(tmp[pos:pos+6].strip())
	91
	92	self.data = cxrayutilities.cbfread(tmp,self.xdim,self.ydim).reshape((self.ydim,self.xdim))
	93	if wasclosed:
	94	self.fid.close()
	95
	96	def Save2HDF5(self,h5,group="/",comp=True):
	97	"""
	98	Saves the data stored in the EDF file in a HDF5 file as a HDF5 array.
	99	By default the data is stored in the root group of the HDF5 file - this
	100	can be changed by passing the name of a target group or a path to the
	101	target group via the "group" keyword argument.
	102
	103	required arguments.
	104	h5 ................... a HDF5 file object
	105
	106	optional keyword arguments:
	107	group ................ group where to store the data (default to the root of the file)
	108	comp ................. activate compression - true by default
	109	"""
	110
	111	if isinstance(group,str):
	112	g = h5.getNode(group)
	113	else:
	114	g = group
	115
	116	#create the array name
	117	ca_name = os.path.split(self.filename)[-1]
	118	ca_name = os.path.splitext(ca_name)[0]
	119	ca_name = os.path.splitext(ca_name)[0] # perform a second time for case of .cbf.gz files
	120	ca_name = ca_name.replace("-","_")
	121	if cbf_name_start_num.match(ca_name):
	122	ca_name = "ccd_"+ca_name
	123	if config.VERBOSITY >= config.INFO_ALL:
	124	print(ca_name)
	125	ca_name = ca_name.replace(" ","_")
	126
	127	#create the array description
	128	ca_desc = "CBF CCD data from file %s " %(self.filename)
	129
	130	#create the Atom for the array
	131	a = tables.Atom.from_dtype(self.data.dtype)
	132	f = tables.Filters(complevel=7,complib="zlib",fletcher32=True)
	133	if comp:
	134	try:
	135	ca = h5.createCArray(g,ca_name,a,self.data.shape,ca_desc,filters=f)
	136	except:
	137	h5.removeNode(g,ca_name,recursive=True)
	138	ca = h5.createCArray(g,ca_name,a,self.data.shape,ca_desc,filters=f)
	139	else:
	140	try:
	141	ca = h5.createCArray(g,ca_name,a,self.data.shape,ca_desc)
	142	except:
	143	h5.removeNode(g,ca_name,recursive=True)
	144	ca = h5.createCArray(g,ca_name,a,self.data.shape,ca_desc)
	145
	146	#write the data
	147	ca[...] = self.data[...]
	148
	149
	150	class CBFDirectory(object):
	151	"""
	152	Parses a directory for CBF files, which can be stored to a HDF5 file for further usage
	153	"""
	154	def __init__(self,datapath,ext="cbf",**keyargs):
	155	"""
	156
	157	required arguments:
	158	datapath ............. directory of the CBF file
	159
	160	optional keyword arguments:
	161	ext .................. extension of the ccd files in the datapath (default: "cbf")
	162
	163	further keyword arguments are passed to CBFFile
	164	"""
	165
	166	self.datapath = os.path.normpath(datapath)
	167	self.extension = ext
	168
	169	#create list of files to read
	170	self.files = glob.glob( os.path.join(self.datapath, '*.%s' %(self.extension)))
	171
	172	if len(self.files) == 0:
	173	print("XU.io.CBFDirectory: no files found in %s" %(self.datapath))
	174	return
	175
	176	if config.VERBOSITY >= config.INFO_ALL:
	177	print("XU.io.CBFDirectory: %d files found in %s" %(len(self.files),self.datapath))
	178
	179	self.init_keyargs = keyargs
	180
	181
	182	def Save2HDF5(self,h5,group="",comp=True):
	183	"""
	184	Saves the data stored in the CBF files in the specified directory
	185	in a HDF5 file as a HDF5 arrays in a subgroup.
	186	By default the data is stored in a group given by the foldername - this
	187	can be changed by passing the name of a target group or a path to the
	188	target group via the "group" keyword argument.
	189
	190	required arguments.
	191	h5 ................... a HDF5 file object
	192
	193	optional keyword arguments:
	194	group ................ group where to store the data (defaults to pathname if group is empty string)
	195	comp ................. activate compression - true by default
	196	"""
	197
	198	if isinstance(group,str):
	199	if group == "":
	200	group = os.path.split(self.datapath)[1]
	201	try:
	202	g = h5.getNode(h5.root,group)
	203	except:
	204	g = h5.createGroup(h5.root,group)
	205	else:
	206	g = group
	207
	208	if "comp" in keyargs:
	209	compflag = keyargs["comp"]
	210	else:
	211	compflag = True
	212
	213	for infile in self.files:
	214	# read EDFFile and save to hdf5
	215	filename = os.path.split(infile)[1]
	216	e = CBFFile(filename,path=self.datapath,**self.init_keyargs)
	217	#e.ReadData()
	218	e.Save2HDF5(h5,group=g)
	219

+208

-0

xrayutilities/io/desy_tty08.py less more

	0	# This file is part of xrayutilities.
	1	#
	2	# xrayutilities is free software; you can redistribute it and/or modify
	3	# it under the terms of the GNU General Public License as published by
	4	# the Free Software Foundation; either version 2 of the License, or
	5	# (at your option) any later version.
	6	#
	7	# This program is distributed in the hope that it will be useful,
	8	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	9	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	10	# GNU General Public License for more details.
	11	#
	12	# You should have received a copy of the GNU General Public License
	13	# along with this program; if not, see <http://www.gnu.org/licenses/>.
	14	#
	15	# Copyright (C) 2013 Dominik Kriegner <dominik.kriegner@gmail.com>
	16
	17
	18	"""
	19	class for reading data+header information from tty08 data files
	20
	21	tty08 is system used at beamline P08 at Hasylab Hamburg and creates simple ASCII files to save the data. Information is easily read from the multicolumn data file.
	22	the functions below enable also to parse the information of the header
	23	"""
	24
	25	import re
	26	import numpy
	27	import os
	28	import matplotlib
	29
	30	# relative imports from xrayutilities
	31	from .helper import xu_open
	32	from .. import config
	33	from ..exception import InputError
	34
	35	re_columns= re.compile(r"/\*H")
	36	re_command = re.compile(r"^/\*C command")
	37	re_comment = re.compile(r"^/\*")
	38	re_date = re.compile(r"^/\*D date")
	39	re_epoch = re.compile(r"^/\*T epoch")
	40	re_initmopo = re.compile(r"^/\*M")
	41
	42	class tty08File(object):
	43	"""
	44	Represents a tty08 data file. The file is read during the
	45	Constructor call. This class should work for data stored at
	46	beamline P08 using the tty08 acquisition system.
	47
	48	Required constructor arguments:
	49	------------------------------
	50	filename: a string with the name of the tty08-file
	51
	52	Optional keyword arguments:
	53	--------------------------
	54	mcadir .................. directory name of MCA files
	55
	56	"""
	57
	58	def __init__(self,filename,path=None,mcadir=None):
	59	self.filename = filename
	60	if path == None:
	61	self.full_filename = self.filename
	62	else:
	63	self.full_filename = os.path.join(path,self.filename)
	64
	65	self.Read()
	66
	67	if mcadir!=None:
	68	self.mca_directory= mcadir
	69	self.mca_files = sorted(glob.glob(os.path.join(self.mca_directory,'*')))
	70
	71	if len(self.mca_files):
	72	self.ReadMCA()
	73
	74	def ReadMCA(self):
	75
	76	mca = numpy.empty((len(raws),numpy.loadtxt(raws[0]).shape[0]),dtype=numpy.float)
	77	for i in range(len(raws)):
	78	mca[i,:] = numpy.loadtxt(self.mca_files[i])[:,1]
	79
	80	fname = self.mca_file_template %i
	81	data = numpy.loadtxt(fname)
	82
	83	if i==self.mca_start_index:
	84	if len(data.shape)==2:
	85	self.mca_channels = data[:,0]
	86	else:
	87	self.mca_channels = numpy.arange(0,data.shape[0])
	88
	89	if len(data.shape)==2:
	90	dlist.append(data[:,1].tolist())
	91	else:
	92	dlist.append(data.tolist())
	93
	94	self.mca= mca
	95	self.data = matplotlib.mlab.rec_append_fields(self.data,'MCA',self.mca,dtypes=[(numpy.double,self.mca.shape[1])])
	96
	97	def Read(self):
	98	"""
	99	Read the data from the file
	100	"""
	101
	102	with xu_open(self.full_filename) as fid:
	103	# read header
	104	self.init_mopo = {}
	105	while True:
	106	line = fid.readline()
	107	#if DEGUG: print line
	108	if not line:
	109	break
	110
	111	if re_command.match(line):
	112	m = line.split(':')
	113	self.scan_command = m[1].strip()
	114	if re_date.match(line):
	115	m = line.split(':',1)
	116	self.scan_date = m[1].strip()
	117	if re_epoch.match(line):
	118	m = line.split(':',1)
	119	self.epoch = float(m[1])
	120	if re_initmopo.match(line):
	121	m = line[3:]
	122	m = m.split(';')
	123	for e in m:
	124	e= e.split('=')
	125	self.init_mopo[e[0].strip()] = float(e[1])
	126
	127	if re_columns.match(line):
	128	self.columns = tuple(line.split()[1:])
	129	break # here all necessary information is read and we can start reading the data
	130	self.data = numpy.loadtxt(fid,comments="/")
	131
	132	self.data = numpy.rec.fromrecords(self.data,names=self.columns)
	133
	134	def gettty08_scan(scanname,scannumbers,*args):
	135	"""
	136	function to obtain the angular cooridinates as well as intensity values
	137	saved in TTY08 datafiles. Especially usefull for reciprocal space map measurements,
	138	and to combine date from several scans
	139
	140	further more it is possible to obtain even more positions from
	141	the data file if more than two string arguments with its names are given
	142
	143	Parameters
	144	----------
	145	scanname: name of the scans, for multiple scans this needs to be a template string
	146	scannumber: number of the scans of the reciprocal space map (int,tuple or list)
	147
	148	*args: names of the motors (optional) (strings)
	149	to read reciprocal space maps measured in coplanar diffraction give:
	150	omname: e.g. name of the omega motor (or its equivalent)
	151	ttname: e.g. name of the two theta motor (or its equivalent)
	152
	153	Returns
	154	-------
	155	MAP
	156
	157	or
	158
	159	[ang1,ang2,...],MAP:
	160	angular positions of the center channel of the position
	161	sensitive detector (numpy.ndarray 1D) together with all the
	162	data values as stored in the data file (includes the
	163	intensities e.g. MAP['MCA']).
	164
	165	Example
	166	-------
	167	>>> [om,tt],MAP = xu.io.gettty08_scan('text%05d.dat',36,'omega','gamma')
	168	"""
	169
	170	if isinstance(scannumbers,(list,tuple)):
	171	scanlist = scannumbers
	172	else:
	173	scanlist = list([scannumbers])
	174
	175	angles = dict.fromkeys(args)
	176	for key in angles.keys():
	177	if not isinstance(key,str):
	178	raise InputError("*arg values need to be strings with motornames")
	179	angles[key] = numpy.zeros(0)
	180	buf=numpy.zeros(0)
	181	MAP = numpy.zeros(0)
	182
	183	for nr in scanlist:
	184	scan = tty08File(scanname%nr)
	185	sdata = scan.data
	186	if MAP.dtype == numpy.float64: MAP.dtype = sdata.dtype
	187	# append scan data to MAP, where all data are stored
	188	MAP = numpy.append(MAP,sdata)
	189	#check type of scan
	190	for i in range(len(args)):
	191	motname = args[i]
	192	scanlength = len(sdata)
	193	try:
	194	buf = sdata[motname]
	195	except:
	196	buf = scan.init_mopos[motname]*numpy.ones(scanlength)
	197	angles[motname] =numpy.concatenate((angles[motname],buf))
	198
	199	retval = []
	200	for motname in args:
	201	#create return values in correct order
	202	retval.append(angles[motname])
	203
	204	if len(args)==0:
	205	return MAP
	206	else: return retval,MAP
	207

+59

-0

xrayutilities/io/helper.py less more

	0	# This file is part of xrayutilities.
	1	#
	2	# xrayutilities is free software; you can redistribute it and/or modify
	3	# it under the terms of the GNU General Public License as published by
	4	# the Free Software Foundation; either version 2 of the License, or
	5	# (at your option) any later version.
	6	#
	7	# This program is distributed in the hope that it will be useful,
	8	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	9	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	10	# GNU General Public License for more details.
	11	#
	12	# You should have received a copy of the GNU General Public License
	13	# along with this program; if not, see <http://www.gnu.org/licenses/>.
	14	#
	15	# Copyright (C) 2013 Dominik Kriegner <dominik.kriegner@gmail.com>
	16
	17
	18	"""
	19	convenience functions to open files for various data file reader
	20
	21	these functions should be used in new parsers since they transparently allow to open gzipped and bzipped files
	22	"""
	23
	24	import os
	25	import gzip
	26	import bz2
	27
	28	from .. import config
	29	from ..exception import InputError
	30
	31	def xu_open(filename,mode='rb'):
	32	"""
	33	function to open a file no matter if zipped or not. Files with extension
	34	'.gz' or '.bz2' are assumed to be compressed and transparently opened to read like
	35	usual files.
	36
	37	Parameters
	38	----------
	39	filename: filename of the file to open (full including path)
	40	mode: mode in which the file should be opened
	41
	42	Returns
	43	-------
	44	file handle of the opened file
	45
	46	If the file does not exist an IOError is raised by the open routine, which is not
	47	caught within the function
	48	"""
	49
	50	if os.path.splitext(filename)[-1] == '.gz':
	51	fid = gzip.open(filename,mode)
	52	elif os.path.splitext(filename)[-1] == '.bz2':
	53	fid = bz2.BZ2File(filename,mode)
	54	else:
	55	fid = open(filename,mode)
	56
	57	return fid
	58

+14

-4

xrayutilities/io/spec.py less more

13	13	# along with this program; if not, see <http://www.gnu.org/licenses/>.
14	14	#
15	15	# Copyright (C) 2009-2010 Eugen Wintersberger <eugen.wintersberger@desy.de>
16		# Copyright (C) 2009-2012 Dominik Kriegner <dominik.kriegner@gmail.com>
	16	# Copyright (C) 2009-2013 Dominik Kriegner <dominik.kriegner@gmail.com>
17	17
18	18	"""
19	19	a threaded class for observing a SPEC data file

31	31	import os
32	32	import time
33	33	import tables
	34	import gzip
34	35
35	36	# relative imports from xrayutilities
36	37	from .. import config

550	551	self.scan_list = []
551	552	#open the file for reading
552	553	try:
553		self.fid = open(self.full_filename,"rb")
	554	if os.path.splitext(self.full_filename)[-1] == '.gz':
	555	self.fid = gzip.open(self.full_filename,"rb")
	556	else :
	557	self.fid = open(self.full_filename,"rb")
554	558	self.last_offset = self.fid.tell()
555	559	except:
556	560	self.fid = None

627	631
628	632	self.fid.close()
629	633	try:
630		self.fid = open(self.full_filename,"r")
	634	if os.path.splitext(self.full_filename)[-1] == '.gz':
	635	self.fid = gzip.open(self.full_filename,"rb")
	636	else :
	637	self.fid = open(self.full_filename,"rb")
631	638	except:
632	639	self.fid = None
633	640	raise IOError("error opening SPEC file %s" %(self.full_filename))

871	878	self.full_filename = os.path.join(path,self.filename)
872	879
873	880	try:
874		self.fid = open(self.full_filename,"r")
	881	if os.path.splitext(self.full_filename)[-1] == '.gz':
	882	self.fid = gzip.open(self.full_filename,"r")
	883	else :
	884	self.fid = open(self.full_filename,"r")
875	885	except:
876	886	raise IOError("cannot open log file %s" %(self.full_filename))
877	887

+68

-0

xrayutilities/materials/lattice.py less more

807	807	l = Lattice(a1,a2,a3,base=lb)
808	808
809	809	return l
	810
	811	def MagnetiteLattice( aa , ab , ac , a ,x = 0.255):
	812	lb = LatticeBase()
	813	#Fe1
	814	lb.append(aa,[0.125,0.125,0.125])
	815	lb.append(aa,[0.875,0.375,0.375])
	816	lb.append(aa,[0.375,.875,.375])
	817	lb.append(aa,[0.375,0.375,.875])
	818	lb.append(aa,[.875,.875,.875])
	819	lb.append(aa,[.125,.625,.625])
	820	lb.append(aa,[.625,.125,.625])
	821	lb.append(aa,[.625,.625,.125])
	822	#Fe2
	823	lb.append(ab,[0.5,0.5,0.5])
	824	lb.append(ab,[.5,0.75,0.75])
	825	lb.append(ab,[0.75,0.75,0.5])
	826	lb.append(ab,[0.75,0.5,0.75])
	827	lb.append(ab,[0.5,0.25,.25])
	828	lb.append(ab,[.25,.25,.5])
	829	lb.append(ab,[.25,.5,.25])
	830	lb.append(ab,[.5,0,0])
	831	lb.append(ab,[.75,.25,0])
	832	lb.append(ab,[.75,0,.25])
	833	lb.append(ab,[.25,.75,0])
	834	lb.append(ab,[.25,0,.75])
	835	lb.append(ab,[0,.5,0])
	836	lb.append(ab,[0,.75,.25])
	837	lb.append(ab,[0,.25,.75])
	838	lb.append(ab,[0,0,.5])
	839	#O
	840	lb.append(ac,[x,x,x])
	841	lb.append(ac,[1-x,x+0.25,x+0.25])
	842	lb.append(ac,[1-x+0.25,1-x+0.25,x])
	843	lb.append(ac,[x+0.25,1-x,x+0.25])
	844	lb.append(ac,[x,1-x+0.25,1-x+0.25])
	845	lb.append(ac,[x+0.25,x+0.25,1-x])
	846	lb.append(ac,[1-x+0.25,x,1-x+0.25])
	847	lb.append(ac,[1-x,1-x,1-x])
	848	lb.append(ac,[x,.75-x,.75-x])
	849	lb.append(ac,[x-.25,x-.25,1-x])
	850	lb.append(ac,[.75-x,x,.75-x])
	851	lb.append(ac,[1-x,x-.25,x-.25])
	852	lb.append(ac,[.75-x,.75-x,x])
	853	lb.append(ac,[x-.25,1-x,x-.25])
	854	lb.append(ac,[x,.5+x,.5+x])
	855	lb.append(ac,[1-x+.25,.25+x,.5+x])
	856	lb.append(ac,[.25+x,.5-x,x-.25])
	857	lb.append(ac,[x+.25,x-.25,.5-x])
	858	lb.append(ac,[1-x+.25,.5+x,1-.25-x])
	859	lb.append(ac,[1-x,.5-x,.5-x])
	860	lb.append(ac,[x-.25,.25+x,.5-x])
	861	lb.append(ac,[.75-x,.5+x,1-x+.25])
	862	lb.append(ac,[.75-x,1-x+.25,.5+x])
	863	lb.append(ac,[x-.25,.5-x,.25+x])
	864	lb.append(ac,[.5+x,x,.5+x])
	865	lb.append(ac,[.5-x,.25+x,x-.25])
	866	lb.append(ac,[.5+x,1-x+.25,.75-x])
	867	lb.append(ac,[.5-x,1-x,.5-x])
	868	lb.append(ac,[.5+x,.75-x,1-x+.25])
	869	lb.append(ac,[.5-x,x-.25,x+.25])
	870	lb.append(ac,[.5+x,.5+x,x])
	871	lb.append(ac,[.5-x,.5-x,1-x])
	872
	873	a1=[ a , 0 , 0 ]
	874	a2=[ 0 , a , 0 ]
	875	a3=[ 0 , 0 , a ]
	876	l = Lattice(a1,a2,a3,base=lb)
	877	return l⏎

-0

xrayutilities/materials/predefined_materials.py less more

54	54	CubicElasticTensor(93.6,7.7,13.4))
55	55	PbSe = Material("PbSe",lattice.RockSalt_Cubic_Lattice(elements.Pb,elements.Se,6.128),
56	56	CubicElasticTensor(123.7,19.3,15.9))
	57	NaCl = Material("NaCl",lattice.RockSalt_Cubic_Lattice(elements.Na,elements.Cl,5.6402))
57	58	GaN = Material("GaN",lattice.WurtziteLattice(elements.Ga,elements.N,3.189,5.186),
58	59	HexagonalElasticTensor(390.e9,145.e9,106.e9,398.e9,105.e9),thetaDebye=600)
59	60	BaF2 = Material("BaF2",lattice.CubicFm3mBaF2(elements.Ba,elements.F,6.2001))
	61	SrF2 = Material("SrF2",lattice.CubicFm3mBaF2(elements.Sr,elements.F,5.8007))
60	62	MnTe = Material("MnTe",lattice.NiAsLattice(elements.Mn,elements.Te,4.1429,6.7031))
61	63	GeTe = Material("GeTe",lattice.GeTeRhombohedral(elements.Ge,elements.Te,5.996,88.18,0.237))
62	64	Al = Material("Al",lattice.FCCLattice(elements.Al,4.04958))

75	77	CuMnAs = Material("CuMnAs",lattice.CuMnAsLattice(elements.Cu,elements.Mn,elements.As,3.82,3.82,6.30))
76	78	CaTiO3 = Material("CaTiO3",lattice.PerovskiteTypeRhombohedral(elements.Ca,elements.Ti,elements.O,3.795,90))
77	79	BiFeO3 = Material("BiFeO3",lattice.PerovskiteTypeRhombohedral(elements.Bi,elements.Fe,elements.O,3.965,89.3))
	80	FeO = Material("FeO",lattice.RockSalt_Cubic_Lattice(elements.Fe,elements.O, 4.332))
	81	CoO = Material("CoO",lattice.RockSalt_Cubic_Lattice(elements.Co,elements.O, 4.214))
	82	Fe3O4 = Material("Fe3O4",lattice.MagnetiteLattice(elements.Fe,elements.Fe,elements.O, 8.3958))
	83	Co3O4 = Material("Co3O4",lattice.MagnetiteLattice(elements.Co,elements.Co,elements.O, 8.0821))
78	84
79	85	# materials defined from cif file
80	86	try:

-0

xrayutilities/math/__init__.py less more

46	46	from .functions import Lorentz1d_der_p
47	47	from .functions import Lorentz2d
48	48	from .functions import PseudoVoigt1d
	49	from .functions import PseudoVoigt1dArea
49	50
50	51	from .fit import fit_peak2d
51	52	from .fit import gauss_fit
	53	from .fit import peak_fit
52	54	from .fit import multPeakFit
53	55	from .fit import multPeakPlot
54	56	from .fit import multGaussFit

+59

-18

xrayutilities/math/fit.py less more

27	27	from scipy.odr import models
28	28
29	29	from .. import config
	30	from .. exception import InputError
30	31	from .functions import Gauss1d,Gauss1d_der_x,Gauss1d_der_p
31	32	from .functions import Lorentz1d,Lorentz1d_der_x,Lorentz1d_der_p
	33	from .functions import PseudoVoigt1d
32	34
33	35	try:
34	36	from matplotlib import pyplot as plt

36	38	if config.VERBOSITY >= config.INFO_ALL:
37	39	print("XU.analysis.sample_align: warning; plotting functionality not available")
38	40
39		def gauss_fit(xdata,ydata,iparams=[],maxit=200):
40		"""
41		Gauss fit function using odr-pack wrapper in scipy similar to
	41	def peak_fit(xdata,ydata,iparams=[],peaktype='Gauss',maxit=200):
	42	"""
	43	fit function using odr-pack wrapper in scipy similar to
42	44	https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fitting
	45	for Gauss, Lorentz or Pseudovoigt-functions
43	46
44	47	Parameters
45	48	----------

47	50	ydata: ycoordinates of the data which should be fit
48	51
49	52	keyword parameters:
50		iparams: initial paramters for the fit (determined automatically if nothing is given
	53	iparams: initial paramters for the fit (determined automatically if nothing is given)
	54	peaktype: type of peak to fit ('Gauss','Lorentz','PseudoVoigt')
51	55	maxit: maximal iteration number of the fit
52	56
53	57	Returns
54	58	-------
55	59	params,sd_params,itlim
56	60
57		the Gauss parameters as defined in function Gauss1d(x, *param)
	61	the parameters as defined in function Gauss1d/Lorentz1d or PseudoVoigt1d(x, *param)
58	62	and their errors of the fit, as well as a boolean flag which is false in the case of a
59	63	successful fit
60	64	"""
61	65
62		gfunc = lambda param,x: Gauss1d(x, *param)
63		gfunc_dx = lambda param,x: Gauss1d_der_x(x, *param)
64		gfunc_dp = lambda param,x: Gauss1d_der_p(x, *param)
	66	if peaktype=='Gauss':
	67	gfunc = lambda param,x: Gauss1d(x, *param)
	68	gfunc_dx = lambda param,x: Gauss1d_der_x(x, *param)
	69	gfunc_dp = lambda param,x: Gauss1d_der_p(x, *param)
	70	elif peaktype=='Lorentz':
	71	gfunc = lambda param,x: Lorentz1d(x, *param)
	72	gfunc_dx = lambda param,x: Lorentz1d_der_x(x, *param)
	73	gfunc_dp = lambda param,x: Lorentz1d_der_p(x, *param)
	74	elif peaktype=='PseudoVoigt':
	75	gfunc = lambda param,x: PseudoVoigt1d(x, *param)
	76	gfunc_dx = None
	77	gfunc_dp = None
	78	else:
	79	raise InputError("keyword rgument peaktype takes invalid value!")
65	80
66	81	if not any(iparams):
67	82	cen = numpy.sum(xdata*ydata)/numpy.sum(ydata)
68		iparams = numpy.array([cen,\
	83	iparams = [cen,\
69	84	numpy.sqrt(numpy.abs(numpy.sum((xdata-cen)*2ydata)/numpy.sum(ydata))),\
70	85	numpy.max(ydata),\
71		numpy.min(ydata)])
	86	numpy.min(ydata)]
	87	if peaktype=='PseudoVoigt':
	88	iparams.append(0.5) # set ETA parameter to be between Gauss and Lorentz shape
72	89
73	90	if config.VERBOSITY >= config.DEBUG:
74		print("XU.math.gauss_fit: iparams: [%f %f %f %f]" %tuple(iparams))
75
76		gauss = odr.Model(gfunc, fjacd=gfunc_dx, fjacb=gfunc_dp)
	91	print("XU.math.peak_fit: iparams: %s"%str(tuple(iparams)))
	92
	93	peak = odr.Model(gfunc, fjacd=gfunc_dx, fjacb=gfunc_dp)
77	94
78	95	sy = numpy.sqrt(ydata)
79	96	sy[sy==0] = 1
80	97	mydata = odr.RealData(xdata, ydata, sy=sy)
81	98
82		myodr = odr.ODR(mydata, gauss, beta0=iparams,maxit=maxit)
	99	myodr = odr.ODR(mydata, peak, beta0=iparams,maxit=maxit)
83	100
84	101	# use least-square fit
85	102	myodr.set_job(fit_type=2)

93	110	#fit.pprint() # prints final message from odrpack
94	111
95	112	if config.VERBOSITY >= config.DEBUG:
96		print("XU.math.gauss_fit: params: [%f %f %f %f]" %tuple(fit.beta))
97		print("XU.math.gauss_fit: params std: [%f %f %f %f]" %tuple(fit.sd_beta))
98		print("XU.math.gauss_fit: %s" %fit.stopreason[0])
	113	print("XU.math.peak_fit: params: %s" %str(tuple(fit.beta)))
	114	print("XU.math.peak_fit: params std: %s" %str(tuple(fit.sd_beta)))
	115	print("XU.math.peak_fit: %s" %fit.stopreason[0])
99	116
100	117	itlim = False
101	118	if fit.stopreason[0] == 'Iteration limit reached':
102	119	itlim = True
103	120	if config.VERBOSITY >= config.INFO_LOW:
104		print("XU.math.gauss_fit: Iteration limit reached, do not trust the result!")
	121	print("XU.math.peak_fit: Iteration limit reached, do not trust the result!")
105	122
106	123	return fit.beta, fit.sd_beta, itlim
107	124
	125	def gauss_fit(xdata,ydata,iparams=[],maxit=200):
	126	"""
	127	Gauss fit function using odr-pack wrapper in scipy similar to
	128	https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fitting
	129
	130	Parameters
	131	----------
	132	xdata: xcoordinates of the data to be fitted
	133	ydata: ycoordinates of the data which should be fit
	134
	135	keyword parameters:
	136	iparams: initial paramters for the fit (determined automatically if nothing is given
	137	maxit: maximal iteration number of the fit
	138
	139	Returns
	140	-------
	141	params,sd_params,itlim
	142
	143	the Gauss parameters as defined in function Gauss1d(x, *param)
	144	and their errors of the fit, as well as a boolean flag which is false in the case of a
	145	successful fit
	146	"""
	147
	148	return peak_fit(xdata,ydata,iparams=iparams,peaktype='Gauss',maxit=maxit)
108	149
109	150	def fit_peak2d(x,y,data,start,drange,fit_function,maxfev=2000):
110	151	"""

+22

-2

xrayutilities/math/functions.py less more

214	214	----------
215	215	p: list of parameters of the Lorentz-function
216	216	[XCEN,FWHM,AMP,BACKGROUND]
217		x,y: coordinate(s) where the function should be evaluated
	217	x: coordinate(s) where the function should be evaluated
218	218
219	219	Returns
220	220	-------

287	287	p: list of parameters of the Lorentz-function
288	288	[XCEN,FWHM,AMP,BACKGROUND,ETA]
289	289	ETA: 0 ...1 0 means pure Gauss and 1 means pure Lorentz
290		x,y: coordinate(s) where the function should be evaluated
	290	x: coordinate(s) where the function should be evaluated
291	291
292	292	Returns
293	293	-------

300	300
301	301	return f
302	302
	303	def PseudoVoigt1dArea(*p):
	304	"""
	305	function to calculate the area of a pseudo Voigt function with neglected background
	306
	307	Parameters
	308	----------
	309	p: list of parameters of the Lorentz-function
	310	[XCEN,FWHM,AMP,BACKGROUND,ETA]
	311	ETA: 0 ...1 0 means pure Gauss and 1 means pure Lorentz
	312
	313	Returns
	314	-------
	315	the value of the PseudoVoigt described by the parameters p
	316	at position (x,y)
	317
	318	"""
	319
	320	f = p[2] * ( p[4]numpy.pi/(4./(p[1])) + (1.-p[4])numpy.sqrt(numpy.pi)/(numpy.sqrt(numpy.log(2))*2./(p[1])) )
	321
	322	return f
303	323
304	324	def Debye1(x):
305	325	"""

+17

-14

xrayutilities/src/block_average.c less more

20	20
21	21	#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
22	22	#define PY_ARRAY_UNIQUE_SYMBOL XU_UNIQUE_SYMBOL
	23	#define NO_IMPORT_ARRAY
23	24	#include <numpy/arrayobject.h>
24	25	#include <math.h>
25	26	#ifdef __OPENMP__

58	59	PyArrayObject input=NULL, outarr=NULL;
59	60	double cin,cout;
60	61	double buf;
	62	npy_intp nout;
61	63
62	64	// Python argument conversion code
63		if (!PyArg_ParseTuple(args, "O!i",&PyArray_Type, &input, &Nav)) return NULL;
64
	65	if (!PyArg_ParseTuple(args, "O!i",&PyArray_Type, &input, &Nav)) return NULL;
	66
65	67	PYARRAY_CHECK(input,1,NPY_DOUBLE,"input must be a 1D double array!");
66	68	N = PyArray_SIZE(input);
67	69	cin = (double *) PyArray_DATA(input);
68	70
69	71	// create output ndarray
70		npy_intp nout;
71	72	nout = ((int)ceil(N/(float)Nav));
72	73	outarr = (PyArrayObject *) PyArray_SimpleNew(1, &nout, NPY_DOUBLE);
73	74	cout = (double *) PyArray_DATA(outarr);
74
	75
75	76	// c-code following is performing the block averaging
76	77	for(i=0; i<N; i=i+Nav) {
77	78	buf=0;

81	82	}
82	83	cout[i/Nav] = buf/(float)(j-i); //save average to output array
83	84	}
84
	85
85	86	// return output array
86	87	return PyArray_Return(outarr);
87	88	}

92	93	* Parameters
93	94	* ----------
94	95	* ccd: input array/CCD frame
95		* size = (Nch2, Nch1)
	96	* size = (Nch2, Nch1)
96	97	* Nch1 is the fast varying index
97	98	* Nav1,2: number of channels to average in each dimension
98	99	* in total a block of Nav1 x Nav2 is averaged

101	102	* Returns
102	103	* -------
103	104	* block_av: block averaged output array
104		* size = (ceil(Nch2/Nav2) , ceil(Nch1/Nav1))
	105	* size = (ceil(Nch2/Nav2) , ceil(Nch1/Nav1))
105	106	*
106	107	*/
107	108

112	113	PyArrayObject input=NULL, outarr=NULL;
113	114	double cin,cout;
114	115	double buf;
	116	npy_intp nout[2];
115	117
116	118	// Python argument conversion code
117	119	if (!PyArg_ParseTuple(args, "O!iiI",&PyArray_Type, &input, &Nav2, &Nav1, &nthreads)) return NULL;
118
	120
119	121	PYARRAY_CHECK(input,2,NPY_DOUBLE,"input must be a 2D double array!");
120	122	Nch2 = PyArray_DIMS(input)[0];
121	123	Nch1 = PyArray_DIMS(input)[1];
122	124	cin = (double *) PyArray_DATA(input);
123	125
124	126	// create output ndarray
125		npy_intp nout[2];
126	127	nout[0] = ((int)ceil(Nch2/(float)Nav2));
127	128	nout[1] = ((int)ceil(Nch1/(float)Nav1));
128	129	outarr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);

170	171	unsigned int nthreads; //number of threads to use
171	172	PyArrayObject input=NULL, outarr=NULL;
172	173	double cin,cout;
173		double buf;
	174	double buf;
	175	npy_intp nout[2];
174	176
175	177	// Python argument conversion code
176	178	if (!PyArg_ParseTuple(args, "O!iI",&PyArray_Type, &input, &Nav, &nthreads)) return NULL;
177
	179
178	180	PYARRAY_CHECK(input,2,NPY_DOUBLE,"input must be a 2D double array!");
179	181	Nspec = PyArray_DIMS(input)[0];
180	182	Nch = PyArray_DIMS(input)[1];
181	183	cin = (double *) PyArray_DATA(input);
182	184
183	185	// create output ndarray
184		npy_intp nout[2];
185	186	nout[0] = Nspec;
186	187	nout[1] = ((int)ceil(Nch/(float)Nav));
187	188	outarr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
188		cout = (double *) PyArray_DATA(outarr);
189
	189	cout = (double *) PyArray_DATA(outarr);
	190
190	191	#ifdef __OPENMP__
191	192	//set openmp thread numbers dynamically
192	193	OMPSETNUMTHREADS(nthreads);

208	209	// return output array
209	210	return PyArray_Return(outarr);
210	211	}
	212
	213	#undef PY_ARRAY_UNIQUE_SYMBOL

+23

-10

xrayutilities/src/cxrayutilities.c less more

42	42	extern PyObject* ang2q_conversion_area_pixel(PyObject self, PyObject args);
43	43	extern PyObject* ang2q_conversion_area_pixel2(PyObject self, PyObject args);
44	44
	45	/* functions from file_io.c */
	46	extern PyObject* cbfread(PyObject self, PyObject args);
	47
45	48	static PyMethodDef XRU_Methods[] = {
46	49	{"block_average1d", (PyCFunction)block_average1d, METH_VARARGS,
47	50	"block average for one-dimensional numpy double array\n\n"

52	55	"Returns\n"
53	56	"-------\n"
54	57	" block_av: block averaged output array\n"
55		" size = ceil(N/Nav) \n"
	58	" size = ceil(N/Nav) \n"
56	59	},
57	60	{"block_average2d", block_average2d, METH_VARARGS,
58	61	"2D block average for one CCD frame\n\n"

67	70	"Returns\n"
68	71	"-------\n"
69	72	" block_av: block averaged output array\n"
70		" size = (ceil(Nch2/Nav2) , ceil(Nch1/Nav1))\n"
	73	" size = (ceil(Nch2/Nav2) , ceil(Nch1/Nav1))\n"
71	74	},
72	75	{"block_average_PSD", block_average_PSD, METH_VARARGS,
73	76	"1D block average for a bunch of PSD spectra in a 2D array\n"

83	86	" block_av: block averaged output array\n"
84	87	" size = (Nspec , ceil(Nch/Nav)) (out)\n"
85	88	},
86		{"gridder2d",pygridder2d,METH_VARARGS,
87		"Function performs 2D gridding on 1D input data. \n\n"
	89	{"gridder2d",pygridder2d,METH_VARARGS,
	90	"Function performs 2D gridding on 1D input data. \n\n"
88	91	"Parameters\n"
89	92	"----------\n"
90	93	" x ...... input x-values (1D numpy array - float64)\n"

98	101	" ymax ... minimum y-value of the grid\n"
99	102	" out .... output data\n"
100	103	},
101		{"gridder3d",pygridder3d,METH_VARARGS,
102		"Function performs 2D gridding on 1D input data. \n\n"
	104	{"gridder3d",pygridder3d,METH_VARARGS,
	105	"Function performs 2D gridding on 1D input data. \n\n"
103	106	"Parameters\n"
104	107	"----------\n"
105	108	" x ...... input x-values (1D numpy array - float64)\n"

160	163	"Returns\n"
161	164	"-------\n"
162	165	" qpos ............ momentum transfer (Npoints*Nch,3)\n"
163		" \n"
	166	" \n"
164	167	},
165	168	{"ang2q_conversion_area", ang2q_conversion_area, METH_VARARGS,
166	169	"conversion of Npoints of goniometer positions to reciprocal space\n"

191	194	"Returns\n"
192	195	"-------\n"
193	196	" qpos ............ momentum transfer (NpointsNpix1Npix2,3)\n"
194		"\n"
	197	"\n"
195	198	},
196	199	{"ang2q_conversion_area_pixel", ang2q_conversion_area_pixel, METH_VARARGS,
197	200	"conversion of Npoints of detector positions to Q\n"

257	260	"-------\n"
258	261	" qpos ............ momentum transfer (Npoints,3)\n"
259	262	},
	263	{"cbfread", cbfread, METH_VARARGS,
	264	"parser for cbf data arrays from Pilatus detector images\n\n"
	265	" Parameters\n"
	266	" ----------\n"
	267	" data: data stream (character array)\n"
	268	" nx,ny: number of entries of the two dimensional image\n\n"
	269	" Returns\n"
	270	" -------\n"
	271	" the parsed data values as float ndarray\n"
	272	},
260	273	{NULL, NULL, 0, NULL} /* Sentinel */
261	274	};
262	275

279	292	#if PY_MAJOR_VERSION >= 3
280	293	PyInit_cxrayutilities(void)
281	294	#else
282		initcxrayutilities(void)
	295	initcxrayutilities(void)
283	296	#endif
284	297	{
285	298	PyObject *m;

287	300	#if PY_MAJOR_VERSION >= 3
288	301	m = PyModule_Create(&moduledef);
289	302	#else
290		m = Py_InitModule3("cxrayutilities", XRU_Methods,
	303	m = Py_InitModule3("cxrayutilities", XRU_Methods,
291	304	"Python C extension including performance critical parts\n"
292	305	"of xrayutilities (gridder, qconversion, block-averageing)\n");
293	306	#endif

+128

-0

xrayutilities/src/file_io.c less more

	0	/*
	1	* This file is part of xrayutilities.
	2	*
	3	* xrayutilities is free software; you can redistribute it and/or modify
	4	* it under the terms of the GNU General Public License as published by
	5	* the Free Software Foundation; either version 2 of the License, or
	6	* (at your option) any later version.
	7	*
	8	* This program is distributed in the hope that it will be useful,
	9	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	10	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	11	* GNU General Public License for more details.
	12	*
	13	* You should have received a copy of the GNU General Public License
	14	* along with this program; if not, see <http://www.gnu.org/licenses/>.
	15	*
	16	* Copyright (C) 2013 Dominik Kriegner <dominik.kriegner@gmail.com>
	17	*/
	18
	19	#include <Python.h>
	20
	21	#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
	22	#define PY_ARRAY_UNIQUE_SYMBOL XU_UNIQUE_SYMBOL
	23	#define NO_IMPORT_ARRAY
	24	#include <numpy/arrayobject.h>
	25
	26	#include <stdio.h>
	27
	28	PyObject* cbfread(PyObject self, PyObject args) {
	29	/* parser for cbf data arrays from Pilatus detector images
	30	*
	31	* Parameters
	32	* ----------
	33	* data: data stream (character array)
	34	* nx,ny: number of entries of the two dimensional image
	35	*
	36	* Returns
	37	* -------
	38	* the parsed data values as float ndarray
	39	*/
	40
	41	unsigned int i,start=0,nx,ny,len;
	42	unsigned int parsed = 0;
	43	PyArrayObject input=NULL, outarr=NULL;
	44	unsigned char *cin;
	45	float *cout;
	46	npy_intp nout;
	47
	48	int cur = 0;
	49	int diff = 0;
	50	unsigned int np = 0;
	51
	52	union {
	53	const unsigned char* uint8;
	54	const unsigned short* uint16;
	55	const unsigned int* uint32;
	56	const char* int8;
	57	const short* int16;
	58	const int* int32;
	59	} parser;
	60
	61	// Python argument conversion code
	62	if (!PyArg_ParseTuple(args, "s#ii", &cin, &len, &nx, &ny)) return NULL;
	63	/*printf("stream length: %d\n",len);
	64	printf("entries: %d %d\n",nx,ny);*/
	65
	66	// create output ndarray
	67	nout = nx*ny;
	68	outarr = (PyArrayObject *) PyArray_SimpleNew(1, &nout, NPY_FLOAT);
	69	cout = (float *) PyArray_DATA(outarr);
	70
	71	i = 0;
	72	while (i<len-10) { // find the start of the array
	73	if ((cin[i]==0x0c)&&(cin[i+1]==0x1a)&&(cin[i+2]==0x04)&&(cin[i+3]==0xd5)) {
	74	start = i+4;
	75	i = len+10;
	76	}
	77	i++;
	78	}
	79	if(i==len-10) {
	80	PyErr_SetString(PyExc_ValueError,"start of data in stream not found!\n");
	81	return NULL;
	82	}
	83	/*else {
	84	printf("found start at %d\n",start);
	85	}*/
	86
	87	// next while part was taken from pilatus code and adapted by O. Seeck and D. Kriegner
	88	parser.uint8 = (const unsigned char*) cin+start;
	89
	90	while (parsed<(len-start)) {
	91	//printf("%d ",parsed);
	92	if (*parser.uint8 != 0x80) {
	93	diff = (int) *parser.int8;
	94	parser.int8++;
	95	parsed += 1;
	96	}
	97	else {
	98	parser.uint8++;
	99	parsed += 1;
	100	if (*parser.uint16 != 0x8000) {
	101	diff = (int) *parser.int16;
	102	parser.int16++;
	103	parsed += 2;
	104	}
	105	else {
	106	parser.uint16++;
	107	parsed += 2;
	108	diff = (int) *parser.int32;
	109	parser.int32++;
	110	parsed += 4;
	111	}
	112	}
	113	cur += diff;
	114	*cout++ = (float) cur;
	115	np++;
	116	// check if we already have all data (file might be longer)
	117	if(np==nout) {
	118	//printf("all data read (%d,%d)\n",np,parsed);
	119	break;
	120	}
	121	}
	122
	123	// return output array
	124	return PyArray_Return(outarr);
	125	}
	126
	127	#undef PY_ARRAY_UNIQUE_SYMBOL

-2

xrayutilities/src/gridder.h less more

41	41	\brief python interface function
42	42
43	43	Python interface function for gridder2d. This function is virtually doing all
44		the Python related stuff to run gridder2d function.
	44	the Python related stuff to run gridder2d function.
45	45	\param self reference to the module
46	46	\param args function arguments
47	47	\return return value of the function

77	77	\brief 3D gridder python interface function
78	78
79	79	Python interface function for gridder3d. This function is virtually doing all
80		the Python related stuff to run gridder2d function.
	80	the Python related stuff to run gridder2d function.
81	81	\param self reference to the module
82	82	\param args function arguments
83	83	\return return value of the function

+13

-7

xrayutilities/src/gridder2d.c less more

24	24
25	25	#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
26	26	#define PY_ARRAY_UNIQUE_SYMBOL XU_UNIQUE_SYMBOL
	27	#define NO_IMPORT_ARRAY
27	28	#include "gridder.h"
28	29	#include "gridder_utils.h"
29	30

39	40	double xmin,xmax,ymin,ymax;
40	41	unsigned int nx,ny;
41	42	int flags;
	43	int n,result;
42	44
43	45	if(!PyArg_ParseTuple(args,"O!O!O!IIddddO!\|O!i",
44	46	&PyArray_Type,&py_x,

66	68	if(norm!=NULL) norm = (double *)PyArray_DATA(py_norm);
67	69
68	70	//get the total number of points
69		int n = PyArray_SIZE(py_x);
	71	n = PyArray_SIZE(py_x);
70	72
71	73	//call the actual gridder routine
72		int result = gridder2d(x,y,data,n,nx,ny,xmin,xmax,ymin,ymax,odata,norm,flags);
	74	result = gridder2d(x,y,data,n,nx,ny,xmin,xmax,ymin,ymax,odata,norm,flags);
73	75	return Py_BuildValue("i",&result);
74	76
75	77	}

85	87
86	88	double dx = delta(xmin,xmax,nx);
87	89	double dy = delta(ymin,ymax,ny);
88
	90
	91	unsigned int i; //loop index
	92
89	93	/check if normalization array is passed/
90	94	if(norm==NULL)
91	95	{

100	104	}
101	105	else
102	106	{
103		if(flags&VERBOSE)
	107	if(flags&VERBOSE)
104	108	{
105	109	fprintf(stdout,"XU.Gridder2D(c): use user provided buffer for normalization data\n");
106	110	}

108	112	}
109	113
110	114	/the master loop over all data points/
111		for(unsigned int i=0;i<n;i++)
	115	for(i=0;i<n;i++)
112	116	{
113	117	//if the x and y values are outside the grids boundaries continue with
114	118	//the next point

125	129	/perform normalization/
126	130	if(!(flags&NO_NORMALIZATION))
127	131	{
128		if(flags&VERBOSE)
	132	if(flags&VERBOSE)
129	133	fprintf(stdout,"XU.Gridder2D(c): perform normalization ...\n");
130	134
131		for(unsigned int i=0;i<nx*ny;i++)
	135	for(i=0;i<nx*ny;i++)
132	136	if(gnorm[i]>1.e-16) odata[i] = odata[i]/gnorm[i];
133	137
134	138	}

138	142
139	143	return(0);
140	144	}
	145
	146	#undef PY_ARRAY_UNIQUE_SYMBOL

-4

xrayutilities/src/gridder3d.c less more

24	24
25	25	#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
26	26	#define PY_ARRAY_UNIQUE_SYMBOL XU_UNIQUE_SYMBOL
	27	#define NO_IMPORT_ARRAY
27	28	#include "gridder.h"
28	29	#include "gridder_utils.h"
29	30	#include <numpy/arrayobject.h>

37	38	double xmin,xmax,ymin,ymax,zmin,zmax;
38	39	unsigned int nx,ny,nz;
39	40	int flags;
	41	int n,result;
40	42
41	43	if(!PyArg_ParseTuple(args,"O!O!O!O!IIIddddddO!\|O!i",
42	44	&PyArray_Type,&py_x,

68	70	if(norm!=NULL) norm = (double *)PyArray_DATA(py_norm);
69	71
70	72	//get the total number of points
71		int n = PyArray_SIZE(py_x);
	73	n = PyArray_SIZE(py_x);
72	74
73	75	//call the actual gridder routine
74		int result = gridder3d(x,y,z,data,n,nx,ny,nz,
	76	result = gridder3d(x,y,z,data,n,nx,ny,nz,
75	77	xmin,xmax,ymin,ymax,zmin,zmax,odata,norm,flags);
76	78	return Py_BuildValue("i",&result);
77	79

87	89	double *gnorm; //pointer to normalization data
88	90	unsigned int offset; //linear offset for the grid data
89	91	unsigned int ntot = nxnynz; //total number of points on the grid
	92	unsigned int i; //loop index variable
90	93
91	94
92	95	//compute step width for the grid

113	116	gnorm = norm;
114	117
115	118	/the master loop over all data points/
116		for(unsigned int i=0;i<n;i++)
	119	for(i=0;i<n;i++)
117	120	{
118	121	//check if the current point is within the bounds of the grid
119	122	if((x[i]<xmin)\|\|(x[i]>xmax)) continue;

132	135	/perform normalization/
133	136	if(!(flags&NO_NORMALIZATION))
134	137	{
135		for(unsigned int i=0;i<ntot;i++)
	138	for(i=0;i<ntot;i++)
136	139	if(gnorm[i]>1.e-16) odata[i] = odata[i]/gnorm[i];
137	140	}
138	141

141	144
142	145	return(0);
143	146	}
	147
	148	#undef PY_ARRAY_UNIQUE_SYMBOL

+24

-3

xrayutilities/src/gridder_utils.c less more

28	28	double get_min(double *a,unsigned int n)
29	29	{
30	30	double m = a[0];
	31	unsigned int i;
31	32
32		for(unsigned int i=0;i<n;i++) if(m<a[i]) m = a[i];
	33	for(i=0;i<n;i++) {
	34	if(m<a[i]) {
	35	m = a[i];
	36	}
	37	}
33	38
34	39	return(m);
35	40	}

38	43	double get_max(double *a,unsigned int n)
39	44	{
40	45	double m=a[0];
	46	unsigned int i;
41	47
42		for(unsigned int i=0;i<n;i++) if(m>a[i]) m = a[i];
	48	for(i=0;i<n;i++) {
	49	if(m>a[i]) {
	50	m = a[i];
	51	}
	52	}
43	53
44	54	return(m);
45	55	}

47	57	//-----------------------------------------------------------------------------
48	58	void set_array(double *a,unsigned int n,double value)
49	59	{
50		for(unsigned int i=0;i<n;++i) a[i] = value;
	60	unsigned int i;
	61
	62	for(i=0;i<n;++i) {
	63	a[i] = value;
	64	}
51	65	}
52	66
53	67	//-----------------------------------------------------------------------------

61	75	{
62	76	return (unsigned int)rint((x-min)/d);
63	77	}
	78	//-----------------------------------------------------------------------------
	79	#ifdef _WIN32
	80	double rint(double x)
	81	{
	82	return x < 0.0 ? ceil(x-0.5) : floor(x+0.5);
	83	}
	84	#endif⏎

-4

xrayutilities/src/gridder_utils.h less more

37	37	}
38	38
39	39	/*!
40		\brief find minimum
	40	\brief find minimum
41	41
42		Finds the minimum in an array.
	42	Finds the minimum in an array.
43	43	\param a input data
44	44	\param n number of elements
45	45	\return minimum value

48	48
49	49	//-----------------------------------------------------------------------------
50	50	/*!
51		\brief find maximum
	51	\brief find maximum
52	52
53	53	Finds the maximum value in an array.
54	54	\param a input data

72	72	/*!
73	73	\brief compute step width
74	74
75		Computes the stepwidth of a grid.
	75	Computes the stepwidth of a grid.
76	76	\param min minimum value
77	77	\param max maximum value
78	78	\param n number of steps

86	86
87	87	*/
88	88	unsigned int gindex(double x,double min,double d);
	89
	90	#ifdef _WIN32
	91	double rint(double x);
	92	#endif⏎

+194

-164

xrayutilities/src/qconversion.c less more

28	28
29	29	#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
30	30	#define PY_ARRAY_UNIQUE_SYMBOL XU_UNIQUE_SYMBOL
	31	#define NO_IMPORT_ARRAY
31	32	#include <numpy/arrayobject.h>
32	33
33	34	#include "qconversion.h"

49	50	#define OMPSETNUMTHREADS(nth) \
50	51	if(nth==0) omp_set_num_threads(omp_get_max_threads());\
51	52	else omp_set_num_threads(nth);
52
	53
53	54	/* ###################################
54	55	* matrix vector operations for
55	56	* 3x3 matrices and vectors of length

63	64	}
64	65
65	66	INLINE void sumvec(double RESTRICT v1,double RESTRICT v2) {
66		for(int i=0; i<3; ++i)
	67	unsigned int i;
	68	for(i=0; i<3; ++i)
67	69	v1[i] += v2[i];
68	70	}
69	71
70	72	INLINE void diffvec(double RESTRICT v1,double RESTRICT v2) {
71		for(int i=0; i<3; ++i)
	73	unsigned int i;
	74	for(i=0; i<3; ++i)
72	75	v1[i] -= v2[i];
73	76	}
74	77
75	78	INLINE double norm(double *v) {
76	79	double n=0.;
77		for(int i=0; i<3; ++i)
	80	unsigned int i;
	81	for(i=0; i<3; ++i)
78	82	n += v[i]*v[i];
79	83	return sqrt(n);
80	84	}
81	85
82	86	INLINE void normalize(double *v) {
83	87	double n=norm(v);
84		for(int i=0; i<3; ++i)
	88	unsigned int i;
	89	for(i=0; i<3; ++i)
85	90	v[i] /= n;
86	91	}
87	92
88	93	INLINE void veccopy(double RESTRICT v1, double RESTRICT v2) {
89		for(int i=0; i<3; ++i)
	94	unsigned int i;
	95	for(i=0; i<3; ++i)
90	96	v1[i] = v2[i];
91	97	}
92	98
93	99	INLINE void vecmul(double *RESTRICT r, double a) {
94		for(int i=0; i<3; ++i)
	100	unsigned int i;
	101	for(i=0; i<3; ++i)
95	102	r[i] *= a;
96	103	}
97	104

102	109	}
103	110
104	111	INLINE void vecmatcross(double RESTRICT v, double RESTRICT m, double *RESTRICT mr) {
105		for (int i=0; i<9; i=i+3) {
	112	unsigned int i;
	113	for (i=0; i<9; i=i+3) {
106	114	mr[0+i] = v[1]m[2+i] - v[2]m[1+i];
107	115	mr[1+i] = -v[0]m[2+i] + v[2]m[0+i];
108	116	mr[2+i] = v[0]m[1+i] - v[1]m[0+i];

110	118	}
111	119
112	120	INLINE void matmulc(double *RESTRICT m, double c) {
113		for (int i=0; i<9; i=i+1) {
	121	unsigned int i;
	122	for (i=0; i<9; i=i+1) {
114	123	m[i] *= c;
115	124	}
116	125	}

123	132
124	133	INLINE void matmul(double RESTRICT m1, double RESTRICT m2) {
125	134	double a,b,c;
126		for(int i=0; i<9; i=i+3) {
	135
	136	unsigned int i;
	137	for(i=0; i<9; i=i+3) {
127	138	a = m1[i]m2[0] + m1[i+1]m2[3] + m1[i+2]*m2[6];
128	139	b = m1[i]m2[1] + m1[i+1]m2[4] + m1[i+2]*m2[7];
129	140	c = m1[i]m2[2] + m1[i+1]m2[5] + m1[i+2]*m2[8];

134	145	}
135	146
136	147	INLINE void tensorprod(double RESTRICT v1, double RESTRICT v2, double *RESTRICT m) {
137		for(int i=0; i<3; i=i+1) {
138		for(int j=0; j<3; j=j+1) {
	148	unsigned int i,j;
	149	for(i=0; i<3; i=i+1) {
	150	for(j=0; j<3; j=j+1) {
139	151	m[i3+j] = v1[i]v2[j];
140	152	}
141	153	}
142	154	}
143	155
144	156	INLINE void summat(double RESTRICT m1,double RESTRICT m2) {
145		for(int i=0; i<9; ++i)
	157	unsigned int i;
	158	for(i=0; i<9; ++i)
146	159	m1[i] += m2[i];
147	160	}
148	161
149	162	INLINE void diffmat(double RESTRICT m1,double RESTRICT m2) {
150		for(int i=0; i<9; ++i)
	163	unsigned int i;
	164	for(i=0; i<9; ++i)
151	165	m1[i] -= m2[i];
152	166	}
153	167
154	168	INLINE void inversemat(double RESTRICT m, double RESTRICT i) {
155	169	double det;
156	170	double h1,h2,h3,h4,h5,h6;
	171	unsigned int j;
157	172
158	173	h1 = m[4]m[8]; // m11m22
159	174	h2 = m[5]m[6]; // m12m20

173	188	i[7] = (m[1]m[6] - m[0]m[7]);
174	189	i[8] = (m[0]m[4] - m[1]m[3]);
175	190
176		for(int j=0; j<9; ++j)
	191	for(j=0; j<9; ++j)
177	192	i[j] /= det;
178	193	}
179	194

279	294	* #######################################*/
280	295
281	296	int print_matrix(double *m) {
282		for(int i=0;i<9;i+=3) {
	297	unsigned int i;
	298	for(i=0;i<9;i+=3) {
283	299	printf("%8.5g %8.5g %8.5g\n",m[i],m[i+1],m[i+2]);
284	300	}
285	301	printf("\n");

304	320	* */
305	321
306	322	double tiltaxis[3], tiltmat[9];
307
308		for(int i=0; i<3; ++i)
	323	unsigned int i;
	324
	325	for(i=0; i<3; ++i)
309	326	rpixel[i] = 0.;
310	327
311	328	switch(tolower(dir[0])) {

372	389	/* feed the function pointer array with the correct
373	390	* rotation matrix generating functions
374	391	* */
375
376		for(int i=0; i<n; ++i) {
	392	unsigned int i;
	393
	394	for(i=0; i<n; ++i) {
377	395	switch(tolower(stringAxis[2*i])) {
378	396	case 'x':
379	397	switch(stringAxis[2*i+1]) {

437	455	*
438	456	* Parameters
439	457	* ----------
440		* sampleAngles .. angular positions of the sample goniometer (Npoints,Ns)
441		* detectorAngles. angular positions of the detector goniometer (Npoints,Nd)
442		* ri ............ direction of primary beam (length irrelevant) (angles zero)
443		* sampleAxis .... string with sample axis directions
444		* detectorAxis .. string with detector axis directions
445		* kappadir ...... rotation axis of a possible kappa circle
446		* UB ............ orientation matrix and reciprocal space conversion of investigated crystal (3,3)
447		* lambda ........ wavelength of the used x-rays (Angstreom)
	458	* sampleAngles .. angular positions of the sample goniometer (Npoints,Ns)
	459	* detectorAngles. angular positions of the detector goniometer (Npoints,Nd)
	460	* ri ............ direction of primary beam (length irrelevant) (angles zero)
	461	* sampleAxis .... string with sample axis directions
	462	* detectorAxis .. string with detector axis directions
	463	* kappadir ...... rotation axis of a possible kappa circle
	464	* UB ............ orientation matrix and reciprocal space conversion of investigated crystal (3,3)
	465	* lambda ........ wavelength of the used x-rays (Angstreom)
448	466	* nthreads ...... number of threads to use in parallel section of the code
449	467	*
450	468	* Returns
451	469	* -------
452		* qpos .......... momentum transfer (Npoints,3)
453		*
	470	* qpos .......... momentum transfer (Npoints,3)
	471	*
454	472	* */
455	473	{
456	474	double mtemp[9],mtemp2[9], ms[9], md[9]; //matrices

462	480	double lambda; // x-ray wavelength
463	481	char sampleAxis,detectorAxis; // string with sample and detector axis
464	482	double sampleAngles,detectorAngles, ri, kappadir, UB, qpos; // c-arrays for further usage
465
466		PyArrayObject sampleAnglesArr=NULL, detectorAnglesArr=NULL, riArr=NULL, kappadirArr=NULL,
	483	npy_intp nout[2];
	484	// arrays with function pointers to rotation matrix functions
	485	fp_rot *sampleRotation;
	486	fp_rot *detectorRotation;
	487
	488	PyArrayObject sampleAnglesArr=NULL, detectorAnglesArr=NULL, riArr=NULL, kappadirArr=NULL,
467	489	UBArr=NULL, qposArr=NULL; // numpy arrays
468	490
469	491	// Python argument conversion code
470		if (!PyArg_ParseTuple(args, "O!O!O!ssO!O!dI",&PyArray_Type, &sampleAnglesArr,
	492	if (!PyArg_ParseTuple(args, "O!O!O!ssO!O!dI",&PyArray_Type, &sampleAnglesArr,
471	493	&PyArray_Type, &detectorAnglesArr,
472	494	&PyArray_Type, &riArr,
473	495	&sampleAxis, &detectorAxis,
474	496	&PyArray_Type, &kappadirArr,
475	497	&PyArray_Type, &UBArr,
476		&lambda, &nthreads)) {
	498	&lambda, &nthreads)) {
477	499	return NULL;
478	500	}
479
	501
480	502	// check Python array dimensions and types
481	503	PYARRAY_CHECK(sampleAnglesArr,2,NPY_DOUBLE,"sampleAngles must be a 2D double array");
482	504	PYARRAY_CHECK(detectorAnglesArr,2,NPY_DOUBLE,"detectorAngles must be a 2D double array");
483	505	PYARRAY_CHECK(riArr,1,NPY_DOUBLE,"r_i must be a 1D double array");
484		if (PyArray_SIZE(riArr) != 3) { PyErr_SetString(PyExc_ValueError,"r_i needs to be of length 3");
	506	if (PyArray_SIZE(riArr) != 3) { PyErr_SetString(PyExc_ValueError,"r_i needs to be of length 3");
485	507	return NULL; }
486	508	PYARRAY_CHECK(kappadirArr,1,NPY_DOUBLE,"kappa_dir must be a 1D double array");
487		if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
	509	if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
488	510	return NULL; }
489	511	PYARRAY_CHECK(UBArr,2,NPY_DOUBLE,"UB must be a 2D double array");
490	512	if (PyArray_DIMS(UBArr)[0] != 3 \|\| PyArray_DIMS(UBArr)[1] != 3) {
491	513	PyErr_SetString(PyExc_ValueError,"UB must be of shape (3,3)"); return NULL; }
492
	514
493	515	Npoints = PyArray_DIMS(sampleAnglesArr)[0];
494	516	Ns = PyArray_DIMS(sampleAnglesArr)[1];
495	517	Nd = PyArray_DIMS(detectorAnglesArr)[1];

501	523	UB = (double *) PyArray_DATA(UBArr);
502	524
503	525	// create output ndarray
504		npy_intp nout[2];
505	526	nout[0] = Npoints;
506	527	nout[1] = 3;
507	528	qposArr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
508		qpos = (double *) PyArray_DATA(qposArr);
	529	qpos = (double *) PyArray_DATA(qposArr);
509	530
510	531	#ifdef __OPENMP__
511	532	//set openmp thread numbers dynamically

513	534	#endif
514	535
515	536	// arrays with function pointers to rotation matrix functions
516		fp_rot sampleRotation[Ns];
517		fp_rot detectorRotation[Nd];
	537	sampleRotation = (fp_rot) malloc(Nssizeof(fp_rot));
	538	detectorRotation = (fp_rot) malloc(Ndsizeof(fp_rot));
518	539
519	540	// determine axes directions
520	541	if(determine_axes_directions(sampleRotation,sampleAxis,Ns) != 0) { return NULL; }

572	593	*
573	594	* Parameters
574	595	* ----------
575		* sampleAngles .... angular positions of the goniometer (Npoints,Ns)
576		* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
577		* rcch ............ direction + distance of center channel (angles zero)
578		* sampleAxis ...... string with sample axis directions
579		* detectorAxis .... string with detector axis directions
580		* kappadir ........ rotation axis of a possible kappa circle
581		* cch ............. center channel of the detector
582		* dpixel .......... width of one pixel, same unit as distance rcch
583		* roi ............. region of interest of the detector
584		* dir ............. direction of the detector, e.g.: "x+"
585		* tilt ............ tilt of the detector direction from dir
586		* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
587		* lambda .......... wavelength of the used x-rays in Angstroem
	596	* sampleAngles .... angular positions of the goniometer (Npoints,Ns)
	597	* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
	598	* rcch ............ direction + distance of center channel (angles zero)
	599	* sampleAxis ...... string with sample axis directions
	600	* detectorAxis .... string with detector axis directions
	601	* kappadir ........ rotation axis of a possible kappa circle
	602	* cch ............. center channel of the detector
	603	* dpixel .......... width of one pixel, same unit as distance rcch
	604	* roi ............. region of interest of the detector
	605	* dir ............. direction of the detector, e.g.: "x+"
	606	* tilt ............ tilt of the detector direction from dir
	607	* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
	608	* lambda .......... wavelength of the used x-rays in Angstroem
588	609	* nthreads ........ number of threads to use in parallel section of the code
589		*
	610	*
590	611	* Returns
591	612	* -------
592		* qpos ............ momentum transfer (Npoints*Nch,3)
	613	* qpos ............ momentum transfer (Npoints*Nch,3)
593	614	* */
594	615	{
595	616	double mtemp[9],mtemp2[9], ms[9], md[9]; //matrices

604	625	char sampleAxis,detectorAxis,*dir; // string with sample and detector axis, and detector direction
605	626	double sampleAngles,detectorAngles, rcch, kappadir, UB, qpos; // c-arrays for further usage
606	627	int *roi; // region of interest integer array
607
	628	npy_intp nout[2];
	629	fp_rot *sampleRotation;
	630	fp_rot *detectorRotation;
	631
608	632	PyArrayObject sampleAnglesArr=NULL, detectorAnglesArr=NULL, *rcchArr=NULL,
609	633	kappadirArr=NULL, roiArr=NULL, UBArr=NULL, qposArr=NULL; // numpy arrays
610	634
611	635	// Python argument conversion code
612		if (!PyArg_ParseTuple(args, "O!O!O!ssO!ddO!sdO!dI",&PyArray_Type, &sampleAnglesArr,
	636	if (!PyArg_ParseTuple(args, "O!O!O!ssO!ddO!sdO!dI",&PyArray_Type, &sampleAnglesArr,
613	637	&PyArray_Type, &detectorAnglesArr,
614	638	&PyArray_Type, &rcchArr,
615	639	&sampleAxis, &detectorAxis,
616	640	&PyArray_Type, &kappadirArr,
617		&cch, &dpixel, &PyArray_Type, &roiArr,
	641	&cch, &dpixel, &PyArray_Type, &roiArr,
618	642	&dir, &tilt,
619	643	&PyArray_Type, &UBArr,
620		&lambda, &nthreads)) {
	644	&lambda, &nthreads)) {
621	645	return NULL;
622	646	}
623
	647
624	648	// check Python array dimensions and types
625	649	PYARRAY_CHECK(sampleAnglesArr,2,NPY_DOUBLE,"sampleAngles must be a 2D double array");
626	650	PYARRAY_CHECK(detectorAnglesArr,2,NPY_DOUBLE,"detectorAngles must be a 2D double array");
627	651	PYARRAY_CHECK(rcchArr,1,NPY_DOUBLE,"rcch must be a 1D double array");
628		if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
	652	if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
629	653	return NULL; }
630	654	PYARRAY_CHECK(kappadirArr,1,NPY_DOUBLE,"kappa_dir must be a 1D double array");
631		if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
	655	if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
632	656	return NULL; }
633	657	PYARRAY_CHECK(UBArr,2,NPY_DOUBLE,"UB must be a 2D double array");
634	658	if (PyArray_DIMS(UBArr)[0] != 3 \|\| PyArray_DIMS(UBArr)[1] != 3) {
635	659	PyErr_SetString(PyExc_ValueError,"UB must be of shape (3,3)"); return NULL; }
636	660	PYARRAY_CHECK(roiArr,1,NPY_INT32,"roi must be a 1D int array");
637		if (PyArray_SIZE(roiArr) != 2) {
	661	if (PyArray_SIZE(roiArr) != 2) {
638	662	PyErr_SetString(PyExc_ValueError,"roi must be of length 2"); return NULL; }
639
	663
640	664	Npoints = PyArray_DIMS(sampleAnglesArr)[0];
641	665	Ns = PyArray_DIMS(sampleAnglesArr)[1];
642	666	Nd = PyArray_DIMS(detectorAnglesArr)[1];

647	671	kappadir = (double *) PyArray_DATA(kappadirArr);
648	672	UB = (double *) PyArray_DATA(UBArr);
649	673	roi = (int *) PyArray_DATA(roiArr);
650
	674
651	675	// derived values from input parameters
652	676	Nch = roi[1]-roi[0]; // number of channels
653	677	f=M_2PI/lambda;
654
	678
655	679	// create output ndarray
656		npy_intp nout[2];
657	680	nout[0] = Npoints*Nch;
658	681	nout[1] = 3;
659	682	qposArr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
660		qpos = (double *) PyArray_DATA(qposArr);
	683	qpos = (double *) PyArray_DATA(qposArr);
661	684
662	685	#ifdef __OPENMP__
663	686	//set openmp thread numbers dynamically

665	688	#endif
666	689
667	690	// arrays with function pointers to rotation matrix functions
668		fp_rot sampleRotation[Ns];
669		fp_rot detectorRotation[Nd];
	691	sampleRotation = (fp_rot) malloc(Nssizeof(fp_rot));
	692	detectorRotation = (fp_rot) malloc(Ndsizeof(fp_rot));
670	693
671	694	// determine axes directions
672	695	if(determine_axes_directions(sampleRotation,sampleAxis,Ns) != 0) { return NULL; }

676	699	normalize(r_i);
677	700	// determine detector pixel vector
678	701	if(determine_detector_pixel(rpixel, dir, dpixel, r_i, tilt) != 0) { return NULL; }
679		for(int k=0; k<3; ++k)
	702	for(k=0; k<3; ++k)
680	703	rcchp[k] = rpixel[k]*cch;
681	704
682	705	// calculate rotation matices and perform rotations

732	755	*
733	756	* Parameters
734	757	* ----------
735		* sampleAngles .... angular positions of the sample goniometer (Npoints,Ns)
736		* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
737		* rcch ............ direction + distance of center pixel (angles zero)
738		* sampleAxis ...... string with sample axis directions
739		* detectorAxis .... string with detector axis directions
740		* kappadir ...... rotation axis of a possible kappa circle
741		* cch1 ............ center channel of the detector
742		* cch2 ............ center channel of the detector
743		* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
744		* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
	758	* sampleAngles .... angular positions of the sample goniometer (Npoints,Ns)
	759	* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
	760	* rcch ............ direction + distance of center pixel (angles zero)
	761	* sampleAxis ...... string with sample axis directions
	762	* detectorAxis .... string with detector axis directions
	763	* kappadir ...... rotation axis of a possible kappa circle
	764	* cch1 ............ center channel of the detector
	765	* cch2 ............ center channel of the detector
	766	* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
	767	* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
745	768	* roi ............. region of interest for the area detector [dir1min,dir1max,dir2min,dir2max]
746		* dir1 ............ first direction of the detector, e.g.: "x+"
747		* dir2 ............ second direction of the detector, e.g.: "z+"
748		* tiltazimuth ..... azimuth of the tilt
	769	* dir1 ............ first direction of the detector, e.g.: "x+"
	770	* dir2 ............ second direction of the detector, e.g.: "z+"
	771	* tiltazimuth ..... azimuth of the tilt
749	772	* tilt ............ tilt of the detector plane (rotation around axis normal to the direction
750		* given by the tiltazimuth
751		* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
752		* lambda .......... wavelength of the used x-rays
	773	* given by the tiltazimuth
	774	* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
	775	* lambda .......... wavelength of the used x-rays
753	776	* nthreads ........ number of threads to use in parallelization
754	777	*
755	778	* Returns
756	779	* -------
757		* qpos ............ momentum transfer (NpointsNpix1Npix2,3)
	780	* qpos ............ momentum transfer (NpointsNpix1Npix2,3)
758	781	* */
759	782	{
760	783	double mtemp[9],mtemp2[9], ms[9], md[9]; //matrices

769	792	char sampleAxis,detectorAxis,dir1,dir2; // string with sample and detector axis, and detector direction
770	793	double sampleAngles,detectorAngles, rcch, kappadir, UB, qpos; // c-arrays for further usage
771	794	int *roi; // region of interest integer array
772
	795	fp_rot *sampleRotation;
	796	fp_rot *detectorRotation;
	797	npy_intp nout[2];
	798
773	799	PyArrayObject sampleAnglesArr=NULL, detectorAnglesArr=NULL, *rcchArr=NULL,
774	800	kappadirArr=NULL, roiArr=NULL, UBArr=NULL, qposArr=NULL; // numpy arrays
775	801
776	802	// Python argument conversion code
777		if (!PyArg_ParseTuple(args, "O!O!O!ssO!ddddO!ssddO!dI",&PyArray_Type, &sampleAnglesArr,
	803	if (!PyArg_ParseTuple(args, "O!O!O!ssO!ddddO!ssddO!dI",&PyArray_Type, &sampleAnglesArr,
778	804	&PyArray_Type, &detectorAnglesArr,
779	805	&PyArray_Type, &rcchArr,
780	806	&sampleAxis, &detectorAxis,
781	807	&PyArray_Type, &kappadirArr,
782		&cch1, &cch2, &dpixel1, &dpixel2,
783		&PyArray_Type, &roiArr,
	808	&cch1, &cch2, &dpixel1, &dpixel2,
	809	&PyArray_Type, &roiArr,
784	810	&dir1, &dir2, &tiltazimuth, &tilt,
785	811	&PyArray_Type, &UBArr,
786		&lambda, &nthreads)) {
	812	&lambda, &nthreads)) {
787	813	return NULL;
788	814	}
789
	815
790	816	// check Python array dimensions and types
791	817	PYARRAY_CHECK(sampleAnglesArr,2,NPY_DOUBLE,"sampleAngles must be a 2D double array");
792	818	PYARRAY_CHECK(detectorAnglesArr,2,NPY_DOUBLE,"detectorAngles must be a 2D double array");
793	819	PYARRAY_CHECK(rcchArr,1,NPY_DOUBLE,"rcch must be a 1D double array");
794		if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
	820	if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
795	821	return NULL; }
796	822	PYARRAY_CHECK(kappadirArr,1,NPY_DOUBLE,"kappa_dir must be a 1D double array");
797		if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
	823	if (PyArray_SIZE(kappadirArr) != 3) { PyErr_SetString(PyExc_ValueError,"kappa_dir needs to be of length 3");
798	824	return NULL; }
799	825	PYARRAY_CHECK(UBArr,2,NPY_DOUBLE,"UB must be a 2D double array");
800	826	if (PyArray_DIMS(UBArr)[0] != 3 \|\| PyArray_DIMS(UBArr)[1] != 3) {
801	827	PyErr_SetString(PyExc_ValueError,"UB must be of shape (3,3)"); return NULL; }
802	828	PYARRAY_CHECK(roiArr,1,NPY_INT32,"roi must be a 1D int array");
803		if (PyArray_SIZE(roiArr) != 4) {
	829	if (PyArray_SIZE(roiArr) != 4) {
804	830	PyErr_SetString(PyExc_ValueError,"roi must be of length 4"); return NULL; }
805
	831
806	832	Npoints = PyArray_DIMS(sampleAnglesArr)[0];
807	833	Ns = PyArray_DIMS(sampleAnglesArr)[1];
808	834	Nd = PyArray_DIMS(detectorAnglesArr)[1];

813	839	kappadir = (double *) PyArray_DATA(kappadirArr);
814	840	UB = (double *) PyArray_DATA(UBArr);
815	841	roi = (int *) PyArray_DATA(roiArr);
816
	842
817	843	// derived values from input parameters
818	844	f=M_2PI/lambda;
819	845	// calculate some index shortcuts
820	846	idxh1 = (roi[1]-roi[0])*(roi[3]-roi[2]);
821	847	idxh2 = roi[3]-roi[2];
822
	848
823	849	// create output ndarray
824		npy_intp nout[2];
825	850	nout[0] = Npoints*idxh1;
826	851	nout[1] = 3;
827	852	qposArr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
828		qpos = (double *) PyArray_DATA(qposArr);
	853	qpos = (double *) PyArray_DATA(qposArr);
829	854
830	855	#ifdef __OPENMP__
831	856	//set openmp thread numbers dynamically

833	858	#endif
834	859
835	860	// arrays with function pointers to rotation matrix functions
836		fp_rot sampleRotation[Ns];
837		fp_rot detectorRotation[Nd];
	861	sampleRotation = (fp_rot) malloc(Nssizeof(fp_rot));
	862	detectorRotation = (fp_rot) malloc(Ndsizeof(fp_rot));
838	863
839	864	// determine axes directions
840	865	if(determine_axes_directions(sampleRotation,sampleAxis,Ns) != 0) { return NULL; }

872	897	print_vector(rpixel2);*/
873	898
874	899	// calculate center channel position in detector plane
875		for(int k=0; k<3; ++k)
	900	for(k=0; k<3; ++k)
876	901	rcchp[k] = rpixel1[k]cch1 + rpixel2[k]cch2;
877	902
878	903	// calculate rotation matices and perform rotations

935	960	*
936	961	* Parameters
937	962	* ----------
938		* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
	963	* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
939	964	* n1 .............. detector pixel numbers dim1 (Npoints)
940		* n2 .............. detector pixel numbers dim2 (Npoints)
941		* rcch ............ direction + distance of center pixel (angles zero)
942		* detectorAxis .... string with detector axis directions
	965	* n2 .............. detector pixel numbers dim2 (Npoints)
	966	* rcch ............ direction + distance of center pixel (angles zero)
	967	* detectorAxis .... string with detector axis directions
943	968	* cch1 ............ center channel of the detector
944		* cch2 ............ center channel of the detector
945		* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
946		* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
947		* dir1 ............ first direction of the detector, e.g.: "x+"
948		* dir2 ............ second direction of the detector, e.g.: "z+"
949		* tiltazimuth ..... azimuth of the tilt
	969	* cch2 ............ center channel of the detector
	970	* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
	971	* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
	972	* dir1 ............ first direction of the detector, e.g.: "x+"
	973	* dir2 ............ second direction of the detector, e.g.: "z+"
	974	* tiltazimuth ..... azimuth of the tilt
950	975	* tilt ............ tilt of the detector plane (rotation around axis normal to the direction
951		* given by the tiltazimuth
952		* lambda .......... wavelength of the used x-rays
	976	* given by the tiltazimuth
	977	* lambda .......... wavelength of the used x-rays
953	978	* nthreads ........ number of threads to use in parallel section of the code
954	979	*
955	980	* Returns
956	981	* -------
957		* qpos ............ momentum transfer (Npoints,3)
	982	* qpos ............ momentum transfer (Npoints,3)
958	983	* */
959	984	{
960	985	double mtemp[9], md[9]; //matrices

967	992	double f,lambda,cch1,cch2,dpixel1,dpixel2,tilt,tiltazimuth; // x-ray wavelength, f=M_2PI/lambda and detector parameters
968	993	char detectorAxis,dir1,*dir2; // string with detector axis, and detector direction
969	994	double detectorAngles, n1, n2, rcch, *qpos; // c-arrays for further usage
970
	995	fp_rot *detectorRotation;
	996	npy_intp nout[2];
	997
971	998	PyArrayObject detectorAnglesArr=NULL, n1Arr=NULL, n2Arr=NULL, rcchArr=NULL, *qposArr=NULL; // numpy arrays
972	999
973	1000	// Python argument conversion code
974	1001	if (!PyArg_ParseTuple(args, "O!O!O!O!sddddssdddI", &PyArray_Type, &detectorAnglesArr,
975	1002	&PyArray_Type, &n1Arr, &PyArray_Type, &n2Arr, &PyArray_Type, &rcchArr,
976		&detectorAxis, &cch1, &cch2, &dpixel1, &dpixel2,
	1003	&detectorAxis, &cch1, &cch2, &dpixel1, &dpixel2,
977	1004	&dir1, &dir2, &tiltazimuth, &tilt,
978		&lambda, &nthreads)) {
	1005	&lambda, &nthreads)) {
979	1006	return NULL;
980	1007	}
981
	1008
982	1009	// check Python array dimensions and types
983	1010	PYARRAY_CHECK(detectorAnglesArr,2,NPY_DOUBLE,"detectorAngles must be a 2D double array");
984	1011	PYARRAY_CHECK(rcchArr,1,NPY_DOUBLE,"rcch must be a 1D double array");
985		if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
	1012	if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
986	1013	return NULL; }
987	1014	PYARRAY_CHECK(n1Arr,1,NPY_DOUBLE,"n1 must be a 1D double array");
988	1015	PYARRAY_CHECK(n2Arr,1,NPY_DOUBLE,"n2 must be a 1D double array");
989
	1016
990	1017	Npoints = PyArray_DIMS(detectorAnglesArr)[0];
991	1018	if (PyArray_SIZE(n1Arr) != Npoints \|\| PyArray_SIZE(n2Arr) != Npoints) {
992	1019	PyErr_SetString(PyExc_ValueError,"n1,n2 must be of length Npoints");

995	1022
996	1023	detectorAngles = (double *) PyArray_DATA(detectorAnglesArr);
997	1024	rcch = (double *) PyArray_DATA(rcchArr);
998		n1 = (double *) PyArray_DATA(n1Arr);
999		n2 = (double *) PyArray_DATA(n2Arr);
1000
	1025	n1 = (double *) PyArray_DATA(n1Arr);
	1026	n2 = (double *) PyArray_DATA(n2Arr);
	1027
1001	1028	// derived values from input parameters
1002	1029	f=M_2PI/lambda;
1003
	1030
1004	1031	// create output ndarray
1005		npy_intp nout[2];
1006	1032	nout[0] = Npoints;
1007	1033	nout[1] = 3;
1008	1034	qposArr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
1009		qpos = (double *) PyArray_DATA(qposArr);
	1035	qpos = (double *) PyArray_DATA(qposArr);
1010	1036
1011	1037	#ifdef __OPENMP__
1012	1038	//set openmp thread numbers dynamically

1014	1040	#endif
1015	1041
1016	1042	// arrays with function pointers to rotation matrix functions
1017		fp_rot detectorRotation[Nd];
	1043	detectorRotation = (fp_rot) malloc(Ndsizeof(fp_rot));
1018	1044
1019	1045	// determine axes directions
1020	1046	if(determine_axes_directions(detectorRotation,detectorAxis,Nd) != 0) { return NULL; }

1046	1072	matvec(mtemp,rtemp,rpixel2);
1047	1073
1048	1074	// calculate center channel position in detector plane
1049		for(int k=0; k<3; ++k)
	1075	for(k=0; k<3; ++k)
1050	1076	rcchp[k] = rpixel1[k]cch1 + rpixel2[k]cch2;
1051	1077
1052	1078	// calculate rotation matices and perform rotations

1089	1115	* as input to allow for a simultaneous fit of reference samples orientation
1090	1116	*
1091	1117	* Interface:
1092		* sampleAngles .... angular positions of the sample goniometer (Npoints,Ns)
1093		* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
1094		* n1 .............. detector pixel numbers dim1 (Npoints)
1095		* n2 .............. detector pixel numbers dim2 (Npoints)
1096		* rcch ............ direction + distance of center pixel (angles zero)
1097		* sampleAxis ...... string with sample axis directions
1098		* detectorAxis .... string with detector axis directions
1099		* cch1 ............ center channel of the detector
1100		* cch2 ............ center channel of the detector
1101		* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
1102		* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
1103		* dir1 ............ first direction of the detector, e.g.: "x+"
1104		* dir2 ............ second direction of the detector, e.g.: "z+"
1105		* tiltazimuth ..... azimuth of the tilt
	1118	* sampleAngles .... angular positions of the sample goniometer (Npoints,Ns)
	1119	* detectorAngles .. angular positions of the detector goniometer (Npoints,Nd)
	1120	* n1 .............. detector pixel numbers dim1 (Npoints)
	1121	* n2 .............. detector pixel numbers dim2 (Npoints)
	1122	* rcch ............ direction + distance of center pixel (angles zero)
	1123	* sampleAxis ...... string with sample axis directions
	1124	* detectorAxis .... string with detector axis directions
	1125	* cch1 ............ center channel of the detector
	1126	* cch2 ............ center channel of the detector
	1127	* dpixel1 ......... width of one pixel in first direction, same unit as distance rcch
	1128	* dpixel2 ......... width of one pixel in second direction, same unit as distance rcch
	1129	* dir1 ............ first direction of the detector, e.g.: "x+"
	1130	* dir2 ............ second direction of the detector, e.g.: "z+"
	1131	* tiltazimuth ..... azimuth of the tilt
1106	1132	* tilt ............ tilt of the detector plane (rotation around axis normal to the direction
1107		* given by the tiltazimuth
1108		* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
1109		* lambda .......... wavelength of the used x-rays
	1133	* given by the tiltazimuth
	1134	* UB .............. orientation matrix and reciprocal space conversion of investigated crystal (3,3)
	1135	* lambda .......... wavelength of the used x-rays
1110	1136	* nthreads ........ number of threads to use in parallel section of the code
1111	1137	*
1112	1138	* Returns
1113	1139	* -------
1114		* qpos ............ momentum transfer (Npoints,3)
	1140	* qpos ............ momentum transfer (Npoints,3)
1115	1141	* */
1116	1142	{
1117	1143	double mtemp[9],mtemp2[9], ms[9], md[9]; //matrices

1124	1150	double f,lambda,cch1,cch2,dpixel1,dpixel2,tilt,tiltazimuth; // x-ray wavelength, f=M_2PI/lambda and detector parameters
1125	1151	char sampleAxis,detectorAxis,dir1,dir2; // string with sample and detector axis, and detector direction
1126	1152	double sampleAngles, detectorAngles, n1, n2, rcch, UB, *qpos; // c-arrays for further usage
	1153	fp_rot *sampleRotation;
	1154	fp_rot *detectorRotation;
	1155	npy_intp nout[2];
1127	1156
1128	1157	PyArrayObject sampleAnglesArr=NULL, detectorAnglesArr=NULL, n1Arr=NULL, n2Arr=NULL, rcchArr=NULL, UBArr=NULL, *qposArr=NULL; // numpy arrays
1129	1158
1130	1159	// Python argument conversion code
1131	1160	if (!PyArg_ParseTuple(args, "O!O!O!O!O!ssddddssddO!dI", &PyArray_Type, &sampleAnglesArr, &PyArray_Type, &detectorAnglesArr,
1132	1161	&PyArray_Type, &n1Arr, &PyArray_Type, &n2Arr, &PyArray_Type, &rcchArr,
1133		&sampleAxis, &detectorAxis, &cch1, &cch2, &dpixel1, &dpixel2,
	1162	&sampleAxis, &detectorAxis, &cch1, &cch2, &dpixel1, &dpixel2,
1134	1163	&dir1, &dir2, &tiltazimuth, &tilt, &PyArray_Type, &UBArr,
1135		&lambda, &nthreads)) {
	1164	&lambda, &nthreads)) {
1136	1165	return NULL;
1137	1166	}
1138	1167

1140	1169	PYARRAY_CHECK(sampleAnglesArr,2,NPY_DOUBLE,"sampleAngles must be a 2D double array");
1141	1170	PYARRAY_CHECK(detectorAnglesArr,2,NPY_DOUBLE,"detectorAngles must be a 2D double array");
1142	1171	PYARRAY_CHECK(rcchArr,1,NPY_DOUBLE,"rcch must be a 1D double array");
1143		if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
	1172	if (PyArray_SIZE(rcchArr) != 3) { PyErr_SetString(PyExc_ValueError,"rcch needs to be of length 3");
1144	1173	return NULL; }
1145	1174	PYARRAY_CHECK(UBArr,2,NPY_DOUBLE,"UB must be a 2D double array");
1146	1175	if (PyArray_DIMS(UBArr)[0] != 3 \|\| PyArray_DIMS(UBArr)[1] != 3) {
1147	1176	PyErr_SetString(PyExc_ValueError,"UB must be of shape (3,3)"); return NULL; }
1148	1177	PYARRAY_CHECK(n1Arr,1,NPY_DOUBLE,"n1 must be a 1D double array");
1149	1178	PYARRAY_CHECK(n2Arr,1,NPY_DOUBLE,"n2 must be a 1D double array");
1150
	1179
1151	1180	Npoints = PyArray_DIMS(detectorAnglesArr)[0];
1152	1181	if (PyArray_SIZE(n1Arr) != Npoints \|\| PyArray_SIZE(n2Arr) != Npoints) {
1153	1182	PyErr_SetString(PyExc_ValueError,"n1,n2 must be of length Npoints");

1159	1188	sampleAngles = (double *) PyArray_DATA(sampleAnglesArr);
1160	1189	rcch = (double *) PyArray_DATA(rcchArr);
1161	1190	UB = (double *) PyArray_DATA(UBArr);
1162		n1 = (double *) PyArray_DATA(n1Arr);
1163		n2 = (double *) PyArray_DATA(n2Arr);
1164
	1191	n1 = (double *) PyArray_DATA(n1Arr);
	1192	n2 = (double *) PyArray_DATA(n2Arr);
	1193
1165	1194	// derived values from input parameters
1166	1195	f=M_2PI/lambda;
1167	1196
1168	1197	// create output ndarray
1169		npy_intp nout[2];
1170	1198	nout[0] = Npoints;
1171	1199	nout[1] = 3;
1172	1200	qposArr = (PyArrayObject *) PyArray_SimpleNew(2, nout, NPY_DOUBLE);
1173		qpos = (double *) PyArray_DATA(qposArr);
	1201	qpos = (double *) PyArray_DATA(qposArr);
1174	1202
1175	1203	#ifdef __OPENMP__
1176	1204	//set openmp thread numbers dynamically

1178	1206	#endif
1179	1207
1180	1208	// arrays with function pointers to rotation matrix functions
1181		fp_rot sampleRotation[Ns];
1182		fp_rot detectorRotation[Nd];
	1209	sampleRotation = (fp_rot) malloc(Nssizeof(fp_rot));
	1210	detectorRotation = (fp_rot) malloc(Ndsizeof(fp_rot));
	1211
1183	1212
1184	1213	// determine axes directions
1185	1214	if(determine_axes_directions(sampleRotation,sampleAxis,Ns) != 0) { return NULL; }

1212	1241	matvec(mtemp,rtemp,rpixel2);
1213	1242
1214	1243	// calculate center channel position in detector plane
1215		for(int k=0; k<3; ++k)
	1244	for(k=0; k<3; ++k)
1216	1245	rcchp[k] = rpixel1[k]cch1 + rpixel2[k]cch2;
1217	1246
1218	1247	// calculate rotation matices and perform rotations

1257	1286	return PyArray_Return(qposArr);
1258	1287	}
1259	1288
	1289	#undef PY_ARRAY_UNIQUE_SYMBOL

-3

xrayutilities/src/qconversion.h less more

18	18
19	19	#pragma once
20	20
21		#define M_PI 3.14159265358979323846
	21	#ifndef M_PI
	22	# define M_PI 3.14159265358979323846
	23	#endif
22	24	#define M_2PI (2*M_PI)
23	25
24	26	#define cdeg2rad (M_PI/180.)

27	29	#define deg2rad(ang) (ang*cdeg2rad)
28	30	#define rad2deg(rad) (rad*crad2deg)
29	31
30		/* 'extern inline' seems to work only on newer version of gcc (>4.6 tested)
	32	/* 'extern inline' seems to work only on newer version of gcc (>4.6 tested)
31	33	* gcc 4.1 seems to need this empty, i am not sure if there is a speed gain
32	34	* by inlining since the calls to those functions are anyhow built dynamically
33	35	* for compatibility keep this empty unless you can test with several compilers */
34		#define INLINE
	36	#define INLINE
	37	#ifdef _WIN32
	38	#define RESTRICT
	39	#else
35	40	#define RESTRICT restrict
	41	#endif
36	42
37	43	/* ###################################
38	44	* matrix vector operations for

xrayutilities.pdf less more

Binary diff not shown