Context Navigation

source: trunk/GSASIIlattice.py @ 4434

Last change on this file since 4434 was 4434, checked in by vondreele, 5 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

Rev	Line
[762]	1	# -- coding: utf-8 --
[939]	2	'''
	3	GSASIIlattice: Unit cells
	4	---------------------------
	5
[989]	6	Perform lattice-related computations
	7
[3720]	8	Note that G is the reciprocal lattice tensor, and g is its inverse,
[989]	9	:math:`G = g^{-1}`, where
	10
	11	.. math::
	12
[3720]	13	g = \\left( \\begin{matrix}
[989]	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:
[4421]	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
[989]	22	d is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing,
[4421]	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.
[989]	33	'''
[762]	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 ###################
[3136]	41	from __future__ import division, print_function
[762]	42	import math
[4248]	43	import time
[2212]	44	import copy
	45	import sys
	46	import random as ran
[762]	47	import numpy as np
	48	import numpy.linalg as nl
[1046]	49	import GSASIIpath
[1372]	50	import GSASIImath as G2mth
[1594]	51	import GSASIIspc as G2spc
[2481]	52	import GSASIIElem as G2elem
[1046]	53	GSASIIpath.SetVersionNumber("$Revision: 4434 $")
[762]	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)
[4213]	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
[2126]	71	nxs = np.newaxis
[762]	72
	73	def sec2HMS(sec):
	74	"""Convert time in sec to H:M:S string
	75
	76	:param sec: time in seconds
[939]	77	:return: H:M:S string (to nearest 100th second)
[762]	78
	79	"""
[3136]	80	H = int(sec//3600)
	81	M = int(sec//60-H*60)
[762]	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):
[989]	137	"""Fill real & reciprocal metric tensor (G) from A.
[762]	138
[939]	139	:param A: reciprocal metric tensor elements as [G11,G22,G33,2G12,2G13,2*G23]
[949]	140	:param bool inverse: if True return both G and g; else just G
[762]	141	:return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
	142
	143	"""
[2560]	144	G = np.array([
[762]	145	[A[0], A[3]/2., A[4]/2.],
	146	[A[3]/2.,A[1], A[5]/2.],
[2560]	147	[A[4]/2.,A[5]/2., A[2]]])
[762]	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
[3103]	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
[2154]	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	'''
[3617]	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
[2154]	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):
[3406]	289	Uij = np.inner(Trans,np.inner(U6toUij(U6),Trans).T)/nl.det(Trans)
[2154]	290	return UijtoU6(Uij)
[3646]	291
	292	def ExpandCell(Atoms,atCodes,cx,Trans):
[4195]	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
[2212]	308
[4415]	309	def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag,Force=True):
[3406]	310	'''Transform atoms from oldPhase to newPhase
	311	M' is inv(M)
[3617]	312	does X' = M(X-U)+V transformation for coordinates and U' = MUM/det(M)
[3406]	313	for anisotropic thermal parameters
[2212]	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
[3406]	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
[2481]	321	:param ifMag: bool True if convert to magnetic phase;
	322	if True all nonmagnetic atoms will be removed
[2482]	323
[3654]	324	:return: newPhase dict modified G2 phase info
	325	:return: atCodes list atom transformation codes
	326
[2212]	327	'''
	328	cx,ct,cs,cia = oldPhase['General']['AtomPtrs']
[2484]	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])
[2212]	334	SGData = newPhase['General']['SGData']
	335	invTrans = nl.inv(Trans)
[4415]	336	newAtoms,atCodes = FillUnitCell(oldPhase,Force)
[3646]	337	newAtoms,atCodes = ExpandCell(newAtoms,atCodes,cx,Trans)
[2480]	338	if ifMag:
	339	cia += 3
	340	cs += 3
	341	newPhase['General']['Type'] = 'magnetic'
	342	newPhase['General']['AtomPtrs'] = [cx,ct,cs,cia]
	343	magAtoms = []
[2512]	344	magatCodes = []
[2481]	345	Landeg = 2.0
[2512]	346	for iat,atom in enumerate(newAtoms):
[2481]	347	if len(G2elem.GetMFtable([atom[ct],],[Landeg,])):
	348	magAtoms.append(atom[:cx+4]+[0.,0.,0.]+atom[cx+4:])
[2512]	349	magatCodes.append(atCodes[iat])
[2480]	350	newAtoms = magAtoms
[2512]	351	atCodes = magatCodes
[2480]	352	newPhase['Draw Atoms'] = []
[2212]	353	for atom in newAtoms:
[4415]	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
[2212]	358	if atom[cia] == 'A':
[3403]	359	atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],Trans)
[2470]	360	atom[cs:cs+2] = G2spc.SytSym(atom[cx:cx+3],SGData)[:2]
[3136]	361	atom[cia+8] = ran.randint(0,sys.maxsize)
[2484]	362	if cm:
	363	mag = np.sqrt(np.sum(np.array(atom[cm:cm+3])**2))
[2486]	364	if mag:
	365	mom = np.inner(np.array(atom[cm:cm+3]),oBmat)
[3402]	366	mom = np.inner(mom,invTrans)
[2486]	367	mom = np.inner(mom,nAmat)
	368	mom /= np.sqrt(np.sum(mom**2))
	369	atom[cm:cm+3] = mom*mag
[2212]	370	newPhase['Atoms'] = newAtoms
[4248]	371	if SGData['SpGrp'] != 'P 1':
	372	newPhase['Atoms'],atCodes = GetUnique(newPhase,atCodes)
[3092]	373	newPhase['Drawing'] = []
[3136]	374	newPhase['ranId'] = ran.randint(0,sys.maxsize)
[2512]	375	return newPhase,atCodes
[2212]	376
[3654]	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	'''
[3583]	386	abc = np.eye(3)
	387	cba = np.rot90(np.eye(3))
[3617]	388	cba[1,1] *= -1 #makes c-ba
[3621]	389	Mats = {'abc':abc,'cab':np.roll(abc,2,1),'bca':np.roll(abc,1,1),
[3591]	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'},
[3583]	392	'B':{'abc':'B','cab':'A','bca':'C','acb':'C','bac':'A','cba':'B'},
[3591]	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
[3583]	395	'b':{'abc':'b','cab':'a','bca':'c','acb':'c','bac':'a','cba':'b'},
[3591]	396	'c':{'abc':'c','cab':'b','bca':'a','acb':'b','bac':'c','cba':'a'},
[3583]	397	'S':{'abc':'S','cab':'S','bca':'S','acb':'S','bac':'S','cba':'S'},
[3636]	398	'I':{'abc':'I','cab':'I','bca':'I','acb':'I','bac':'I','cba':'I'},
[3583]	399	}
	400	Trans = Phase['Trans']
	401	Uvec = Phase['Uvec']
	402	SGData = Phase['SGData']
[3686]	403	MSG = SGData.get('MagSpGrp',SGData['SpGrp']).split(' ',1)
[3657]	404	MSG[0] += ' '
[3583]	405	bns = ''
	406	if '_' in MSG[0]:
	407	bns = MSG[0][2]
[3686]	408	spn = SGData.get('SGSpin',[])
[3583]	409	if 'ortho' in SGData['SGSys']:
[3591]	410	lattSym = G2spc.getlattSym(Trans)
[3583]	411	SpGrp = SGData['SpGrp']
[3621]	412	NTrans = np.inner(Mats[lattSym].T,Trans.T) #ok
[3686]	413	if len(spn): spn[1:4] = np.inner(np.abs(nl.inv(Mats[lattSym])),spn[1:4]) #ok
[3621]	414	SGsym = G2spc.getlattSym(nl.inv(Mats[lattSym]))
[3591]	415
	416	if lattSym != 'abc':
[3594]	417	NSG = G2spc.altSettingOrtho[SpGrp][SGsym].replace("'",'').split(' ')
[3622]	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
[3583]	424	Bns = ''
	425	if bns:
	426	Bns = BNS[bns][lattSym]
	427	NSG[0] += '_'+Bns+' '
[3686]	428	elif len(spn):
[3583]	429	for ifld in [1,2,3]:
	430	if spn[ifld] < 0:
	431	NSG[ifld] += "'"
	432	Nresult = [''.join(NSG)+' ',Bns]
[3594]	433	return Nresult,Uvec,NTrans
[3583]	434	else:
	435	return None
[3659]	436	elif 'mono' in SGData['SGSys']: # and not 'P_A' in Phase['Name']: #skip the one that doesn't work
[3654]	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')
[3659]	442	MSG[0] = MSG[0].replace('C_c','C_B').replace('P_A','P ')
[3654]	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:
[3658]	449	MSG[1] = MSG[1].replace('n','c')
	450	MSG[0] = MSG[0].replace('C ','I ')
[3654]	451	Mats = MatsA.T
	452	elif newcell[4] < 70.:
	453	if newcell[2] > newcell[0]:
	454	Mats = MatsB
	455	else:
[3658]	456	MSG[1] = MSG[1].replace('n','c')
[3657]	457	MSG[0] = MSG[0].replace('C ','I ')
[3654]	458	Mats = MatsB.T
[3658]	459	Nresult = [' '.join(MSG)+' ',bns]
[3654]	460	NTrans = np.inner(Mats,Trans.T)
	461	return Nresult,Uvec,NTrans
[3583]	462	return None
[3736]	463
[3686]	464	def makeBilbaoPhase(result,uvec,trans,ifMag=False):
[3583]	465	phase = {}
	466	phase['Name'] = result[0].strip()
	467	phase['Uvec'] = uvec
	468	phase['Trans'] = trans
	469	phase['Keep'] = False
	470	phase['Use'] = False
[3586]	471	phase['aType'] = ''
[3583]	472	SpGp = result[0].replace("'",'')
	473	SpGrp = G2spc.StandardizeSpcName(SpGp)
	474	phase['SGData'] = G2spc.SpcGroup(SpGrp)[1]
[3686]	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]
[3716]	482	if 'P_S' in BNSlatt: BNSlatt = 'P_c' #triclinic fix
[3686]	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'])
[3583]	487	return phase
	488
[4415]	489	def FillUnitCell(Phase,Force=True):
[3097]	490	Atoms = copy.deepcopy(Phase['Atoms'])
[2218]	491	atomData = []
[2512]	492	atCodes = []
[2212]	493	SGData = Phase['General']['SGData']
[2484]	494	SpnFlp = SGData.get('SpnFlp',[])
	495	Amat,Bmat = cell2AB(Phase['General']['Cell'][1:7])
[2212]	496	cx,ct,cs,cia = Phase['General']['AtomPtrs']
[2484]	497	cm = 0
	498	if Phase['General']['Type'] == 'magnetic':
	499	cm = cx+4
[2512]	500	for iat,atom in enumerate(Atoms):
[2212]	501	XYZ = np.array(atom[cx:cx+3])
[4415]	502	xyz = XYZ
	503	if Force:
	504	xyz %= 1.
[2212]	505	if atom[cia] == 'A':
	506	Uij = atom[cia+2:cia+8]
[4434]	507	result = G2spc.GenAtom(xyz,SGData,False,Uij,True)
[2212]	508	for item in result:
[4415]	509	# if item[0][2] >= .95: item[0][2] -= 1.
[2228]	510	atom[cx:cx+3] = item[0]
[2212]	511	atom[cia+2:cia+8] = item[1]
[2484]	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]
[2512]	518	atCodes.append('%d:%s'%(iat,str(item[2])))
[2212]	519	atomData.append(atom[:cia+9]) #not SS stuff
	520	else:
[2228]	521	result = G2spc.GenAtom(xyz,SGData,False,Move=True)
[2212]	522	for item in result:
[2218]	523	if item[0][2] >= .95: item[0][2] -= 1.
[2228]	524	atom[cx:cx+3] = item[0]
[2484]	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]
[2512]	531	atCodes.append('%d:%s'%(iat,str(item[1])))
[2212]	532	atomData.append(atom[:cia+9]) #not SS stuff
[2484]	533
[2512]	534	return atomData,atCodes
[2212]	535
[2512]	536	def GetUnique(Phase,atCodes):
[2218]	537
[3568]	538	def noDuplicate(xyzA,XYZ):
[3685]	539	if True in [np.allclose(xyzA%1.,xyzB%1.,atol=0.0002) for xyzB in XYZ]:
[2218]	540	return False
	541	return True
[2212]	542
[3638]	543	cx,ct = Phase['General']['AtomPtrs'][:2]
[2218]	544	SGData = Phase['General']['SGData']
	545	Atoms = Phase['Atoms']
	546	Ind = len(Atoms)
	547	newAtoms = []
[2512]	548	newAtCodes = []
[2218]	549	Indx = {}
	550	XYZ = {}
	551	for ind in range(Ind):
[2228]	552	XYZ[ind] = np.array(Atoms[ind][cx:cx+3])%1.
[2218]	553	Indx[ind] = True
	554	for ind in range(Ind):
	555	if Indx[ind]:
	556	xyz = XYZ[ind]
	557	for jnd in range(Ind):
[3568]	558	if Atoms[ind][ct-1] == Atoms[jnd][ct-1]:
[3474]	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])
[3568]	562	Indx[jnd] = noDuplicate(xyz,xyzs)
[2218]	563	Ind = []
	564	for ind in Indx:
	565	if Indx[ind]:
	566	newAtoms.append(Atoms[ind])
[2512]	567	newAtCodes.append(atCodes[ind])
	568	return newAtoms,newAtCodes
[2212]	569
[762]	570	def calc_rVsq(A):
[939]	571	"""Compute the square of the reciprocal lattice volume (1/V**2) from A'
[762]	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
[939]	599
	600	:returns: tuple of two 3x3 numpy arrays (A,B)
	601
[3000]	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)
[939]	604
[762]	605	"""
[3801]	606	# cellstar = Gmat2cell(G)
[762]	607	g = nl.inv(G)
	608	cell = Gmat2cell(g)
[3801]	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
[762]	620
[4268]	621	def cell2AB(cell,alt=False):
[762]	622	"""Computes orthogonalization matrix from unit cell constants
[939]	623
	624	:param tuple cell: a,b,c, alpha, beta, gamma (degrees)
	625	:returns: tuple of two 3x3 numpy arrays (A,B)
[762]	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))
[4268]	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
[762]	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)
[3136]	647	A[2][2] = 1./cellstar[2] # 1/c*
[762]	648	B = nl.inv(A)
	649	return A,B
	650
[2367]	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)
[2369]	660	xH = np.inner(B.T,H)
[2367]	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
[762]	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
[939]	671
	672	:param list U6: 6 terms of u11,u22,...
	673	:returns:
[762]	674	Uij - numpy [3][3] array of uij
	675	"""
[885]	676	U = np.array([
[1999]	677	[U6[0], U6[3], U6[4]],
	678	[U6[3], U6[1], U6[5]],
	679	[U6[4], U6[5], U6[2]]])
[762]	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	"""
[1999]	686	U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
[762]	687	return U6
	688
[2038]	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
[762]	701	def Uij2betaij(Uij,G):
	702	"""
[939]	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]
[762]	708	"""
	709	pass
	710
[959]	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)
[961]	725	return np.sum(E)/3.
[959]	726
[1075]	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
[762]	742	def CosSinAngle(U,V,G):
[949]	743	""" calculate sin & cos of angle between U & V in generalized coordinates
[762]	744	defined by metric tensor G
[939]	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:
[762]	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))
[3136]	780	cellGen = np.array([A//N2-1,A//N%N-1,A%N-1]).T
[762]	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):
[939]	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	'''
[762]	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
[1367]	836	def Pos2dsp(Inst,pos):
	837	''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
	838	'''
[1820]	839	if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
[1367]	840	wave = G2mth.getWave(Inst)
[1475]	841	return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
	842	else: #'T'OF - ignore difB
[1587]	843	return TOF2dsp(Inst,pos)
[1585]	844
	845	def TOF2dsp(Inst,Pos):
[1587]	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)
[2120]	852	N = 0
[1587]	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
[2120]	858	N += 1
	859	if N > 10:
	860	return dsp
[1367]	861
	862	def Dsp2pos(Inst,dsp):
	863	''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
	864	'''
[1820]	865	if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
[1367]	866	wave = G2mth.getWave(Inst)
[2430]	867	val = min(0.995,wave/(2.*dsp)) #set max at 168deg
	868	pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]
[1439]	869	else: #'T'OF
[1475]	870	pos = Inst['difC'][1]dsp+Inst['Zero'][1]+Inst['difA'][1]dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
[1443]	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
[1475]	879	pos = parmdict['difC']dsp+parmdict['difA']dsp**2+parmdict['difB']/dsp+parmdict['Zero']
[1443]	880	return pos
	881
[762]	882	def calc_rDsq(H,A):
[939]	883	'needs doc string'
[762]	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):
[939]	888	'needs doc string'
[762]	889	return np.inner(H,np.inner(G,H))
	890
[1597]	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
[762]	896	def calc_rDsqZ(H,A,Z,tth,lam):
[939]	897	'needs doc string'
[762]	898	rdsq = calc_rDsq(H,A)+Zsind(tth)2.0rpd/lam*2
	899	return rdsq
	900
[1578]	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
[1445]	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
[1578]	911	def calc_rDsqTSS(H,A,vec,Z,tof,difC):
	912	'needs doc string'
[2022]	913	rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
[1578]	914	return rdsq
[2476]	915
[3888]	916	def PlaneIntercepts(Amat,H,phase,stack):
[2476]	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]
[3888]	931	intcpt = ((phase)/360.+step-H[h]Ux[j,0]-H[k]Ux[j,1])/H[l]
[2476]	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
[1578]	951
[762]	952	def MaxIndex(dmin,A):
[939]	953	'needs doc string'
[762]	954	Hmax = [0,0,0]
	955	try:
	956	cell = A2cell(A)
	957	except:
[3346]	958	cell = [1.,1.,1.,90.,90.,90.]
[762]	959	for i in range(3):
	960	Hmax[i] = int(round(cell[i]/dmin))
	961	return Hmax
	962
[2129]	963	def transposeHKLF(transMat,Super,refList):
[2147]	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	'''
[2129]	971	newRefs = np.copy(refList)
[2147]	972	badRefs = []
[2129]	973	for H in newRefs:
[2146]	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):
[2147]	977	badRefs.append(newH)
	978	return newRefs,badRefs
[2129]	979
[1578]	980	def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
[2129]	981	'''sort reflection list on d-spacing; can sort in either order
[939]	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
[2802]	986	:return: sorted reflection list
[939]	987	'''
[762]	988	T = []
[1578]	989	N = 3
	990	if ifSS:
	991	N = 4
[762]	992	for i,H in enumerate(HKLd):
	993	if ifdup:
[1578]	994	T.append((H[N],i))
[762]	995	else:
[1578]	996	T.append(H[N])
[762]	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):
[939]	1009	'needs doc string'
[762]	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]]
[4195]	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
[762]	1025
	1026	def Rh2Hx(Rh):
[939]	1027	'needs doc string'
[762]	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):
[939]	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:
[3136]	1041	Rh[1] = itk//3
[939]	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
[762]	1048
	1049	def CentCheck(Cent,H):
[939]	1050	'needs doc string'
[762]	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)
[939]	1072	'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
[762]	1073	:return: a number between 0 and 13
[939]	1074	or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
	1075
[762]	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
[3633]	1095	elif center.upper() == 'A' and system.lower() == 'orthorhombic':
	1096	return 9
	1097	elif center.upper() == 'B' and system.lower() == 'orthorhombic':
	1098	return 10
[762]	1099	elif center.upper() == 'C' and system.lower() == 'orthorhombic':
[3633]	1100	return 11
[762]	1101	elif center.upper() == 'P' and system.lower() == 'orthorhombic':
[3633]	1102	return 12
[762]	1103	elif center.upper() == 'C' and system.lower() == 'monoclinic':
[3633]	1104	return 13
[762]	1105	elif center.upper() == 'P' and system.lower() == 'monoclinic':
[3633]	1106	return 14
[762]	1107	elif center.upper() == 'P' and system.lower() == 'triclinic':
[3633]	1108	return 15
[3136]	1109	raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
[762]	1110
	1111	def GenHBravais(dmin,Bravais,A):
	1112	"""Generate the positionally unique powder diffraction reflections
	1113
	1114	:param dmin: minimum d-spacing in A
[2802]	1115	:param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
[2764]	1116
[2802]	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
[3633]	1126	* 9 A orthorhombic
	1127	* 10 B orthorhombic
	1128	* 11 C orthorhombic
	1129	* 12 P orthorhombic
[3783]	1130	* 13 I monoclinic
	1131	* 14 C monoclinic
	1132	* 15 P monoclinic
	1133	* 16 P triclinic
[939]	1134
[762]	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	"""
[3633]	1139	if Bravais in [9,]:
	1140	Cent = 'A'
	1141	elif Bravais in [10,]:
	1142	Cent = 'B'
[3783]	1143	elif Bravais in [11,14]:
[762]	1144	Cent = 'C'
[3783]	1145	elif Bravais in [1,5,8,13]:
[762]	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 = []
[3783]	1156	if Bravais == 16: #triclinic
[762]	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])
[3783]	1167	elif Bravais in [13,14,15]: #monoclinic - b unique
[762]	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])
[3633]	1182	elif Bravais in [7,8,9,10,11,12]: #orthorhombic
[762]	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):
[939]	1230	'finds maximum allowed hkl for given A within dmin'
[762]	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
[939]	1249	:param SGData: space group dictionary with at least
[762]	1250
[939]	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
[762]	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]
[812]	1277	if H:
	1278	rdsq = calc_rDsq(H,A)
	1279	if 0 < rdsq <= dminsq:
[3136]	1280	HKL.append([h,k,l,1./math.sqrt(rdsq)])
[762]	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:
[3136]	1295	HKL.append([h,k,l,1./math.sqrt(rdsq)])
[762]	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:
[3136]	1308	HKL.append([h,k,l,1./math.sqrt(rdsq)])
[762]	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:
[3136]	1318	HKL.append([h,k,l,1./math.sqrt(rdsq)])
[762]	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:
[3136]	1340	HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
[762]	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:
[3136]	1355	HKL.append([h,k,l,1./math.sqrt(rdsq)])
[762]	1356	return sortHKLd(HKL,True,True)
[1773]	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]]
[762]	1376
[1594]	1377	def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
	1378	'needs a doc string'
[3774]	1379	ifMag = False
	1380	if 'MagSpGrp' in SGData:
	1381	ifMag = True
[1594]	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:
[1598]	1390	ext = G2spc.GenHKLf([h,k,l],SGData)[0] #h,k,l must be integral values here
[1594]	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]]
[3136]	1397	d = 1./np.sqrt(calc_rDsq(H,A))
[1594]	1398	if d >= dmin:
	1399	HKLM = np.array([h,k,l,dH])
[3774]	1400	if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
[1594]	1401	HKLs.append([h,k,l,dH,d])
[1596]	1402	return HKLs
[2126]	1403
[2136]	1404	def LaueUnique2(SGData,refList):
	1405	''' Impose Laue symmetry on hkl
[2764]	1406
[2136]	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
[2126]	1419	def LaueUnique(Laue,HKLF):
	1420	''' Impose Laue symmetry on hkl
[2764]	1421
[2775]	1422	:param str Laue: Laue symbol, as below
	1423
[3000]	1424	centrosymmetric Laue groups::
[2764]	1425
[2802]	1426	['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
	1427	'-31m','-3m1','6/m','6/mmm','m3','m3m']
[2764]	1428
[3000]	1429	noncentrosymmetric Laue groups::
[2764]	1430
[2802]	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']
[2764]	1434
[2126]	1435	:param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
	1436
[2802]	1437	:returns: HKLF new reflection array with imposed Laue symmetry
[2126]	1438	'''
[2136]	1439
[2126]	1440	HKLFT = HKLF.T
[2136]	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
[2140]	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
[2143]	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
[2140]	1449	mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65
[2136]	1450	matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l
[2140]	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
[2136]	1453	matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l
[2143]	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
[2140]	1456	matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l
[2126]	1457	#triclinic
[2129]	1458	if Laue == '1': #ok
[2126]	1459	pass
[2129]	1460	elif Laue == '-1': #ok
[2126]	1461	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
[2139]	1462	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
[2143]	1463	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
[2136]	1464	#monoclinic
[2143]	1465	#noncentrosymmetric - all ok
[2136]	1466	elif Laue == '2':
[2143]	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])
[2136]	1469	elif Laue == '1 1 2':
[2143]	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])
[2136]	1472	elif Laue == '2 1 1':
[2143]	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])
[2136]	1475	elif Laue == 'm':
[2129]	1476	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
[2143]	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':
[2129]	1480	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2143]	1481	#centrosymmetric - all ok
[2136]	1482	elif Laue == '2/m 1 1':
[2143]	1483	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
[2129]	1484	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
[2143]	1485	HKLFT[:3] = np.where((HKLFT[2]HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]np.array([1,-1,1])[:,nxs],HKLFT[:3])
[2136]	1486	elif Laue == '2/m':
[2126]	1487	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
[2129]	1488	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
[2143]	1489	HKLFT[:3] = np.where((HKLFT[0]HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2136]	1490	elif Laue == '1 1 2/m':
[2129]	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])
[2143]	1493	HKLFT[:3] = np.where((HKLFT[1]HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]np.array([-1,1,1])[:,nxs],HKLFT[:3])
[2136]	1494	#orthorhombic
	1495	#noncentrosymmetric - all OK
	1496	elif Laue == '2 2 2':
[2126]	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])
[2139]	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])
[2143]	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])
[2136]	1504	elif Laue == '2 m m':
[2129]	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])
[2136]	1507	elif Laue == 'm 2 m':
[2126]	1508	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
[2129]	1509	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2136]	1510	#centrosymmetric - all ok
	1511	elif Laue == 'm m m':
[2126]	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])
[2135]	1515	#tetragonal
[2143]	1516	#noncentrosymmetric - all ok
[2136]	1517	elif Laue == '4':
[2126]	1518	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2135]	1519	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
[2139]	1520	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
[2143]	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])
[2136]	1527	elif Laue == '4 2 2':
[2139]	1528	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
[2136]	1529	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2135]	1530	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
[2139]	1531	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
[2143]	1532	HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold
[2136]	1533	elif Laue == '4 m m':
[2131]	1534	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2136]	1535	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2135]	1536	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
[2136]	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':
[2143]	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])
[2139]	1551	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
[2143]	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])
[2136]	1554	#centrosymmetric - all ok
[2135]	1555	elif Laue == '4/m':
[2136]	1556	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2135]	1557	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2136]	1558	HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
[2139]	1559	HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
[2136]	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])
[2135]	1564	HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
[2136]	1565	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
[2138]	1566	#trigonal - all hex cell
[2143]	1567	#noncentrosymmetric - all ok
[2139]	1568	elif Laue == '3':
[2140]	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])
[2143]	1572	elif Laue == '3 1 2':
[2140]	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])
[2143]	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
[2144]	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])
[2143]	1585	elif Laue == '3 1 m':
[2140]	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])
[2143]	1589	HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
	1590	elif Laue == '3 m 1':
[2140]	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])
[2138]	1594	#centrosymmetric
[2140]	1595	elif Laue == '-3': #ok
[2138]	1596	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
[2140]	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])
[2143]	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])
[2140]	1606	HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
[2143]	1607	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
[2140]	1608	elif Laue == '-3 1 m': #ok
[2143]	1609	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
[2140]	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])
[2143]	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])
[2135]	1615	#hexagonal
	1616	#noncentrosymmetric
[2139]	1617	elif Laue == '6': #ok
[2140]	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])
[2139]	1622	elif Laue == '-6': #ok
	1623	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2140]	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])
[2139]	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])
[2140]	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])
[2139]	1632	HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
[2143]	1633	HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
[2139]	1634	elif Laue == '6 m m': #ok
[2140]	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])
[2139]	1639	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
[2143]	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])
[2140]	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])
[2143]	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])
[2144]	1647	elif Laue == '-6 2 m': #ok
[2143]	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
[2144]	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])
[2135]	1656	#centrosymmetric
[2139]	1657	elif Laue == '6/m': #ok
	1658	HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
[2140]	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])
[2139]	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])
[2140]	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])
[2139]	1668	HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
[2143]	1669	#cubic - all ok
	1670	#noncentrosymmetric -
	1671	elif Laue == '2 3':
[2139]	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])
[2143]	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])
[2139]	1683	HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
[2143]	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])
[2139]	1707	HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
[2143]	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])
[2126]	1718	return HKLFT.T
	1719
[1594]	1720
[762]	1721	#Spherical harmonics routines
	1722	def OdfChk(SGLaue,L,M):
[939]	1723	'needs doc string'
[762]	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:
[3136]	1748	if M <= L//12: return True
[762]	1749	else:
[3136]	1750	if M <= L//12+1: return True
[762]	1751	elif SGLaue == 'm3m':
	1752	if M > 0:
	1753	if L%12 == 2:
[3136]	1754	if M <= L//12: return True
[762]	1755	else:
[3136]	1756	if M <= L//12+1: return True
[762]	1757	return False
	1758
	1759	def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
[939]	1760	'needs doc string'
[762]	1761	coeffNames = []
[3136]	1762	for iord in [2*i+2 for i in range(L//2)]:
[762]	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
[825]	1772
[762]	1773	def CrsAng(H,cell,SGData):
[939]	1774	'needs doc string'
[762]	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']
[1792]	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])
[762]	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)
[780]	1805	phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
[762]	1806	if Laue == '-1':
[1792]	1807	BA = H.T[1]a/(b-H.T[0]cosd(ga))
	1808	BB = H.T[0]sind(ga)*2
[762]	1809	elif Laue == '2/m':
	1810	if Naxis == 'a':
[1792]	1811	BA = H.T[2]b/(c-H.T[1]cosd(al))
	1812	BB = H.T[1]sind(al)*2
[762]	1813	elif Naxis == 'b':
[1792]	1814	BA = H.T[0]c/(a-H.T[2]cosd(be))
	1815	BB = H.T[2]sind(be)*2
[762]	1816	else:
[1792]	1817	BA = H.T[1]a/(b-H.T[0]cosd(ga))
	1818	BB = H.T[0]sind(ga)*2
[762]	1819	elif Laue in ['mmm','4/m','4/mmm']:
[1792]	1820	BA = H.T[1]*a
	1821	BB = H.T[0]*b
[762]	1822	elif Laue in ['3R','3mR']:
[1792]	1823	BA = H.T[0]+H.T[1]-2.0*H.T[2]
	1824	BB = SQ3*(H.T[0]-H.T[1])
[762]	1825	elif Laue in ['m3','m3m']:
[1792]	1826	BA = H.T[1]
	1827	BB = H.T[0]
[762]	1828	else:
[1792]	1829	BA = H.T[0]+2.0*H.T[1]
	1830	BB = SQ3*H.T[0]
[762]	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
[939]	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:
[762]	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
[1792]	1877	C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
[762]	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
[1792]	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)])
[762]	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):
[939]	1940	'needs doc string'
[762]	1941	import pytexture as ptx
	1942	if SGLaue in ['m3','m3m']:
[1947]	1943	if 'array' in str(type(phi)) and np.any(phi.shape):
[1792]	1944	Kcl = np.zeros_like(phi)
	1945	else:
	1946	Kcl = 0.
[762]	1947	for j in range(0,L+1,4):
[3136]	1948	im = j//4
[1947]	1949	if 'array' in str(type(phi)) and np.any(phi.shape):
[1945]	1950	pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
[1792]	1951	else:
[1947]	1952	pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
	1953	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
[762]	1954	else:
[1947]	1955	if 'array' in str(type(phi)) and np.any(phi.shape):
[1945]	1956	pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
[1792]	1957	else:
[1947]	1958	pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
[762]	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):
[939]	1975	'needs doc string'
[762]	1976	import pytexture as ptx
[1947]	1977	if 'array' in str(type(psi)) and np.any(psi.shape):
	1978	psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
	1979	else:
[1793]	1980	psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
[1766]	1981	psrs *= RSQ2PI
	1982	dpdps *= RSQ2PI
	1983	if M:
	1984	psrs *= SQ2
	1985	dpdps *= SQ2
[762]	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):
[3136]	2009	im = j//4
[762]	2010	pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
[1947]	2011	Kcl += BOH['L=%d'%(L)][N-1][im]pcrscosd(j*beta)
[762]	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):
[939]	2030	'needs doc string'
[762]	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):
[939]	2052	'needs doc string'
[762]	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):
[3136]	2064	im = j//4
[762]	2065	pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
[1947]	2066	Kcl += BOH['L='+str(l)][n-1][im]pcrscosd(j*beta)
[762]	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):
[1452]	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	'''
[762]	2091	import pytexture as ptx
[1452]	2092	PolVal = np.ones_like(psi)
[762]	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']:
[1766]	2098	if m:
[762]	2099	Ksl = RSQPIpsrs(cosd(mgam)+sind(mgam))
	2100	else:
	2101	Ksl = RSQPI*psrs/SQ2
	2102	else:
[1766]	2103	if m:
[762]	2104	Ksl = RSQPIpsrscosd(m*gam)
	2105	else:
	2106	Ksl = RSQPI*psrs/SQ2
[1452]	2107	PolVal += ODFln[term][1]*Ksl
[762]	2108	return PolVal
	2109
	2110	def invpolfcal(ODFln,SGData,phi,beta):
[939]	2111	'needs doc string'
[762]	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):
[3136]	2121	im = j//4
[762]	2122	pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
[1947]	2123	Kcl += BOH['L=%d'%(l)][n-1][im]pcrscosd(j*beta)
[762]	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):
[939]	2144	'needs doc string'
[762]	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
[939]	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)
[762]	2165	NeedTestData = True
	2166	def TestData():
	2167	array = np.array
	2168	global NeedTestData
	2169	NeedTestData = False
	2170	global CellTestData
[939]	2171	# output from uctbx computed on platform darwin on 2010-05-28
[762]	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
[2802]	2287	if __name__ == '__main__': selftestlist.append(test0)
[762]	2288
	2289	def test1():
[939]	2290	'test cell2A and A2Gmat'
	2291	_ReportTest()
[762]	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
[2802]	2298	if __name__ == '__main__': selftestlist.append(test1)
[762]	2299
	2300	def test2():
[939]	2301	'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
	2302	_ReportTest()
[762]	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
[2802]	2309	if __name__ == '__main__': selftestlist.append(test2)
[762]	2310
	2311	def test3():
[939]	2312	'test invcell2Gmat'
	2313	_ReportTest()
[762]	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
[2802]	2320	if __name__ == '__main__': selftestlist.append(test3)
[762]	2321
	2322	def test4():
[939]	2323	'test calc_rVsq, calc_rV, calc_V'
	2324	_ReportTest()
[762]	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
[2802]	2330	if __name__ == '__main__': selftestlist.append(test4)
[762]	2331
	2332	def test5():
[939]	2333	'test A2invcell'
	2334	_ReportTest()
[762]	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
[2802]	2340	if __name__ == '__main__': selftestlist.append(test5)
[762]	2341
	2342	def test6():
[939]	2343	'test cell2AB'
	2344	_ReportTest()
[762]	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
[2802]	2358	if __name__ == '__main__': selftestlist.append(test6)
[762]	2359
	2360	def test7():
[939]	2361	'test GetBraviasNum(...) and GenHBravais(...)'
	2362	_ReportTest()
[762]	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)
[2802]	2408	if __name__ == '__main__': selftestlist.append(test7)
[762]	2409
	2410	def test8():
[939]	2411	'test GenHLaue'
	2412	_ReportTest()
[762]	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
[2802]	2470	if __name__ == '__main__': selftestlist.append(test8)
[762]	2471
	2472	def test9():
[939]	2473	'test GenHLaue'
	2474	_ReportTest()
[762]	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
[939]	2488	#print spc,hklO.shape,hklN.shape
	2489	err = True
[762]	2490	for H in hklO:
	2491	if H not in hklN:
[3136]	2492	print ('%d %s'%(H,' missing from hkl from GSASII'))
[939]	2493	err = False
	2494	assert(err)
[2802]	2495	if __name__ == '__main__': selftestlist.append(test9)
[762]	2496
	2497
	2498
	2499
	2500	if __name__ == '__main__':
[939]	2501	# run self-tests
	2502	selftestquiet = False
	2503	for test in selftestlist:
	2504	test()
[3136]	2505	print ("OK")

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

Download in other formats: