Context Navigation

source: trunk/GSASIIlattice.py @ 2480

Last change on this file since 2480 was 2480, checked in by vondreele, 7 years ago
add make magnetic phase to General/Transform? option trap missing rigid bodies to define torsion seq option fix lighting issues for polygons fix xye importer to stop on trailing blank lines rather than crashing fix error in powder structure factor calc add magnetic structure factor calc. (some error still & no derivatives yet)
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 93.4 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIlattice: Unit cells
4	---------------------------
5
6	Perform lattice-related computations
7
8	Note that g is the reciprocal lattice tensor, and G is its inverse,
9	:math:`G = g^{-1}`, where
10
11	.. math::
12
13	G = \\left( \\begin{matrix}
14	a^2 & a b\\cos\gamma & a c\\cos\\beta \\\\
15	a b\\cos\\gamma & b^2 & b c \cos\\alpha \\\\
16	a c\\cos\\beta & b c \\cos\\alpha & c^2
17	\\end{matrix}\\right)
18
19	The "A tensor" terms are defined as
20	:math:`A = (\\begin{matrix} G_{11} & G_{22} & G_{33} & 2G_{12} & 2G_{13} & 2G_{23}\\end{matrix})` and A can be used in this fashion:
21	:math:`d^* = \sqrt {A_1 h^2 + A_2 k^2 + A_3 l^2 + A_4 hk + A_5 hl + A_6 kl}`, where
22	d is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing,
23	:math:`Q = 2 \\pi d^* = 2 \\pi / d`
24	'''
25	########### SVN repository information ###################
26	# $Date: 2016-09-28 17:16:00 +0000 (Wed, 28 Sep 2016) $
27	# $Author: vondreele $
28	# $Revision: 2480 $
29	# $URL: trunk/GSASIIlattice.py $
30	# $Id: GSASIIlattice.py 2480 2016-09-28 17:16:00Z vondreele $
31	########### SVN repository information ###################
32	import math
33	import copy
34	import sys
35	import random as ran
36	import numpy as np
37	import numpy.linalg as nl
38	import GSASIIpath
39	import GSASIImath as G2mth
40	import GSASIIspc as G2spc
41	GSASIIpath.SetVersionNumber("$Revision: 2480 $")
42	# trig functions in degrees
43	sind = lambda x: np.sin(x*np.pi/180.)
44	asind = lambda x: 180.*np.arcsin(x)/np.pi
45	tand = lambda x: np.tan(x*np.pi/180.)
46	atand = lambda x: 180.*np.arctan(x)/np.pi
47	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
48	cosd = lambda x: np.cos(x*np.pi/180.)
49	acosd = lambda x: 180.*np.arccos(x)/np.pi
50	rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
51	rpd = np.pi/180.
52	RSQ2PI = 1./np.sqrt(2.*np.pi)
53	SQ2 = np.sqrt(2.)
54	RSQPI = 1./np.sqrt(np.pi)
55	R2pisq = 1./(2.np.pi*2)
56	nxs = np.newaxis
57
58	def sec2HMS(sec):
59	"""Convert time in sec to H:M:S string
60
61	:param sec: time in seconds
62	:return: H:M:S string (to nearest 100th second)
63
64	"""
65	H = int(sec/3600)
66	M = int(sec/60-H*60)
67	S = sec-3600H-60M
68	return '%d:%2d:%.2f'%(H,M,S)
69
70	def rotdMat(angle,axis=0):
71	"""Prepare rotation matrix for angle in degrees about axis(=0,1,2)
72
73	:param angle: angle in degrees
74	:param axis: axis (0,1,2 = x,y,z) about which for the rotation
75	:return: rotation matrix - 3x3 numpy array
76
77	"""
78	if axis == 2:
79	return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
80	elif axis == 1:
81	return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
82	else:
83	return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
84
85	def rotdMat4(angle,axis=0):
86	"""Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
87
88	:param angle: angle in degrees
89	:param axis: axis (0,1,2 = x,y,z) about which for the rotation
90	:return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
91
92	"""
93	Mat = rotdMat(angle,axis)
94	return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
95
96	def fillgmat(cell):
97	"""Compute lattice metric tensor from unit cell constants
98
99	:param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
100	:return: 3x3 numpy array
101
102	"""
103	a,b,c,alp,bet,gam = cell
104	g = np.array([
105	[aa, abcosd(gam), ac*cosd(bet)],
106	[abcosd(gam), bb, bc*cosd(alp)],
107	[accosd(bet) ,bccosd(alp), c*c]])
108	return g
109
110	def cell2Gmat(cell):
111	"""Compute real and reciprocal lattice metric tensor from unit cell constants
112
113	:param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
114	:return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
115
116	"""
117	g = fillgmat(cell)
118	G = nl.inv(g)
119	return G,g
120
121	def A2Gmat(A,inverse=True):
122	"""Fill real & reciprocal metric tensor (G) from A.
123
124	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
125	:param bool inverse: if True return both G and g; else just G
126	:return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
127
128	"""
129	G = np.zeros(shape=(3,3))
130	G = [
131	[A[0], A[3]/2., A[4]/2.],
132	[A[3]/2.,A[1], A[5]/2.],
133	[A[4]/2.,A[5]/2., A[2]]]
134	if inverse:
135	g = nl.inv(G)
136	return G,g
137	else:
138	return G
139
140	def Gmat2A(G):
141	"""Extract A from reciprocal metric tensor (G)
142
143	:param G: reciprocal maetric tensor (3x3 numpy array
144	:return: A = [G11,G22,G33,2G12,2G13,2*G23]
145
146	"""
147	return [G[0][0],G[1][1],G[2][2],2.G[0][1],2.G[0][2],2.*G[1][2]]
148
149	def cell2A(cell):
150	"""Obtain A = [G11,G22,G33,2G12,2G13,2*G23] from lattice parameters
151
152	:param cell: [a,b,c,alpha,beta,gamma] (degrees)
153	:return: G reciprocal metric tensor as 3x3 numpy array
154
155	"""
156	G,g = cell2Gmat(cell)
157	return Gmat2A(G)
158
159	def A2cell(A):
160	"""Compute unit cell constants from A
161
162	:param A: [G11,G22,G33,2G12,2G13,2*G23] G - reciprocal metric tensor
163	:return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
164
165	"""
166	G,g = A2Gmat(A)
167	return Gmat2cell(g)
168
169	def Gmat2cell(g):
170	"""Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
171	The math works the same either way.
172
173	:param g (or G): real (or reciprocal) metric tensor 3x3 array
174	:return: a,b,c,alpha, beta, gamma (degrees) (or a,b,c,alpha,beta,gamma degrees)
175
176	"""
177	oldset = np.seterr('raise')
178	a = np.sqrt(max(0,g[0][0]))
179	b = np.sqrt(max(0,g[1][1]))
180	c = np.sqrt(max(0,g[2][2]))
181	alp = acosd(g[2][1]/(b*c))
182	bet = acosd(g[2][0]/(a*c))
183	gam = acosd(g[0][1]/(a*b))
184	np.seterr(**oldset)
185	return a,b,c,alp,bet,gam
186
187	def invcell2Gmat(invcell):
188	"""Compute real and reciprocal lattice metric tensor from reciprocal
189	unit cell constants
190
191	:param invcell: [a,b,c,alpha, beta, gamma] (degrees)
192	:return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
193
194	"""
195	G = fillgmat(invcell)
196	g = nl.inv(G)
197	return G,g
198
199	def prodMGMT(G,Mat):
200	'''Transform metric tensor by matrix
201
202	:param G: array metric tensor
203	:param Mat: array transformation matrix
204	:return: array new metric tensor
205
206	'''
207	return np.inner(Mat,np.inner(G,Mat).T)
208
209	def TransformCell(cell,Trans):
210	'''Transform lattice parameters by matrix
211
212	:param cell: list a,b,c,alpha,beta,gamma,(volume)
213	:param Trans: array transformation matrix
214	:return: array transformed a,b,c,alpha,beta,gamma,volume
215
216	'''
217	newCell = np.zeros(7)
218	g = cell2Gmat(cell)[1]
219	newg = prodMGMT(g,Trans)
220	newCell[:6] = Gmat2cell(newg)
221	newCell[6] = calc_V(cell2A(newCell[:6]))
222	return newCell
223
224	def TransformXYZ(XYZ,Trans,Vec):
225	return np.inner(XYZ,Trans)+Vec
226
227	def TransformU6(U6,Trans):
228	Uij = np.inner(Trans,np.inner(U6toUij(U6),Trans))
229	return UijtoU6(Uij)
230
231	def TransformPhase(oldPhase,newPhase,Trans,Vec,ifMag):
232	'''Transform atoms from oldPhase to newPhase by Trans & Vec
233	NB: doesnt transform moments correctly - TBD
234
235	:param oldPhase: dict G2 phase info for old phase
236	:param newPhase: dict G2 phase info for new phase; with new cell & space group
237	atoms are from oldPhase & will be transformed
238	:param Trans: array transformation matrix
239	:param Vec: array transformation vector
240	:param ifMag: bool True if convert to magnetic phase
241	'''
242
243	cx,ct,cs,cia = oldPhase['General']['AtomPtrs']
244	SGData = newPhase['General']['SGData']
245	invTrans = nl.inv(Trans)
246	newAtoms = FillUnitCell(oldPhase)
247	Unit =[abs(int(max(unit))-1) for unit in Trans]
248	for i,unit in enumerate(Unit):
249	if unit > 0:
250	for j in range(unit):
251	moreAtoms = copy.deepcopy(newAtoms)
252	for atom in moreAtoms:
253	atom[cx+i] += 1.
254	newAtoms += moreAtoms
255	if ifMag:
256	cia += 3
257	cs += 3
258	newPhase['General']['Type'] = 'magnetic'
259	newPhase['General']['AtomPtrs'] = [cx,ct,cs,cia]
260	magAtoms = []
261	for atom in newAtoms:
262	magAtoms.append(atom[:cx+4]+[0.,0.,0.]+atom[cx+4:])
263	newAtoms = magAtoms
264	newPhase['Draw Atoms'] = []
265	for atom in newAtoms:
266	atom[cx:cx+3] = TransformXYZ(atom[cx:cx+3],invTrans.T,Vec)%1.
267	if atom[cia] == 'A':
268	atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],invTrans)
269	atom[cs:cs+2] = G2spc.SytSym(atom[cx:cx+3],SGData)[:2]
270	atom[cia+8] = ran.randint(0,sys.maxint)
271	newPhase['Atoms'] = newAtoms
272	newPhase['Atoms'] = GetUnique(newPhase)
273	newPhase['Drawing']['Atoms'] = []
274	return newPhase
275
276	def FillUnitCell(Phase):
277	Atoms = Phase['Atoms']
278	atomData = []
279	SGData = Phase['General']['SGData']
280	cx,ct,cs,cia = Phase['General']['AtomPtrs']
281	unit = np.zeros(3)
282	for atom in Atoms:
283	XYZ = np.array(atom[cx:cx+3])
284	xyz = XYZ%1.
285	unit = XYZ-xyz
286	if atom[cia] == 'A':
287	Uij = atom[cia+2:cia+8]
288	result = G2spc.GenAtom(xyz,SGData,False,Uij,True)
289	for item in result:
290	if item[0][2] >= .95: item[0][2] -= 1.
291	atom[cx:cx+3] = item[0]
292	atom[cia+2:cia+8] = item[1]
293	atomData.append(atom[:cia+9]) #not SS stuff
294	else:
295	result = G2spc.GenAtom(xyz,SGData,False,Move=True)
296	for item in result:
297	if item[0][2] >= .95: item[0][2] -= 1.
298	atom[cx:cx+3] = item[0]
299	atomData.append(atom[:cia+9]) #not SS stuff
300	return atomData
301
302	def GetUnique(Phase):
303
304	def noDuplicate(xyzA,XYZ,Amat):
305	if True in [np.allclose(np.inner(Amat,xyzA),np.inner(Amat,xyzB),atol=0.05) for xyzB in XYZ]:
306	return False
307	return True
308
309	cx,ct,cs,cia = Phase['General']['AtomPtrs']
310	cell = Phase['General']['Cell'][1:7]
311	Amat,Bmat = cell2AB(cell)
312	SGData = Phase['General']['SGData']
313	Atoms = Phase['Atoms']
314	Ind = len(Atoms)
315	newAtoms = []
316	Indx = {}
317	XYZ = {}
318	for ind in range(Ind):
319	XYZ[ind] = np.array(Atoms[ind][cx:cx+3])%1.
320	Indx[ind] = True
321	for ind in range(Ind):
322	if Indx[ind]:
323	xyz = XYZ[ind]
324	for jnd in range(Ind):
325	if ind != jnd and Indx[jnd]:
326	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
327	xyzs = np.array([equiv[0] for equiv in Equiv])
328	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
329	Ind = []
330	for ind in Indx:
331	if Indx[ind]:
332	newAtoms.append(Atoms[ind])
333	return newAtoms
334
335	def calc_rVsq(A):
336	"""Compute the square of the reciprocal lattice volume (1/V**2) from A'
337
338	"""
339	G,g = A2Gmat(A)
340	rVsq = nl.det(G)
341	if rVsq < 0:
342	return 1
343	return rVsq
344
345	def calc_rV(A):
346	"""Compute the reciprocal lattice volume (V*) from A
347	"""
348	return np.sqrt(calc_rVsq(A))
349
350	def calc_V(A):
351	"""Compute the real lattice volume (V) from A
352	"""
353	return 1./calc_rV(A)
354
355	def A2invcell(A):
356	"""Compute reciprocal unit cell constants from A
357	returns tuple with a,b,c,alpha, beta, gamma (degrees)
358	"""
359	G,g = A2Gmat(A)
360	return Gmat2cell(G)
361
362	def Gmat2AB(G):
363	"""Computes orthogonalization matrix from reciprocal metric tensor G
364
365	:returns: tuple of two 3x3 numpy arrays (A,B)
366
367	* A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
368	* B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
369
370	"""
371	cellstar = Gmat2cell(G)
372	g = nl.inv(G)
373	cell = Gmat2cell(g)
374	A = np.zeros(shape=(3,3))
375	# from Giacovazzo (Fundamentals 2nd Ed.) p.75
376	A[0][0] = cell[0] # a
377	A[0][1] = cell[1]*cosd(cell[5]) # b cos(gamma)
378	A[0][2] = cell[2]*cosd(cell[4]) # c cos(beta)
379	A[1][1] = cell[1]*sind(cell[5]) # b sin(gamma)
380	A[1][2] = -cell[2]cosd(cellstar[3])sind(cell[4]) # - c cos(alpha*) sin(beta)
381	A[2][2] = 1/cellstar[2] # 1/c*
382	B = nl.inv(A)
383	return A,B
384
385	def cell2AB(cell):
386	"""Computes orthogonalization matrix from unit cell constants
387
388	:param tuple cell: a,b,c, alpha, beta, gamma (degrees)
389	:returns: tuple of two 3x3 numpy arrays (A,B)
390	A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
391	B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
392	"""
393	G,g = cell2Gmat(cell)
394	cellstar = Gmat2cell(G)
395	A = np.zeros(shape=(3,3))
396	# from Giacovazzo (Fundamentals 2nd Ed.) p.75
397	A[0][0] = cell[0] # a
398	A[0][1] = cell[1]*cosd(cell[5]) # b cos(gamma)
399	A[0][2] = cell[2]*cosd(cell[4]) # c cos(beta)
400	A[1][1] = cell[1]*sind(cell[5]) # b sin(gamma)
401	A[1][2] = -cell[2]cosd(cellstar[3])sind(cell[4]) # - c cos(alpha*) sin(beta)
402	A[2][2] = 1/cellstar[2] # 1/c*
403	B = nl.inv(A)
404	return A,B
405
406	def HKL2SpAng(H,cell,SGData):
407	"""Computes spherical coords for hkls; view along 001
408
409	:param array H: arrays of hkl
410	:param tuple cell: a,b,c, alpha, beta, gamma (degrees)
411	:param dict SGData: space group dictionary
412	:returns: arrays of r,phi,psi (radius,inclination,azimuth) about 001
413	"""
414	A,B = cell2AB(cell)
415	xH = np.inner(B.T,H)
416	r = np.sqrt(np.sum(xH**2,axis=0))
417	phi = acosd(xH[2]/r)
418	psi = atan2d(xH[1],xH[0])
419	phi = np.where(phi>90.,180.-phi,phi)
420	# GSASIIpath.IPyBreak()
421	return r,phi,psi
422
423	def U6toUij(U6):
424	"""Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
425	NB: there is a non numpy version in GSASIIspc: U2Uij
426
427	:param list U6: 6 terms of u11,u22,...
428	:returns:
429	Uij - numpy [3][3] array of uij
430	"""
431	U = np.array([
432	[U6[0], U6[3], U6[4]],
433	[U6[3], U6[1], U6[5]],
434	[U6[4], U6[5], U6[2]]])
435	return U
436
437	def UijtoU6(U):
438	"""Fill vector [U11,U22,U33,U12,U13,U23] from Uij
439	NB: there is a non numpy version in GSASIIspc: Uij2U
440	"""
441	U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
442	return U6
443
444	def betaij2Uij(betaij,G):
445	"""
446	Convert beta-ij to Uij tensors
447
448	:param beta-ij - numpy array [beta-ij]
449	:param G: reciprocal metric tensor
450	:returns: Uij: numpy array [Uij]
451	"""
452	ast = np.sqrt(np.diag(G)) #a, b, c*
453	Mast = np.multiply.outer(ast,ast)
454	return R2pisq*UijtoU6(U6toUij(betaij)/Mast)
455
456	def Uij2betaij(Uij,G):
457	"""
458	Convert Uij to beta-ij tensors -- stub for eventual completion
459
460	:param Uij: numpy array [Uij]
461	:param G: reciprocal metric tensor
462	:returns: beta-ij - numpy array [beta-ij]
463	"""
464	pass
465
466	def cell2GS(cell):
467	''' returns Uij to betaij conversion matrix'''
468	G,g = cell2Gmat(cell)
469	GS = G
470	GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
471	GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
472	GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
473	return GS
474
475	def Uij2Ueqv(Uij,GS,Amat):
476	''' returns 1/3 trace of diagonalized U matrix'''
477	U = np.multiply(U6toUij(Uij),GS)
478	U = np.inner(Amat,np.inner(U,Amat).T)
479	E,R = nl.eigh(U)
480	return np.sum(E)/3.
481
482	def CosAngle(U,V,G):
483	""" calculate cos of angle between U & V in generalized coordinates
484	defined by metric tensor G
485
486	:param U: 3-vectors assume numpy arrays, can be multiple reflections as (N,3) array
487	:param V: 3-vectors assume numpy arrays, only as (3) vector
488	:param G: metric tensor for U & V defined space assume numpy array
489	:returns:
490	cos(phi)
491	"""
492	u = (U.T/np.sqrt(np.sum(np.inner(U,G)*U,axis=1))).T
493	v = V/np.sqrt(np.inner(V,np.inner(G,V)))
494	cosP = np.inner(u,np.inner(G,v))
495	return cosP
496
497	def CosSinAngle(U,V,G):
498	""" calculate sin & cos of angle between U & V in generalized coordinates
499	defined by metric tensor G
500
501	:param U: 3-vectors assume numpy arrays
502	:param V: 3-vectors assume numpy arrays
503	:param G: metric tensor for U & V defined space assume numpy array
504	:returns:
505	cos(phi) & sin(phi)
506	"""
507	u = U/np.sqrt(np.inner(U,np.inner(G,U)))
508	v = V/np.sqrt(np.inner(V,np.inner(G,V)))
509	cosP = np.inner(u,np.inner(G,v))
510	sinP = np.sqrt(max(0.0,1.0-cosP**2))
511	return cosP,sinP
512
513	def criticalEllipse(prob):
514	"""
515	Calculate critical values for probability ellipsoids from probability
516	"""
517	if not ( 0.01 <= prob < 1.0):
518	return 1.54
519	coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
520	0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
521	llpr = math.log(-math.log(prob))
522	return np.polyval(coeff,llpr)
523
524	def CellBlock(nCells):
525	"""
526	Generate block of unit cells nnn on a side; [0,0,0] centered, n = 2*nCells+1
527	currently only works for nCells = 0 or 1 (not >1)
528	"""
529	if nCells:
530	N = 2*nCells+1
531	N2 = N*N
532	N3 = NNN
533	cellArray = []
534	A = np.array(range(N3))
535	cellGen = np.array([A/N2-1,A/N%N-1,A%N-1]).T
536	for cell in cellGen:
537	cellArray.append(cell)
538	return cellArray
539	else:
540	return [0,0,0]
541
542	def CellAbsorption(ElList,Volume):
543	'''Compute unit cell absorption
544
545	:param dict ElList: dictionary of element contents including mu and
546	number of atoms be cell
547	:param float Volume: unit cell volume
548	:returns: mu-total/Volume
549	'''
550	muT = 0
551	for El in ElList:
552	muT += ElList[El]['mu']*ElList[El]['FormulaNo']
553	return muT/Volume
554
555	#Permutations and Combinations
556	# Four routines: combinations,uniqueCombinations, selections & permutations
557	#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
558	#
559	def _combinators(_handle, items, n):
560	""" factored-out common structure of all following combinators """
561	if n==0:
562	yield [ ]
563	return
564	for i, item in enumerate(items):
565	this_one = [ item ]
566	for cc in _combinators(_handle, _handle(items, i), n-1):
567	yield this_one + cc
568	def combinations(items, n):
569	""" take n distinct items, order matters """
570	def skipIthItem(items, i):
571	return items[:i] + items[i+1:]
572	return _combinators(skipIthItem, items, n)
573	def uniqueCombinations(items, n):
574	""" take n distinct items, order is irrelevant """
575	def afterIthItem(items, i):
576	return items[i+1:]
577	return _combinators(afterIthItem, items, n)
578	def selections(items, n):
579	""" take n (not necessarily distinct) items, order matters """
580	def keepAllItems(items, i):
581	return items
582	return _combinators(keepAllItems, items, n)
583	def permutations(items):
584	""" take all items, order matters """
585	return combinations(items, len(items))
586
587	#reflection generation routines
588	#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
589	# cell - a,b,c,alp,bet,gam in A & deg
590
591	def Pos2dsp(Inst,pos):
592	''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
593	'''
594	if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
595	wave = G2mth.getWave(Inst)
596	return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
597	else: #'T'OF - ignore difB
598	return TOF2dsp(Inst,pos)
599
600	def TOF2dsp(Inst,Pos):
601	''' convert powder pattern TOF, musec to d-spacing by successive approximation
602	Pos can be numpy array
603	'''
604	def func(d,pos,Inst):
605	return (pos-Inst['difA'][1]d*2-Inst['Zero'][1]-Inst['difB'][1]/d)/Inst['difC'][1]
606	dsp0 = np.ones_like(Pos)
607	N = 0
608	while True: #successive approximations
609	dsp = func(dsp0,Pos,Inst)
610	if np.allclose(dsp,dsp0,atol=0.000001):
611	return dsp
612	dsp0 = dsp
613	N += 1
614	if N > 10:
615	return dsp
616
617	def Dsp2pos(Inst,dsp):
618	''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
619	'''
620	if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
621	wave = G2mth.getWave(Inst)
622	val = min(0.995,wave/(2.*dsp)) #set max at 168deg
623	pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]
624	else: #'T'OF
625	pos = Inst['difC'][1]dsp+Inst['Zero'][1]+Inst['difA'][1]dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
626	return pos
627
628	def getPeakPos(dataType,parmdict,dsp):
629	''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
630	'''
631	if 'C' in dataType:
632	pos = 2.0asind(parmdict['Lam']/(2.dsp))+parmdict['Zero']
633	else: #'T'OF
634	pos = parmdict['difC']dsp+parmdict['difA']dsp**2+parmdict['difB']/dsp+parmdict['Zero']
635	return pos
636
637	def calc_rDsq(H,A):
638	'needs doc string'
639	rdsq = H[0]H[0]A[0]+H[1]H[1]A[1]+H[2]H[2]A[2]+H[0]H[1]A[3]+H[0]H[2]A[4]+H[1]H[2]A[5]
640	return rdsq
641
642	def calc_rDsq2(H,G):
643	'needs doc string'
644	return np.inner(H,np.inner(G,H))
645
646	def calc_rDsqSS(H,A,vec):
647	'needs doc string'
648	rdsq = calc_rDsq(H[:3]+(H[3]*vec).T,A)
649	return rdsq
650
651	def calc_rDsqZ(H,A,Z,tth,lam):
652	'needs doc string'
653	rdsq = calc_rDsq(H,A)+Zsind(tth)2.0rpd/lam*2
654	return rdsq
655
656	def calc_rDsqZSS(H,A,vec,Z,tth,lam):
657	'needs doc string'
658	rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]vec).T,A)+Zsind(tth)2.0rpd/lam**2
659	return rdsq
660
661	def calc_rDsqT(H,A,Z,tof,difC):
662	'needs doc string'
663	rdsq = calc_rDsq(H,A)+Z/difC
664	return rdsq
665
666	def calc_rDsqTSS(H,A,vec,Z,tof,difC):
667	'needs doc string'
668	rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
669	return rdsq
670
671	def PlaneIntercepts(Amat,Bmat,H,phase,stack):
672	''' find unit cell intercepts for a stack of hkl planes
673	'''
674	Steps = range(-1,2,2)
675	if stack:
676	Steps = range(-10,10,1)
677	Stack = []
678	Ux = np.array([[0,0],[1,0],[1,1],[0,1]])
679	for step in Steps:
680	HX = []
681	for i in [0,1,2]:
682	if H[i]:
683	h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3]
684	for j in [0,1,2,3]:
685	hx = [0,0,0]
686	intcpt = (phase/360.+step-H[h]Ux[j,0]-H[k]Ux[j,1])/H[l]
687	if 0. <= intcpt <= 1.:
688	hx[h] = Ux[j,0]
689	hx[k] = Ux[j,1]
690	hx[l] = intcpt
691	HX.append(hx)
692	if len(HX)> 2:
693	HX = np.array(HX)
694	DX = np.inner(HX-HX[0],Amat)
695	D = np.sqrt(np.sum(DX**2,axis=1))
696	Dsort = np.argsort(D)
697	HX = HX[Dsort]
698	DX = DX[Dsort]
699	D = D[Dsort]
700	DX[1:,:] = DX[1:,:]/D[1:,nxs]
701	A = 2.*np.ones(HX.shape[0])
702	A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
703	HX = HX[np.argsort(A)]
704	# GSASIIpath.IPyBreak()
705	Stack.append(HX)
706	return Stack
707
708	def MaxIndex(dmin,A):
709	'needs doc string'
710	Hmax = [0,0,0]
711	try:
712	cell = A2cell(A)
713	except:
714	cell = [1,1,1,90,90,90]
715	for i in range(3):
716	Hmax[i] = int(round(cell[i]/dmin))
717	return Hmax
718
719	def transposeHKLF(transMat,Super,refList):
720	''' Apply transformation matrix to hkl(m)
721	param: transmat: 3x3 or 4x4 array
722	param: Super: 0 or 1 for extra index
723	param: refList list of h,k,l,....
724	return: newRefs transformed list of h',k',l',,,
725	return: badRefs list of noninteger h',k',l'...
726	'''
727	newRefs = np.copy(refList)
728	badRefs = []
729	for H in newRefs:
730	newH = np.inner(transMat,H[:3+Super])
731	H[:3+Super] = np.rint(newH)
732	if not np.allclose(newH,H[:3+Super],atol=0.01):
733	badRefs.append(newH)
734	return newRefs,badRefs
735
736	def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
737	'''sort reflection list on d-spacing; can sort in either order
738
739	:param HKLd: a list of [h,k,l,d,...];
740	:param ifreverse: True for largest d first
741	:param ifdup: True if duplicate d-spacings allowed
742	:return sorted reflection list
743	'''
744	T = []
745	N = 3
746	if ifSS:
747	N = 4
748	for i,H in enumerate(HKLd):
749	if ifdup:
750	T.append((H[N],i))
751	else:
752	T.append(H[N])
753	D = dict(zip(T,HKLd))
754	T.sort()
755	if ifreverse:
756	T.reverse()
757	X = []
758	okey = ''
759	for key in T:
760	if key != okey: X.append(D[key]) #remove duplicate d-spacings
761	okey = key
762	return X
763
764	def SwapIndx(Axis,H):
765	'needs doc string'
766	if Axis in [1,-1]:
767	return H
768	elif Axis in [2,-3]:
769	return [H[1],H[2],H[0]]
770	else:
771	return [H[2],H[0],H[1]]
772
773	def Rh2Hx(Rh):
774	'needs doc string'
775	Hx = [0,0,0]
776	Hx[0] = Rh[0]-Rh[1]
777	Hx[1] = Rh[1]-Rh[2]
778	Hx[2] = np.sum(Rh)
779	return Hx
780
781	def Hx2Rh(Hx):
782	'needs doc string'
783	Rh = [0,0,0]
784	itk = -Hx[0]+Hx[1]+Hx[2]
785	if itk%3 != 0:
786	return 0 #error - not rhombohedral reflection
787	else:
788	Rh[1] = itk/3
789	Rh[0] = Rh[1]+Hx[0]
790	Rh[2] = Rh[1]-Hx[1]
791	if Rh[0] < 0:
792	for i in range(3):
793	Rh[i] = -Rh[i]
794	return Rh
795
796	def CentCheck(Cent,H):
797	'needs doc string'
798	h,k,l = H
799	if Cent == 'A' and (k+l)%2:
800	return False
801	elif Cent == 'B' and (h+l)%2:
802	return False
803	elif Cent == 'C' and (h+k)%2:
804	return False
805	elif Cent == 'I' and (h+k+l)%2:
806	return False
807	elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
808	return False
809	elif Cent == 'R' and (-h+k+l)%3:
810	return False
811	else:
812	return True
813
814	def GetBraviasNum(center,system):
815	"""Determine the Bravais lattice number, as used in GenHBravais
816
817	:param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
818	:param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
819	'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
820	:return: a number between 0 and 13
821	or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
822
823	"""
824	if center.upper() == 'F' and system.lower() == 'cubic':
825	return 0
826	elif center.upper() == 'I' and system.lower() == 'cubic':
827	return 1
828	elif center.upper() == 'P' and system.lower() == 'cubic':
829	return 2
830	elif center.upper() == 'R' and system.lower() == 'trigonal':
831	return 3
832	elif center.upper() == 'P' and system.lower() == 'hexagonal':
833	return 4
834	elif center.upper() == 'I' and system.lower() == 'tetragonal':
835	return 5
836	elif center.upper() == 'P' and system.lower() == 'tetragonal':
837	return 6
838	elif center.upper() == 'F' and system.lower() == 'orthorhombic':
839	return 7
840	elif center.upper() == 'I' and system.lower() == 'orthorhombic':
841	return 8
842	elif center.upper() == 'C' and system.lower() == 'orthorhombic':
843	return 9
844	elif center.upper() == 'P' and system.lower() == 'orthorhombic':
845	return 10
846	elif center.upper() == 'C' and system.lower() == 'monoclinic':
847	return 11
848	elif center.upper() == 'P' and system.lower() == 'monoclinic':
849	return 12
850	elif center.upper() == 'P' and system.lower() == 'triclinic':
851	return 13
852	raise ValueError,'non-standard Bravais lattice center=%s, cell=%s' % (center,system)
853
854	def GenHBravais(dmin,Bravais,A):
855	"""Generate the positionally unique powder diffraction reflections
856
857	:param dmin: minimum d-spacing in A
858	:param Bravais: lattice type (see GetBraviasNum). Bravais is one of::
859	0 F cubic
860	1 I cubic
861	2 P cubic
862	3 R hexagonal (trigonal not rhombohedral)
863	4 P hexagonal
864	5 I tetragonal
865	6 P tetragonal
866	7 F orthorhombic
867	8 I orthorhombic
868	9 C orthorhombic
869	10 P orthorhombic
870	11 C monoclinic
871	12 P monoclinic
872	13 P triclinic
873
874	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
875	:return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
876
877	"""
878	import math
879	if Bravais in [9,11]:
880	Cent = 'C'
881	elif Bravais in [1,5,8]:
882	Cent = 'I'
883	elif Bravais in [0,7]:
884	Cent = 'F'
885	elif Bravais in [3]:
886	Cent = 'R'
887	else:
888	Cent = 'P'
889	Hmax = MaxIndex(dmin,A)
890	dminsq = 1./(dmin**2)
891	HKL = []
892	if Bravais == 13: #triclinic
893	for l in range(-Hmax[2],Hmax[2]+1):
894	for k in range(-Hmax[1],Hmax[1]+1):
895	hmin = 0
896	if (k < 0): hmin = 1
897	if (k ==0 and l < 0): hmin = 1
898	for h in range(hmin,Hmax[0]+1):
899	H=[h,k,l]
900	rdsq = calc_rDsq(H,A)
901	if 0 < rdsq <= dminsq:
902	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
903	elif Bravais in [11,12]: #monoclinic - b unique
904	Hmax = SwapIndx(2,Hmax)
905	for h in range(Hmax[0]+1):
906	for k in range(-Hmax[1],Hmax[1]+1):
907	lmin = 0
908	if k < 0:lmin = 1
909	for l in range(lmin,Hmax[2]+1):
910	[h,k,l] = SwapIndx(-2,[h,k,l])
911	H = []
912	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
913	if H:
914	rdsq = calc_rDsq(H,A)
915	if 0 < rdsq <= dminsq:
916	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
917	[h,k,l] = SwapIndx(2,[h,k,l])
918	elif Bravais in [7,8,9,10]: #orthorhombic
919	for h in range(Hmax[0]+1):
920	for k in range(Hmax[1]+1):
921	for l in range(Hmax[2]+1):
922	H = []
923	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
924	if H:
925	rdsq = calc_rDsq(H,A)
926	if 0 < rdsq <= dminsq:
927	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
928	elif Bravais in [5,6]: #tetragonal
929	for l in range(Hmax[2]+1):
930	for k in range(Hmax[1]+1):
931	for h in range(k,Hmax[0]+1):
932	H = []
933	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
934	if H:
935	rdsq = calc_rDsq(H,A)
936	if 0 < rdsq <= dminsq:
937	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
938	elif Bravais in [3,4]:
939	lmin = 0
940	if Bravais == 3: lmin = -Hmax[2] #hexagonal/trigonal
941	for l in range(lmin,Hmax[2]+1):
942	for k in range(Hmax[1]+1):
943	hmin = k
944	if l < 0: hmin += 1
945	for h in range(hmin,Hmax[0]+1):
946	H = []
947	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
948	if H:
949	rdsq = calc_rDsq(H,A)
950	if 0 < rdsq <= dminsq:
951	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
952
953	else: #cubic
954	for l in range(Hmax[2]+1):
955	for k in range(l,Hmax[1]+1):
956	for h in range(k,Hmax[0]+1):
957	H = []
958	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
959	if H:
960	rdsq = calc_rDsq(H,A)
961	if 0 < rdsq <= dminsq:
962	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
963	return sortHKLd(HKL,True,False)
964
965	def getHKLmax(dmin,SGData,A):
966	'finds maximum allowed hkl for given A within dmin'
967	SGLaue = SGData['SGLaue']
968	if SGLaue in ['3R','3mR']: #Rhombohedral axes
969	Hmax = [0,0,0]
970	cell = A2cell(A)
971	aHx = cell[0]math.sqrt(2.0(1.0-cosd(cell[3])))
972	cHx = cell[0]math.sqrt(3.0(1.0+2.0*cosd(cell[3])))
973	Hmax[0] = Hmax[1] = int(round(aHx/dmin))
974	Hmax[2] = int(round(cHx/dmin))
975	#print Hmax,aHx,cHx
976	else: # all others
977	Hmax = MaxIndex(dmin,A)
978	return Hmax
979
980	def GenHLaue(dmin,SGData,A):
981	"""Generate the crystallographically unique powder diffraction reflections
982	for a lattice and Bravais type
983
984	:param dmin: minimum d-spacing
985	:param SGData: space group dictionary with at least
986
987	* 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
988	* 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
989	* 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
990
991	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
992	:return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
993	part of reciprocal space ignoring anomalous dispersion
994
995	"""
996	import math
997	SGLaue = SGData['SGLaue']
998	SGLatt = SGData['SGLatt']
999	SGUniq = SGData['SGUniq']
1000	#finds maximum allowed hkl for given A within dmin
1001	Hmax = getHKLmax(dmin,SGData,A)
1002
1003	dminsq = 1./(dmin**2)
1004	HKL = []
1005	if SGLaue == '-1': #triclinic
1006	for l in range(-Hmax[2],Hmax[2]+1):
1007	for k in range(-Hmax[1],Hmax[1]+1):
1008	hmin = 0
1009	if (k < 0) or (k ==0 and l < 0): hmin = 1
1010	for h in range(hmin,Hmax[0]+1):
1011	H = []
1012	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1013	if H:
1014	rdsq = calc_rDsq(H,A)
1015	if 0 < rdsq <= dminsq:
1016	HKL.append([h,k,l,1/math.sqrt(rdsq)])
1017	elif SGLaue == '2/m': #monoclinic
1018	axisnum = 1 + ['a','b','c'].index(SGUniq)
1019	Hmax = SwapIndx(axisnum,Hmax)
1020	for h in range(Hmax[0]+1):
1021	for k in range(-Hmax[1],Hmax[1]+1):
1022	lmin = 0
1023	if k < 0:lmin = 1
1024	for l in range(lmin,Hmax[2]+1):
1025	[h,k,l] = SwapIndx(-axisnum,[h,k,l])
1026	H = []
1027	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1028	if H:
1029	rdsq = calc_rDsq(H,A)
1030	if 0 < rdsq <= dminsq:
1031	HKL.append([h,k,l,1/math.sqrt(rdsq)])
1032	[h,k,l] = SwapIndx(axisnum,[h,k,l])
1033	elif SGLaue in ['mmm','4/m','6/m']: #orthorhombic
1034	for l in range(Hmax[2]+1):
1035	for h in range(Hmax[0]+1):
1036	kmin = 1
1037	if SGLaue == 'mmm' or h ==0: kmin = 0
1038	for k in range(kmin,Hmax[1]+1):
1039	H = []
1040	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1041	if H:
1042	rdsq = calc_rDsq(H,A)
1043	if 0 < rdsq <= dminsq:
1044	HKL.append([h,k,l,1/math.sqrt(rdsq)])
1045	elif SGLaue in ['4/mmm','6/mmm']: #tetragonal & hexagonal
1046	for l in range(Hmax[2]+1):
1047	for h in range(Hmax[0]+1):
1048	for k in range(h+1):
1049	H = []
1050	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1051	if H:
1052	rdsq = calc_rDsq(H,A)
1053	if 0 < rdsq <= dminsq:
1054	HKL.append([h,k,l,1/math.sqrt(rdsq)])
1055	elif SGLaue in ['3m1','31m','3','3R','3mR']: #trigonals
1056	for l in range(-Hmax[2],Hmax[2]+1):
1057	hmin = 0
1058	if l < 0: hmin = 1
1059	for h in range(hmin,Hmax[0]+1):
1060	if SGLaue in ['3R','3']:
1061	kmax = h
1062	kmin = -int((h-1.)/2.)
1063	else:
1064	kmin = 0
1065	kmax = h
1066	if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1067	if SGLaue == '31m' and l < 0: kmin = 1
1068	for k in range(kmin,kmax+1):
1069	H = []
1070	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1071	if SGLaue in ['3R','3mR']:
1072	H = Hx2Rh(H)
1073	if H:
1074	rdsq = calc_rDsq(H,A)
1075	if 0 < rdsq <= dminsq:
1076	HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)])
1077	else: #cubic
1078	for h in range(Hmax[0]+1):
1079	for k in range(h+1):
1080	lmin = 0
1081	lmax = k
1082	if SGLaue =='m3':
1083	lmax = h-1
1084	if h == k: lmax += 1
1085	for l in range(lmin,lmax+1):
1086	H = []
1087	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1088	if H:
1089	rdsq = calc_rDsq(H,A)
1090	if 0 < rdsq <= dminsq:
1091	HKL.append([h,k,l,1/math.sqrt(rdsq)])
1092	return sortHKLd(HKL,True,True)
1093
1094	def GenPfHKLs(nMax,SGData,A):
1095	"""Generate the unique pole figure reflections for a lattice and Bravais type.
1096	Min d-spacing=1.0A & no more than nMax returned
1097
1098	:param nMax: maximum number of hkls returned
1099	:param SGData: space group dictionary with at least
1100
1101	* 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1102	* 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1103	* 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1104
1105	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
1106	:return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1107
1108	"""
1109	HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T #strip d-spacings
1110	N = min(nMax,len(HKL))
1111	return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]
1112
1113	def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1114	'needs a doc string'
1115	HKLs = []
1116	vec = np.array(Vec)
1117	vstar = np.sqrt(calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1118	dvec = 1./(maxH*vstar+1./dmin)
1119	HKL = GenHLaue(dvec,SGData,A)
1120	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1121	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1122	for h,k,l,d in HKL:
1123	ext = G2spc.GenHKLf([h,k,l],SGData)[0] #h,k,l must be integral values here
1124	if not ext and d >= dmin:
1125	HKLs.append([h,k,l,0,d])
1126	for dH in SSdH:
1127	if dH:
1128	DH = SSdH[dH]
1129	H = [h+DH[0],k+DH[1],l+DH[2]]
1130	d = 1/np.sqrt(calc_rDsq(H,A))
1131	if d >= dmin:
1132	HKLM = np.array([h,k,l,dH])
1133	if G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData):
1134	HKLs.append([h,k,l,dH,d])
1135	return HKLs
1136
1137	def LaueUnique2(SGData,refList):
1138	''' Impose Laue symmetry on hkl
1139	:param SGData: space group data from 'P '+Laue
1140	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1141
1142	:return: HKLF new reflection array with imposed Laue symmetry
1143	'''
1144	for ref in refList:
1145	H = ref[:3]
1146	Uniq = G2spc.GenHKLf(H,SGData)[2]
1147	Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1148	ref[:3] = Uniq[-1]
1149	return refList
1150
1151	def LaueUnique(Laue,HKLF):
1152	''' Impose Laue symmetry on hkl
1153	:param Laue: str Laue symbol
1154	centrosymmetric Laue groups
1155	['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3','-31m','-3m1',
1156	'6/m','6/mmm','m3','m3m']
1157	noncentrosymmetric Laue groups
1158	['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1159	'4','-4','422','4mm','3','312','321','31m','3m1',
1160	'6','-6','622','6mm','-62m','-6m2','23','432','-43m']
1161	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1162
1163	:return: HKLF new reflection array with imposed Laue symmetry
1164	'''
1165
1166	HKLFT = HKLF.T
1167	mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]]) #hkl -> k,-h,l
1168	mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]]) #hkl -> -k,h,l
1169	mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1170	mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]]) #hkl -> ihl = -h-k,h,l
1171	mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]]) #hkl -> kil = k,-h-k,l
1172	matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]]) #hkl -> k,l,h
1173	matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]]) #hkl -> -l,-h,k
1174	matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]]) #hkl -> -l,h,-k
1175	matd3p = np.array([[0,1,0],[0,0,-1],[-1,0,0]]) #hkl -> k,-l,-h
1176	mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65
1177	matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l
1178	matdmt = np.array([[0,-1,0],[-1,0,0],[0,0,1]]) #hkl -> -k,-h,l
1179	matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]]) #hkl -> -h-k,k,l
1180	matdmq = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l
1181	matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l
1182	matkmp = np.array([[1,0,0],[-1,-1,0],[0,0,1]]) #hkl -> h,-h-k,l
1183	matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l
1184	matd2p = np.array([[-1,-1,0],[0,1,0],[0,0,-1]]) #hkl -> -h-k,k,-l
1185	matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]]) #hkl -> h,l,k
1186	mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,-h,l
1187	math2 = np.array([[0,-1,0],[-1,0,0],[0,0,-1]]) #hkl -> -k,-h,-l
1188	matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l
1189	#triclinic
1190	if Laue == '1': #ok
1191	pass
1192	elif Laue == '-1': #ok
1193	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1194	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1195	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1196	#monoclinic
1197	#noncentrosymmetric - all ok
1198	elif Laue == '2':
1199	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1200	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1201	elif Laue == '1 1 2':
1202	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1203	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1204	elif Laue == '2 1 1':
1205	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1206	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1207	elif Laue == 'm':
1208	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1209	elif Laue == 'm 1 1':
1210	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1211	elif Laue == '1 1 m':
1212	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1213	#centrosymmetric - all ok
1214	elif Laue == '2/m 1 1':
1215	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1216	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1217	HKLFT[:3] = np.where((HKLFT[2]HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]np.array([1,-1,1])[:,nxs],HKLFT[:3])
1218	elif Laue == '2/m':
1219	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1220	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1221	HKLFT[:3] = np.where((HKLFT[0]HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1222	elif Laue == '1 1 2/m':
1223	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1224	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1225	HKLFT[:3] = np.where((HKLFT[1]HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]np.array([-1,1,1])[:,nxs],HKLFT[:3])
1226	#orthorhombic
1227	#noncentrosymmetric - all OK
1228	elif Laue == '2 2 2':
1229	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1230	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1231	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1232	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1233	elif Laue == 'm m 2':
1234	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1235	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1236	elif Laue == '2 m m':
1237	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1238	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1239	elif Laue == 'm 2 m':
1240	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1241	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1242	#centrosymmetric - all ok
1243	elif Laue == 'm m m':
1244	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1245	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1246	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1247	#tetragonal
1248	#noncentrosymmetric - all ok
1249	elif Laue == '4':
1250	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1251	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1252	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1253	elif Laue == '-4':
1254	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1255	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1256	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1257	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1258	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1259	elif Laue == '4 2 2':
1260	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1261	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1262	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1263	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1264	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1265	elif Laue == '4 m m':
1266	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1267	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1268	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1269	HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1270	elif Laue == '-4 2 m':
1271	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1272	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1273	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1274	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1275	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1276	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1277	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1278	elif Laue == '-4 m 2':
1279	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1280	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1281	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1282	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1283	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1284	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1285	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1286	#centrosymmetric - all ok
1287	elif Laue == '4/m':
1288	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1289	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1290	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1291	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1292	elif Laue == '4/m m m':
1293	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1294	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1295	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1296	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1297	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1298	#trigonal - all hex cell
1299	#noncentrosymmetric - all ok
1300	elif Laue == '3':
1301	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1302	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1303	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1304	elif Laue == '3 1 2':
1305	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1306	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1307	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1308	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1309	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1310	elif Laue == '3 2 1':
1311	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1312	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1313	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1314	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1315	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1316	HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2HKLFT[1]),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1317	elif Laue == '3 1 m':
1318	HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1319	HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1320	HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1321	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1322	elif Laue == '3 m 1':
1323	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1324	HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1325	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1326	#centrosymmetric
1327	elif Laue == '-3': #ok
1328	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1329	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1330	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1331	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1332	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1333	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1334	elif Laue == '-3 m 1': #ok
1335	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1336	HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1337	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1338	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1339	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1340	elif Laue == '-3 1 m': #ok
1341	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1342	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1343	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1344	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1345	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1346	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1347	#hexagonal
1348	#noncentrosymmetric
1349	elif Laue == '6': #ok
1350	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1351	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1352	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1353	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1354	elif Laue == '-6': #ok
1355	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1356	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1357	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1358	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1359	elif Laue == '6 2 2': #ok
1360	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1361	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1362	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1363	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1364	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1365	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1366	elif Laue == '6 m m': #ok
1367	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1368	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1369	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1370	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1371	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1372	elif Laue == '-6 m 2': #ok
1373	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1374	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1375	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1376	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1377	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1378	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1379	elif Laue == '-6 2 m': #ok
1380	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1381	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1382	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1383	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1384	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1385	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1386	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1387	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1388	#centrosymmetric
1389	elif Laue == '6/m': #ok
1390	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1391	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1392	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1393	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1394	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1395	elif Laue == '6/m m m': #ok
1396	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1397	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1398	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1399	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1400	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1401	#cubic - all ok
1402	#noncentrosymmetric -
1403	elif Laue == '2 3':
1404	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1405	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1406	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1407	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1408	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1409	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1410	HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])\|(HKLFT[1]>-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1411	HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])\|(HKLFT[1]>=-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1412	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1413	elif Laue == '4 3 2':
1414	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1415	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1416	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1417	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1418	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1419	HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1420	HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1421	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1422	elif Laue == '-4 3 m':
1423	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1424	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1425	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1426	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1427	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1428	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1429	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1430	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1431	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1432	HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1433	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1434	HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]<-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3q[nxs,:,:])).T,HKLFT[:3])
1435	HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]>=-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matdm3[nxs,:,:])).T,HKLFT[:3])
1436	#centrosymmetric
1437	elif Laue == 'm 3':
1438	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1439	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1440	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1441	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1442	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1443	elif Laue == 'm 3 m':
1444	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1445	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1446	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1447	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1448	HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1449	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1450	return HKLFT.T
1451
1452
1453	#Spherical harmonics routines
1454	def OdfChk(SGLaue,L,M):
1455	'needs doc string'
1456	if not L%2 and abs(M) <= L:
1457	if SGLaue == '0': #cylindrical symmetry
1458	if M == 0: return True
1459	elif SGLaue == '-1':
1460	return True
1461	elif SGLaue == '2/m':
1462	if not abs(M)%2: return True
1463	elif SGLaue == 'mmm':
1464	if not abs(M)%2 and M >= 0: return True
1465	elif SGLaue == '4/m':
1466	if not abs(M)%4: return True
1467	elif SGLaue == '4/mmm':
1468	if not abs(M)%4 and M >= 0: return True
1469	elif SGLaue in ['3R','3']:
1470	if not abs(M)%3: return True
1471	elif SGLaue in ['3mR','3m1','31m']:
1472	if not abs(M)%3 and M >= 0: return True
1473	elif SGLaue == '6/m':
1474	if not abs(M)%6: return True
1475	elif SGLaue == '6/mmm':
1476	if not abs(M)%6 and M >= 0: return True
1477	elif SGLaue == 'm3':
1478	if M > 0:
1479	if L%12 == 2:
1480	if M <= L/12: return True
1481	else:
1482	if M <= L/12+1: return True
1483	elif SGLaue == 'm3m':
1484	if M > 0:
1485	if L%12 == 2:
1486	if M <= L/12: return True
1487	else:
1488	if M <= L/12+1: return True
1489	return False
1490
1491	def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1492	'needs doc string'
1493	coeffNames = []
1494	for iord in [2*i+2 for i in range(L/2)]:
1495	for m in [i-iord for i in range(2*iord+1)]:
1496	if OdfChk(SamSym,iord,m):
1497	for n in [i-iord for i in range(2*iord+1)]:
1498	if OdfChk(SGLaue,iord,n):
1499	if IfLMN:
1500	coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1501	else:
1502	coeffNames.append('C(%d,%d)'%(iord,n))
1503	return coeffNames
1504
1505	def CrsAng(H,cell,SGData):
1506	'needs doc string'
1507	a,b,c,al,be,ga = cell
1508	SQ3 = 1.732050807569
1509	H1 = np.array([1,0,0])
1510	H2 = np.array([0,1,0])
1511	H3 = np.array([0,0,1])
1512	H4 = np.array([1,1,1])
1513	G,g = cell2Gmat(cell)
1514	Laue = SGData['SGLaue']
1515	Naxis = SGData['SGUniq']
1516	if len(H.shape) == 1:
1517	DH = np.inner(H,np.inner(G,H))
1518	else:
1519	DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1520	if Laue == '2/m':
1521	if Naxis == 'a':
1522	DR = np.inner(H1,np.inner(G,H1))
1523	DHR = np.inner(H,np.inner(G,H1))
1524	elif Naxis == 'b':
1525	DR = np.inner(H2,np.inner(G,H2))
1526	DHR = np.inner(H,np.inner(G,H2))
1527	else:
1528	DR = np.inner(H3,np.inner(G,H3))
1529	DHR = np.inner(H,np.inner(G,H3))
1530	elif Laue in ['R3','R3m']:
1531	DR = np.inner(H4,np.inner(G,H4))
1532	DHR = np.inner(H,np.inner(G,H4))
1533	else:
1534	DR = np.inner(H3,np.inner(G,H3))
1535	DHR = np.inner(H,np.inner(G,H3))
1536	DHR /= np.sqrt(DR*DH)
1537	phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1538	if Laue == '-1':
1539	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1540	BB = H.T[0]sind(ga)*2
1541	elif Laue == '2/m':
1542	if Naxis == 'a':
1543	BA = H.T[2]b/(c-H.T[1]cosd(al))
1544	BB = H.T[1]sind(al)*2
1545	elif Naxis == 'b':
1546	BA = H.T[0]c/(a-H.T[2]cosd(be))
1547	BB = H.T[2]sind(be)*2
1548	else:
1549	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1550	BB = H.T[0]sind(ga)*2
1551	elif Laue in ['mmm','4/m','4/mmm']:
1552	BA = H.T[1]*a
1553	BB = H.T[0]*b
1554	elif Laue in ['3R','3mR']:
1555	BA = H.T[0]+H.T[1]-2.0*H.T[2]
1556	BB = SQ3*(H.T[0]-H.T[1])
1557	elif Laue in ['m3','m3m']:
1558	BA = H.T[1]
1559	BB = H.T[0]
1560	else:
1561	BA = H.T[0]+2.0*H.T[1]
1562	BB = SQ3*H.T[0]
1563	beta = atan2d(BA,BB)
1564	return phi,beta
1565
1566	def SamAng(Tth,Gangls,Sangl,IFCoup):
1567	"""Compute sample orientation angles vs laboratory coord. system
1568
1569	:param Tth: Signed theta
1570	:param Gangls: Sample goniometer angles phi,chi,omega,azmuth
1571	:param Sangl: Sample angle zeros om-0, chi-0, phi-0
1572	:param IFCoup: True if omega & 2-theta coupled in CW scan
1573	:returns:
1574	psi,gam: Sample odf angles
1575	dPSdA,dGMdA: Angle zero derivatives
1576	"""
1577
1578	if IFCoup:
1579	GSomeg = sind(Gangls[2]+Tth)
1580	GComeg = cosd(Gangls[2]+Tth)
1581	else:
1582	GSomeg = sind(Gangls[2])
1583	GComeg = cosd(Gangls[2])
1584	GSTth = sind(Tth)
1585	GCTth = cosd(Tth)
1586	GSazm = sind(Gangls[3])
1587	GCazm = cosd(Gangls[3])
1588	GSchi = sind(Gangls[1])
1589	GCchi = cosd(Gangls[1])
1590	GSphi = sind(Gangls[0]+Sangl[2])
1591	GCphi = cosd(Gangls[0]+Sangl[2])
1592	SSomeg = sind(Sangl[0])
1593	SComeg = cosd(Sangl[0])
1594	SSchi = sind(Sangl[1])
1595	SCchi = cosd(Sangl[1])
1596	AT = -GSTthGComeg+GCTthGCazm*GSomeg
1597	BT = GSTthGSomeg+GCTthGCazm*GComeg
1598	CT = -GCTthGSazmGSchi
1599	DT = -GCTthGSazmGCchi
1600
1601	BC1 = -ATGSphi+(CT+BTGCchi)*GCphi
1602	BC2 = DT-BT*GSchi
1603	BC3 = ATGCphi+(CT+BTGCchi)*GSphi
1604
1605	BC = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1606	psi = acosd(BC)
1607
1608	BD = 1.0-BC**2
1609	C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1610	dPSdA = [-C(-BC1SSomegSCchi-BC2SSomegSSchi-BC3SComeg),
1611	-C(-BC1SComegSSchi+BC2SComeg*SCchi),
1612	-C(-BC1SSomeg-BC3SComegSCchi)]
1613
1614	BA = -BC1SSchi+BC2SCchi
1615	BB = BC1SSomegSCchi+BC2SSomegSSchi+BC3*SComeg
1616	gam = atan2d(BB,BA)
1617
1618	BD = (BA2+BB2)/rpd
1619
1620	dBAdO = 0
1621	dBAdC = -BC1SCchi-BC2SSchi
1622	dBAdF = BC3*SSchi
1623
1624	dBBdO = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1625	dBBdC = -BC1SSomegSSchi+BC2SSomegSCchi
1626	dBBdF = BC1SComeg-BC3SSomeg*SCchi
1627
1628	dGMdA = np.where(BD > 1.e-6,[(BAdBBdO-BBdBAdO)/BD,(BAdBBdC-BBdBAdC)/BD, \
1629	(BAdBBdF-BBdBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1630
1631	return psi,gam,dPSdA,dGMdA
1632
1633	BOH = {
1634	'L=2':[[],[],[]],
1635	'L=4':[[0.30469720,0.36418281],[],[]],
1636	'L=6':[[-0.14104740,0.52775103],[],[]],
1637	'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1638	'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1639	'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1640	[0.09298802,-0.23773812,0.49446631,0.0],[]],
1641	'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1642	'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1643	[0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1644	'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1645	[-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1646	'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1647	[0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1648	'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1649	[-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1650	'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1651	[0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1652	[0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1653	'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1654	[-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1655	'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1656	[0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1657	[0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1658	'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1659	[-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1660	[-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1661	'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1662	[0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1663	[0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1664	'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1665	[-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1666	[-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1667	}
1668
1669	Lnorm = lambda L: 4.np.pi/(2.0L+1.)
1670
1671	def GetKcl(L,N,SGLaue,phi,beta):
1672	'needs doc string'
1673	import pytexture as ptx
1674	if SGLaue in ['m3','m3m']:
1675	if 'array' in str(type(phi)) and np.any(phi.shape):
1676	Kcl = np.zeros_like(phi)
1677	else:
1678	Kcl = 0.
1679	for j in range(0,L+1,4):
1680	im = j/4
1681	if 'array' in str(type(phi)) and np.any(phi.shape):
1682	pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1683	else:
1684	pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1685	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
1686	else:
1687	if 'array' in str(type(phi)) and np.any(phi.shape):
1688	pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1689	else:
1690	pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1691	pcrs *= RSQ2PI
1692	if N:
1693	pcrs *= SQ2
1694	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1695	if SGLaue in ['3mR','3m1','31m']:
1696	if N%6 == 3:
1697	Kcl = pcrssind(Nbeta)
1698	else:
1699	Kcl = pcrscosd(Nbeta)
1700	else:
1701	Kcl = pcrscosd(Nbeta)
1702	else:
1703	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
1704	return Kcl
1705
1706	def GetKsl(L,M,SamSym,psi,gam):
1707	'needs doc string'
1708	import pytexture as ptx
1709	if 'array' in str(type(psi)) and np.any(psi.shape):
1710	psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1711	else:
1712	psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1713	psrs *= RSQ2PI
1714	dpdps *= RSQ2PI
1715	if M:
1716	psrs *= SQ2
1717	dpdps *= SQ2
1718	if SamSym in ['mmm',]:
1719	dum = cosd(M*gam)
1720	Ksl = psrs*dum
1721	dKsdp = dpdps*dum
1722	dKsdg = -psrsMsind(M*gam)
1723	else:
1724	dum = cosd(Mgam)+sind(Mgam)
1725	Ksl = psrs*dum
1726	dKsdp = dpdps*dum
1727	dKsdg = psrsM(-sind(Mgam)+cosd(Mgam))
1728	return Ksl,dKsdp,dKsdg
1729
1730	def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1731	"""
1732	This is used for spherical harmonics description of preferred orientation;
1733	cylindrical symmetry only (M=0) and no sample angle derivatives returned
1734	"""
1735	import pytexture as ptx
1736	Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1737	Ksl *= RSQ2PI
1738	if SGLaue in ['m3','m3m']:
1739	Kcl = 0.0
1740	for j in range(0,L+1,4):
1741	im = j/4
1742	pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1743	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
1744	else:
1745	pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1746	pcrs *= RSQ2PI
1747	if N:
1748	pcrs *= SQ2
1749	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1750	if SGLaue in ['3mR','3m1','31m']:
1751	if N%6 == 3:
1752	Kcl = pcrssind(Nbeta)
1753	else:
1754	Kcl = pcrscosd(Nbeta)
1755	else:
1756	Kcl = pcrscosd(Nbeta)
1757	else:
1758	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
1759	return Kcl*Ksl,Lnorm(L)
1760
1761	def Glnh(Start,SHCoef,psi,gam,SamSym):
1762	'needs doc string'
1763	import pytexture as ptx
1764
1765	if Start:
1766	ptx.pyqlmninit()
1767	Start = False
1768	Fln = np.zeros(len(SHCoef))
1769	for i,term in enumerate(SHCoef):
1770	l,m,n = eval(term.strip('C'))
1771	pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1772	pcrs *= RSQPI
1773	if m == 0:
1774	pcrs /= SQ2
1775	if SamSym in ['mmm',]:
1776	Ksl = pcrscosd(mgam)
1777	else:
1778	Ksl = pcrs(cosd(mgam)+sind(m*gam))
1779	Fln[i] = SHCoef[term]KslLnorm(l)
1780	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1781	return ODFln
1782
1783	def Flnh(Start,SHCoef,phi,beta,SGData):
1784	'needs doc string'
1785	import pytexture as ptx
1786
1787	if Start:
1788	ptx.pyqlmninit()
1789	Start = False
1790	Fln = np.zeros(len(SHCoef))
1791	for i,term in enumerate(SHCoef):
1792	l,m,n = eval(term.strip('C'))
1793	if SGData['SGLaue'] in ['m3','m3m']:
1794	Kcl = 0.0
1795	for j in range(0,l+1,4):
1796	im = j/4
1797	pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1798	Kcl += BOH['L='+str(l)][n-1][im]pcrscosd(j*beta)
1799	else: #all but cubic
1800	pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1801	pcrs *= RSQPI
1802	if n == 0:
1803	pcrs /= SQ2
1804	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1805	if SGData['SGLaue'] in ['3mR','3m1','31m']:
1806	if n%6 == 3:
1807	Kcl = pcrssind(nbeta)
1808	else:
1809	Kcl = pcrscosd(nbeta)
1810	else:
1811	Kcl = pcrscosd(nbeta)
1812	else:
1813	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
1814	Fln[i] = SHCoef[term]KclLnorm(l)
1815	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1816	return ODFln
1817
1818	def polfcal(ODFln,SamSym,psi,gam):
1819	'''Perform a pole figure computation.
1820	Note that the the number of gam values must either be 1 or must
1821	match psi. Updated for numpy 1.8.0
1822	'''
1823	import pytexture as ptx
1824	PolVal = np.ones_like(psi)
1825	for term in ODFln:
1826	if abs(ODFln[term][1]) > 1.e-3:
1827	l,m,n = eval(term.strip('C'))
1828	psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1829	if SamSym in ['-1','2/m']:
1830	if m:
1831	Ksl = RSQPIpsrs(cosd(mgam)+sind(mgam))
1832	else:
1833	Ksl = RSQPI*psrs/SQ2
1834	else:
1835	if m:
1836	Ksl = RSQPIpsrscosd(m*gam)
1837	else:
1838	Ksl = RSQPI*psrs/SQ2
1839	PolVal += ODFln[term][1]*Ksl
1840	return PolVal
1841
1842	def invpolfcal(ODFln,SGData,phi,beta):
1843	'needs doc string'
1844	import pytexture as ptx
1845
1846	invPolVal = np.ones_like(beta)
1847	for term in ODFln:
1848	if abs(ODFln[term][1]) > 1.e-3:
1849	l,m,n = eval(term.strip('C'))
1850	if SGData['SGLaue'] in ['m3','m3m']:
1851	Kcl = 0.0
1852	for j in range(0,l+1,4):
1853	im = j/4
1854	pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1855	Kcl += BOH['L=%d'%(l)][n-1][im]pcrscosd(j*beta)
1856	else: #all but cubic
1857	pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1858	pcrs *= RSQPI
1859	if n == 0:
1860	pcrs /= SQ2
1861	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1862	if SGData['SGLaue'] in ['3mR','3m1','31m']:
1863	if n%6 == 3:
1864	Kcl = pcrssind(nbeta)
1865	else:
1866	Kcl = pcrscosd(nbeta)
1867	else:
1868	Kcl = pcrscosd(nbeta)
1869	else:
1870	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
1871	invPolVal += ODFln[term][1]*Kcl
1872	return invPolVal
1873
1874
1875	def textureIndex(SHCoef):
1876	'needs doc string'
1877	Tindx = 1.0
1878	for term in SHCoef:
1879	l = eval(term.strip('C'))[0]
1880	Tindx += SHCoef[term]*2/(2.0l+1.)
1881	return Tindx
1882
1883	# self-test materials follow.
1884	selftestlist = []
1885	'''Defines a list of self-tests'''
1886	selftestquiet = True
1887	def _ReportTest():
1888	'Report name and doc string of current routine when ``selftestquiet`` is False'
1889	if not selftestquiet:
1890	import inspect
1891	caller = inspect.stack()[1][3]
1892	doc = eval(caller).__doc__
1893	if doc is not None:
1894	print('testing '+__file__+' with '+caller+' ('+doc+')')
1895	else:
1896	print('testing '+__file__()+" with "+caller)
1897	NeedTestData = True
1898	def TestData():
1899	array = np.array
1900	global NeedTestData
1901	NeedTestData = False
1902	global CellTestData
1903	# output from uctbx computed on platform darwin on 2010-05-28
1904	CellTestData = [
1905	# cell, g, G, cell, V, V
1906	[(4, 4, 4, 90, 90, 90),
1907	array([[ 1.60000000e+01, 9.79717439e-16, 9.79717439e-16],
1908	[ 9.79717439e-16, 1.60000000e+01, 9.79717439e-16],
1909	[ 9.79717439e-16, 9.79717439e-16, 1.60000000e+01]]), array([[ 6.25000000e-02, 3.82702125e-18, 3.82702125e-18],
1910	[ 3.82702125e-18, 6.25000000e-02, 3.82702125e-18],
1911	[ 3.82702125e-18, 3.82702125e-18, 6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625],
1912	# cell, g, G, cell, V, V
1913	[(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130),
1914	array([[ 16.81 , -13.70423184, 4.48533243],
1915	[-13.70423184, 27.04 , -5.6887143 ],
1916	[ 4.48533243, -5.6887143 , 39.69 ]]), array([[ 0.10206349, 0.05083339, -0.00424823],
1917	[ 0.05083339, 0.06344997, 0.00334956],
1918	[-0.00424823, 0.00334956, 0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
1919	# cell, g, G, cell, V, V
1920	[(3.5, 3.5, 6, 90, 90, 120),
1921	array([[ 1.22500000e+01, -6.12500000e+00, 1.28587914e-15],
1922	[ -6.12500000e+00, 1.22500000e+01, 1.28587914e-15],
1923	[ 1.28587914e-15, 1.28587914e-15, 3.60000000e+01]]), array([[ 1.08843537e-01, 5.44217687e-02, 3.36690552e-18],
1924	[ 5.44217687e-02, 1.08843537e-01, 3.36690552e-18],
1925	[ 3.36690552e-18, 3.36690552e-18, 2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
1926	]
1927	global CoordTestData
1928	CoordTestData = [
1929	# cell, ((frac, ortho),...)
1930	((4,4,4,90,90,90,), [
1931	((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
1932	((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
1933	((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
1934	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
1935	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
1936	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
1937	((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
1938	]),
1939	# cell, ((frac, ortho),...)
1940	((4.1,5.2,6.3,100,80,130,), [
1941	((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
1942	((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
1943	((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
1944	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
1945	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
1946	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
1947	((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
1948	]),
1949	# cell, ((frac, ortho),...)
1950	((3.5,3.5,6,90,90,120,), [
1951	((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
1952	((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
1953	((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
1954	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
1955	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
1956	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
1957	((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
1958	]),
1959	]
1960	global LaueTestData #generated by GSAS
1961	LaueTestData = {
1962	'R 3 m':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),
1963	(1,1,3,12),(1,0,4,6),(2,1,1,12),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),(3,0,-3,6),(3,0,3,6),
1964	(0,0,6,2),(2,2,0,6),(2,1,4,12),(2,0,5,6),(3,1,-1,12),(3,1,2,12),(1,1,6,12),(2,2,3,12),(2,1,-5,12))],
1965	'R 3':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),(1,1,3,6),
1966	(1,1,-3,6),(1,0,4,6),(3,-1,1,6),(2,1,1,6),(3,-1,-2,6),(2,1,-2,6),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),
1967	(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(3,-1,4,6),(2,0,5,6),(2,1,4,6),(4,-1,-1,6),(3,1,-1,6),
1968	(3,1,2,6),(4,-1,2,6),(2,2,-3,6),(1,1,-6,6),(1,1,6,6),(2,2,3,6),(2,1,-5,6),(3,-1,-5,6))],
1969	'P 3':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,1,6),(0,0,2,2),(1,0,-1,6),(1,0,2,6),(1,0,-2,6),
1970	(1,1,0,6),(0,0,3,2),(1,1,1,6),(1,1,-1,6),(1,0,3,6),(1,0,-3,6),(2,0,0,6),(2,0,-1,6),(1,1,-2,6),
1971	(1,1,2,6),(2,0,1,6),(2,0,-2,6),(2,0,2,6),(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,-4,6),(1,0,4,6),
1972	(2,0,-3,6),(2,1,0,6),(2,0,3,6),(3,-1,0,6),(2,1,1,6),(3,-1,-1,6),(2,1,-1,6),(3,-1,1,6),(1,1,4,6),
1973	(3,-1,2,6),(3,-1,-2,6),(1,1,-4,6),(0,0,5,2),(2,1,2,6),(2,1,-2,6),(3,0,0,6),(3,0,1,6),(2,0,4,6),
1974	(2,0,-4,6),(3,0,-1,6),(1,0,-5,6),(1,0,5,6),(3,-1,-3,6),(2,1,-3,6),(2,1,3,6),(3,-1,3,6),(3,0,-2,6),
1975	(3,0,2,6),(1,1,5,6),(1,1,-5,6),(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(2,0,-5,6),(2,1,-4,6),
1976	(2,2,-1,6),(3,-1,-4,6),(2,2,1,6),(3,-1,4,6),(2,1,4,6),(2,0,5,6),(1,0,-6,6),(1,0,6,6),(4,-1,0,6),
1977	(3,1,0,6),(3,1,-1,6),(3,1,1,6),(4,-1,-1,6),(2,2,2,6),(4,-1,1,6),(2,2,-2,6),(3,1,2,6),(3,1,-2,6),
1978	(3,0,4,6),(3,0,-4,6),(4,-1,-2,6),(4,-1,2,6),(2,2,-3,6),(1,1,6,6),(1,1,-6,6),(2,2,3,6),(3,-1,5,6),
1979	(2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
1980	'P 3 m 1':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,-1,6),(1,0,1,6),(0,0,2,2),(1,0,-2,6),
1981	(1,0,2,6),(1,1,0,6),(0,0,3,2),(1,1,1,12),(1,0,-3,6),(1,0,3,6),(2,0,0,6),(1,1,2,12),(2,0,1,6),
1982	(2,0,-1,6),(0,0,4,2),(2,0,-2,6),(2,0,2,6),(1,1,3,12),(1,0,-4,6),(1,0,4,6),(2,0,3,6),(2,1,0,12),
1983	(2,0,-3,6),(2,1,1,12),(2,1,-1,12),(1,1,4,12),(2,1,2,12),(0,0,5,2),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),
1984	(3,0,1,6),(3,0,-1,6),(1,0,5,6),(2,0,4,6),(2,0,-4,6),(2,1,3,12),(2,1,-3,12),(3,0,-2,6),(3,0,2,6),
1985	(1,1,5,12),(3,0,-3,6),(0,0,6,2),(2,2,0,6),(3,0,3,6),(2,1,4,12),(2,2,1,12),(2,0,5,6),(2,1,-4,12),
1986	(2,0,-5,6),(1,0,-6,6),(1,0,6,6),(3,1,0,12),(3,1,-1,12),(3,1,1,12),(2,2,2,12),(3,1,2,12),
1987	(3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
1988	'P 3 1 m':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(0,0,2,2),(1,0,1,12),(1,0,2,12),(1,1,0,6),
1989	(0,0,3,2),(1,1,-1,6),(1,1,1,6),(1,0,3,12),(2,0,0,6),(2,0,1,12),(1,1,2,6),(1,1,-2,6),(2,0,2,12),
1990	(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,4,12),(2,1,0,12),(2,0,3,12),(2,1,1,12),(2,1,-1,12),(1,1,-4,6),
1991	(1,1,4,6),(0,0,5,2),(2,1,-2,12),(2,1,2,12),(3,0,0,6),(1,0,5,12),(2,0,4,12),(3,0,1,12),(2,1,-3,12),
1992	(2,1,3,12),(3,0,2,12),(1,1,5,6),(1,1,-5,6),(3,0,3,12),(0,0,6,2),(2,2,0,6),(2,1,-4,12),(2,0,5,12),
1993	(2,2,-1,6),(2,2,1,6),(2,1,4,12),(3,1,0,12),(1,0,6,12),(2,2,2,6),(3,1,-1,12),(2,2,-2,6),(3,1,1,12),
1994	(3,1,-2,12),(3,0,4,12),(3,1,2,12),(1,1,-6,6),(2,2,3,6),(2,2,-3,6),(1,1,6,6))],
1995	}
1996
1997	global FLnhTestData
1998	FLnhTestData = [{
1999	'C(4,0,0)': (0.965, 0.42760447),
2000	'C(2,0,0)': (1.0122, -0.80233610),
2001	'C(2,0,2)': (0.0061, 8.37491546E-03),
2002	'C(6,0,4)': (-0.0898, 4.37985696E-02),
2003	'C(6,0,6)': (-0.1369, -9.04081762E-02),
2004	'C(6,0,0)': (0.5935, -0.18234928),
2005	'C(4,0,4)': (0.1872, 0.16358127),
2006	'C(6,0,2)': (0.6193, 0.27573633),
2007	'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2008	def test0():
2009	if NeedTestData: TestData()
2010	msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2011	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2012	G, g = cell2Gmat(cell)
2013	assert np.allclose(G,tG),msg
2014	assert np.allclose(g,tg),msg
2015	tcell = Gmat2cell(g)
2016	assert np.allclose(cell,tcell),msg
2017	tcell = Gmat2cell(G)
2018	assert np.allclose(tcell,trcell),msg
2019	selftestlist.append(test0)
2020
2021	def test1():
2022	'test cell2A and A2Gmat'
2023	_ReportTest()
2024	if NeedTestData: TestData()
2025	msg = 'test cell2A and A2Gmat'
2026	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2027	G, g = A2Gmat(cell2A(cell))
2028	assert np.allclose(G,tG),msg
2029	assert np.allclose(g,tg),msg
2030	selftestlist.append(test1)
2031
2032	def test2():
2033	'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2034	_ReportTest()
2035	if NeedTestData: TestData()
2036	msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2037	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2038	G, g = cell2Gmat(cell)
2039	tcell = A2cell(Gmat2A(G))
2040	assert np.allclose(cell,tcell),msg
2041	selftestlist.append(test2)
2042
2043	def test3():
2044	'test invcell2Gmat'
2045	_ReportTest()
2046	if NeedTestData: TestData()
2047	msg = 'test invcell2Gmat'
2048	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2049	G, g = invcell2Gmat(trcell)
2050	assert np.allclose(G,tG),msg
2051	assert np.allclose(g,tg),msg
2052	selftestlist.append(test3)
2053
2054	def test4():
2055	'test calc_rVsq, calc_rV, calc_V'
2056	_ReportTest()
2057	if NeedTestData: TestData()
2058	msg = 'test calc_rVsq, calc_rV, calc_V'
2059	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2060	assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2061	assert np.allclose(calc_V(cell2A(cell)),tV), msg
2062	selftestlist.append(test4)
2063
2064	def test5():
2065	'test A2invcell'
2066	_ReportTest()
2067	if NeedTestData: TestData()
2068	msg = 'test A2invcell'
2069	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2070	rcell = A2invcell(cell2A(cell))
2071	assert np.allclose(rcell,trcell),msg
2072	selftestlist.append(test5)
2073
2074	def test6():
2075	'test cell2AB'
2076	_ReportTest()
2077	if NeedTestData: TestData()
2078	msg = 'test cell2AB'
2079	for (cell,coordlist) in CoordTestData:
2080	A,B = cell2AB(cell)
2081	for (frac,ortho) in coordlist:
2082	to = np.inner(A,frac)
2083	tf = np.inner(B,to)
2084	assert np.allclose(ortho,to), msg
2085	assert np.allclose(frac,tf), msg
2086	to = np.sum(A*frac,axis=1)
2087	tf = np.sum(B*to,axis=1)
2088	assert np.allclose(ortho,to), msg
2089	assert np.allclose(frac,tf), msg
2090	selftestlist.append(test6)
2091
2092	def test7():
2093	'test GetBraviasNum(...) and GenHBravais(...)'
2094	_ReportTest()
2095	import os.path
2096	import sys
2097	import GSASIIspc as spc
2098	testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2099	if os.path.exists(testdir):
2100	if testdir not in sys.path: sys.path.insert(0,testdir)
2101	import sgtbxlattinp
2102	derror = 1e-4
2103	def indexmatch(hklin, hkllist, system):
2104	for hklref in hkllist:
2105	hklref = list(hklref)
2106	# these permutations are far from complete, but are sufficient to
2107	# allow the test to complete
2108	if system == 'cubic':
2109	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2110	elif system == 'monoclinic':
2111	permlist = [(1,2,3),(-1,2,-3)]
2112	else:
2113	permlist = [(1,2,3)]
2114
2115	for perm in permlist:
2116	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2117	if hkl == hklref: return True
2118	if [-i for i in hkl] == hklref: return True
2119	else:
2120	return False
2121
2122	for key in sgtbxlattinp.sgtbx7:
2123	spdict = spc.SpcGroup(key)
2124	cell = sgtbxlattinp.sgtbx7[key][0]
2125	system = spdict[1]['SGSys']
2126	center = spdict[1]['SGLatt']
2127
2128	bravcode = GetBraviasNum(center, system)
2129
2130	g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2131
2132	assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2133	for h,k,l,d,num in g2list:
2134	for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]:
2135	if abs(d-dref) < derror:
2136	if indexmatch((h,k,l,), hkllist, system):
2137	break
2138	else:
2139	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2140	selftestlist.append(test7)
2141
2142	def test8():
2143	'test GenHLaue'
2144	_ReportTest()
2145	import GSASIIspc as spc
2146	import sgtbxlattinp
2147	derror = 1e-4
2148	dmin = sgtbxlattinp.dmin
2149
2150	def indexmatch(hklin, hklref, system, axis):
2151	# these permutations are far from complete, but are sufficient to
2152	# allow the test to complete
2153	if system == 'cubic':
2154	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2155	elif system == 'monoclinic' and axis=='b':
2156	permlist = [(1,2,3),(-1,2,-3)]
2157	elif system == 'monoclinic' and axis=='a':
2158	permlist = [(1,2,3),(1,-2,-3)]
2159	elif system == 'monoclinic' and axis=='c':
2160	permlist = [(1,2,3),(-1,-2,3)]
2161	elif system == 'trigonal':
2162	permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2163	elif system == 'rhombohedral':
2164	permlist = [(1,2,3),(2,3,1),(3,1,2)]
2165	else:
2166	permlist = [(1,2,3)]
2167
2168	hklref = list(hklref)
2169	for perm in permlist:
2170	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2171	if hkl == hklref: return True
2172	if [-i for i in hkl] == hklref: return True
2173	return False
2174
2175	for key in sgtbxlattinp.sgtbx8:
2176	spdict = spc.SpcGroup(key)[1]
2177	cell = sgtbxlattinp.sgtbx8[key][0]
2178	center = spdict['SGLatt']
2179	Laue = spdict['SGLaue']
2180	Axis = spdict['SGUniq']
2181	system = spdict['SGSys']
2182
2183	g2list = GenHLaue(dmin,spdict,cell2A(cell))
2184	#if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2185	# print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2186	# print 'GSAS-II:'
2187	# for h,k,l,d in g2list: print ' ',(h,k,l),d
2188	# print 'SGTBX:'
2189	# for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2190	assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2191	'Reflection lists differ for %s' % key
2192	)
2193	#match = True
2194	for h,k,l,d in g2list:
2195	for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]:
2196	if abs(d-dref) < derror:
2197	if indexmatch((h,k,l,), hkllist, system, Axis): break
2198	else:
2199	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2200	#match = False
2201	#if not match:
2202	#for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2203	#print center, Laue, Axis, system
2204	selftestlist.append(test8)
2205
2206	def test9():
2207	'test GenHLaue'
2208	_ReportTest()
2209	import GSASIIspc as G2spc
2210	if NeedTestData: TestData()
2211	for spc in LaueTestData:
2212	data = LaueTestData[spc]
2213	cell = data[0]
2214	hklm = np.array(data[1])
2215	H = hklm[-1][:3]
2216	hklO = hklm.T[:3].T
2217	A = cell2A(cell)
2218	dmin = 1./np.sqrt(calc_rDsq(H,A))
2219	SGData = G2spc.SpcGroup(spc)[1]
2220	hkls = np.array(GenHLaue(dmin,SGData,A))
2221	hklN = hkls.T[:3].T
2222	#print spc,hklO.shape,hklN.shape
2223	err = True
2224	for H in hklO:
2225	if H not in hklN:
2226	print H,' missing from hkl from GSASII'
2227	err = False
2228	assert(err)
2229	selftestlist.append(test9)
2230
2231
2232
2233
2234	if __name__ == '__main__':
2235	# run self-tests
2236	selftestquiet = False
2237	for test in selftestlist:
2238	test()
2239	print "OK"

Note: See TracBrowser for help on using the repository browser.

Download in other formats: