Context Navigation

source: trunk/GSASIImath.py @ 936

Last change on this file since 936 was 936, checked in by vondreele, 12 years ago
more MC/SA stuff further drawing optimization
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 53.0 KB

Rev	Line
[762]	1	# -- coding: utf-8 --
	2	#GSASIImath - major mathematics routines
	3	########### SVN repository information ###################
	4	# $Date: 2013-05-30 14:05:36 +0000 (Thu, 30 May 2013) $
	5	# $Author: vondreele $
	6	# $Revision: 936 $
	7	# $URL: trunk/GSASIImath.py $
	8	# $Id: GSASIImath.py 936 2013-05-30 14:05:36Z vondreele $
	9	########### SVN repository information ###################
	10	import sys
	11	import os
	12	import os.path as ospath
	13	import random as rn
	14	import numpy as np
	15	import numpy.linalg as nl
	16	import numpy.ma as ma
	17	import cPickle
	18	import time
	19	import math
	20	import copy
	21	import GSASIIpath
	22	GSASIIpath.SetVersionNumber("$Revision: 936 $")
	23	import GSASIIElem as G2el
	24	import GSASIIlattice as G2lat
	25	import GSASIIspc as G2spc
[763]	26	import numpy.fft as fft
	27
	28	sind = lambda x: np.sin(x*np.pi/180.)
	29	cosd = lambda x: np.cos(x*np.pi/180.)
	30	tand = lambda x: np.tan(x*np.pi/180.)
	31	asind = lambda x: 180.*np.arcsin(x)/np.pi
	32	acosd = lambda x: 180.*np.arccos(x)/np.pi
[823]	33	atand = lambda x: 180.*np.arctan(x)/np.pi
[763]	34	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
	35
	36	def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
	37
	38	"""
	39	Minimize the sum of squares of a set of equations.
	40
	41	::
	42
	43	Nobs
	44	x = arg min(sum(func(y)**2,axis=0))
	45	y=0
	46
	47	Parameters
	48	----------
	49	func : callable
	50	should take at least one (possibly length N vector) argument and
	51	returns M floating point numbers.
	52	x0 : ndarray
	53	The starting estimate for the minimization of length N
	54	Hess : callable
	55	A required function or method to compute the weighted vector and Hessian for func.
	56	It must be a symmetric NxN array
	57	args : tuple
	58	Any extra arguments to func are placed in this tuple.
	59	ftol : float
	60	Relative error desired in the sum of squares.
	61	xtol : float
	62	Relative error desired in the approximate solution.
	63	maxcyc : int
	64	The maximum number of cycles of refinement to execute, if -1 refine
	65	until other limits are met (ftol, xtol)
	66
	67	Returns
	68	-------
	69	x : ndarray
	70	The solution (or the result of the last iteration for an unsuccessful
	71	call).
	72	cov_x : ndarray
	73	Uses the fjac and ipvt optional outputs to construct an
	74	estimate of the jacobian around the solution. ``None`` if a
	75	singular matrix encountered (indicates very flat curvature in
	76	some direction). This matrix must be multiplied by the
	77	residual standard deviation to get the covariance of the
	78	parameter estimates -- see curve_fit.
	79	infodict : dict
	80	a dictionary of optional outputs with the key s::
	81
	82	- 'fvec' : the function evaluated at the output
	83
	84
	85	Notes
	86	-----
	87
	88	"""
	89
	90	x0 = np.array(x0, ndmin=1) #might be redundant?
	91	n = len(x0)
	92	if type(args) != type(()):
	93	args = (args,)
	94
	95	icycle = 0
	96	One = np.ones((n,n))
	97	lam = 0.001
	98	lamMax = lam
	99	nfev = 0
	100	while icycle < maxcyc:
	101	lamMax = max(lamMax,lam)
	102	M = func(x0,*args)
	103	nfev += 1
	104	chisq0 = np.sum(M**2)
	105	Yvec,Amat = Hess(x0,*args)
	106	Adiag = np.sqrt(np.diag(Amat))
	107	psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
	108	if np.any(psing): #hard singularity in matrix
	109	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
	110	Anorm = np.outer(Adiag,Adiag)
	111	Yvec /= Adiag
	112	Amat /= Anorm
	113	while True:
	114	Lam = np.eye(Amat.shape[0])*lam
	115	Amatlam = Amat*(One+Lam)
	116	try:
	117	Xvec = nl.solve(Amatlam,Yvec)
	118	except nl.LinAlgError:
	119	print 'ouch #1'
	120	psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
	121	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
	122	Xvec /= Adiag
	123	M2 = func(x0+Xvec,*args)
	124	nfev += 1
	125	chisq1 = np.sum(M2**2)
	126	if chisq1 > chisq0:
	127	lam *= 10.
	128	else:
	129	x0 += Xvec
	130	lam /= 10.
	131	break
	132	if (chisq0-chisq1)/chisq0 < ftol:
	133	break
	134	icycle += 1
	135	M = func(x0,*args)
	136	nfev += 1
	137	Yvec,Amat = Hess(x0,*args)
	138	try:
	139	Bmat = nl.inv(Amat)
	140	return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
	141	except nl.LinAlgError:
	142	print 'ouch #2 linear algebra error in LS'
	143	psing = []
	144	if maxcyc:
	145	psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
	146	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
	147
	148	def getVCov(varyNames,varyList,covMatrix):
	149	vcov = np.zeros((len(varyNames),len(varyNames)))
	150	for i1,name1 in enumerate(varyNames):
	151	for i2,name2 in enumerate(varyNames):
	152	try:
	153	vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
	154	except ValueError:
	155	vcov[i1][i2] = 0.0
	156	return vcov
	157
[808]	158	def FindAtomIndexByIDs(atomData,IDs,Draw=True):
	159	indx = []
	160	for i,atom in enumerate(atomData):
	161	if Draw and atom[-3] in IDs:
	162	indx.append(i)
	163	elif atom[-1] in IDs:
	164	indx.append(i)
	165	return indx
	166
	167	def FillAtomLookUp(atomData):
	168	atomLookUp = {}
	169	for iatm,atom in enumerate(atomData):
	170	atomLookUp[atom[-1]] = iatm
	171	return atomLookUp
	172
	173	def GetAtomsById(atomData,atomLookUp,IdList):
	174	atoms = []
	175	for id in IdList:
	176	atoms.append(atomData[atomLookUp[id]])
	177	return atoms
	178
	179	def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
	180	Items = []
[811]	181	if not isinstance(IdList,list):
	182	IdList = [IdList,]
[808]	183	for id in IdList:
	184	if numItems == 1:
	185	Items.append(atomData[atomLookUp[id]][itemLoc])
	186	else:
	187	Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
	188	return Items
[820]	189
	190	def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
	191	pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
	192	dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
	193	XYZ = []
	194	for ind in indx:
	195	names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
	196	dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
	197	XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
	198	return XYZ
[826]	199
[881]	200	def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used
[859]	201	for atom in atomData:
	202	XYZ = np.inner(Amat,atom[cx:cx+3])
	203	if atom[cia] == 'A':
	204	UIJ = atom[cia+2:cia+8]
[902]	205
[910]	206	def TLS2Uij(xyz,g,Amat,rbObj): #not used anywhere, but could be?
[902]	207	TLStype,TLS = rbObj['ThermalMotion'][:2]
	208	Tmat = np.zeros((3,3))
	209	Lmat = np.zeros((3,3))
	210	Smat = np.zeros((3,3))
	211	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
	212	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
	213	if 'T' in TLStype:
	214	Tmat = G2lat.U6toUij(TLS[:6])
	215	if 'L' in TLStype:
	216	Lmat = G2lat.U6toUij(TLS[6:12])
	217	if 'S' in TLStype:
	218	Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
	219	XYZ = np.inner(Amat,xyz)
	220	Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
	221	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
	222	beta = np.inner(np.inner(g,Umat),g)
	223	return G2lat.UijtoU6(beta)*gvec
[859]	224
[905]	225	def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj): #not used anywhere, but could be?
[859]	226	cx,ct,cs,cia = atPtrs
	227	TLStype,TLS = rbObj['ThermalMotion'][:2]
	228	Tmat = np.zeros((3,3))
	229	Lmat = np.zeros((3,3))
	230	Smat = np.zeros((3,3))
	231	G,g = G2lat.A2Gmat(Amat)
	232	gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
	233	if 'T' in TLStype:
	234	Tmat = G2lat.U6toUij(TLS[:6])
	235	if 'L' in TLStype:
	236	Lmat = G2lat.U6toUij(TLS[6:12])
	237	if 'S' in TLStype:
[902]	238	Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
[859]	239	for atom in atomData:
	240	XYZ = np.inner(Amat,atom[cx:cx+3])
[902]	241	Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
[859]	242	if 'U' in TSLtype:
	243	atom[cia+1] = TLS[0]
	244	atom[cia] = 'I'
	245	else:
	246	atom[cia] = 'A'
	247	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
	248	beta = np.inner(np.inner(g,Umat),g)
[881]	249	atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
[859]	250
[855]	251	def GetXYZDist(xyz,XYZ,Amat):
	252	''' gets distance from position xyz to all XYZ, xyz & XYZ are np.array
	253	and are in crystal coordinates; Amat is crystal to Cart matrix
	254	'''
	255	return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
	256
	257	def getAtomXYZ(atoms,cx):
	258	XYZ = []
	259	for atom in atoms:
	260	XYZ.append(atom[cx:cx+3])
	261	return np.array(XYZ)
	262
[881]	263	def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
	264	''' Returns crystal coordinates for atoms described by RBObj
	265	'''
	266	RBRes = RBData[RBType][RBObj['RBId']]
[851]	267	if RBType == 'Vector':
	268	vecs = RBRes['rbVect']
	269	mags = RBRes['VectMag']
[881]	270	Cart = np.zeros_like(vecs[0])
[851]	271	for vec,mag in zip(vecs,mags):
[881]	272	Cart += vec*mag
[851]	273	elif RBType == 'Residue':
[881]	274	Cart = np.array(RBRes['rbXYZ'])
[851]	275	for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
[881]	276	QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
[851]	277	for ride in seq[3]:
[882]	278	Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
[881]	279	XYZ = np.zeros_like(Cart)
	280	for i,xyz in enumerate(Cart):
[882]	281	X = prodQVQ(RBObj['Orient'][0],xyz)
	282	XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
[881]	283	return XYZ,Cart
[934]	284
[936]	285	def UpdateMCSAxyz(Bmat,MCSA):
[934]	286	xyz = []
	287	atTypes = []
	288	iatm = 0
[936]	289	for model in MCSA['Models'][1:]: #skip the MD model
[934]	290	if model['Type'] == 'Atom':
	291	xyz.append(model['Pos'][0])
	292	atTypes.append(model['atType'])
	293	iatm += 1
	294	else:
[936]	295	RBRes = MCSA['rbData'][model['Type']][model['RBId']]
[934]	296	Pos = np.array(model['Pos'][0])
	297	Qori = np.array(model['Ori'][0])
	298	if model['Type'] == 'Vector':
	299	vecs = RBRes['rbVect']
	300	mags = RBRes['VectMag']
	301	Cart = np.zeros_like(vecs[0])
	302	for vec,mag in zip(vecs,mags):
	303	Cart += vec*mag
	304	elif model['Type'] == 'Residue':
	305	Cart = np.array(RBRes['rbXYZ'])
	306	for itor,seq in enumerate(RBRes['rbSeq']):
	307	QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
	308	for ride in seq[3]:
	309	Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
[936]	310	if model['MolCent'][1]:
	311	Cart -= model['MolCent'][0]
[934]	312	for i,x in enumerate(Cart):
	313	xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
	314	atType = RBRes['rbTypes'][i]
	315	atTypes.append(atType)
	316	iatm += 1
	317	return np.array(xyz),atTypes
[851]	318
[936]	319	def SetMolCent(model,RBData):
	320	rideList = []
	321	RBRes = RBData[model['Type']][model['RBId']]
	322	if model['Type'] == 'Vector':
	323	vecs = RBRes['rbVect']
	324	mags = RBRes['VectMag']
	325	Cart = np.zeros_like(vecs[0])
	326	for vec,mag in zip(vecs,mags):
	327	Cart += vec*mag
	328	elif model['Type'] == 'Residue':
	329	Cart = np.array(RBRes['rbXYZ'])
	330	for seq in RBRes['rbSeq']:
	331	rideList += seq[3]
	332	centList = set(range(len(Cart)))-set(rideList)
	333	cent = np.zeros(3)
	334	for i in centList:
	335	cent += Cart[i]
	336	model['MolCent'][0] = cent/len(centList)
	337
[881]	338	def UpdateRBUIJ(Bmat,Cart,RBObj):
	339	''' Returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
	340	'''
	341	TLStype,TLS = RBObj['ThermalMotion'][:2]
[885]	342	T = np.zeros(6)
	343	L = np.zeros(6)
	344	S = np.zeros(8)
[881]	345	if 'T' in TLStype:
[882]	346	T = TLS[:6]
[881]	347	if 'L' in TLStype:
[882]	348	L = np.array(TLS[6:12])(np.pi/180.)*2
[881]	349	if 'S' in TLStype:
[882]	350	S = np.array(TLS[12:])*(np.pi/180.)
[915]	351	g = nl.inv(np.inner(Bmat,Bmat))
	352	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
[881]	353	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
	354	Uout = []
[882]	355	Q = RBObj['Orient'][0]
[881]	356	for X in Cart:
[882]	357	X = prodQVQ(Q,X)
[881]	358	if 'U' in TLStype:
	359	Uout.append(['I',TLS[0],0,0,0,0,0,0])
[882]	360	elif not 'N' in TLStype:
	361	U = [0,0,0,0,0,0]
	362	U[0] = T[0]+L[1]X[2]2+L[2]X[1]*2-2.0L[5]X[1]X[2]+2.0(S[2]X[2]-S[4]*X[1])
	363	U[1] = T[1]+L[0]X[2]2+L[2]X[0]*2-2.0L[4]X[0]X[2]+2.0(S[5]X[0]-S[0]*X[2])
	364	U[2] = T[2]+L[1]X[0]2+L[0]X[1]*2-2.0L[3]X[1]X[0]+2.0(S[1]X[1]-S[3]*X[0])
	365	U[3] = T[3]+L[4]X[1]X[2]+L[5]X[0]X[2]-L[3]X[2]2-L[2]X[0]*X[1]+ \
	366	S[4]X[0]-S[5]X[1]-(S[6]+S[7])*X[2]
	367	U[4] = T[4]+L[3]X[1]X[2]+L[5]X[0]X[1]-L[4]X[1]2-L[1]X[0]*X[2]+ \
	368	S[3]X[2]-S[2]X[0]+S[6]*X[1]
	369	U[5] = T[5]+L[3]X[0]X[2]+L[4]X[0]X[1]-L[5]X[0]2-L[0]X[2]*X[1]+ \
	370	S[0]X[1]-S[1]X[2]+S[7]*X[0]
	371	Umat = G2lat.U6toUij(U)
[910]	372	beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
[881]	373	Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
[885]	374	else:
	375	Uout.append(['N',])
[881]	376	return Uout
	377
[826]	378	def GetSHCoeff(pId,parmDict,SHkeys):
	379	SHCoeff = {}
	380	for shkey in SHkeys:
	381	shname = str(pId)+'::'+shkey
	382	SHCoeff[shkey] = parmDict[shname]
	383	return SHCoeff
[808]	384
[763]	385	def getMass(generalData):
	386	mass = 0.
	387	for i,elem in enumerate(generalData['AtomTypes']):
	388	mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
	389	return mass
	390
	391	def getDensity(generalData):
	392
	393	mass = getMass(generalData)
	394	Volume = generalData['Cell'][7]
	395	density = mass/(0.6022137*Volume)
	396	return density,Volume/mass
	397
[815]	398	def getSyXYZ(XYZ,ops,SGData):
	399	XYZout = np.zeros_like(XYZ)
	400	for i,[xyz,op] in enumerate(zip(XYZ,ops)):
	401	if op == '1':
	402	XYZout[i] = xyz
	403	else:
	404	oprs = op.split('+')
	405	unit = [0,0,0]
	406	if oprs[1]:
	407	unit = np.array(list(eval(oprs[1])))
	408	syop =int(oprs[0])
	409	inv = syop/abs(syop)
	410	syop *= inv
	411	cent = syop/100
	412	syop %= 100
	413	syop -= 1
	414	M,T = SGData['SGOps'][syop]
	415	XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
	416	return XYZout
	417
[763]	418	def getRestDist(XYZ,Amat):
	419	return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
[815]	420
[820]	421	def getRestDeriv(Func,XYZ,Amat,ops,SGData):
	422	deriv = np.zeros((len(XYZ),3))
[815]	423	dx = 0.00001
	424	for j,xyz in enumerate(XYZ):
[820]	425	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
[902]	426	XYZ[j] -= x
[820]	427	d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
[902]	428	XYZ[j] += 2*x
[820]	429	d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
[902]	430	XYZ[j] -= x
[820]	431	deriv[j][i] = (d1-d2)/(2*dx)
	432	return deriv.flatten()
[763]	433
	434	def getRestAngle(XYZ,Amat):
	435
	436	def calcVec(Ox,Tx,Amat):
	437	return np.inner(Amat,(Tx-Ox))
	438
	439	VecA = calcVec(XYZ[1],XYZ[0],Amat)
	440	VecA /= np.sqrt(np.sum(VecA**2))
	441	VecB = calcVec(XYZ[1],XYZ[2],Amat)
	442	VecB /= np.sqrt(np.sum(VecB**2))
	443	edge = VecB-VecA
	444	edge = np.sum(edge**2)
	445	angle = (2.-edge)/2.
	446	angle = max(angle,-1.)
	447	return acosd(angle)
	448
	449	def getRestPlane(XYZ,Amat):
	450	sumXYZ = np.zeros(3)
	451	for xyz in XYZ:
	452	sumXYZ += xyz
	453	sumXYZ /= len(XYZ)
	454	XYZ = np.array(XYZ)-sumXYZ
	455	XYZ = np.inner(Amat,XYZ).T
	456	Zmat = np.zeros((3,3))
	457	for i,xyz in enumerate(XYZ):
	458	Zmat += np.outer(xyz.T,xyz)
	459	Evec,Emat = nl.eig(Zmat)
	460	Evec = np.sqrt(Evec)/(len(XYZ)-3)
	461	Order = np.argsort(Evec)
	462	return Evec[Order[0]]
	463
[815]	464	def getRestChiral(XYZ,Amat):
[763]	465	VecA = np.empty((3,3))
	466	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
	467	VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
	468	VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
	469	return nl.det(VecA)
[815]	470
[818]	471	def getRestTorsion(XYZ,Amat):
[815]	472	VecA = np.empty((3,3))
	473	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
	474	VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
	475	VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
	476	D = nl.det(VecA)
	477	Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
	478	P12 = np.sum(VecA[0]VecA[1])/(Mag[0]Mag[1])
	479	P13 = np.sum(VecA[0]VecA[2])/(Mag[0]Mag[2])
	480	P23 = np.sum(VecA[1]VecA[2])/(Mag[1]Mag[2])
	481	Ang = 1.0
	482	if abs(P12) < 1.0 and abs(P23) < 1.0:
	483	Ang = (P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23**2))
	484	TOR = (acosd(Ang)*D/abs(D)+720.)%360.
[818]	485	return TOR
	486
	487	def calcTorsionEnergy(TOR,Coeff=[]):
[815]	488	sum = 0.
	489	if len(Coeff):
	490	cof = np.reshape(Coeff,(3,3)).T
	491	delt = TOR-cof[1]
	492	delt = np.where(delt<-180.,delt+360.,delt)
	493	delt = np.where(delt>180.,delt-360.,delt)
	494	term = -cof[2]delt*2
[817]	495	val = cof[0]*np.exp(term/1000.0)
	496	pMax = cof[0][np.argmin(val)]
[818]	497	Eval = np.sum(val)
	498	sum = Eval-pMax
	499	return sum,Eval
[817]	500
[822]	501	def getTorsionDeriv(XYZ,Amat,Coeff):
	502	deriv = np.zeros((len(XYZ),3))
	503	dx = 0.00001
	504	for j,xyz in enumerate(XYZ):
	505	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
[902]	506	XYZ[j] -= x
[822]	507	tor = getRestTorsion(XYZ,Amat)
	508	p,d1 = calcTorsionEnergy(tor,Coeff)
[902]	509	XYZ[j] += 2*x
[822]	510	tor = getRestTorsion(XYZ,Amat)
	511	p,d2 = calcTorsionEnergy(tor,Coeff)
[902]	512	XYZ[j] -= x
	513	deriv[j][i] = (d2-d1)/(2*dx)
[822]	514	return deriv.flatten()
	515
	516	def getRestRama(XYZ,Amat):
[818]	517	phi = getRestTorsion(XYZ[:5],Amat)
	518	psi = getRestTorsion(XYZ[1:],Amat)
	519	return phi,psi
	520
	521	def calcRamaEnergy(phi,psi,Coeff=[]):
	522	sum = 0.
[817]	523	if len(Coeff):
	524	cof = Coeff.T
	525	dPhi = phi-cof[1]
	526	dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
	527	dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
	528	dPsi = psi-cof[2]
	529	dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
	530	dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
	531	val = -cof[3]dPhi2-cof[4]dPsi*2-2.0cof[5]dPhidPsi
	532	val = cof[0]*np.exp(val/1000.)
	533	pMax = cof[0][np.argmin(val)]
[818]	534	Eval = np.sum(val)
	535	sum = Eval-pMax
	536	return sum,Eval
[822]	537
	538	def getRamaDeriv(XYZ,Amat,Coeff):
	539	deriv = np.zeros((len(XYZ),3))
	540	dx = 0.00001
	541	for j,xyz in enumerate(XYZ):
	542	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
[902]	543	XYZ[j] -= x
[822]	544	phi,psi = getRestRama(XYZ,Amat)
	545	p,d1 = calcRamaEnergy(phi,psi,Coeff)
[902]	546	XYZ[j] += 2*x
[822]	547	phi,psi = getRestRama(XYZ,Amat)
	548	p,d2 = calcRamaEnergy(phi,psi,Coeff)
[902]	549	XYZ[j] -= x
	550	deriv[j][i] = (d2-d1)/(2*dx)
[822]	551	return deriv.flatten()
[823]	552
	553	def getRestPolefig(ODFln,SamSym,Grid):
	554	X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
[824]	555	R,P = np.sqrt(X2+Y2).flatten(),atan2d(Y,X).flatten()
[823]	556	R = np.where(R <= 1.,2.*atand(R),0.0)
	557	Z = np.zeros_like(R)
	558	Z = G2lat.polfcal(ODFln,SamSym,R,P)
	559	Z = np.reshape(Z,(Grid,Grid))
	560	return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
	561
	562	def getRestPolefigDerv(HKL,Grid,SHCoeff):
	563	pass
[818]	564
[763]	565	def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
	566
	567	def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
	568	TxT = inv*(np.inner(M,Tx)+T)+C+U
	569	return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
	570
	571	inv = Top/abs(Top)
	572	cent = abs(Top)/100
	573	op = abs(Top)%100-1
	574	M,T = SGData['SGOps'][op]
	575	C = SGData['SGCen'][cent]
	576	dx = .00001
	577	deriv = np.zeros(6)
	578	for i in [0,1,2]:
[902]	579	Oxyz[i] -= dx
[763]	580	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
[902]	581	Oxyz[i] += 2*dx
[763]	582	deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
[902]	583	Oxyz[i] -= dx
	584	Txyz[i] -= dx
[763]	585	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
[902]	586	Txyz[i] += 2*dx
[763]	587	deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
[902]	588	Txyz[i] -= dx
[763]	589	return deriv
	590
	591	def getAngSig(VA,VB,Amat,SGData,covData={}):
	592
	593	def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
	594	TxT = inv*(np.inner(M,Tx)+T)+C
	595	TxT = G2spc.MoveToUnitCell(TxT)+U
	596	return np.inner(Amat,(TxT-Ox))
	597
	598	def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
	599	VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
	600	VecA /= np.sqrt(np.sum(VecA**2))
	601	VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
	602	VecB /= np.sqrt(np.sum(VecB**2))
	603	edge = VecB-VecA
	604	edge = np.sum(edge**2)
	605	angle = (2.-edge)/2.
	606	angle = max(angle,-1.)
	607	return acosd(angle)
	608
	609	OxAN,OxA,TxAN,TxA,unitA,TopA = VA
	610	OxBN,OxB,TxBN,TxB,unitB,TopB = VB
	611	invA = invB = 1
	612	invA = TopA/abs(TopA)
	613	invB = TopB/abs(TopB)
	614	centA = abs(TopA)/100
	615	centB = abs(TopB)/100
	616	opA = abs(TopA)%100-1
	617	opB = abs(TopB)%100-1
	618	MA,TA = SGData['SGOps'][opA]
	619	MB,TB = SGData['SGOps'][opB]
	620	CA = SGData['SGCen'][centA]
	621	CB = SGData['SGCen'][centB]
	622	if 'covMatrix' in covData:
	623	covMatrix = covData['covMatrix']
	624	varyList = covData['varyList']
	625	AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
	626	dx = .00001
	627	dadx = np.zeros(9)
	628	Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
	629	for i in [0,1,2]:
[902]	630	OxA[i] -= dx
[763]	631	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
[902]	632	OxA[i] += 2*dx
	633	dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
	634	OxA[i] -= dx
[763]	635
[902]	636	TxA[i] -= dx
[763]	637	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
[902]	638	TxA[i] += 2*dx
	639	dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
	640	TxA[i] -= dx
[763]	641
[902]	642	TxB[i] -= dx
[763]	643	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
[902]	644	TxB[i] += 2*dx
	645	dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
	646	TxB[i] -= dx
[763]	647
	648	sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
	649	if sigAng < 0.01:
	650	sigAng = 0.0
	651	return Ang,sigAng
	652	else:
	653	return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
	654
	655	def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
	656
	657	def calcDist(Atoms,SyOps,Amat):
	658	XYZ = []
	659	for i,atom in enumerate(Atoms):
	660	Inv,M,T,C,U = SyOps[i]
	661	XYZ.append(np.array(atom[1:4]))
	662	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	663	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	664	V1 = XYZ[1]-XYZ[0]
	665	return np.sqrt(np.sum(V1**2))
	666
	667	Inv = []
	668	SyOps = []
	669	names = []
	670	for i,atom in enumerate(Oatoms):
	671	names += atom[-1]
	672	Op,unit = Atoms[i][-1]
	673	inv = Op/abs(Op)
	674	m,t = SGData['SGOps'][abs(Op)%100-1]
	675	c = SGData['SGCen'][abs(Op)/100]
	676	SyOps.append([inv,m,t,c,unit])
	677	Dist = calcDist(Oatoms,SyOps,Amat)
	678
	679	sig = -0.001
	680	if 'covMatrix' in covData:
	681	parmNames = []
	682	dx = .00001
	683	dadx = np.zeros(6)
	684	for i in range(6):
	685	ia = i/3
	686	ix = i%3
	687	Oatoms[ia][ix+1] += dx
	688	a0 = calcDist(Oatoms,SyOps,Amat)
	689	Oatoms[ia][ix+1] -= 2*dx
	690	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
	691	covMatrix = covData['covMatrix']
	692	varyList = covData['varyList']
	693	DistVcov = getVCov(names,varyList,covMatrix)
	694	sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
	695	if sig < 0.001:
	696	sig = -0.001
	697
	698	return Dist,sig
	699
	700	def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
	701
	702	def calcAngle(Atoms,SyOps,Amat):
	703	XYZ = []
	704	for i,atom in enumerate(Atoms):
	705	Inv,M,T,C,U = SyOps[i]
	706	XYZ.append(np.array(atom[1:4]))
	707	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	708	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	709	V1 = XYZ[1]-XYZ[0]
	710	V1 /= np.sqrt(np.sum(V1**2))
	711	V2 = XYZ[1]-XYZ[2]
	712	V2 /= np.sqrt(np.sum(V2**2))
	713	V3 = V2-V1
	714	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
	715	return acosd(cang)
	716
	717	Inv = []
	718	SyOps = []
	719	names = []
	720	for i,atom in enumerate(Oatoms):
	721	names += atom[-1]
	722	Op,unit = Atoms[i][-1]
	723	inv = Op/abs(Op)
	724	m,t = SGData['SGOps'][abs(Op)%100-1]
	725	c = SGData['SGCen'][abs(Op)/100]
	726	SyOps.append([inv,m,t,c,unit])
	727	Angle = calcAngle(Oatoms,SyOps,Amat)
	728
	729	sig = -0.01
	730	if 'covMatrix' in covData:
	731	parmNames = []
	732	dx = .00001
	733	dadx = np.zeros(9)
	734	for i in range(9):
	735	ia = i/3
	736	ix = i%3
	737	Oatoms[ia][ix+1] += dx
	738	a0 = calcAngle(Oatoms,SyOps,Amat)
	739	Oatoms[ia][ix+1] -= 2*dx
	740	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
	741	covMatrix = covData['covMatrix']
	742	varyList = covData['varyList']
	743	AngVcov = getVCov(names,varyList,covMatrix)
	744	sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
	745	if sig < 0.01:
	746	sig = -0.01
	747
	748	return Angle,sig
	749
	750	def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
	751
	752	def calcTorsion(Atoms,SyOps,Amat):
	753
	754	XYZ = []
	755	for i,atom in enumerate(Atoms):
	756	Inv,M,T,C,U = SyOps[i]
	757	XYZ.append(np.array(atom[1:4]))
	758	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	759	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	760	V1 = XYZ[1]-XYZ[0]
	761	V2 = XYZ[2]-XYZ[1]
	762	V3 = XYZ[3]-XYZ[2]
	763	V1 /= np.sqrt(np.sum(V1**2))
	764	V2 /= np.sqrt(np.sum(V2**2))
	765	V3 /= np.sqrt(np.sum(V3**2))
	766	M = np.array([V1,V2,V3])
	767	D = nl.det(M)
	768	Ang = 1.0
	769	P12 = np.dot(V1,V2)
	770	P13 = np.dot(V1,V3)
	771	P23 = np.dot(V2,V3)
	772	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
	773	return Tors
	774
	775	Inv = []
	776	SyOps = []
	777	names = []
	778	for i,atom in enumerate(Oatoms):
	779	names += atom[-1]
	780	Op,unit = Atoms[i][-1]
	781	inv = Op/abs(Op)
	782	m,t = SGData['SGOps'][abs(Op)%100-1]
	783	c = SGData['SGCen'][abs(Op)/100]
	784	SyOps.append([inv,m,t,c,unit])
	785	Tors = calcTorsion(Oatoms,SyOps,Amat)
	786
	787	sig = -0.01
	788	if 'covMatrix' in covData:
	789	parmNames = []
	790	dx = .00001
	791	dadx = np.zeros(12)
	792	for i in range(12):
	793	ia = i/3
	794	ix = i%3
[902]	795	Oatoms[ia][ix+1] -= dx
[763]	796	a0 = calcTorsion(Oatoms,SyOps,Amat)
[902]	797	Oatoms[ia][ix+1] += 2*dx
[763]	798	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
[902]	799	Oatoms[ia][ix+1] -= dx
[763]	800	covMatrix = covData['covMatrix']
	801	varyList = covData['varyList']
	802	TorVcov = getVCov(names,varyList,covMatrix)
	803	sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
	804	if sig < 0.01:
	805	sig = -0.01
	806
	807	return Tors,sig
	808
	809	def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
	810
	811	def calcDist(Atoms,SyOps,Amat):
	812	XYZ = []
	813	for i,atom in enumerate(Atoms):
	814	Inv,M,T,C,U = SyOps[i]
	815	XYZ.append(np.array(atom[1:4]))
	816	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	817	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	818	V1 = XYZ[1]-XYZ[0]
	819	return np.sqrt(np.sum(V1**2))
	820
	821	def calcAngle(Atoms,SyOps,Amat):
	822	XYZ = []
	823	for i,atom in enumerate(Atoms):
	824	Inv,M,T,C,U = SyOps[i]
	825	XYZ.append(np.array(atom[1:4]))
	826	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	827	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	828	V1 = XYZ[1]-XYZ[0]
	829	V1 /= np.sqrt(np.sum(V1**2))
	830	V2 = XYZ[1]-XYZ[2]
	831	V2 /= np.sqrt(np.sum(V2**2))
	832	V3 = V2-V1
	833	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
	834	return acosd(cang)
	835
	836	def calcTorsion(Atoms,SyOps,Amat):
	837
	838	XYZ = []
	839	for i,atom in enumerate(Atoms):
	840	Inv,M,T,C,U = SyOps[i]
	841	XYZ.append(np.array(atom[1:4]))
	842	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
	843	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
	844	V1 = XYZ[1]-XYZ[0]
	845	V2 = XYZ[2]-XYZ[1]
	846	V3 = XYZ[3]-XYZ[2]
	847	V1 /= np.sqrt(np.sum(V1**2))
	848	V2 /= np.sqrt(np.sum(V2**2))
	849	V3 /= np.sqrt(np.sum(V3**2))
	850	M = np.array([V1,V2,V3])
	851	D = nl.det(M)
	852	Ang = 1.0
	853	P12 = np.dot(V1,V2)
	854	P13 = np.dot(V1,V3)
	855	P23 = np.dot(V2,V3)
	856	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
	857	return Tors
	858
	859	Inv = []
	860	SyOps = []
	861	names = []
	862	for i,atom in enumerate(Oatoms):
	863	names += atom[-1]
	864	Op,unit = Atoms[i][-1]
	865	inv = Op/abs(Op)
	866	m,t = SGData['SGOps'][abs(Op)%100-1]
	867	c = SGData['SGCen'][abs(Op)/100]
	868	SyOps.append([inv,m,t,c,unit])
	869	M = len(Oatoms)
	870	if M == 2:
	871	Val = calcDist(Oatoms,SyOps,Amat)
	872	elif M == 3:
	873	Val = calcAngle(Oatoms,SyOps,Amat)
	874	else:
	875	Val = calcTorsion(Oatoms,SyOps,Amat)
	876
	877	sigVals = [-0.001,-0.01,-0.01]
	878	sig = sigVals[M-3]
	879	if 'covMatrix' in covData:
	880	parmNames = []
	881	dx = .00001
	882	N = M*3
	883	dadx = np.zeros(N)
	884	for i in range(N):
	885	ia = i/3
	886	ix = i%3
	887	Oatoms[ia][ix+1] += dx
	888	if M == 2:
	889	a0 = calcDist(Oatoms,SyOps,Amat)
	890	elif M == 3:
	891	a0 = calcAngle(Oatoms,SyOps,Amat)
	892	else:
	893	a0 = calcTorsion(Oatoms,SyOps,Amat)
	894	Oatoms[ia][ix+1] -= 2*dx
	895	if M == 2:
	896	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
	897	elif M == 3:
	898	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
	899	else:
	900	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
	901	covMatrix = covData['covMatrix']
	902	varyList = covData['varyList']
	903	Vcov = getVCov(names,varyList,covMatrix)
	904	sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
	905	if sig < sigVals[M-3]:
	906	sig = sigVals[M-3]
	907
	908	return Val,sig
	909
	910
	911	def ValEsd(value,esd=0,nTZ=False): #NOT complete - don't use
	912	# returns value(esd) string; nTZ=True for no trailing zeros
	913	# use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
	914	#get the 2 significant digits in the esd
	915	edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
	916	#get the number of digits to represent them
	917	epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
	918
	919	mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
	920	ndec = lambda esd: int(1.545-math.log10(abs(esd)))
	921	if esd > 0:
	922	fmt = '"%.'+str(ndec(esd))+'f(%d)"'
	923	return str(fmt%(value,int(round(esd10*(mdec(esd)))))).strip('"')
	924	elif esd < 0:
	925	return str(round(value,mdec(esd)-1))
	926	else:
	927	text = str("%f"%(value))
	928	if nTZ:
	929	return text.rstrip('0')
	930	else:
	931	return text
	932
	933	def adjHKLmax(SGData,Hmax):
	934	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
	935	Hmax[0] = ((Hmax[0]+3)/6)*6
	936	Hmax[1] = ((Hmax[1]+3)/6)*6
	937	Hmax[2] = ((Hmax[2]+1)/4)*4
	938	else:
	939	Hmax[0] = ((Hmax[0]+2)/4)*4
	940	Hmax[1] = ((Hmax[1]+2)/4)*4
	941	Hmax[2] = ((Hmax[2]+1)/4)*4
	942
[881]	943	def OmitMap(data,reflData):
	944	generalData = data['General']
	945	if not generalData['Map']['MapType']:
	946	print '**** ERROR - Fourier map not defined'
	947	return
	948	mapData = generalData['Map']
	949	dmin = mapData['Resolution']
	950	SGData = generalData['SGData']
	951	cell = generalData['Cell'][1:8]
	952	A = G2lat.cell2A(cell[:6])
	953	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
	954	adjHKLmax(SGData,Hmax)
	955	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
	956	time0 = time.time()
	957	for ref in reflData:
	958	if ref[4] >= dmin:
	959	Fosq,Fcsq,ph = ref[8:11]
	960	for i,hkl in enumerate(ref[11]):
	961	hkl = np.asarray(hkl,dtype='i')
	962	dp = 360.*ref[12][i]
	963	a = cosd(ph+dp)
	964	b = sind(ph+dp)
	965	phasep = complex(a,b)
	966	phasem = complex(a,-b)
	967	F = np.sqrt(Fosq)
	968	h,k,l = hkl+Hmax
	969	Fhkl[h,k,l] = F*phasep
	970	h,k,l = -hkl+Hmax
	971	Fhkl[h,k,l] = F*phasem
	972	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
	973	print 'NB: this is just an Fobs map for now - under development'
	974	print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
	975	print rho.shape
	976	mapData['rho'] = np.real(rho)
	977	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
	978	return mapData
	979
[763]	980	def FourierMap(data,reflData):
	981
	982	generalData = data['General']
	983	if not generalData['Map']['MapType']:
	984	print '**** ERROR - Fourier map not defined'
	985	return
	986	mapData = generalData['Map']
	987	dmin = mapData['Resolution']
	988	SGData = generalData['SGData']
	989	cell = generalData['Cell'][1:8]
	990	A = G2lat.cell2A(cell[:6])
	991	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
	992	adjHKLmax(SGData,Hmax)
	993	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
	994	# Fhkl[0,0,0] = generalData['F000X']
	995	time0 = time.time()
	996	for ref in reflData:
	997	if ref[4] >= dmin:
	998	Fosq,Fcsq,ph = ref[8:11]
	999	for i,hkl in enumerate(ref[11]):
	1000	hkl = np.asarray(hkl,dtype='i')
	1001	dp = 360.*ref[12][i]
	1002	a = cosd(ph+dp)
	1003	b = sind(ph+dp)
	1004	phasep = complex(a,b)
	1005	phasem = complex(a,-b)
	1006	if 'Fobs' in mapData['MapType']:
	1007	F = np.sqrt(Fosq)
	1008	h,k,l = hkl+Hmax
	1009	Fhkl[h,k,l] = F*phasep
	1010	h,k,l = -hkl+Hmax
	1011	Fhkl[h,k,l] = F*phasem
	1012	elif 'Fcalc' in mapData['MapType']:
	1013	F = np.sqrt(Fcsq)
	1014	h,k,l = hkl+Hmax
	1015	Fhkl[h,k,l] = F*phasep
	1016	h,k,l = -hkl+Hmax
	1017	Fhkl[h,k,l] = F*phasem
	1018	elif 'delt-F' in mapData['MapType']:
	1019	dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
	1020	h,k,l = hkl+Hmax
	1021	Fhkl[h,k,l] = dF*phasep
	1022	h,k,l = -hkl+Hmax
	1023	Fhkl[h,k,l] = dF*phasem
	1024	elif '2*Fo-Fc' in mapData['MapType']:
	1025	F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
	1026	h,k,l = hkl+Hmax
	1027	Fhkl[h,k,l] = F*phasep
	1028	h,k,l = -hkl+Hmax
	1029	Fhkl[h,k,l] = F*phasem
	1030	elif 'Patterson' in mapData['MapType']:
	1031	h,k,l = hkl+Hmax
	1032	Fhkl[h,k,l] = complex(Fosq,0.)
	1033	h,k,l = -hkl+Hmax
	1034	Fhkl[h,k,l] = complex(Fosq,0.)
	1035	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
	1036	print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
	1037	mapData['rho'] = np.real(rho)
	1038	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
	1039	return mapData
	1040
	1041	# map printing for testing purposes
	1042	def printRho(SGLaue,rho,rhoMax):
	1043	dim = len(rho.shape)
	1044	if dim == 2:
	1045	ix,jy = rho.shape
	1046	for j in range(jy):
	1047	line = ''
	1048	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
	1049	line += (jy-j)*' '
	1050	for i in range(ix):
	1051	r = int(100*rho[i,j]/rhoMax)
	1052	line += '%4d'%(r)
	1053	print line+'\n'
	1054	else:
	1055	ix,jy,kz = rho.shape
	1056	for k in range(kz):
	1057	print 'k = ',k
	1058	for j in range(jy):
	1059	line = ''
	1060	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
	1061	line += (jy-j)*' '
	1062	for i in range(ix):
	1063	r = int(100*rho[i,j,k]/rhoMax)
	1064	line += '%4d'%(r)
	1065	print line+'\n'
	1066	## keep this
	1067
	1068	def findOffset(SGData,A,Fhkl):
	1069	if SGData['SpGrp'] == 'P 1':
	1070	return [0,0,0]
	1071	hklShape = Fhkl.shape
	1072	hklHalf = np.array(hklShape)/2
	1073	sortHKL = np.argsort(Fhkl.flatten())
	1074	Fdict = {}
	1075	for hkl in sortHKL:
	1076	HKL = np.unravel_index(hkl,hklShape)
	1077	F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
	1078	if F == 0.:
	1079	break
	1080	Fdict['%.6f'%(np.absolute(F))] = hkl
	1081	Flist = np.flipud(np.sort(Fdict.keys()))
	1082	F = str(1.e6)
	1083	i = 0
	1084	DH = []
	1085	Dphi = []
[805]	1086	Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
	1087	for F in Flist:
[763]	1088	hkl = np.unravel_index(Fdict[F],hklShape)
	1089	iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
	1090	Uniq = np.array(Uniq,dtype='i')
	1091	Phi = np.array(Phi)
	1092	Uniq = np.concatenate((Uniq,-Uniq))+hklHalf # put in Friedel pairs & make as index to Farray
	1093	Phi = np.concatenate((Phi,-Phi)) # and their phase shifts
	1094	Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
	1095	ang0 = np.angle(Fh0,deg=True)/360.
[803]	1096	for H,phi in zip(Uniq,Phi)[1:]:
	1097	ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
[763]	1098	dH = H-hkl
	1099	dang = ang-ang0
	1100	if np.any(np.abs(dH)-Hmax > 0): #keep low order DHs
	1101	continue
	1102	DH.append(dH)
[803]	1103	Dphi.append((dang+.5) % 1.0)
[805]	1104	if i > 20 or len(DH) > 30:
	1105	break
[763]	1106	i += 1
	1107	DH = np.array(DH)
[805]	1108	print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
[763]	1109	Dphi = np.array(Dphi)
[805]	1110	steps = np.array(hklShape)
[763]	1111	X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
	1112	XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
[805]	1113	Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
	1114	Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
	1115	hist,bins = np.histogram(Mmap,bins=1000)
[803]	1116	# for i,item in enumerate(hist[:10]):
	1117	# print item,bins[i]
[763]	1118	chisq = np.min(Mmap)
	1119	DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
	1120	print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
[805]	1121	# print (np.dot(DX,DH.T)+.5)%1.-Dphi
[763]	1122	return DX
	1123
	1124	def ChargeFlip(data,reflData,pgbar):
	1125	generalData = data['General']
	1126	mapData = generalData['Map']
	1127	flipData = generalData['Flip']
	1128	FFtable = {}
	1129	if 'None' not in flipData['Norm element']:
	1130	normElem = flipData['Norm element'].upper()
	1131	FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
	1132	for ff in FFs:
	1133	if ff['Symbol'] == normElem:
	1134	FFtable.update(ff)
	1135	dmin = flipData['Resolution']
	1136	SGData = generalData['SGData']
	1137	cell = generalData['Cell'][1:8]
	1138	A = G2lat.cell2A(cell[:6])
	1139	Vol = cell[6]
	1140	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
	1141	adjHKLmax(SGData,Hmax)
	1142	Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no.
	1143	time0 = time.time()
	1144	for ref in reflData:
	1145	dsp = ref[4]
	1146	if dsp >= dmin:
	1147	ff = 0.1Vol #est. no. atoms for ~10A*3/atom
	1148	if FFtable:
	1149	SQ = 0.25/dsp**2
	1150	ff *= G2el.ScatFac(FFtable,SQ)[0]
[803]	1151	if ref[8] > 0.: #use only +ve Fobs**2
[763]	1152	E = np.sqrt(ref[8])/ff
	1153	else:
	1154	E = 0.
	1155	ph = ref[10]
	1156	ph = rn.uniform(0.,360.)
	1157	for i,hkl in enumerate(ref[11]):
	1158	hkl = np.asarray(hkl,dtype='i')
	1159	dp = 360.*ref[12][i]
	1160	a = cosd(ph+dp)
	1161	b = sind(ph+dp)
	1162	phasep = complex(a,b)
	1163	phasem = complex(a,-b)
	1164	h,k,l = hkl+Hmax
	1165	Ehkl[h,k,l] = E*phasep
	1166	h,k,l = -hkl+Hmax #Friedel pair refl.
	1167	Ehkl[h,k,l] = E*phasem
	1168	# Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0]
	1169	CEhkl = copy.copy(Ehkl)
	1170	MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
	1171	Emask = ma.getmask(MEhkl)
	1172	sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
	1173	Ncyc = 0
	1174	old = np.seterr(all='raise')
	1175	while True:
	1176	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
	1177	CEsig = np.std(CErho)
	1178	CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
[795]	1179	CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho) #solves U atom problem! make 20. adjustible
[763]	1180	CFhkl = fft.ifftshift(fft.ifftn(CFrho))
	1181	CFhkl = np.where(CFhkl,CFhkl,1.0) #avoid divide by zero
	1182	phase = CFhkl/np.absolute(CFhkl)
	1183	CEhkl = np.absolute(Ehkl)*phase
	1184	Ncyc += 1
	1185	sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
	1186	DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
	1187	Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
	1188	if Rcf < 5.:
	1189	break
	1190	GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
	1191	if not GoOn or Ncyc > 10000:
	1192	break
	1193	np.seterr(**old)
	1194	print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
	1195	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
	1196	print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
	1197	roll = findOffset(SGData,A,CEhkl)
	1198
	1199	mapData['Rcf'] = Rcf
	1200	mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
	1201	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
	1202	mapData['rollMap'] = [0,0,0]
	1203	return mapData
	1204
	1205	def SearchMap(data):
	1206	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
	1207
	1208	norm = 1./(np.sqrt(3.)np.sqrt(2.np.pi)**3)
	1209
	1210	def noDuplicate(xyz,peaks,Amat):
	1211	XYZ = np.inner(Amat,xyz)
	1212	if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
	1213	print ' Peak',xyz,' <0.5A from another peak'
	1214	return False
	1215	return True
	1216
	1217	def fixSpecialPos(xyz,SGData,Amat):
	1218	equivs = G2spc.GenAtom(xyz,SGData,Move=True)
	1219	X = []
	1220	xyzs = [equiv[0] for equiv in equivs]
	1221	for x in xyzs:
	1222	if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
	1223	X.append(x)
	1224	if len(X) > 1:
	1225	return np.average(X,axis=0)
	1226	else:
	1227	return xyz
	1228
	1229	def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
	1230	Mag,x0,y0,z0,sig = parms
[779]	1231	z = -((x0-rX)2+(y0-rY)2+(z0-rZ)*2)/(2.sig**2)
	1232	# return normMagnp.exp(z)/(sigres*3) #not slower but some faults in LS
	1233	return normMag(1.+z+z*2/2.)/(sigres**3)
[763]	1234
	1235	def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
	1236	Mag,x0,y0,z0,sig = parms
	1237	M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
	1238	return M
	1239
	1240	def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
	1241	Mag,x0,y0,z0,sig = parms
	1242	dMdv = np.zeros(([5,]+list(rX.shape)))
	1243	delt = .01
	1244	for i in range(5):
	1245	parms[i] -= delt
	1246	rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
	1247	parms[i] += 2.*delt
	1248	rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
	1249	parms[i] -= delt
	1250	dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
	1251	rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
	1252	Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
	1253	dMdv = np.reshape(dMdv,(5,rX.size))
	1254	Hess = np.inner(dMdv,dMdv)
	1255
	1256	return Vec,Hess
	1257
	1258	generalData = data['General']
	1259	phaseName = generalData['Name']
	1260	SGData = generalData['SGData']
	1261	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
	1262	drawingData = data['Drawing']
	1263	peaks = []
	1264	mags = []
	1265	dzeros = []
	1266	try:
	1267	mapData = generalData['Map']
	1268	contLevel = mapData['cutOff']*mapData['rhoMax']/100.
	1269	rho = copy.copy(mapData['rho']) #don't mess up original
	1270	mapHalf = np.array(rho.shape)/2
	1271	res = mapData['Resolution']
	1272	incre = np.array(rho.shape,dtype=np.float)
	1273	step = max(1.0,1./res)+1
	1274	steps = np.array(3*[step,])
	1275	except KeyError:
	1276	print '**** ERROR - Fourier map not defined'
	1277	return peaks,mags
	1278	rhoMask = ma.array(rho,mask=(rho<contLevel))
	1279	indices = (-1,0,1)
	1280	rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
	1281	for roll in rolls:
	1282	if np.any(roll):
	1283	rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
	1284	indx = np.transpose(rhoMask.nonzero())
	1285	peaks = indx/incre
	1286	mags = rhoMask[rhoMask.nonzero()]
	1287	for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
	1288	rho = rollMap(rho,ind)
	1289	rMM = mapHalf-steps
	1290	rMP = mapHalf+steps+1
	1291	rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
	1292	peakInt = np.sum(rhoPeak)res*3
	1293	rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
	1294	x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig)
	1295	result = HessianLSQ(peakFunc,x0,Hess=peakHess,
	1296	args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
	1297	x1 = result[0]
	1298	if not np.any(x1 < 0):
	1299	mag = x1[0]
	1300	peak = (np.array(x1[1:4])-ind)/incre
	1301	peak = fixSpecialPos(peak,SGData,Amat)
	1302	rho = rollMap(rho,-ind)
	1303	dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
	1304	return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
	1305
	1306	def sortArray(data,pos,reverse=False):
	1307	#data is a list of items
	1308	#sort by pos in list; reverse if True
	1309	T = []
	1310	for i,M in enumerate(data):
	1311	T.append((M[pos],i))
	1312	D = dict(zip(T,data))
	1313	T.sort()
	1314	if reverse:
	1315	T.reverse()
	1316	X = []
	1317	for key in T:
	1318	X.append(D[key])
	1319	return X
[774]	1320
	1321	def PeaksEquiv(data,Ind):
	1322
	1323	def Duplicate(xyz,peaks,Amat):
	1324	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
	1325	return True
	1326	return False
	1327
	1328	generalData = data['General']
	1329	cell = generalData['Cell'][1:7]
	1330	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
	1331	A = G2lat.cell2A(cell)
	1332	SGData = generalData['SGData']
	1333	mapPeaks = data['Map Peaks']
	1334	XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
	1335	Indx = {}
	1336	for ind in Ind:
	1337	xyz = np.array(mapPeaks[ind][1:4])
[805]	1338	xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
	1339	# for x in xyzs: print x
[774]	1340	for jnd,xyz in enumerate(XYZ):
	1341	Indx[jnd] = Duplicate(xyz,xyzs,Amat)
	1342	Ind = []
	1343	for ind in Indx:
	1344	if Indx[ind]:
	1345	Ind.append(ind)
	1346	return Ind
[763]	1347
	1348	def PeaksUnique(data,Ind):
[774]	1349	# XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
[763]	1350
	1351	def noDuplicate(xyz,peaks,Amat):
	1352	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
	1353	return False
	1354	return True
	1355
	1356	generalData = data['General']
	1357	cell = generalData['Cell'][1:7]
	1358	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
	1359	A = G2lat.cell2A(cell)
	1360	SGData = generalData['SGData']
	1361	mapPeaks = data['Map Peaks']
	1362	Indx = {}
	1363	XYZ = {}
	1364	for ind in Ind:
	1365	XYZ[ind] = np.array(mapPeaks[ind][1:4])
	1366	Indx[ind] = True
	1367	for ind in Ind:
	1368	if Indx[ind]:
	1369	xyz = XYZ[ind]
	1370	for jnd in Ind:
	1371	if ind != jnd and Indx[jnd]:
	1372	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
	1373	xyzs = np.array([equiv[0] for equiv in Equiv])
	1374	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
	1375	Ind = []
	1376	for ind in Indx:
	1377	if Indx[ind]:
	1378	Ind.append(ind)
	1379	return Ind
	1380
[803]	1381	def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
	1382	ind = 0
	1383	if useFit:
	1384	ind = 1
[795]	1385	ins = {}
	1386	if 'C' in Parms['Type'][0]: #CW data - TOF later in an elif
	1387	for x in ['U','V','W','X','Y']:
[803]	1388	ins[x] = Parms[x][ind]
[795]	1389	if ifQ: #qplot - convert back to 2-theta
	1390	pos = 2.0asind(poswave/(4*math.pi))
	1391	sig = ins['U']tand(pos/2.0)2+ins['V']tand(pos/2.0)+ins['W']
	1392	gam = ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
	1393	XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st
	1394	else:
[798]	1395	if ifQ:
	1396	dsp = 2.*np.pi/pos
	1397	pos = Parms['difC']*dsp
	1398	else:
	1399	dsp = pos/Parms['difC'][1]
[796]	1400	if 'Pdabc' in Parms2:
[802]	1401	for x in ['sig-0','sig-1','X','Y']:
[803]	1402	ins[x] = Parms[x][ind]
[796]	1403	Pdabc = Parms2['Pdabc'].T
[795]	1404	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
	1405	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
	1406	else:
[798]	1407	for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
[803]	1408	ins[x] = Parms[x][ind]
[797]	1409	alp = ins['alpha']/dsp
[798]	1410	bet = ins['beta-0']+ins['beta-1']/dsp**4
	1411	sig = ins['sig-0']+ins['sig-1']dsp*2
[795]	1412	gam = ins['X']dsp+ins['Y']dsp**2
	1413	XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
	1414	return XY
[934]	1415
	1416	def mcsaSearch(data,reflType,reflData,covData,pgbar):
	1417	generalData = data['General']
	1418	mcsaControls = generalData['MCSA controls']
	1419	reflName = mcsaData['Data source']
	1420	phaseName = generalData['Name']
	1421	MCSAObjs = data['MCSAModels']
	1422
	1423	# = {'':'','Annealing':[50.0,0.001,0.90,1000],
	1424	# 'dmin':2.0,'Algolrithm':'Normal','Jump coeff':[0.95,0.5]} #'Normal','Random jump','Tremayne jump'
	1425
	1426	return {}
	1427
[795]	1428
	1429	def getWave(Parms):
	1430	try:
	1431	return Parms['Lam'][1]
	1432	except KeyError:
	1433	return Parms['Lam1'][1]
	1434
[763]	1435	def prodQQ(QA,QB):
	1436	''' Grassman quaternion product
	1437	QA,QB quaternions; q=r+ai+bj+ck
	1438	'''
	1439	D = np.zeros(4)
	1440	D[0] = QA[0]QB[0]-QA[1]QB[1]-QA[2]QB[2]-QA[3]QB[3]
	1441	D[1] = QA[0]QB[1]+QA[1]QB[0]+QA[2]QB[3]-QA[3]QB[2]
	1442	D[2] = QA[0]QB[2]-QA[1]QB[3]+QA[2]QB[0]+QA[3]QB[1]
	1443	D[3] = QA[0]QB[3]+QA[1]QB[2]-QA[2]QB[1]+QA[3]QB[0]
	1444	return D
	1445
	1446	def normQ(QA):
	1447	''' get length of quaternion & normalize it
	1448	q=r+ai+bj+ck
	1449	'''
	1450	n = np.sqrt(np.sum(np.array(QA)**2))
	1451	return QA/n
	1452
	1453	def invQ(Q):
	1454	'''
	1455	get inverse of quaternion
	1456	q=r+ai+bj+ck; q* = r-ai-bj-ck
	1457	'''
	1458	return Q*np.array([1,-1,-1,-1])
	1459
	1460	def prodQVQ(Q,V):
	1461	''' compute the quaternion vector rotation qvq-1 = v'
	1462	q=r+ai+bj+ck
	1463	'''
	1464	VP = np.zeros(3)
	1465	T2 = Q[0]*Q[1]
	1466	T3 = Q[0]*Q[2]
	1467	T4 = Q[0]*Q[3]
	1468	T5 = -Q[1]*Q[1]
	1469	T6 = Q[1]*Q[2]
	1470	T7 = Q[1]*Q[3]
	1471	T8 = -Q[2]*Q[2]
	1472	T9 = Q[2]*Q[3]
	1473	T10 = -Q[3]*Q[3]
	1474	VP[0] = 2.((T8+T10)V[0]+(T6-T4)V[1]+(T3+T7)V[2])+V[0]
	1475	VP[1] = 2.((T4+T6)V[0]+(T5+T10)V[1]+(T9-T2)V[2])+V[1]
	1476	VP[2] = 2.((T7-T3)V[0]+(T2+T9)V[1]+(T5+T8)V[2])+V[2]
	1477	return VP
	1478
	1479	def Q2Mat(Q):
	1480	''' make rotation matrix from quaternion
	1481	q=r+ai+bj+ck
	1482	'''
[915]	1483	QN = normQ(Q)
	1484	aa = QN[0]**2
	1485	ab = QN[0]*QN[1]
	1486	ac = QN[0]*QN[2]
	1487	ad = QN[0]*QN[3]
	1488	bb = QN[1]**2
	1489	bc = QN[1]*QN[2]
	1490	bd = QN[1]*QN[3]
	1491	cc = QN[2]**2
	1492	cd = QN[2]*QN[3]
	1493	dd = QN[3]**2
[763]	1494	M = [[aa+bb-cc-dd, 2.(bc-ad), 2.(ac+bd)],
	1495	[2(ad+bc), aa-bb+cc-dd, 2.(cd-ab)],
	1496	[2(bd-ac), 2.(ab+cd), aa-bb-cc+dd]]
	1497	return np.array(M)
	1498
	1499	def AV2Q(A,V):
[859]	1500	''' convert angle (radians) & vector to quaternion
[763]	1501	q=r+ai+bj+ck
	1502	'''
	1503	Q = np.zeros(4)
	1504	d = np.sqrt(np.sum(np.array(V)**2))
	1505	if d:
	1506	V /= d
	1507	else:
[934]	1508	V = np.array([0.,0.,1.])
[763]	1509	p = A/2.
	1510	Q[0] = np.cos(p)
[849]	1511	Q[1:4] = V*np.sin(p)
[763]	1512	return Q
	1513
	1514	def AVdeg2Q(A,V):
[859]	1515	''' convert angle (degrees) & vector to quaternion
[763]	1516	q=r+ai+bj+ck
	1517	'''
	1518	Q = np.zeros(4)
	1519	d = np.sqrt(np.sum(np.array(V)**2))
	1520	if d:
	1521	V /= d
	1522	else:
[934]	1523	V = np.array([0.,0.,1.])
[763]	1524	p = A/2.
	1525	Q[0] = cosd(p)
[849]	1526	Q[1:4] = V*sind(p)
[763]	1527	return Q
[859]	1528
	1529	def Q2AVdeg(Q):
	1530	''' convert quaternion to angle (degrees 0-360) & normalized vector
	1531	q=r+ai+bj+ck
	1532	'''
	1533	A = 2.*acosd(Q[0])
[885]	1534	V = np.array(Q[1:])
[934]	1535	d = np.sqrt(np.sum(V**2))
	1536	if d:
	1537	V /= d
[889]	1538	else:
	1539	A = 0.
[934]	1540	V = np.array([0.,0.,0.])
[859]	1541	return A,V
	1542
	1543	def Q2AV(Q):
	1544	''' convert quaternion to angle (radians 0-2pi) & normalized vector
	1545	q=r+ai+bj+ck
	1546	'''
	1547	A = 2.*np.arccos(Q[0])
[885]	1548	V = np.array(Q[1:])
[934]	1549	d = np.sqrt(np.sum(V**2))
	1550	if d:
	1551	V /= d
[889]	1552	else:
	1553	A = 0.
[934]	1554	V = np.array([0.,0.,0.])
[859]	1555	return A,V
	1556
[848]	1557	def makeQuat(A,B,C):
	1558	''' Make quaternion from rotation of A vector to B vector about C axis
	1559	A,B,C are np.array Cartesian 3-vectors
	1560	Returns quaternion & rotation angle in radians
[849]	1561	q=r+ai+bj+ck
[848]	1562	'''
	1563
	1564	V1 = np.cross(A,C)
	1565	V2 = np.cross(B,C)
	1566	if nl.norm(V1)*nl.norm(V2):
	1567	V1 /= nl.norm(V1)
	1568	V2 /= nl.norm(V2)
	1569	V3 = np.cross(V1,V2)
	1570	else:
[853]	1571	V3 = np.zeros(3)
[860]	1572	Q = np.array([0.,0.,0.,1.])
[848]	1573	D = 0.
	1574	if nl.norm(V3):
	1575	V3 /= nl.norm(V3)
	1576	D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
	1577	D = np.arccos(D1)/2.0
	1578	V1 = C-V3
	1579	V2 = C+V3
	1580	DM = nl.norm(V1)
	1581	DP = nl.norm(V2)
	1582	S = np.sin(D)
	1583	Q[0] = np.cos(D)
	1584	Q[1:] = V3*S
	1585	D *= 2.
	1586	if DM > DP:
	1587	D *= -1.
	1588	return Q,D
[763]	1589

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

Download in other formats: