Context Navigation

source: trunk/GSASIIlattice.py @ 4862

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

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

Download in other formats: