Context Navigation

source: trunk/GSASIIlattice.py @ 4434

Last change on this file since 4434 was 4434, checked in by vondreele, 3 years ago
Change Force=True for anisotropic atms in call to GenAtom? in FillUnitCell? - otherwise breaks RMCProfile data prep fix a 'Va' atom check in RMCProfile setup
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 105.3 KB

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

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

Download in other formats: