Context Navigation

source: trunk/GSASIIlattice.py @ 2825

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

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

Download in other formats: