source: trunk/GSASIImath.py @ 936

Last change on this file since 936 was 936, checked in by vondreele, 10 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
Line 
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 ###################
10import sys
11import os
12import os.path as ospath
13import random as rn
14import numpy as np
15import numpy.linalg as nl
16import numpy.ma as ma
17import cPickle
18import time
19import math
20import copy
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 936 $")
23import GSASIIElem as G2el
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import numpy.fft as fft
27
28sind = lambda x: np.sin(x*np.pi/180.)
29cosd = lambda x: np.cos(x*np.pi/180.)
30tand = lambda x: np.tan(x*np.pi/180.)
31asind = lambda x: 180.*np.arcsin(x)/np.pi
32acosd = lambda x: 180.*np.arccos(x)/np.pi
33atand = lambda x: 180.*np.arctan(x)/np.pi
34atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35
36def 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   
148def 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
158def 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
167def FillAtomLookUp(atomData):
168    atomLookUp = {}
169    for iatm,atom in enumerate(atomData):
170        atomLookUp[atom[-1]] = iatm
171    return atomLookUp
172
173def GetAtomsById(atomData,atomLookUp,IdList):
174    atoms = []
175    for id in IdList:
176        atoms.append(atomData[atomLookUp[id]])
177    return atoms
178   
179def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
180    Items = []
181    if not isinstance(IdList,list):
182        IdList = [IdList,]
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
189   
190def 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
199
200def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
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]
205               
206def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
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   
224       
225def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
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:
238        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
239    for atom in atomData:
240        XYZ = np.inner(Amat,atom[cx:cx+3])
241        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
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)
249            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
250
251def 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
257def getAtomXYZ(atoms,cx):
258    XYZ = []
259    for atom in atoms:
260        XYZ.append(atom[cx:cx+3])
261    return np.array(XYZ)
262
263def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
264    ''' Returns crystal coordinates for atoms described by RBObj
265    '''
266    RBRes = RBData[RBType][RBObj['RBId']]
267    if RBType == 'Vector':
268        vecs = RBRes['rbVect']
269        mags = RBRes['VectMag']
270        Cart = np.zeros_like(vecs[0])
271        for vec,mag in zip(vecs,mags):
272            Cart += vec*mag
273    elif RBType == 'Residue':
274        Cart = np.array(RBRes['rbXYZ'])
275        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
276            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
277            for ride in seq[3]:
278                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
279    XYZ = np.zeros_like(Cart)
280    for i,xyz in enumerate(Cart):
281        X = prodQVQ(RBObj['Orient'][0],xyz)
282        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
283    return XYZ,Cart
284
285def UpdateMCSAxyz(Bmat,MCSA):
286    xyz = []
287    atTypes = []
288    iatm = 0
289    for model in MCSA['Models'][1:]:        #skip the MD model
290        if model['Type'] == 'Atom':
291            xyz.append(model['Pos'][0])
292            atTypes.append(model['atType'])
293            iatm += 1
294        else:
295            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
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]]
310            if model['MolCent'][1]:
311                Cart -= model['MolCent'][0]
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
318   
319def 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   
338def 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]
342    T = np.zeros(6)
343    L = np.zeros(6)
344    S = np.zeros(8)
345    if 'T' in TLStype:
346        T = TLS[:6]
347    if 'L' in TLStype:
348        L = np.array(TLS[6:12])*(np.pi/180.)**2
349    if 'S' in TLStype:
350        S = np.array(TLS[12:])*(np.pi/180.)
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,
353        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
354    Uout = []
355    Q = RBObj['Orient'][0]
356    for X in Cart:
357        X = prodQVQ(Q,X)
358        if 'U' in TLStype:
359            Uout.append(['I',TLS[0],0,0,0,0,0,0])
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.0*L[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.0*L[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.0*L[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)
372            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
373            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
374        else:
375            Uout.append(['N',])
376    return Uout
377   
378def GetSHCoeff(pId,parmDict,SHkeys):
379    SHCoeff = {}
380    for shkey in SHkeys:
381        shname = str(pId)+'::'+shkey
382        SHCoeff[shkey] = parmDict[shname]
383    return SHCoeff
384       
385def 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
391def getDensity(generalData):
392   
393    mass = getMass(generalData)
394    Volume = generalData['Cell'][7]
395    density = mass/(0.6022137*Volume)
396    return density,Volume/mass
397   
398def 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   
418def getRestDist(XYZ,Amat):
419    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
420   
421def getRestDeriv(Func,XYZ,Amat,ops,SGData):
422    deriv = np.zeros((len(XYZ),3))
423    dx = 0.00001
424    for j,xyz in enumerate(XYZ):
425        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
426            XYZ[j] -= x
427            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
428            XYZ[j] += 2*x
429            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
430            XYZ[j] -= x
431            deriv[j][i] = (d1-d2)/(2*dx)
432    return deriv.flatten()
433
434def 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   
449def 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   
464def getRestChiral(XYZ,Amat):   
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)
470   
471def getRestTorsion(XYZ,Amat):
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 = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
484    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
485    return TOR
486   
487def calcTorsionEnergy(TOR,Coeff=[]):
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
495        val = cof[0]*np.exp(term/1000.0)
496        pMax = cof[0][np.argmin(val)]
497        Eval = np.sum(val)
498        sum = Eval-pMax
499    return sum,Eval
500
501def 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]])):
506            XYZ[j] -= x
507            tor = getRestTorsion(XYZ,Amat)
508            p,d1 = calcTorsionEnergy(tor,Coeff)
509            XYZ[j] += 2*x
510            tor = getRestTorsion(XYZ,Amat)
511            p,d2 = calcTorsionEnergy(tor,Coeff)           
512            XYZ[j] -= x
513            deriv[j][i] = (d2-d1)/(2*dx)
514    return deriv.flatten()
515
516def getRestRama(XYZ,Amat):
517    phi = getRestTorsion(XYZ[:5],Amat)
518    psi = getRestTorsion(XYZ[1:],Amat)
519    return phi,psi
520   
521def calcRamaEnergy(phi,psi,Coeff=[]):
522    sum = 0.
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]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
532        val = cof[0]*np.exp(val/1000.)
533        pMax = cof[0][np.argmin(val)]
534        Eval = np.sum(val)
535        sum = Eval-pMax
536    return sum,Eval
537
538def 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]])):
543            XYZ[j] -= x
544            phi,psi = getRestRama(XYZ,Amat)
545            p,d1 = calcRamaEnergy(phi,psi,Coeff)
546            XYZ[j] += 2*x
547            phi,psi = getRestRama(XYZ,Amat)
548            p,d2 = calcRamaEnergy(phi,psi,Coeff)
549            XYZ[j] -= x
550            deriv[j][i] = (d2-d1)/(2*dx)
551    return deriv.flatten()
552
553def getRestPolefig(ODFln,SamSym,Grid):
554    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
555    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
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
562def getRestPolefigDerv(HKL,Grid,SHCoeff):
563    pass
564       
565def 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]:
579        Oxyz[i] -= dx
580        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
581        Oxyz[i] += 2*dx
582        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
583        Oxyz[i] -= dx
584        Txyz[i] -= dx
585        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
586        Txyz[i] += 2*dx
587        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
588        Txyz[i] -= dx
589    return deriv
590   
591def 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]:
630            OxA[i] -= dx
631            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
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
635           
636            TxA[i] -= dx
637            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
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
641           
642            TxB[i] -= dx
643            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
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
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
655def 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
700def 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
750def 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((P12*P23-P13)/(np.sqrt(1.-P12**2)*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
795            Oatoms[ia][ix+1] -= dx
796            a0 = calcTorsion(Oatoms,SyOps,Amat)
797            Oatoms[ia][ix+1] += 2*dx
798            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
799            Oatoms[ia][ix+1] -= dx           
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       
809def 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((P12*P23-P13)/(np.sqrt(1.-P12**2)*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   
911def 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(esd*10**(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           
933def 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
943def 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
980def 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
1042def 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               
1068def 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 = []
1086    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1087    for F in Flist:
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.
1096        for H,phi in zip(Uniq,Phi)[1:]:
1097            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
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)
1103            Dphi.append((dang+.5) % 1.0)
1104        if i > 20 or len(DH) > 30:
1105            break
1106        i += 1
1107    DH = np.array(DH)
1108    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1109    Dphi = np.array(Dphi)
1110    steps = np.array(hklShape)
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()))
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)
1116#    for i,item in enumerate(hist[:10]):
1117#        print item,bins[i]
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])
1121#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1122    return DX
1123   
1124def 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.1*Vol    #est. no. atoms for ~10A**3/atom
1148            if FFtable:
1149                SQ = 0.25/dsp**2
1150                ff *= G2el.ScatFac(FFtable,SQ)[0]
1151            if ref[8] > 0.:         #use only +ve Fobs**2
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)
1179        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
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   
1205def 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
1231        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1232#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1233        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
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   
1306def 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
1320
1321def 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])
1338        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1339#        for x in xyzs: print x
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
1347               
1348def PeaksUnique(data,Ind):
1349#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
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   
1381def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1382    ind = 0
1383    if useFit:
1384        ind = 1
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']:
1388            ins[x] = Parms[x][ind]
1389        if ifQ:                              #qplot - convert back to 2-theta
1390            pos = 2.0*asind(pos*wave/(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:
1395        if ifQ:
1396            dsp = 2.*np.pi/pos
1397            pos = Parms['difC']*dsp
1398        else:
1399            dsp = pos/Parms['difC'][1]
1400        if 'Pdabc' in Parms2:
1401            for x in ['sig-0','sig-1','X','Y']:
1402                ins[x] = Parms[x][ind]
1403            Pdabc = Parms2['Pdabc'].T
1404            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1405            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1406        else:
1407            for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
1408                ins[x] = Parms[x][ind]
1409            alp = ins['alpha']/dsp
1410            bet = ins['beta-0']+ins['beta-1']/dsp**4
1411        sig = ins['sig-0']+ins['sig-1']*dsp**2
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
1415
1416def 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
1428       
1429def getWave(Parms):
1430    try:
1431        return Parms['Lam'][1]
1432    except KeyError:
1433        return Parms['Lam1'][1]
1434   
1435def 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   
1446def 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   
1453def 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   
1460def 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   
1479def Q2Mat(Q):
1480    ''' make rotation matrix from quaternion
1481        q=r+ai+bj+ck
1482    '''
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
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   
1499def AV2Q(A,V):
1500    ''' convert angle (radians) & vector to quaternion
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:
1508        V = np.array([0.,0.,1.])
1509    p = A/2.
1510    Q[0] = np.cos(p)
1511    Q[1:4] = V*np.sin(p)
1512    return Q
1513   
1514def AVdeg2Q(A,V):
1515    ''' convert angle (degrees) & vector to quaternion
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:
1523        V = np.array([0.,0.,1.])
1524    p = A/2.
1525    Q[0] = cosd(p)
1526    Q[1:4] = V*sind(p)
1527    return Q
1528   
1529def 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])
1534    V = np.array(Q[1:])
1535    d = np.sqrt(np.sum(V**2))
1536    if d:
1537        V /= d
1538    else:
1539        A = 0.
1540        V = np.array([0.,0.,0.])
1541    return A,V
1542   
1543def 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])
1548    V = np.array(Q[1:])
1549    d = np.sqrt(np.sum(V**2))
1550    if d:
1551        V /= d
1552    else:
1553        A = 0.
1554        V = np.array([0.,0.,0.])
1555    return A,V
1556   
1557def 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
1561        q=r+ai+bj+ck
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:
1571        V3 = np.zeros(3)
1572    Q = np.array([0.,0.,0.,1.])
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
1589   
Note: See TracBrowser for help on using the repository browser.