Context Navigation

source: trunk/GSASIIlattice.py @ 4804

Last change on this file since 4804 was 4792, checked in by toby, 4 years ago
check in changes to allow cctbx to use DoIndexPeaks and !refinePeaks
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 106.9 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-02-02 17:01:51 +0000 (Tue, 02 Feb 2021) $
36	# $Author: toby $
37	# $Revision: 4792 $
38	# $URL: trunk/GSASIIlattice.py $
39	# $Id: GSASIIlattice.py 4792 2021-02-02 17:01:51Z toby $
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: 4792 $")
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] or 'PKS' in Inst['Type'][0]:
845	return TOF2dsp(Inst,pos)
846	else: #'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] or 'PKS' 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(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 GetBraviasNum(center,system):
1073	"""Determine the Bravais lattice number, as used in GenHBravais
1074
1075	:param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
1076	:param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
1077	'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
1078	:return: a number between 0 and 13
1079	or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
1080
1081	"""
1082	if center.upper() == 'F' and system.lower() == 'cubic':
1083	return 0
1084	elif center.upper() == 'I' and system.lower() == 'cubic':
1085	return 1
1086	elif center.upper() == 'P' and system.lower() == 'cubic':
1087	return 2
1088	elif center.upper() == 'R' and system.lower() == 'trigonal':
1089	return 3
1090	elif center.upper() == 'P' and system.lower() == 'hexagonal':
1091	return 4
1092	elif center.upper() == 'I' and system.lower() == 'tetragonal':
1093	return 5
1094	elif center.upper() == 'P' and system.lower() == 'tetragonal':
1095	return 6
1096	elif center.upper() == 'F' and system.lower() == 'orthorhombic':
1097	return 7
1098	elif center.upper() == 'I' and system.lower() == 'orthorhombic':
1099	return 8
1100	elif center.upper() == 'A' and system.lower() == 'orthorhombic':
1101	return 9
1102	elif center.upper() == 'B' and system.lower() == 'orthorhombic':
1103	return 10
1104	elif center.upper() == 'C' and system.lower() == 'orthorhombic':
1105	return 11
1106	elif center.upper() == 'P' and system.lower() == 'orthorhombic':
1107	return 12
1108	elif center.upper() == 'C' and system.lower() == 'monoclinic':
1109	return 13
1110	elif center.upper() == 'P' and system.lower() == 'monoclinic':
1111	return 14
1112	elif center.upper() == 'P' and system.lower() == 'triclinic':
1113	return 15
1114	raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
1115
1116	def _GenHBravais_cctbx(dmin, Bravais, A, sg_type, uctbx_unit_cell, miller_index_generator):
1117	'''Alternate form of :func:`GenHBravais` that uses CCTBX internals
1118	'''
1119	g_inv = np.array([[A[0], A[3]/2, A[4]/2],
1120	[A[3]/2, A[1], A[5]/2],
1121	[A[4]/2, A[5]/2, A[2]]])
1122	g = np.linalg.inv(g_inv)
1123	g_elems = (g[0][0], g[1][1], g[2][2], g[0][1], g[0][2], g[1][2])
1124	try:
1125	uc = uctbx_unit_cell(metrical_matrix=g_elems)
1126	except ValueError: # this function sometimes receives an A matrix that gives
1127	# numbers <0 in the diagonal elems of g. Not sure why.
1128	return []
1129	#if sg_type is None:
1130	# sg_type = make_sgtype(Bravais)
1131	mig = miller_index_generator(uc, sg_type, 0, dmin)
1132	result = []
1133	for h,k,l in mig:
1134	d = uc.d((h,k,l))
1135	result.append([h, k, l, d, -1])
1136	result.sort(key=lambda l: l[3], reverse=True)
1137	return result
1138
1139	def GenHBravais(dmin, Bravais, A, cctbx_args=None):
1140	"""Generate the positionally unique powder diffraction reflections
1141
1142	:param dmin: minimum d-spacing in A
1143	:param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
1144
1145	* 0 F cubic
1146	* 1 I cubic
1147	* 2 P cubic
1148	* 3 R hexagonal (trigonal not rhombohedral)
1149	* 4 P hexagonal
1150	* 5 I tetragonal
1151	* 6 P tetragonal
1152	* 7 F orthorhombic
1153	* 8 I orthorhombic
1154	* 9 A orthorhombic
1155	* 10 B orthorhombic
1156	* 11 C orthorhombic
1157	* 12 P orthorhombic
1158	* 13 I monoclinic
1159	* 14 C monoclinic
1160	* 15 P monoclinic
1161	* 16 P triclinic
1162
1163	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
1164	:param dict cctbx_args: items defined in CCTBX:
1165
1166	* 'sg_type': value from cctbx.sgtbx.space_group_type(symmorphic_sgs[ibrav])
1167	* 'uctbx_unit_cell': pointer to :meth:`cctbx.uctbx.unit_cell`
1168	* 'miller_index_generator': pointer to :meth:`cctbx.miller.index_generator`
1169
1170	:returns: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
1171
1172	"""
1173	if cctbx_args:
1174	return _GenHBravais_cctbx(dmin, Bravais, A,
1175	cctbx_args['sg_type'], cctbx_args['uctbx_unit_cell'], cctbx_args['miller_index_generator'])
1176
1177	if Bravais in [9,]:
1178	Cent = 'A'
1179	elif Bravais in [10,]:
1180	Cent = 'B'
1181	elif Bravais in [11,14]:
1182	Cent = 'C'
1183	elif Bravais in [1,5,8,13]:
1184	Cent = 'I'
1185	elif Bravais in [0,7]:
1186	Cent = 'F'
1187	elif Bravais in [3]:
1188	Cent = 'R'
1189	else:
1190	Cent = 'P'
1191	Hmax = MaxIndex(dmin,A)
1192	dminsq = 1./(dmin**2)
1193	HKL = []
1194	if Bravais == 16: #triclinic
1195	for l in range(-Hmax[2],Hmax[2]+1):
1196	for k in range(-Hmax[1],Hmax[1]+1):
1197	hmin = 0
1198	if (k < 0): hmin = 1
1199	if (k ==0 and l < 0): hmin = 1
1200	for h in range(hmin,Hmax[0]+1):
1201	H=[h,k,l]
1202	rdsq = calc_rDsq(H,A)
1203	if 0 < rdsq <= dminsq:
1204	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1205	elif Bravais in [13,14,15]: #monoclinic - b unique
1206	Hmax = SwapIndx(2,Hmax)
1207	for h in range(Hmax[0]+1):
1208	for k in range(-Hmax[1],Hmax[1]+1):
1209	lmin = 0
1210	if k < 0:lmin = 1
1211	for l in range(lmin,Hmax[2]+1):
1212	[h,k,l] = SwapIndx(-2,[h,k,l])
1213	H = []
1214	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1215	if H:
1216	rdsq = calc_rDsq(H,A)
1217	if 0 < rdsq <= dminsq:
1218	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1219	[h,k,l] = SwapIndx(2,[h,k,l])
1220	elif Bravais in [7,8,9,10,11,12]: #orthorhombic
1221	for h in range(Hmax[0]+1):
1222	for k in range(Hmax[1]+1):
1223	for l in range(Hmax[2]+1):
1224	H = []
1225	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1226	if H:
1227	rdsq = calc_rDsq(H,A)
1228	if 0 < rdsq <= dminsq:
1229	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1230	elif Bravais in [5,6]: #tetragonal
1231	for l in range(Hmax[2]+1):
1232	for k in range(Hmax[1]+1):
1233	for h in range(k,Hmax[0]+1):
1234	H = []
1235	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1236	if H:
1237	rdsq = calc_rDsq(H,A)
1238	if 0 < rdsq <= dminsq:
1239	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1240	elif Bravais in [3,4]:
1241	lmin = 0
1242	if Bravais == 3: lmin = -Hmax[2] #hexagonal/trigonal
1243	for l in range(lmin,Hmax[2]+1):
1244	for k in range(Hmax[1]+1):
1245	hmin = k
1246	if l < 0: hmin += 1
1247	for h in range(hmin,Hmax[0]+1):
1248	H = []
1249	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1250	if H:
1251	rdsq = calc_rDsq(H,A)
1252	if 0 < rdsq <= dminsq:
1253	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1254
1255	else: #cubic
1256	for l in range(Hmax[2]+1):
1257	for k in range(l,Hmax[1]+1):
1258	for h in range(k,Hmax[0]+1):
1259	H = []
1260	if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1261	if H:
1262	rdsq = calc_rDsq(H,A)
1263	if 0 < rdsq <= dminsq:
1264	HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1265	return sortHKLd(HKL,True,False)
1266
1267	def getHKLmax(dmin,SGData,A):
1268	'finds maximum allowed hkl for given A within dmin'
1269	SGLaue = SGData['SGLaue']
1270	if SGLaue in ['3R','3mR']: #Rhombohedral axes
1271	Hmax = [0,0,0]
1272	cell = A2cell(A)
1273	aHx = cell[0]math.sqrt(2.0(1.0-cosd(cell[3])))
1274	cHx = cell[0]math.sqrt(3.0(1.0+2.0*cosd(cell[3])))
1275	Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1276	Hmax[2] = int(round(cHx/dmin))
1277	#print Hmax,aHx,cHx
1278	else: # all others
1279	Hmax = MaxIndex(dmin,A)
1280	return Hmax
1281
1282	def GenHLaue(dmin,SGData,A):
1283	"""Generate the crystallographically unique powder diffraction reflections
1284	for a lattice and Bravais type
1285
1286	:param dmin: minimum d-spacing
1287	:param SGData: space group dictionary with at least
1288
1289	* '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'
1290	* 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1291	* 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1292
1293	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
1294	:return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1295	part of reciprocal space ignoring anomalous dispersion
1296
1297	"""
1298	import math
1299	SGLaue = SGData['SGLaue']
1300	SGLatt = SGData['SGLatt']
1301	SGUniq = SGData['SGUniq']
1302	#finds maximum allowed hkl for given A within dmin
1303	Hmax = getHKLmax(dmin,SGData,A)
1304
1305	dminsq = 1./(dmin**2)
1306	HKL = []
1307	if SGLaue == '-1': #triclinic
1308	for l in range(-Hmax[2],Hmax[2]+1):
1309	for k in range(-Hmax[1],Hmax[1]+1):
1310	hmin = 0
1311	if (k < 0) or (k ==0 and l < 0): hmin = 1
1312	for h in range(hmin,Hmax[0]+1):
1313	H = []
1314	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1315	if H:
1316	rdsq = calc_rDsq(H,A)
1317	if 0 < rdsq <= dminsq:
1318	HKL.append([h,k,l,1./math.sqrt(rdsq)])
1319	elif SGLaue == '2/m': #monoclinic
1320	axisnum = 1 + ['a','b','c'].index(SGUniq)
1321	Hmax = SwapIndx(axisnum,Hmax)
1322	for h in range(Hmax[0]+1):
1323	for k in range(-Hmax[1],Hmax[1]+1):
1324	lmin = 0
1325	if k < 0:lmin = 1
1326	for l in range(lmin,Hmax[2]+1):
1327	[h,k,l] = SwapIndx(-axisnum,[h,k,l])
1328	H = []
1329	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1330	if H:
1331	rdsq = calc_rDsq(H,A)
1332	if 0 < rdsq <= dminsq:
1333	HKL.append([h,k,l,1./math.sqrt(rdsq)])
1334	[h,k,l] = SwapIndx(axisnum,[h,k,l])
1335	elif SGLaue in ['mmm','4/m','6/m']: #orthorhombic
1336	for l in range(Hmax[2]+1):
1337	for h in range(Hmax[0]+1):
1338	kmin = 1
1339	if SGLaue == 'mmm' or h ==0: kmin = 0
1340	for k in range(kmin,Hmax[1]+1):
1341	H = []
1342	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1343	if H:
1344	rdsq = calc_rDsq(H,A)
1345	if 0 < rdsq <= dminsq:
1346	HKL.append([h,k,l,1./math.sqrt(rdsq)])
1347	elif SGLaue in ['4/mmm','6/mmm']: #tetragonal & hexagonal
1348	for l in range(Hmax[2]+1):
1349	for h in range(Hmax[0]+1):
1350	for k in range(h+1):
1351	H = []
1352	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1353	if H:
1354	rdsq = calc_rDsq(H,A)
1355	if 0 < rdsq <= dminsq:
1356	HKL.append([h,k,l,1./math.sqrt(rdsq)])
1357	elif SGLaue in ['3m1','31m','3','3R','3mR']: #trigonals
1358	for l in range(-Hmax[2],Hmax[2]+1):
1359	hmin = 0
1360	if l < 0: hmin = 1
1361	for h in range(hmin,Hmax[0]+1):
1362	if SGLaue in ['3R','3']:
1363	kmax = h
1364	kmin = -int((h-1.)/2.)
1365	else:
1366	kmin = 0
1367	kmax = h
1368	if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1369	if SGLaue == '31m' and l < 0: kmin = 1
1370	for k in range(kmin,kmax+1):
1371	H = []
1372	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1373	if SGLaue in ['3R','3mR']:
1374	H = Hx2Rh(H)
1375	if H:
1376	rdsq = calc_rDsq(H,A)
1377	if 0 < rdsq <= dminsq:
1378	HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1379	else: #cubic
1380	for h in range(Hmax[0]+1):
1381	for k in range(h+1):
1382	lmin = 0
1383	lmax = k
1384	if SGLaue =='m3':
1385	lmax = h-1
1386	if h == k: lmax += 1
1387	for l in range(lmin,lmax+1):
1388	H = []
1389	if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1390	if H:
1391	rdsq = calc_rDsq(H,A)
1392	if 0 < rdsq <= dminsq:
1393	HKL.append([h,k,l,1./math.sqrt(rdsq)])
1394	return sortHKLd(HKL,True,True)
1395
1396	def GenPfHKLs(nMax,SGData,A):
1397	"""Generate the unique pole figure reflections for a lattice and Bravais type.
1398	Min d-spacing=1.0A & no more than nMax returned
1399
1400	:param nMax: maximum number of hkls returned
1401	:param SGData: space group dictionary with at least
1402
1403	* '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'
1404	* 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1405	* 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1406
1407	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
1408	:return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1409
1410	"""
1411	HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T #strip d-spacings
1412	N = min(nMax,len(HKL))
1413	return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]
1414
1415	def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1416	'needs a doc string'
1417	ifMag = False
1418	if 'MagSpGrp' in SGData:
1419	ifMag = True
1420	HKLs = []
1421	vec = np.array(Vec)
1422	vstar = np.sqrt(calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1423	dvec = 1./(maxH*vstar+1./dmin)
1424	HKL = GenHLaue(dvec,SGData,A)
1425	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1426	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1427	for h,k,l,d in HKL:
1428	ext = G2spc.GenHKLf([h,k,l],SGData)[0] #h,k,l must be integral values here
1429	if not ext and d >= dmin:
1430	HKLs.append([h,k,l,0,d])
1431	for dH in SSdH:
1432	if dH:
1433	DH = SSdH[dH]
1434	H = [h+DH[0],k+DH[1],l+DH[2]]
1435	d = 1./np.sqrt(calc_rDsq(H,A))
1436	if d >= dmin:
1437	HKLM = np.array([h,k,l,dH])
1438	if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
1439	HKLs.append([h,k,l,dH,d])
1440	return HKLs
1441
1442	def LaueUnique2(SGData,refList):
1443	''' Impose Laue symmetry on hkl
1444
1445	:param SGData: space group data from 'P '+Laue
1446	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1447
1448	:return: HKLF new reflection array with imposed Laue symmetry
1449	'''
1450	for ref in refList:
1451	H = ref[:3]
1452	Uniq = G2spc.GenHKLf(H,SGData)[2]
1453	Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1454	ref[:3] = Uniq[-1]
1455	return refList
1456
1457	def LaueUnique(Laue,HKLF):
1458	''' Impose Laue symmetry on hkl
1459
1460	:param str Laue: Laue symbol, as below
1461
1462	centrosymmetric Laue groups::
1463
1464	['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
1465	'-31m','-3m1','6/m','6/mmm','m3','m3m']
1466
1467	noncentrosymmetric Laue groups::
1468
1469	['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1470	'4','-4','422','4mm','3','312','321','31m','3m1','6','-6',
1471	'622','6mm','-62m','-6m2','23','432','-43m']
1472
1473	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1474
1475	:returns: HKLF new reflection array with imposed Laue symmetry
1476	'''
1477
1478	HKLFT = HKLF.T
1479	mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]]) #hkl -> k,-h,l
1480	mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]]) #hkl -> -k,h,l
1481	mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1482	mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]]) #hkl -> ihl = -h-k,h,l
1483	mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]]) #hkl -> kil = k,-h-k,l
1484	matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]]) #hkl -> k,l,h
1485	matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]]) #hkl -> -l,-h,k
1486	matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]]) #hkl -> -l,h,-k
1487	mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65
1488	matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l
1489	matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]]) #hkl -> -h-k,k,l
1490	matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l
1491	matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l
1492	matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]]) #hkl -> h,l,k
1493	mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,-h,l
1494	matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l
1495	#triclinic
1496	if Laue == '1': #ok
1497	pass
1498	elif Laue == '-1': #ok
1499	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1500	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1501	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1502	#monoclinic
1503	#noncentrosymmetric - all ok
1504	elif Laue == '2':
1505	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1506	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1507	elif Laue == '1 1 2':
1508	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1509	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1510	elif Laue == '2 1 1':
1511	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1512	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1513	elif Laue == 'm':
1514	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1515	elif Laue == 'm 1 1':
1516	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1517	elif Laue == '1 1 m':
1518	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1519	#centrosymmetric - all ok
1520	elif Laue == '2/m 1 1':
1521	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1522	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1523	HKLFT[:3] = np.where((HKLFT[2]HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]np.array([1,-1,1])[:,nxs],HKLFT[:3])
1524	elif Laue == '2/m':
1525	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1526	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1527	HKLFT[:3] = np.where((HKLFT[0]HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1528	elif Laue == '1 1 2/m':
1529	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1530	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1531	HKLFT[:3] = np.where((HKLFT[1]HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]np.array([-1,1,1])[:,nxs],HKLFT[:3])
1532	#orthorhombic
1533	#noncentrosymmetric - all OK
1534	elif Laue == '2 2 2':
1535	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1536	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1537	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1538	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1539	elif Laue == 'm m 2':
1540	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1541	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1542	elif Laue == '2 m m':
1543	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1544	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1545	elif Laue == 'm 2 m':
1546	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
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 == 'm m m':
1550	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1551	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1552	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1553	#tetragonal
1554	#noncentrosymmetric - all ok
1555	elif Laue == '4':
1556	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1557	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1558	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1559	elif Laue == '-4':
1560	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1561	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1562	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1563	HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1564	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1565	elif Laue == '4 2 2':
1566	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1567	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1568	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1569	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1570	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1571	elif Laue == '4 m m':
1572	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1573	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1574	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1575	HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1576	elif Laue == '-4 2 m':
1577	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1578	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
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,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1581	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1582	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1583	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1584	elif Laue == '-4 m 2':
1585	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1586	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1587	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1588	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1589	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1590	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1591	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1592	#centrosymmetric - all ok
1593	elif Laue == '4/m':
1594	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1595	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1596	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1597	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1598	elif Laue == '4/m m m':
1599	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1600	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1601	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1602	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1603	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1604	#trigonal - all hex cell
1605	#noncentrosymmetric - all ok
1606	elif Laue == '3':
1607	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1608	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1609	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1610	elif Laue == '3 1 2':
1611	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1612	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1613	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1614	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1615	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1616	elif Laue == '3 2 1':
1617	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1618	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1619	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1620	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1621	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1622	HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2HKLFT[1]),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
1623	elif Laue == '3 1 m':
1624	HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1625	HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1626	HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1627	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1628	elif Laue == '3 m 1':
1629	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1630	HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1631	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1632	#centrosymmetric
1633	elif Laue == '-3': #ok
1634	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1635	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1636	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1637	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1638	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1639	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1640	elif Laue == '-3 m 1': #ok
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[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1643	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1644	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1645	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1646	elif Laue == '-3 1 m': #ok
1647	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1648	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1649	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1650	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1651	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])
1652	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1653	#hexagonal
1654	#noncentrosymmetric
1655	elif Laue == '6': #ok
1656	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1657	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1658	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1659	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1660	elif Laue == '-6': #ok
1661	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1662	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1663	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1664	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1665	elif Laue == '6 2 2': #ok
1666	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1667	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1668	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1669	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1670	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1671	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1672	elif Laue == '6 m m': #ok
1673	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1674	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1675	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1676	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1677	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1678	elif Laue == '-6 m 2': #ok
1679	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1680	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1681	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1682	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1683	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1684	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1685	elif Laue == '-6 2 m': #ok
1686	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1687	HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1688	HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1689	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1690	HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1691	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1692	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1693	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1694	#centrosymmetric
1695	elif Laue == '6/m': #ok
1696	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1697	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1698	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1699	HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1700	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1701	elif Laue == '6/m m m': #ok
1702	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1703	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1704	HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[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.T[nxs,:,:])).T,HKLFT[:3])
1707	#cubic - all ok
1708	#noncentrosymmetric -
1709	elif Laue == '2 3':
1710	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1711	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1712	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1713	HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1714	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])
1715	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])
1716	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])
1717	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])
1718	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1719	elif Laue == '4 3 2':
1720	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1721	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1722	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1723	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1724	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
1725	HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1726	HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])\|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1727	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1728	elif Laue == '-4 3 m':
1729	HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1730	HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1731	HKLFT[:3] = np.where(HKLFT[0]<=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],mat4bar[nxs,:,:])).T,HKLFT[:3])
1733	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1734	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1735	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1736	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])
1737	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])
1738	HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1739	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1740	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])
1741	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])
1742	#centrosymmetric
1743	elif Laue == 'm 3':
1744	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1745	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1746	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1747	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])
1748	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])
1749	elif Laue == 'm 3 m':
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,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1752	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1753	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])
1754	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])
1755	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1756	return HKLFT.T
1757
1758
1759	#Spherical harmonics routines
1760	def OdfChk(SGLaue,L,M):
1761	'needs doc string'
1762	if not L%2 and abs(M) <= L:
1763	if SGLaue == '0': #cylindrical symmetry
1764	if M == 0: return True
1765	elif SGLaue == '-1':
1766	return True
1767	elif SGLaue == '2/m':
1768	if not abs(M)%2: return True
1769	elif SGLaue == 'mmm':
1770	if not abs(M)%2 and M >= 0: return True
1771	elif SGLaue == '4/m':
1772	if not abs(M)%4: return True
1773	elif SGLaue == '4/mmm':
1774	if not abs(M)%4 and M >= 0: return True
1775	elif SGLaue in ['3R','3']:
1776	if not abs(M)%3: return True
1777	elif SGLaue in ['3mR','3m1','31m']:
1778	if not abs(M)%3 and M >= 0: return True
1779	elif SGLaue == '6/m':
1780	if not abs(M)%6: return True
1781	elif SGLaue == '6/mmm':
1782	if not abs(M)%6 and M >= 0: return True
1783	elif SGLaue == 'm3':
1784	if M > 0:
1785	if L%12 == 2:
1786	if M <= L//12: return True
1787	else:
1788	if M <= L//12+1: return True
1789	elif SGLaue == 'm3m':
1790	if M > 0:
1791	if L%12 == 2:
1792	if M <= L//12: return True
1793	else:
1794	if M <= L//12+1: return True
1795	return False
1796
1797	def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1798	'needs doc string'
1799	coeffNames = []
1800	for iord in [2*i+2 for i in range(L//2)]:
1801	for m in [i-iord for i in range(2*iord+1)]:
1802	if OdfChk(SamSym,iord,m):
1803	for n in [i-iord for i in range(2*iord+1)]:
1804	if OdfChk(SGLaue,iord,n):
1805	if IfLMN:
1806	coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1807	else:
1808	coeffNames.append('C(%d,%d)'%(iord,n))
1809	return coeffNames
1810
1811	def CrsAng(H,cell,SGData):
1812	'needs doc string'
1813	a,b,c,al,be,ga = cell
1814	SQ3 = 1.732050807569
1815	H1 = np.array([1,0,0])
1816	H2 = np.array([0,1,0])
1817	H3 = np.array([0,0,1])
1818	H4 = np.array([1,1,1])
1819	G,g = cell2Gmat(cell)
1820	Laue = SGData['SGLaue']
1821	Naxis = SGData['SGUniq']
1822	if len(H.shape) == 1:
1823	DH = np.inner(H,np.inner(G,H))
1824	else:
1825	DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1826	if Laue == '2/m':
1827	if Naxis == 'a':
1828	DR = np.inner(H1,np.inner(G,H1))
1829	DHR = np.inner(H,np.inner(G,H1))
1830	elif Naxis == 'b':
1831	DR = np.inner(H2,np.inner(G,H2))
1832	DHR = np.inner(H,np.inner(G,H2))
1833	else:
1834	DR = np.inner(H3,np.inner(G,H3))
1835	DHR = np.inner(H,np.inner(G,H3))
1836	elif Laue in ['R3','R3m']:
1837	DR = np.inner(H4,np.inner(G,H4))
1838	DHR = np.inner(H,np.inner(G,H4))
1839	else:
1840	DR = np.inner(H3,np.inner(G,H3))
1841	DHR = np.inner(H,np.inner(G,H3))
1842	DHR /= np.sqrt(DR*DH)
1843	phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1844	if Laue == '-1':
1845	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1846	BB = H.T[0]sind(ga)*2
1847	elif Laue == '2/m':
1848	if Naxis == 'a':
1849	BA = H.T[2]b/(c-H.T[1]cosd(al))
1850	BB = H.T[1]sind(al)*2
1851	elif Naxis == 'b':
1852	BA = H.T[0]c/(a-H.T[2]cosd(be))
1853	BB = H.T[2]sind(be)*2
1854	else:
1855	BA = H.T[1]a/(b-H.T[0]cosd(ga))
1856	BB = H.T[0]sind(ga)*2
1857	elif Laue in ['mmm','4/m','4/mmm']:
1858	BA = H.T[1]*a
1859	BB = H.T[0]*b
1860	elif Laue in ['3R','3mR']:
1861	BA = H.T[0]+H.T[1]-2.0*H.T[2]
1862	BB = SQ3*(H.T[0]-H.T[1])
1863	elif Laue in ['m3','m3m']:
1864	BA = H.T[1]
1865	BB = H.T[0]
1866	else:
1867	BA = H.T[0]+2.0*H.T[1]
1868	BB = SQ3*H.T[0]
1869	beta = atan2d(BA,BB)
1870	return phi,beta
1871
1872	def SamAng(Tth,Gangls,Sangl,IFCoup):
1873	"""Compute sample orientation angles vs laboratory coord. system
1874
1875	:param Tth: Signed theta
1876	:param Gangls: Sample goniometer angles phi,chi,omega,azmuth
1877	:param Sangl: Sample angle zeros om-0, chi-0, phi-0
1878	:param IFCoup: True if omega & 2-theta coupled in CW scan
1879	:returns:
1880	psi,gam: Sample odf angles
1881	dPSdA,dGMdA: Angle zero derivatives
1882	"""
1883
1884	if IFCoup:
1885	GSomeg = sind(Gangls[2]+Tth)
1886	GComeg = cosd(Gangls[2]+Tth)
1887	else:
1888	GSomeg = sind(Gangls[2])
1889	GComeg = cosd(Gangls[2])
1890	GSTth = sind(Tth)
1891	GCTth = cosd(Tth)
1892	GSazm = sind(Gangls[3])
1893	GCazm = cosd(Gangls[3])
1894	GSchi = sind(Gangls[1])
1895	GCchi = cosd(Gangls[1])
1896	GSphi = sind(Gangls[0]+Sangl[2])
1897	GCphi = cosd(Gangls[0]+Sangl[2])
1898	SSomeg = sind(Sangl[0])
1899	SComeg = cosd(Sangl[0])
1900	SSchi = sind(Sangl[1])
1901	SCchi = cosd(Sangl[1])
1902	AT = -GSTthGComeg+GCTthGCazm*GSomeg
1903	BT = GSTthGSomeg+GCTthGCazm*GComeg
1904	CT = -GCTthGSazmGSchi
1905	DT = -GCTthGSazmGCchi
1906
1907	BC1 = -ATGSphi+(CT+BTGCchi)*GCphi
1908	BC2 = DT-BT*GSchi
1909	BC3 = ATGCphi+(CT+BTGCchi)*GSphi
1910
1911	BC = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1912	psi = acosd(BC)
1913
1914	BD = 1.0-BC**2
1915	C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1916	dPSdA = [-C(-BC1SSomegSCchi-BC2SSomegSSchi-BC3SComeg),
1917	-C(-BC1SComegSSchi+BC2SComeg*SCchi),
1918	-C(-BC1SSomeg-BC3SComegSCchi)]
1919
1920	BA = -BC1SSchi+BC2SCchi
1921	BB = BC1SSomegSCchi+BC2SSomegSSchi+BC3*SComeg
1922	gam = atan2d(BB,BA)
1923
1924	BD = (BA2+BB2)/rpd
1925
1926	dBAdO = 0
1927	dBAdC = -BC1SCchi-BC2SSchi
1928	dBAdF = BC3*SSchi
1929
1930	dBBdO = BC1SComegSCchi+BC2SComegSSchi-BC3*SSomeg
1931	dBBdC = -BC1SSomegSSchi+BC2SSomegSCchi
1932	dBBdF = BC1SComeg-BC3SSomeg*SCchi
1933
1934	dGMdA = np.where(BD > 1.e-6,[(BAdBBdO-BBdBAdO)/BD,(BAdBBdC-BBdBAdC)/BD, \
1935	(BAdBBdF-BBdBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1936
1937	return psi,gam,dPSdA,dGMdA
1938
1939	BOH = {
1940	'L=2':[[],[],[]],
1941	'L=4':[[0.30469720,0.36418281],[],[]],
1942	'L=6':[[-0.14104740,0.52775103],[],[]],
1943	'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1944	'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1945	'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1946	[0.09298802,-0.23773812,0.49446631,0.0],[]],
1947	'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1948	'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1949	[0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1950	'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1951	[-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1952	'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1953	[0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1954	'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1955	[-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1956	'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1957	[0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1958	[0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1959	'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1960	[-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1961	'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1962	[0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1963	[0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1964	'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1965	[-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1966	[-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1967	'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1968	[0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1969	[0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1970	'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1971	[-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1972	[-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1973	}
1974
1975	Lnorm = lambda L: 4.np.pi/(2.0L+1.)
1976
1977	def GetKcl(L,N,SGLaue,phi,beta):
1978	'needs doc string'
1979	import pytexture as ptx
1980	if SGLaue in ['m3','m3m']:
1981	if 'array' in str(type(phi)) and np.any(phi.shape):
1982	Kcl = np.zeros_like(phi)
1983	else:
1984	Kcl = 0.
1985	for j in range(0,L+1,4):
1986	im = j//4
1987	if 'array' in str(type(phi)) and np.any(phi.shape):
1988	pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1989	else:
1990	pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1991	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
1992	else:
1993	if 'array' in str(type(phi)) and np.any(phi.shape):
1994	pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1995	else:
1996	pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1997	pcrs *= RSQ2PI
1998	if N:
1999	pcrs *= SQ2
2000	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2001	if SGLaue in ['3mR','3m1','31m']:
2002	if N%6 == 3:
2003	Kcl = pcrssind(Nbeta)
2004	else:
2005	Kcl = pcrscosd(Nbeta)
2006	else:
2007	Kcl = pcrscosd(Nbeta)
2008	else:
2009	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
2010	return Kcl
2011
2012	def GetKsl(L,M,SamSym,psi,gam):
2013	'needs doc string'
2014	import pytexture as ptx
2015	if 'array' in str(type(psi)) and np.any(psi.shape):
2016	psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
2017	else:
2018	psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
2019	psrs *= RSQ2PI
2020	dpdps *= RSQ2PI
2021	if M:
2022	psrs *= SQ2
2023	dpdps *= SQ2
2024	if SamSym in ['mmm',]:
2025	dum = cosd(M*gam)
2026	Ksl = psrs*dum
2027	dKsdp = dpdps*dum
2028	dKsdg = -psrsMsind(M*gam)
2029	else:
2030	dum = cosd(Mgam)+sind(Mgam)
2031	Ksl = psrs*dum
2032	dKsdp = dpdps*dum
2033	dKsdg = psrsM(-sind(Mgam)+cosd(Mgam))
2034	return Ksl,dKsdp,dKsdg
2035
2036	def GetKclKsl(L,N,SGLaue,psi,phi,beta):
2037	"""
2038	This is used for spherical harmonics description of preferred orientation;
2039	cylindrical symmetry only (M=0) and no sample angle derivatives returned
2040	"""
2041	import pytexture as ptx
2042	Ksl,x = ptx.pyplmpsi(L,0,1,psi)
2043	Ksl *= RSQ2PI
2044	if SGLaue in ['m3','m3m']:
2045	Kcl = 0.0
2046	for j in range(0,L+1,4):
2047	im = j//4
2048	pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
2049	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
2050	else:
2051	pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
2052	pcrs *= RSQ2PI
2053	if N:
2054	pcrs *= SQ2
2055	if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2056	if SGLaue in ['3mR','3m1','31m']:
2057	if N%6 == 3:
2058	Kcl = pcrssind(Nbeta)
2059	else:
2060	Kcl = pcrscosd(Nbeta)
2061	else:
2062	Kcl = pcrscosd(Nbeta)
2063	else:
2064	Kcl = pcrs(cosd(Nbeta)+sind(N*beta))
2065	return Kcl*Ksl,Lnorm(L)
2066
2067	def Glnh(Start,SHCoef,psi,gam,SamSym):
2068	'needs doc string'
2069	import pytexture as ptx
2070
2071	if Start:
2072	ptx.pyqlmninit()
2073	Start = False
2074	Fln = np.zeros(len(SHCoef))
2075	for i,term in enumerate(SHCoef):
2076	l,m,n = eval(term.strip('C'))
2077	pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
2078	pcrs *= RSQPI
2079	if m == 0:
2080	pcrs /= SQ2
2081	if SamSym in ['mmm',]:
2082	Ksl = pcrscosd(mgam)
2083	else:
2084	Ksl = pcrs(cosd(mgam)+sind(m*gam))
2085	Fln[i] = SHCoef[term]KslLnorm(l)
2086	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2087	return ODFln
2088
2089	def Flnh(Start,SHCoef,phi,beta,SGData):
2090	'needs doc string'
2091	import pytexture as ptx
2092
2093	if Start:
2094	ptx.pyqlmninit()
2095	Start = False
2096	Fln = np.zeros(len(SHCoef))
2097	for i,term in enumerate(SHCoef):
2098	l,m,n = eval(term.strip('C'))
2099	if SGData['SGLaue'] in ['m3','m3m']:
2100	Kcl = 0.0
2101	for j in range(0,l+1,4):
2102	im = j//4
2103	pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
2104	Kcl += BOH['L='+str(l)][n-1][im]pcrscosd(j*beta)
2105	else: #all but cubic
2106	pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
2107	pcrs *= RSQPI
2108	if n == 0:
2109	pcrs /= SQ2
2110	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2111	if SGData['SGLaue'] in ['3mR','3m1','31m']:
2112	if n%6 == 3:
2113	Kcl = pcrssind(nbeta)
2114	else:
2115	Kcl = pcrscosd(nbeta)
2116	else:
2117	Kcl = pcrscosd(nbeta)
2118	else:
2119	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
2120	Fln[i] = SHCoef[term]KclLnorm(l)
2121	ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
2122	return ODFln
2123
2124	def polfcal(ODFln,SamSym,psi,gam):
2125	'''Perform a pole figure computation.
2126	Note that the the number of gam values must either be 1 or must
2127	match psi. Updated for numpy 1.8.0
2128	'''
2129	import pytexture as ptx
2130	PolVal = np.ones_like(psi)
2131	for term in ODFln:
2132	if abs(ODFln[term][1]) > 1.e-3:
2133	l,m,n = eval(term.strip('C'))
2134	psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
2135	if SamSym in ['-1','2/m']:
2136	if m:
2137	Ksl = RSQPIpsrs(cosd(mgam)+sind(mgam))
2138	else:
2139	Ksl = RSQPI*psrs/SQ2
2140	else:
2141	if m:
2142	Ksl = RSQPIpsrscosd(m*gam)
2143	else:
2144	Ksl = RSQPI*psrs/SQ2
2145	PolVal += ODFln[term][1]*Ksl
2146	return PolVal
2147
2148	def invpolfcal(ODFln,SGData,phi,beta):
2149	'needs doc string'
2150	import pytexture as ptx
2151
2152	invPolVal = np.ones_like(beta)
2153	for term in ODFln:
2154	if abs(ODFln[term][1]) > 1.e-3:
2155	l,m,n = eval(term.strip('C'))
2156	if SGData['SGLaue'] in ['m3','m3m']:
2157	Kcl = 0.0
2158	for j in range(0,l+1,4):
2159	im = j//4
2160	pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
2161	Kcl += BOH['L=%d'%(l)][n-1][im]pcrscosd(j*beta)
2162	else: #all but cubic
2163	pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
2164	pcrs *= RSQPI
2165	if n == 0:
2166	pcrs /= SQ2
2167	if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
2168	if SGData['SGLaue'] in ['3mR','3m1','31m']:
2169	if n%6 == 3:
2170	Kcl = pcrssind(nbeta)
2171	else:
2172	Kcl = pcrscosd(nbeta)
2173	else:
2174	Kcl = pcrscosd(nbeta)
2175	else:
2176	Kcl = pcrs(cosd(nbeta)+sind(n*beta))
2177	invPolVal += ODFln[term][1]*Kcl
2178	return invPolVal
2179
2180
2181	def textureIndex(SHCoef):
2182	'needs doc string'
2183	Tindx = 1.0
2184	for term in SHCoef:
2185	l = eval(term.strip('C'))[0]
2186	Tindx += SHCoef[term]*2/(2.0l+1.)
2187	return Tindx
2188
2189	# self-test materials follow.
2190	selftestlist = []
2191	'''Defines a list of self-tests'''
2192	selftestquiet = True
2193	def _ReportTest():
2194	'Report name and doc string of current routine when ``selftestquiet`` is False'
2195	if not selftestquiet:
2196	import inspect
2197	caller = inspect.stack()[1][3]
2198	doc = eval(caller).__doc__
2199	if doc is not None:
2200	print('testing '+__file__+' with '+caller+' ('+doc+')')
2201	else:
2202	print('testing '+__file__()+" with "+caller)
2203	NeedTestData = True
2204	def TestData():
2205	array = np.array
2206	global NeedTestData
2207	NeedTestData = False
2208	global CellTestData
2209	# output from uctbx computed on platform darwin on 2010-05-28
2210	CellTestData = [
2211	# cell, g, G, cell, V, V
2212	[(4, 4, 4, 90, 90, 90),
2213	array([[ 1.60000000e+01, 9.79717439e-16, 9.79717439e-16],
2214	[ 9.79717439e-16, 1.60000000e+01, 9.79717439e-16],
2215	[ 9.79717439e-16, 9.79717439e-16, 1.60000000e+01]]), array([[ 6.25000000e-02, 3.82702125e-18, 3.82702125e-18],
2216	[ 3.82702125e-18, 6.25000000e-02, 3.82702125e-18],
2217	[ 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],
2218	# cell, g, G, cell, V, V
2219	[(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130),
2220	array([[ 16.81 , -13.70423184, 4.48533243],
2221	[-13.70423184, 27.04 , -5.6887143 ],
2222	[ 4.48533243, -5.6887143 , 39.69 ]]), array([[ 0.10206349, 0.05083339, -0.00424823],
2223	[ 0.05083339, 0.06344997, 0.00334956],
2224	[-0.00424823, 0.00334956, 0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2225	# cell, g, G, cell, V, V
2226	[(3.5, 3.5, 6, 90, 90, 120),
2227	array([[ 1.22500000e+01, -6.12500000e+00, 1.28587914e-15],
2228	[ -6.12500000e+00, 1.22500000e+01, 1.28587914e-15],
2229	[ 1.28587914e-15, 1.28587914e-15, 3.60000000e+01]]), array([[ 1.08843537e-01, 5.44217687e-02, 3.36690552e-18],
2230	[ 5.44217687e-02, 1.08843537e-01, 3.36690552e-18],
2231	[ 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],
2232	]
2233	global CoordTestData
2234	CoordTestData = [
2235	# cell, ((frac, ortho),...)
2236	((4,4,4,90,90,90,), [
2237	((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2238	((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2239	((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2240	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2241	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2242	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2243	((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2244	]),
2245	# cell, ((frac, ortho),...)
2246	((4.1,5.2,6.3,100,80,130,), [
2247	((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2248	((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2249	((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2250	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2251	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2252	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2253	((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2254	]),
2255	# cell, ((frac, ortho),...)
2256	((3.5,3.5,6,90,90,120,), [
2257	((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2258	((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2259	((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2260	((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2261	((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2262	((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2263	((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2264	]),
2265	]
2266	global LaueTestData #generated by GSAS
2267	LaueTestData = {
2268	'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),
2269	(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),
2270	(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))],
2271	'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),
2272	(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),
2273	(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),
2274	(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))],
2275	'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),
2276	(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),
2277	(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),
2278	(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),
2279	(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),
2280	(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),
2281	(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),
2282	(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),
2283	(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),
2284	(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),
2285	(2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2286	'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),
2287	(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),
2288	(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),
2289	(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),
2290	(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),
2291	(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),
2292	(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),
2293	(3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2294	'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),
2295	(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),
2296	(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),
2297	(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),
2298	(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),
2299	(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),
2300	(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))],
2301	}
2302
2303	global FLnhTestData
2304	FLnhTestData = [{
2305	'C(4,0,0)': (0.965, 0.42760447),
2306	'C(2,0,0)': (1.0122, -0.80233610),
2307	'C(2,0,2)': (0.0061, 8.37491546E-03),
2308	'C(6,0,4)': (-0.0898, 4.37985696E-02),
2309	'C(6,0,6)': (-0.1369, -9.04081762E-02),
2310	'C(6,0,0)': (0.5935, -0.18234928),
2311	'C(4,0,4)': (0.1872, 0.16358127),
2312	'C(6,0,2)': (0.6193, 0.27573633),
2313	'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2314	def test0():
2315	if NeedTestData: TestData()
2316	msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2317	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2318	G, g = cell2Gmat(cell)
2319	assert np.allclose(G,tG),msg
2320	assert np.allclose(g,tg),msg
2321	tcell = Gmat2cell(g)
2322	assert np.allclose(cell,tcell),msg
2323	tcell = Gmat2cell(G)
2324	assert np.allclose(tcell,trcell),msg
2325	if __name__ == '__main__': selftestlist.append(test0)
2326
2327	def test1():
2328	'test cell2A and A2Gmat'
2329	_ReportTest()
2330	if NeedTestData: TestData()
2331	msg = 'test cell2A and A2Gmat'
2332	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2333	G, g = A2Gmat(cell2A(cell))
2334	assert np.allclose(G,tG),msg
2335	assert np.allclose(g,tg),msg
2336	if __name__ == '__main__': selftestlist.append(test1)
2337
2338	def test2():
2339	'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2340	_ReportTest()
2341	if NeedTestData: TestData()
2342	msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2343	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2344	G, g = cell2Gmat(cell)
2345	tcell = A2cell(Gmat2A(G))
2346	assert np.allclose(cell,tcell),msg
2347	if __name__ == '__main__': selftestlist.append(test2)
2348
2349	def test3():
2350	'test invcell2Gmat'
2351	_ReportTest()
2352	if NeedTestData: TestData()
2353	msg = 'test invcell2Gmat'
2354	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2355	G, g = invcell2Gmat(trcell)
2356	assert np.allclose(G,tG),msg
2357	assert np.allclose(g,tg),msg
2358	if __name__ == '__main__': selftestlist.append(test3)
2359
2360	def test4():
2361	'test calc_rVsq, calc_rV, calc_V'
2362	_ReportTest()
2363	if NeedTestData: TestData()
2364	msg = 'test calc_rVsq, calc_rV, calc_V'
2365	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2366	assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2367	assert np.allclose(calc_V(cell2A(cell)),tV), msg
2368	if __name__ == '__main__': selftestlist.append(test4)
2369
2370	def test5():
2371	'test A2invcell'
2372	_ReportTest()
2373	if NeedTestData: TestData()
2374	msg = 'test A2invcell'
2375	for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2376	rcell = A2invcell(cell2A(cell))
2377	assert np.allclose(rcell,trcell),msg
2378	if __name__ == '__main__': selftestlist.append(test5)
2379
2380	def test6():
2381	'test cell2AB'
2382	_ReportTest()
2383	if NeedTestData: TestData()
2384	msg = 'test cell2AB'
2385	for (cell,coordlist) in CoordTestData:
2386	A,B = cell2AB(cell)
2387	for (frac,ortho) in coordlist:
2388	to = np.inner(A,frac)
2389	tf = np.inner(B,to)
2390	assert np.allclose(ortho,to), msg
2391	assert np.allclose(frac,tf), msg
2392	to = np.sum(A*frac,axis=1)
2393	tf = np.sum(B*to,axis=1)
2394	assert np.allclose(ortho,to), msg
2395	assert np.allclose(frac,tf), msg
2396	if __name__ == '__main__': selftestlist.append(test6)
2397
2398	def test7():
2399	'test GetBraviasNum(...) and GenHBravais(...)'
2400	_ReportTest()
2401	import os.path
2402	import sys
2403	import GSASIIspc as spc
2404	testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2405	if os.path.exists(testdir):
2406	if testdir not in sys.path: sys.path.insert(0,testdir)
2407	import sgtbxlattinp
2408	derror = 1e-4
2409	def indexmatch(hklin, hkllist, system):
2410	for hklref in hkllist:
2411	hklref = list(hklref)
2412	# these permutations are far from complete, but are sufficient to
2413	# allow the test to complete
2414	if system == 'cubic':
2415	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2416	elif system == 'monoclinic':
2417	permlist = [(1,2,3),(-1,2,-3)]
2418	else:
2419	permlist = [(1,2,3)]
2420
2421	for perm in permlist:
2422	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2423	if hkl == hklref: return True
2424	if [-i for i in hkl] == hklref: return True
2425	else:
2426	return False
2427
2428	for key in sgtbxlattinp.sgtbx7:
2429	spdict = spc.SpcGroup(key)
2430	cell = sgtbxlattinp.sgtbx7[key][0]
2431	system = spdict[1]['SGSys']
2432	center = spdict[1]['SGLatt']
2433
2434	bravcode = GetBraviasNum(center, system)
2435
2436	g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2437
2438	assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2439	for h,k,l,d,num in g2list:
2440	for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]:
2441	if abs(d-dref) < derror:
2442	if indexmatch((h,k,l,), hkllist, system):
2443	break
2444	else:
2445	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2446	if __name__ == '__main__': selftestlist.append(test7)
2447
2448	def test8():
2449	'test GenHLaue'
2450	_ReportTest()
2451	import GSASIIspc as spc
2452	import sgtbxlattinp
2453	derror = 1e-4
2454	dmin = sgtbxlattinp.dmin
2455
2456	def indexmatch(hklin, hklref, system, axis):
2457	# these permutations are far from complete, but are sufficient to
2458	# allow the test to complete
2459	if system == 'cubic':
2460	permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2461	elif system == 'monoclinic' and axis=='b':
2462	permlist = [(1,2,3),(-1,2,-3)]
2463	elif system == 'monoclinic' and axis=='a':
2464	permlist = [(1,2,3),(1,-2,-3)]
2465	elif system == 'monoclinic' and axis=='c':
2466	permlist = [(1,2,3),(-1,-2,3)]
2467	elif system == 'trigonal':
2468	permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2469	elif system == 'rhombohedral':
2470	permlist = [(1,2,3),(2,3,1),(3,1,2)]
2471	else:
2472	permlist = [(1,2,3)]
2473
2474	hklref = list(hklref)
2475	for perm in permlist:
2476	hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2477	if hkl == hklref: return True
2478	if [-i for i in hkl] == hklref: return True
2479	return False
2480
2481	for key in sgtbxlattinp.sgtbx8:
2482	spdict = spc.SpcGroup(key)[1]
2483	cell = sgtbxlattinp.sgtbx8[key][0]
2484	Axis = spdict['SGUniq']
2485	system = spdict['SGSys']
2486
2487	g2list = GenHLaue(dmin,spdict,cell2A(cell))
2488	#if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2489	# print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2490	# print 'GSAS-II:'
2491	# for h,k,l,d in g2list: print ' ',(h,k,l),d
2492	# print 'SGTBX:'
2493	# for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2494	assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2495	'Reflection lists differ for %s' % key
2496	)
2497	#match = True
2498	for h,k,l,d in g2list:
2499	for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]:
2500	if abs(d-dref) < derror:
2501	if indexmatch((h,k,l,), hkllist, system, Axis): break
2502	else:
2503	assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2504	#match = False
2505	#if not match:
2506	#for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref
2507	#print center, Laue, Axis, system
2508	if __name__ == '__main__': selftestlist.append(test8)
2509
2510	def test9():
2511	'test GenHLaue'
2512	_ReportTest()
2513	import GSASIIspc as G2spc
2514	if NeedTestData: TestData()
2515	for spc in LaueTestData:
2516	data = LaueTestData[spc]
2517	cell = data[0]
2518	hklm = np.array(data[1])
2519	H = hklm[-1][:3]
2520	hklO = hklm.T[:3].T
2521	A = cell2A(cell)
2522	dmin = 1./np.sqrt(calc_rDsq(H,A))
2523	SGData = G2spc.SpcGroup(spc)[1]
2524	hkls = np.array(GenHLaue(dmin,SGData,A))
2525	hklN = hkls.T[:3].T
2526	#print spc,hklO.shape,hklN.shape
2527	err = True
2528	for H in hklO:
2529	if H not in hklN:
2530	print ('%d %s'%(H,' missing from hkl from GSASII'))
2531	err = False
2532	assert(err)
2533	if __name__ == '__main__': selftestlist.append(test9)
2534
2535
2536
2537
2538	if __name__ == '__main__':
2539	# run self-tests
2540	selftestquiet = False
2541	for test in selftestlist:
2542	test()
2543	print ("OK")

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

Download in other formats: