Context Navigation

source: trunk/GSASIIlattice.py @ 2764

Last change on this file since 2764 was 2764, checked in by toby, 6 years ago
formatting changes for docs; add debug to conf.py
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 95.3 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-03-30 01:28:04 +0000 (Thu, 30 Mar 2017) $
27	# $Author: toby $
28	# $Revision: 2764 $
29	# $URL: trunk/GSASIIlattice.py $
30	# $Id: GSASIIlattice.py 2764 2017-03-30 01:28:04Z toby $
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: 2764 $")
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 Laue: str Laue symbol
1209	centrosymmetric Laue groups::
1210
1211	['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3','-31m','-3m1',
1212	'6/m','6/mmm','m3','m3m']
1213
1214	noncentrosymmetric Laue groups::
1215
1216	['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1217	'4','-4','422','4mm','3','312','321','31m','3m1',
1218	'6','-6','622','6mm','-62m','-6m2','23','432','-43m']
1219
1220	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1221
1222	:return: HKLF new reflection array with imposed Laue symmetry
1223	'''
1224
1225	HKLFT = HKLF.T
1226	mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]]) #hkl -> k,-h,l
1227	mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]]) #hkl -> -k,h,l
1228	mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1229	mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]]) #hkl -> ihl = -h-k,h,l
1230	mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]]) #hkl -> kil = k,-h-k,l
1231	matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]]) #hkl -> k,l,h
1232	matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]]) #hkl -> -l,-h,k
1233	matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]]) #hkl -> -l,h,-k
1234	mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65
1235	matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l
1236	matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]]) #hkl -> -h-k,k,l
1237	matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l
1238	matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l
1239	matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]]) #hkl -> h,l,k
1240	mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,-h,l
1241	matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l
1242	#triclinic
1243	if Laue == '1': #ok
1244	pass
1245	elif Laue == '-1': #ok
1246	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1247	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1248	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1249	#monoclinic
1250	#noncentrosymmetric - all ok
1251	elif Laue == '2':
1252	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1253	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1254	elif Laue == '1 1 2':
1255	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1256	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1257	elif Laue == '2 1 1':
1258	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1259	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1260	elif Laue == 'm':
1261	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1262	elif Laue == 'm 1 1':
1263	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1264	elif Laue == '1 1 m':
1265	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1266	#centrosymmetric - all ok
1267	elif Laue == '2/m 1 1':
1268	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1269	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1270	HKLFT[:3] = np.where((HKLFT[2]HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]np.array([1,-1,1])[:,nxs],HKLFT[:3])
1271	elif Laue == '2/m':
1272	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1273	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1274	HKLFT[:3] = np.where((HKLFT[0]HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1275	elif Laue == '1 1 2/m':
1276	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1277	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1278	HKLFT[:3] = np.where((HKLFT[1]HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]np.array([-1,1,1])[:,nxs],HKLFT[:3])
1279	#orthorhombic
1280	#noncentrosymmetric - all OK
1281	elif Laue == '2 2 2':
1282	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1283	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1284	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1285	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1286	elif Laue == 'm m 2':
1287	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1288	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1289	elif Laue == '2 m m':
1290	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1291	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1292	elif Laue == 'm 2 m':
1293	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1294	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1295	#centrosymmetric - all ok
1296	elif Laue == 'm m m':
1297	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1298	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1299	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1300	#tetragonal
1301	#noncentrosymmetric - all ok
1302	elif Laue == '4':
1303	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1304	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1305	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1306	elif Laue == '-4':
1307	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1308	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1309	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1310	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1311	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1312	elif Laue == '4 2 2':
1313	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1314	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1315	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1316	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1317	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1318	elif Laue == '4 m m':
1319	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1320	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1321	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1322	HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1323	elif Laue == '-4 2 m':
1324	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1325	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1326	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1327	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1328	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1329	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1330	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1331	elif Laue == '-4 m 2':
1332	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1333	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1334	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1335	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1336	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1337	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1338	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1339	#centrosymmetric - all ok
1340	elif Laue == '4/m':
1341	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1342	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1343	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1344	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1345	elif Laue == '4/m m m':
1346	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1347	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1348	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1349	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1350	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1351	#trigonal - all hex cell
1352	#noncentrosymmetric - all ok
1353	elif Laue == '3':
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[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1357	elif Laue == '3 1 2':
1358	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1359	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1360	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1361	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1362	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1363	elif Laue == '3 2 1':
1364	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1365	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1366	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1367	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1368	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1369	HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2HKLFT[1]),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1370	elif Laue == '3 1 m':
1371	HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1372	HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1373	HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1374	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1375	elif Laue == '3 m 1':
1376	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1377	HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1378	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1379	#centrosymmetric
1380	elif Laue == '-3': #ok
1381	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1382	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1383	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1384	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1385	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1386	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1387	elif Laue == '-3 m 1': #ok
1388	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1389	HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1390	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1391	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1392	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1393	elif Laue == '-3 1 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[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1398	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1399	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1400	#hexagonal
1401	#noncentrosymmetric
1402	elif Laue == '6': #ok
1403	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1404	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1405	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1406	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1407	elif Laue == '-6': #ok
1408	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1409	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1410	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1411	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1412	elif Laue == '6 2 2': #ok
1413	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1414	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1415	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1416	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1417	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1418	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1419	elif Laue == '6 m m': #ok
1420	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1421	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1422	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[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]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1425	elif Laue == '-6 m 2': #ok
1426	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1427	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1428	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1429	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1430	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1431	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1432	elif Laue == '-6 2 m': #ok
1433	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1434	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1435	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1436	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1437	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1438	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1439	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1440	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1441	#centrosymmetric
1442	elif Laue == '6/m': #ok
1443	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1444	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1445	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1446	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1447	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1448	elif Laue == '6/m m m': #ok
1449	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1450	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1451	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1452	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1453	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1454	#cubic - all ok
1455	#noncentrosymmetric -
1456	elif Laue == '2 3':
1457	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1458	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1459	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1460	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1461	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])
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],matd3t[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[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1466	elif Laue == '4 3 2':
1467	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1468	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1469	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1470	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1471	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1472	HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
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[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1475	elif Laue == '-4 3 m':
1476	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1477	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1478	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1479	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1480	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1481	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1482	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1483	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])
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[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1486	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1487	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])
1488	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])
1489	#centrosymmetric
1490	elif Laue == 'm 3':
1491	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1492	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1493	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1494	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])
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	elif Laue == 'm 3 m':
1497	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1498	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1499	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1500	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])
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[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1503	return HKLFT.T
1504
1505
1506	#Spherical harmonics routines
1507	def OdfChk(SGLaue,L,M):
1508	'needs doc string'
1509	if not L%2 and abs(M) <= L:
1510	if SGLaue == '0': #cylindrical symmetry
1511	if M == 0: return True
1512	elif SGLaue == '-1':
1513	return True
1514	elif SGLaue == '2/m':
1515	if not abs(M)%2: return True
1516	elif SGLaue == 'mmm':
1517	if not abs(M)%2 and M >= 0: return True
1518	elif SGLaue == '4/m':
1519	if not abs(M)%4: return True
1520	elif SGLaue == '4/mmm':
1521	if not abs(M)%4 and M >= 0: return True
1522	elif SGLaue in ['3R','3']:
1523	if not abs(M)%3: return True
1524	elif SGLaue in ['3mR','3m1','31m']:
1525	if not abs(M)%3 and M >= 0: return True
1526	elif SGLaue == '6/m':
1527	if not abs(M)%6: return True
1528	elif SGLaue == '6/mmm':
1529	if not abs(M)%6 and M >= 0: return True
1530	elif SGLaue == 'm3':
1531	if M > 0:
1532	if L%12 == 2:
1533	if M <= L/12: return True
1534	else:
1535	if M <= L/12+1: return True
1536	elif SGLaue == 'm3m':
1537	if M > 0:
1538	if L%12 == 2:
1539	if M <= L/12: return True
1540	else:
1541	if M <= L/12+1: return True
1542	return False
1543
1544	def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1545	'needs doc string'
1546	coeffNames = []
1547	for iord in [2*i+2 for i in range(L/2)]:
1548	for m in [i-iord for i in range(2*iord+1)]:
1549	if OdfChk(SamSym,iord,m):
1550	for n in [i-iord for i in range(2*iord+1)]:
1551	if OdfChk(SGLaue,iord,n):
1552	if IfLMN:
1553	coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1554	else:
1555	coeffNames.append('C(%d,%d)'%(iord,n))
1556	return coeffNames
1557
1558	def CrsAng(H,cell,SGData):
1559	'needs doc string'
1560	a,b,c,al,be,ga = cell
1561	SQ3 = 1.732050807569
1562	H1 = np.array([1,0,0])
1563	H2 = np.array([0,1,0])
1564	H3 = np.array([0,0,1])
1565	H4 = np.array([1,1,1])
1566	G,g = cell2Gmat(cell)
1567	Laue = SGData['SGLaue']
1568	Naxis = SGData['SGUniq']
1569	if len(H.shape) == 1:
1570	DH = np.inner(H,np.inner(G,H))
1571	else:
1572	DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1573	if Laue == '2/m':
1574	if Naxis == 'a':
1575	DR = np.inner(H1,np.inner(G,H1))
1576	DHR = np.inner(H,np.inner(G,H1))
1577	elif Naxis == 'b':
1578	DR = np.inner(H2,np.inner(G,H2))
1579	DHR = np.inner(H,np.inner(G,H2))
1580	else:
1581	DR = np.inner(H3,np.inner(G,H3))
1582	DHR = np.inner(H,np.inner(G,H3))
1583	elif Laue in ['R3','R3m']:
1584	DR = np.inner(H4,np.inner(G,H4))
1585	DHR = np.inner(H,np.inner(G,H4))
1586	else:
1587	DR = np.inner(H3,np.inner(G,H3))
1588	DHR = np.inner(H,np.inner(G,H3))
1589	DHR /= np.sqrt(DR*DH)
1590	phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1591	if Laue == '-1':
1592	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1593	BB = H.T[0]sind(ga)*2
1594	elif Laue == '2/m':
1595	if Naxis == 'a':
1596	BA = H.T[2]b/(c-H.T[1]cosd(al))
1597	BB = H.T[1]sind(al)*2
1598	elif Naxis == 'b':
1599	BA = H.T[0]c/(a-H.T[2]cosd(be))
1600	BB = H.T[2]sind(be)*2
1601	else:
1602	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1603	BB = H.T[0]sind(ga)*2
1604	elif Laue in ['mmm','4/m','4/mmm']:
1605	BA = H.T[1]*a
1606	BB = H.T[0]*b
1607	elif Laue in ['3R','3mR']:
1608	BA = H.T[0]+H.T[1]-2.0*H.T[2]
1609	BB = SQ3*(H.T[0]-H.T[1])
1610	elif Laue in ['m3','m3m']:
1611	BA = H.T[1]
1612	BB = H.T[0]
1613	else:
1614	BA = H.T[0]+2.0*H.T[1]
1615	BB = SQ3*H.T[0]
1616	beta = atan2d(BA,BB)
1617	return phi,beta
1618
1619	def SamAng(Tth,Gangls,Sangl,IFCoup):
1620	"""Compute sample orientation angles vs laboratory coord. system
1621
1622	:param Tth: Signed theta
1623	:param Gangls: Sample goniometer angles phi,chi,omega,azmuth
1624	:param Sangl: Sample angle zeros om-0, chi-0, phi-0
1625	:param IFCoup: True if omega & 2-theta coupled in CW scan
1626	:returns:
1627	psi,gam: Sample odf angles
1628	dPSdA,dGMdA: Angle zero derivatives
1629	"""
1630
1631	if IFCoup:
1632	GSomeg = sind(Gangls[2]+Tth)
1633	GComeg = cosd(Gangls[2]+Tth)
1634	else:
1635	GSomeg = sind(Gangls[2])
1636	GComeg = cosd(Gangls[2])
1637	GSTth = sind(Tth)
1638	GCTth = cosd(Tth)
1639	GSazm = sind(Gangls[3])
1640	GCazm = cosd(Gangls[3])
1641	GSchi = sind(Gangls[1])
1642	GCchi = cosd(Gangls[1])
1643	GSphi = sind(Gangls[0]+Sangl[2])
1644	GCphi = cosd(Gangls[0]+Sangl[2])
1645	SSomeg = sind(Sangl[0])
1646	SComeg = cosd(Sangl[0])
1647	SSchi = sind(Sangl[1])
1648	SCchi = cosd(Sangl[1])
1649	AT = -GSTthGComeg+GCTthGCazm*GSomeg
1650	BT = GSTthGSomeg+GCTthGCazm*GComeg
1651	CT = -GCTthGSazmGSchi
1652	DT = -GCTthGSazmGCchi
1653
1654	BC1 = -ATGSphi+(CT+BTGCchi)*GCphi
1655	BC2 = DT-BT*GSchi
1656	BC3 = ATGCphi+(CT+BTGCchi)*GSphi
1657
1658	BC = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1659	psi = acosd(BC)
1660
1661	BD = 1.0-BC**2
1662	C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1663	dPSdA = [-C(-BC1SSomegSCchi-BC2SSomegSSchi-BC3SComeg),
1664	-C(-BC1SComegSSchi+BC2SComeg*SCchi),
1665	-C(-BC1SSomeg-BC3SComegSCchi)]
1666
1667	BA = -BC1SSchi+BC2SCchi
1668	BB = BC1SSomegSCchi+BC2SSomegSSchi+BC3*SComeg
1669	gam = atan2d(BB,BA)
1670
1671	BD = (BA2+BB2)/rpd
1672
1673	dBAdO = 0
1674	dBAdC = -BC1SCchi-BC2SSchi
1675	dBAdF = BC3*SSchi
1676
1677	dBBdO = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1678	dBBdC = -BC1SSomegSSchi+BC2SSomegSCchi
1679	dBBdF = BC1SComeg-BC3SSomeg*SCchi
1680
1681	dGMdA = np.where(BD > 1.e-6,[(BAdBBdO-BBdBAdO)/BD,(BAdBBdC-BBdBAdC)/BD, \
1682	(BAdBBdF-BBdBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1683
1684	return psi,gam,dPSdA,dGMdA
1685
1686	BOH = {
1687	'L=2':[[],[],[]],
1688	'L=4':[[0.30469720,0.36418281],[],[]],
1689	'L=6':[[-0.14104740,0.52775103],[],[]],
1690	'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1691	'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1692	'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1693	[0.09298802,-0.23773812,0.49446631,0.0],[]],
1694	'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1695	'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1696	[0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1697	'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1698	[-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1699	'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1700	[0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1701	'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1702	[-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1703	'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1704	[0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1705	[0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1706	'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1707	[-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1708	'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1709	[0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1710	[0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1711	'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1712	[-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1713	[-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1714	'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1715	[0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1716	[0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1717	'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1718	[-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1719	[-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1720	}
1721
1722	Lnorm = lambda L: 4.np.pi/(2.0L+1.)
1723
1724	def GetKcl(L,N,SGLaue,phi,beta):
1725	'needs doc string'
1726	import pytexture as ptx
1727	if SGLaue in ['m3','m3m']:
1728	if 'array' in str(type(phi)) and np.any(phi.shape):
1729	Kcl = np.zeros_like(phi)
1730	else:
1731	Kcl = 0.
1732	for j in range(0,L+1,4):
1733	im = j/4
1734	if 'array' in str(type(phi)) and np.any(phi.shape):
1735	pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1736	else:
1737	pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1738	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
1739	else:
1740	if 'array' in str(type(phi)) and np.any(phi.shape):
1741	pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1742	else:
1743	pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1744	pcrs *= RSQ2PI
1745	if N:
1746	pcrs *= SQ2
1747	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1748	if SGLaue in ['3mR','3m1','31m']:
1749	if N%6 == 3:
1750	Kcl = pcrssind(Nbeta)
1751	else:
1752	Kcl = pcrscosd(Nbeta)
1753	else:
1754	Kcl = pcrscosd(Nbeta)
1755	else:
1756	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
1757	return Kcl
1758
1759	def GetKsl(L,M,SamSym,psi,gam):
1760	'needs doc string'
1761	import pytexture as ptx
1762	if 'array' in str(type(psi)) and np.any(psi.shape):
1763	psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1764	else:
1765	psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1766	psrs *= RSQ2PI
1767	dpdps *= RSQ2PI
1768	if M:
1769	psrs *= SQ2
1770	dpdps *= SQ2
1771	if SamSym in ['mmm',]:
1772	dum = cosd(M*gam)
1773	Ksl = psrs*dum
1774	dKsdp = dpdps*dum
1775	dKsdg = -psrsMsind(M*gam)
1776	else:
1777	dum = cosd(Mgam)+sind(Mgam)
1778	Ksl = psrs*dum
1779	dKsdp = dpdps*dum
1780	dKsdg = psrsM(-sind(Mgam)+cosd(Mgam))
1781	return Ksl,dKsdp,dKsdg
1782
1783	def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1784	"""
1785	This is used for spherical harmonics description of preferred orientation;
1786	cylindrical symmetry only (M=0) and no sample angle derivatives returned
1787	"""
1788	import pytexture as ptx
1789	Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1790	Ksl *= RSQ2PI
1791	if SGLaue in ['m3','m3m']:
1792	Kcl = 0.0
1793	for j in range(0,L+1,4):
1794	im = j/4
1795	pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1796	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
1797	else:
1798	pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1799	pcrs *= RSQ2PI
1800	if N:
1801	pcrs *= SQ2
1802	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1803	if SGLaue in ['3mR','3m1','31m']:
1804	if N%6 == 3:
1805	Kcl = pcrssind(Nbeta)
1806	else:
1807	Kcl = pcrscosd(Nbeta)
1808	else:
1809	Kcl = pcrscosd(Nbeta)
1810	else:
1811	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
1812	return Kcl*Ksl,Lnorm(L)
1813
1814	def Glnh(Start,SHCoef,psi,gam,SamSym):
1815	'needs doc string'
1816	import pytexture as ptx
1817
1818	if Start:
1819	ptx.pyqlmninit()
1820	Start = False
1821	Fln = np.zeros(len(SHCoef))
1822	for i,term in enumerate(SHCoef):
1823	l,m,n = eval(term.strip('C'))
1824	pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1825	pcrs *= RSQPI
1826	if m == 0:
1827	pcrs /= SQ2
1828	if SamSym in ['mmm',]:
1829	Ksl = pcrscosd(mgam)
1830	else:
1831	Ksl = pcrs(cosd(mgam)+sind(m*gam))
1832	Fln[i] = SHCoef[term]KslLnorm(l)
1833	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1834	return ODFln
1835
1836	def Flnh(Start,SHCoef,phi,beta,SGData):
1837	'needs doc string'
1838	import pytexture as ptx
1839
1840	if Start:
1841	ptx.pyqlmninit()
1842	Start = False
1843	Fln = np.zeros(len(SHCoef))
1844	for i,term in enumerate(SHCoef):
1845	l,m,n = eval(term.strip('C'))
1846	if SGData['SGLaue'] in ['m3','m3m']:
1847	Kcl = 0.0
1848	for j in range(0,l+1,4):
1849	im = j/4
1850	pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1851	Kcl += BOH['L='+str(l)][n-1][im]pcrscosd(j*beta)
1852	else: #all but cubic
1853	pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1854	pcrs *= RSQPI
1855	if n == 0:
1856	pcrs /= SQ2
1857	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1858	if SGData['SGLaue'] in ['3mR','3m1','31m']:
1859	if n%6 == 3:
1860	Kcl = pcrssind(nbeta)
1861	else:
1862	Kcl = pcrscosd(nbeta)
1863	else:
1864	Kcl = pcrscosd(nbeta)
1865	else:
1866	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
1867	Fln[i] = SHCoef[term]KclLnorm(l)
1868	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1869	return ODFln
1870
1871	def polfcal(ODFln,SamSym,psi,gam):
1872	'''Perform a pole figure computation.
1873	Note that the the number of gam values must either be 1 or must
1874	match psi. Updated for numpy 1.8.0
1875	'''
1876	import pytexture as ptx
1877	PolVal = np.ones_like(psi)
1878	for term in ODFln:
1879	if abs(ODFln[term][1]) > 1.e-3:
1880	l,m,n = eval(term.strip('C'))
1881	psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1882	if SamSym in ['-1','2/m']:
1883	if m:
1884	Ksl = RSQPIpsrs(cosd(mgam)+sind(mgam))
1885	else:
1886	Ksl = RSQPI*psrs/SQ2
1887	else:
1888	if m:
1889	Ksl = RSQPIpsrscosd(m*gam)
1890	else:
1891	Ksl = RSQPI*psrs/SQ2
1892	PolVal += ODFln[term][1]*Ksl
1893	return PolVal
1894
1895	def invpolfcal(ODFln,SGData,phi,beta):
1896	'needs doc string'
1897	import pytexture as ptx
1898
1899	invPolVal = np.ones_like(beta)
1900	for term in ODFln:
1901	if abs(ODFln[term][1]) > 1.e-3:
1902	l,m,n = eval(term.strip('C'))
1903	if SGData['SGLaue'] in ['m3','m3m']:
1904	Kcl = 0.0
1905	for j in range(0,l+1,4):
1906	im = j/4
1907	pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1908	Kcl += BOH['L=%d'%(l)][n-1][im]pcrscosd(j*beta)
1909	else: #all but cubic
1910	pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1911	pcrs *= RSQPI
1912	if n == 0:
1913	pcrs /= SQ2
1914	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1915	if SGData['SGLaue'] in ['3mR','3m1','31m']:
1916	if n%6 == 3:
1917	Kcl = pcrssind(nbeta)
1918	else:
1919	Kcl = pcrscosd(nbeta)
1920	else:
1921	Kcl = pcrscosd(nbeta)
1922	else:
1923	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
1924	invPolVal += ODFln[term][1]*Kcl
1925	return invPolVal
1926
1927
1928	def textureIndex(SHCoef):
1929	'needs doc string'
1930	Tindx = 1.0
1931	for term in SHCoef:
1932	l = eval(term.strip('C'))[0]
1933	Tindx += SHCoef[term]*2/(2.0l+1.)
1934	return Tindx
1935
1936	# self-test materials follow.
1937	selftestlist = []
1938	'''Defines a list of self-tests'''
1939	selftestquiet = True
1940	def _ReportTest():
1941	'Report name and doc string of current routine when ``selftestquiet`` is False'
1942	if not selftestquiet:
1943	import inspect
1944	caller = inspect.stack()[1][3]
1945	doc = eval(caller).__doc__
1946	if doc is not None:
1947	print('testing '+__file__+' with '+caller+' ('+doc+')')
1948	else:
1949	print('testing '+__file__()+" with "+caller)
1950	NeedTestData = True
1951	def TestData():
1952	array = np.array
1953	global NeedTestData
1954	NeedTestData = False
1955	global CellTestData
1956	# output from uctbx computed on platform darwin on 2010-05-28
1957	CellTestData = [
1958	# cell, g, G, cell, V, V
1959	[(4, 4, 4, 90, 90, 90),
1960	array([[ 1.60000000e+01, 9.79717439e-16, 9.79717439e-16],
1961	[ 9.79717439e-16, 1.60000000e+01, 9.79717439e-16],
1962	[ 9.79717439e-16, 9.79717439e-16, 1.60000000e+01]]), array([[ 6.25000000e-02, 3.82702125e-18, 3.82702125e-18],
1963	[ 3.82702125e-18, 6.25000000e-02, 3.82702125e-18],
1964	[ 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],
1965	# cell, g, G, cell, V, V
1966	[(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130),
1967	array([[ 16.81 , -13.70423184, 4.48533243],
1968	[-13.70423184, 27.04 , -5.6887143 ],
1969	[ 4.48533243, -5.6887143 , 39.69 ]]), array([[ 0.10206349, 0.05083339, -0.00424823],
1970	[ 0.05083339, 0.06344997, 0.00334956],
1971	[-0.00424823, 0.00334956, 0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
1972	# cell, g, G, cell, V, V
1973	[(3.5, 3.5, 6, 90, 90, 120),
1974	array([[ 1.22500000e+01, -6.12500000e+00, 1.28587914e-15],
1975	[ -6.12500000e+00, 1.22500000e+01, 1.28587914e-15],
1976	[ 1.28587914e-15, 1.28587914e-15, 3.60000000e+01]]), array([[ 1.08843537e-01, 5.44217687e-02, 3.36690552e-18],
1977	[ 5.44217687e-02, 1.08843537e-01, 3.36690552e-18],
1978	[ 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],
1979	]
1980	global CoordTestData
1981	CoordTestData = [
1982	# cell, ((frac, ortho),...)
1983	((4,4,4,90,90,90,), [
1984	((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
1985	((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
1986	((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
1987	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
1988	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
1989	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
1990	((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
1991	]),
1992	# cell, ((frac, ortho),...)
1993	((4.1,5.2,6.3,100,80,130,), [
1994	((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
1995	((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
1996	((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
1997	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
1998	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
1999	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2000	((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2001	]),
2002	# cell, ((frac, ortho),...)
2003	((3.5,3.5,6,90,90,120,), [
2004	((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2005	((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2006	((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2007	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2008	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2009	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2010	((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2011	]),
2012	]
2013	global LaueTestData #generated by GSAS
2014	LaueTestData = {
2015	'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),
2016	(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),
2017	(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))],
2018	'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),
2019	(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),
2020	(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),
2021	(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))],
2022	'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),
2023	(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),
2024	(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),
2025	(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),
2026	(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),
2027	(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),
2028	(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),
2029	(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),
2030	(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),
2031	(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),
2032	(2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2033	'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),
2034	(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),
2035	(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),
2036	(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),
2037	(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),
2038	(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),
2039	(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),
2040	(3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2041	'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),
2042	(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),
2043	(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),
2044	(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),
2045	(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),
2046	(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),
2047	(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))],
2048	}
2049
2050	global FLnhTestData
2051	FLnhTestData = [{
2052	'C(4,0,0)': (0.965, 0.42760447),
2053	'C(2,0,0)': (1.0122, -0.80233610),
2054	'C(2,0,2)': (0.0061, 8.37491546E-03),
2055	'C(6,0,4)': (-0.0898, 4.37985696E-02),
2056	'C(6,0,6)': (-0.1369, -9.04081762E-02),
2057	'C(6,0,0)': (0.5935, -0.18234928),
2058	'C(4,0,4)': (0.1872, 0.16358127),
2059	'C(6,0,2)': (0.6193, 0.27573633),
2060	'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2061	def test0():
2062	if NeedTestData: TestData()
2063	msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2064	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2065	G, g = cell2Gmat(cell)
2066	assert np.allclose(G,tG),msg
2067	assert np.allclose(g,tg),msg
2068	tcell = Gmat2cell(g)
2069	assert np.allclose(cell,tcell),msg
2070	tcell = Gmat2cell(G)
2071	assert np.allclose(tcell,trcell),msg
2072	selftestlist.append(test0)
2073
2074	def test1():
2075	'test cell2A and A2Gmat'
2076	_ReportTest()
2077	if NeedTestData: TestData()
2078	msg = 'test cell2A and A2Gmat'
2079	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2080	G, g = A2Gmat(cell2A(cell))
2081	assert np.allclose(G,tG),msg
2082	assert np.allclose(g,tg),msg
2083	selftestlist.append(test1)
2084
2085	def test2():
2086	'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2087	_ReportTest()
2088	if NeedTestData: TestData()
2089	msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2090	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2091	G, g = cell2Gmat(cell)
2092	tcell = A2cell(Gmat2A(G))
2093	assert np.allclose(cell,tcell),msg
2094	selftestlist.append(test2)
2095
2096	def test3():
2097	'test invcell2Gmat'
2098	_ReportTest()
2099	if NeedTestData: TestData()
2100	msg = 'test invcell2Gmat'
2101	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2102	G, g = invcell2Gmat(trcell)
2103	assert np.allclose(G,tG),msg
2104	assert np.allclose(g,tg),msg
2105	selftestlist.append(test3)
2106
2107	def test4():
2108	'test calc_rVsq, calc_rV, calc_V'
2109	_ReportTest()
2110	if NeedTestData: TestData()
2111	msg = 'test calc_rVsq, calc_rV, calc_V'
2112	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2113	assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2114	assert np.allclose(calc_V(cell2A(cell)),tV), msg
2115	selftestlist.append(test4)
2116
2117	def test5():
2118	'test A2invcell'
2119	_ReportTest()
2120	if NeedTestData: TestData()
2121	msg = 'test A2invcell'
2122	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2123	rcell = A2invcell(cell2A(cell))
2124	assert np.allclose(rcell,trcell),msg
2125	selftestlist.append(test5)
2126
2127	def test6():
2128	'test cell2AB'
2129	_ReportTest()
2130	if NeedTestData: TestData()
2131	msg = 'test cell2AB'
2132	for (cell,coordlist) in CoordTestData:
2133	A,B = cell2AB(cell)
2134	for (frac,ortho) in coordlist:
2135	to = np.inner(A,frac)
2136	tf = np.inner(B,to)
2137	assert np.allclose(ortho,to), msg
2138	assert np.allclose(frac,tf), msg
2139	to = np.sum(A*frac,axis=1)
2140	tf = np.sum(B*to,axis=1)
2141	assert np.allclose(ortho,to), msg
2142	assert np.allclose(frac,tf), msg
2143	selftestlist.append(test6)
2144
2145	def test7():
2146	'test GetBraviasNum(...) and GenHBravais(...)'
2147	_ReportTest()
2148	import os.path
2149	import sys
2150	import GSASIIspc as spc
2151	testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2152	if os.path.exists(testdir):
2153	if testdir not in sys.path: sys.path.insert(0,testdir)
2154	import sgtbxlattinp
2155	derror = 1e-4
2156	def indexmatch(hklin, hkllist, system):
2157	for hklref in hkllist:
2158	hklref = list(hklref)
2159	# these permutations are far from complete, but are sufficient to
2160	# allow the test to complete
2161	if system == 'cubic':
2162	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2163	elif system == 'monoclinic':
2164	permlist = [(1,2,3),(-1,2,-3)]
2165	else:
2166	permlist = [(1,2,3)]
2167
2168	for perm in permlist:
2169	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2170	if hkl == hklref: return True
2171	if [-i for i in hkl] == hklref: return True
2172	else:
2173	return False
2174
2175	for key in sgtbxlattinp.sgtbx7:
2176	spdict = spc.SpcGroup(key)
2177	cell = sgtbxlattinp.sgtbx7[key][0]
2178	system = spdict[1]['SGSys']
2179	center = spdict[1]['SGLatt']
2180
2181	bravcode = GetBraviasNum(center, system)
2182
2183	g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2184
2185	assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2186	for h,k,l,d,num in g2list:
2187	for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]:
2188	if abs(d-dref) < derror:
2189	if indexmatch((h,k,l,), hkllist, system):
2190	break
2191	else:
2192	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2193	selftestlist.append(test7)
2194
2195	def test8():
2196	'test GenHLaue'
2197	_ReportTest()
2198	import GSASIIspc as spc
2199	import sgtbxlattinp
2200	derror = 1e-4
2201	dmin = sgtbxlattinp.dmin
2202
2203	def indexmatch(hklin, hklref, system, axis):
2204	# these permutations are far from complete, but are sufficient to
2205	# allow the test to complete
2206	if system == 'cubic':
2207	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2208	elif system == 'monoclinic' and axis=='b':
2209	permlist = [(1,2,3),(-1,2,-3)]
2210	elif system == 'monoclinic' and axis=='a':
2211	permlist = [(1,2,3),(1,-2,-3)]
2212	elif system == 'monoclinic' and axis=='c':
2213	permlist = [(1,2,3),(-1,-2,3)]
2214	elif system == 'trigonal':
2215	permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2216	elif system == 'rhombohedral':
2217	permlist = [(1,2,3),(2,3,1),(3,1,2)]
2218	else:
2219	permlist = [(1,2,3)]
2220
2221	hklref = list(hklref)
2222	for perm in permlist:
2223	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2224	if hkl == hklref: return True
2225	if [-i for i in hkl] == hklref: return True
2226	return False
2227
2228	for key in sgtbxlattinp.sgtbx8:
2229	spdict = spc.SpcGroup(key)[1]
2230	cell = sgtbxlattinp.sgtbx8[key][0]
2231	Axis = spdict['SGUniq']
2232	system = spdict['SGSys']
2233
2234	g2list = GenHLaue(dmin,spdict,cell2A(cell))
2235	#if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2236	# print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2237	# print 'GSAS-II:'
2238	# for h,k,l,d in g2list: print ' ',(h,k,l),d
2239	# print 'SGTBX:'
2240	# for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2241	assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2242	'Reflection lists differ for %s' % key
2243	)
2244	#match = True
2245	for h,k,l,d in g2list:
2246	for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]:
2247	if abs(d-dref) < derror:
2248	if indexmatch((h,k,l,), hkllist, system, Axis): break
2249	else:
2250	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2251	#match = False
2252	#if not match:
2253	#for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2254	#print center, Laue, Axis, system
2255	selftestlist.append(test8)
2256
2257	def test9():
2258	'test GenHLaue'
2259	_ReportTest()
2260	import GSASIIspc as G2spc
2261	if NeedTestData: TestData()
2262	for spc in LaueTestData:
2263	data = LaueTestData[spc]
2264	cell = data[0]
2265	hklm = np.array(data[1])
2266	H = hklm[-1][:3]
2267	hklO = hklm.T[:3].T
2268	A = cell2A(cell)
2269	dmin = 1./np.sqrt(calc_rDsq(H,A))
2270	SGData = G2spc.SpcGroup(spc)[1]
2271	hkls = np.array(GenHLaue(dmin,SGData,A))
2272	hklN = hkls.T[:3].T
2273	#print spc,hklO.shape,hklN.shape
2274	err = True
2275	for H in hklO:
2276	if H not in hklN:
2277	print H,' missing from hkl from GSASII'
2278	err = False
2279	assert(err)
2280	selftestlist.append(test9)
2281
2282
2283
2284
2285	if __name__ == '__main__':
2286	# run self-tests
2287	selftestquiet = False
2288	for test in selftestlist:
2289	test()
2290	print "OK"

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

Download in other formats: