source: trunk/GSASIImath.py @ 934

Last change on this file since 934 was 934, checked in by vondreele, 10 years ago

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