Context Navigation

source: trunk/GSASIIlattice.py @ 2481

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

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

Download in other formats: