Context Navigation

source: trunk/GSASIIlattice.py @ 3048

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

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

Download in other formats: