source: trunk/GSASIImath.py @ 763

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

change peak search algorithm - now much faster & more complete
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 70.6 KB
Line 
1<<<<<<< .mine
2# -*- coding: utf-8 -*-
3#GSASIImath - major mathematics routines
4########### SVN repository information ###################
5# $Date: 2012-09-26 06:34:54 +0000 (Wed, 26 Sep 2012) $
6# $Author: vondreele $
7# $Revision: 763 $
8# $URL: trunk/GSASIImath.py $
9# $Id: GSASIImath.py 763 2012-09-26 06:34:54Z vondreele $
10########### SVN repository information ###################
11import sys
12import os
13import os.path as ospath
14import random as rn
15import numpy as np
16import numpy.linalg as nl
17import numpy.ma as ma
18import cPickle
19import time
20import math
21import copy
22import GSASIIpath
23GSASIIpath.SetVersionNumber("$Revision: 763 $")
24import GSASIIElem as G2el
25import GSASIIlattice as G2lat
26import GSASIIspc as G2spc
27import numpy.fft as fft
28
29sind = lambda x: np.sin(x*np.pi/180.)
30cosd = lambda x: np.cos(x*np.pi/180.)
31tand = lambda x: np.tan(x*np.pi/180.)
32asind = lambda x: 180.*np.arcsin(x)/np.pi
33acosd = lambda x: 180.*np.arccos(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 getMass(generalData):
159    mass = 0.
160    for i,elem in enumerate(generalData['AtomTypes']):
161        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
162    return mass   
163
164def getDensity(generalData):
165   
166    mass = getMass(generalData)
167    Volume = generalData['Cell'][7]
168    density = mass/(0.6022137*Volume)
169    return density,Volume/mass
170   
171def getRestDist(XYZ,Amat):
172    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
173
174def getRestAngle(XYZ,Amat):
175   
176    def calcVec(Ox,Tx,Amat):
177        return np.inner(Amat,(Tx-Ox))
178
179    VecA = calcVec(XYZ[1],XYZ[0],Amat)
180    VecA /= np.sqrt(np.sum(VecA**2))
181    VecB = calcVec(XYZ[1],XYZ[2],Amat)
182    VecB /= np.sqrt(np.sum(VecB**2))
183    edge = VecB-VecA
184    edge = np.sum(edge**2)
185    angle = (2.-edge)/2.
186    angle = max(angle,-1.)
187    return acosd(angle)
188   
189def getRestPlane(XYZ,Amat):
190    sumXYZ = np.zeros(3)
191    for xyz in XYZ:
192        sumXYZ += xyz
193    sumXYZ /= len(XYZ)
194    XYZ = np.array(XYZ)-sumXYZ
195    XYZ = np.inner(Amat,XYZ).T
196    Zmat = np.zeros((3,3))
197    for i,xyz in enumerate(XYZ):
198        Zmat += np.outer(xyz.T,xyz)
199    Evec,Emat = nl.eig(Zmat)
200    Evec = np.sqrt(Evec)/(len(XYZ)-3)
201    Order = np.argsort(Evec)
202    return Evec[Order[0]]
203   
204def getRestChiral(XYZ,Amat):
205   
206    VecA = np.empty((3,3))   
207    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
208    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
209    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
210    return nl.det(VecA)
211       
212def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
213   
214    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
215        TxT = inv*(np.inner(M,Tx)+T)+C+U
216        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
217       
218    inv = Top/abs(Top)
219    cent = abs(Top)/100
220    op = abs(Top)%100-1
221    M,T = SGData['SGOps'][op]
222    C = SGData['SGCen'][cent]
223    dx = .00001
224    deriv = np.zeros(6)
225    for i in [0,1,2]:
226        Oxyz[i] += dx
227        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
228        Oxyz[i] -= 2*dx
229        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
230        Oxyz[i] += dx
231        Txyz[i] += dx
232        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
233        Txyz[i] -= 2*dx
234        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
235        Txyz[i] += dx
236    return deriv
237   
238def getAngSig(VA,VB,Amat,SGData,covData={}):
239   
240    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
241        TxT = inv*(np.inner(M,Tx)+T)+C
242        TxT = G2spc.MoveToUnitCell(TxT)+U
243        return np.inner(Amat,(TxT-Ox))
244       
245    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
246        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
247        VecA /= np.sqrt(np.sum(VecA**2))
248        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
249        VecB /= np.sqrt(np.sum(VecB**2))
250        edge = VecB-VecA
251        edge = np.sum(edge**2)
252        angle = (2.-edge)/2.
253        angle = max(angle,-1.)
254        return acosd(angle)
255       
256    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
257    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
258    invA = invB = 1
259    invA = TopA/abs(TopA)
260    invB = TopB/abs(TopB)
261    centA = abs(TopA)/100
262    centB = abs(TopB)/100
263    opA = abs(TopA)%100-1
264    opB = abs(TopB)%100-1
265    MA,TA = SGData['SGOps'][opA]
266    MB,TB = SGData['SGOps'][opB]
267    CA = SGData['SGCen'][centA]
268    CB = SGData['SGCen'][centB]
269    if 'covMatrix' in covData:
270        covMatrix = covData['covMatrix']
271        varyList = covData['varyList']
272        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
273        dx = .00001
274        dadx = np.zeros(9)
275        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
276        for i in [0,1,2]:
277            OxA[i] += dx
278            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
279            OxA[i] -= 2*dx
280            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
281            OxA[i] += dx
282           
283            TxA[i] += dx
284            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
285            TxA[i] -= 2*dx
286            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
287            TxA[i] += dx
288           
289            TxB[i] += dx
290            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
291            TxB[i] -= 2*dx
292            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
293            TxB[i] += dx
294           
295        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
296        if sigAng < 0.01:
297            sigAng = 0.0
298        return Ang,sigAng
299    else:
300        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
301
302def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
303
304    def calcDist(Atoms,SyOps,Amat):
305        XYZ = []
306        for i,atom in enumerate(Atoms):
307            Inv,M,T,C,U = SyOps[i]
308            XYZ.append(np.array(atom[1:4]))
309            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
310            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
311        V1 = XYZ[1]-XYZ[0]
312        return np.sqrt(np.sum(V1**2))
313       
314    Inv = []
315    SyOps = []
316    names = []
317    for i,atom in enumerate(Oatoms):
318        names += atom[-1]
319        Op,unit = Atoms[i][-1]
320        inv = Op/abs(Op)
321        m,t = SGData['SGOps'][abs(Op)%100-1]
322        c = SGData['SGCen'][abs(Op)/100]
323        SyOps.append([inv,m,t,c,unit])
324    Dist = calcDist(Oatoms,SyOps,Amat)
325   
326    sig = -0.001
327    if 'covMatrix' in covData:
328        parmNames = []
329        dx = .00001
330        dadx = np.zeros(6)
331        for i in range(6):
332            ia = i/3
333            ix = i%3
334            Oatoms[ia][ix+1] += dx
335            a0 = calcDist(Oatoms,SyOps,Amat)
336            Oatoms[ia][ix+1] -= 2*dx
337            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
338        covMatrix = covData['covMatrix']
339        varyList = covData['varyList']
340        DistVcov = getVCov(names,varyList,covMatrix)
341        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
342        if sig < 0.001:
343            sig = -0.001
344   
345    return Dist,sig
346
347def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
348       
349    def calcAngle(Atoms,SyOps,Amat):
350        XYZ = []
351        for i,atom in enumerate(Atoms):
352            Inv,M,T,C,U = SyOps[i]
353            XYZ.append(np.array(atom[1:4]))
354            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
355            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
356        V1 = XYZ[1]-XYZ[0]
357        V1 /= np.sqrt(np.sum(V1**2))
358        V2 = XYZ[1]-XYZ[2]
359        V2 /= np.sqrt(np.sum(V2**2))
360        V3 = V2-V1
361        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
362        return acosd(cang)
363
364    Inv = []
365    SyOps = []
366    names = []
367    for i,atom in enumerate(Oatoms):
368        names += atom[-1]
369        Op,unit = Atoms[i][-1]
370        inv = Op/abs(Op)
371        m,t = SGData['SGOps'][abs(Op)%100-1]
372        c = SGData['SGCen'][abs(Op)/100]
373        SyOps.append([inv,m,t,c,unit])
374    Angle = calcAngle(Oatoms,SyOps,Amat)
375   
376    sig = -0.01
377    if 'covMatrix' in covData:
378        parmNames = []
379        dx = .00001
380        dadx = np.zeros(9)
381        for i in range(9):
382            ia = i/3
383            ix = i%3
384            Oatoms[ia][ix+1] += dx
385            a0 = calcAngle(Oatoms,SyOps,Amat)
386            Oatoms[ia][ix+1] -= 2*dx
387            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
388        covMatrix = covData['covMatrix']
389        varyList = covData['varyList']
390        AngVcov = getVCov(names,varyList,covMatrix)
391        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
392        if sig < 0.01:
393            sig = -0.01
394   
395    return Angle,sig
396
397def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
398   
399    def calcTorsion(Atoms,SyOps,Amat):
400       
401        XYZ = []
402        for i,atom in enumerate(Atoms):
403            Inv,M,T,C,U = SyOps[i]
404            XYZ.append(np.array(atom[1:4]))
405            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
406            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
407        V1 = XYZ[1]-XYZ[0]
408        V2 = XYZ[2]-XYZ[1]
409        V3 = XYZ[3]-XYZ[2]
410        V1 /= np.sqrt(np.sum(V1**2))
411        V2 /= np.sqrt(np.sum(V2**2))
412        V3 /= np.sqrt(np.sum(V3**2))
413        M = np.array([V1,V2,V3])
414        D = nl.det(M)
415        Ang = 1.0
416        P12 = np.dot(V1,V2)
417        P13 = np.dot(V1,V3)
418        P23 = np.dot(V2,V3)
419        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
420        return Tors
421           
422    Inv = []
423    SyOps = []
424    names = []
425    for i,atom in enumerate(Oatoms):
426        names += atom[-1]
427        Op,unit = Atoms[i][-1]
428        inv = Op/abs(Op)
429        m,t = SGData['SGOps'][abs(Op)%100-1]
430        c = SGData['SGCen'][abs(Op)/100]
431        SyOps.append([inv,m,t,c,unit])
432    Tors = calcTorsion(Oatoms,SyOps,Amat)
433   
434    sig = -0.01
435    if 'covMatrix' in covData:
436        parmNames = []
437        dx = .00001
438        dadx = np.zeros(12)
439        for i in range(12):
440            ia = i/3
441            ix = i%3
442            Oatoms[ia][ix+1] += dx
443            a0 = calcTorsion(Oatoms,SyOps,Amat)
444            Oatoms[ia][ix+1] -= 2*dx
445            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
446        covMatrix = covData['covMatrix']
447        varyList = covData['varyList']
448        TorVcov = getVCov(names,varyList,covMatrix)
449        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
450        if sig < 0.01:
451            sig = -0.01
452   
453    return Tors,sig
454       
455def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
456   
457    def calcDist(Atoms,SyOps,Amat):
458        XYZ = []
459        for i,atom in enumerate(Atoms):
460            Inv,M,T,C,U = SyOps[i]
461            XYZ.append(np.array(atom[1:4]))
462            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
463            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
464        V1 = XYZ[1]-XYZ[0]
465        return np.sqrt(np.sum(V1**2))
466       
467    def calcAngle(Atoms,SyOps,Amat):
468        XYZ = []
469        for i,atom in enumerate(Atoms):
470            Inv,M,T,C,U = SyOps[i]
471            XYZ.append(np.array(atom[1:4]))
472            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
473            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
474        V1 = XYZ[1]-XYZ[0]
475        V1 /= np.sqrt(np.sum(V1**2))
476        V2 = XYZ[1]-XYZ[2]
477        V2 /= np.sqrt(np.sum(V2**2))
478        V3 = V2-V1
479        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
480        return acosd(cang)
481
482    def calcTorsion(Atoms,SyOps,Amat):
483       
484        XYZ = []
485        for i,atom in enumerate(Atoms):
486            Inv,M,T,C,U = SyOps[i]
487            XYZ.append(np.array(atom[1:4]))
488            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
489            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
490        V1 = XYZ[1]-XYZ[0]
491        V2 = XYZ[2]-XYZ[1]
492        V3 = XYZ[3]-XYZ[2]
493        V1 /= np.sqrt(np.sum(V1**2))
494        V2 /= np.sqrt(np.sum(V2**2))
495        V3 /= np.sqrt(np.sum(V3**2))
496        M = np.array([V1,V2,V3])
497        D = nl.det(M)
498        Ang = 1.0
499        P12 = np.dot(V1,V2)
500        P13 = np.dot(V1,V3)
501        P23 = np.dot(V2,V3)
502        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
503        return Tors
504           
505    Inv = []
506    SyOps = []
507    names = []
508    for i,atom in enumerate(Oatoms):
509        names += atom[-1]
510        Op,unit = Atoms[i][-1]
511        inv = Op/abs(Op)
512        m,t = SGData['SGOps'][abs(Op)%100-1]
513        c = SGData['SGCen'][abs(Op)/100]
514        SyOps.append([inv,m,t,c,unit])
515    M = len(Oatoms)
516    if M == 2:
517        Val = calcDist(Oatoms,SyOps,Amat)
518    elif M == 3:
519        Val = calcAngle(Oatoms,SyOps,Amat)
520    else:
521        Val = calcTorsion(Oatoms,SyOps,Amat)
522   
523    sigVals = [-0.001,-0.01,-0.01]
524    sig = sigVals[M-3]
525    if 'covMatrix' in covData:
526        parmNames = []
527        dx = .00001
528        N = M*3
529        dadx = np.zeros(N)
530        for i in range(N):
531            ia = i/3
532            ix = i%3
533            Oatoms[ia][ix+1] += dx
534            if M == 2:
535                a0 = calcDist(Oatoms,SyOps,Amat)
536            elif M == 3:
537                a0 = calcAngle(Oatoms,SyOps,Amat)
538            else:
539                a0 = calcTorsion(Oatoms,SyOps,Amat)
540            Oatoms[ia][ix+1] -= 2*dx
541            if M == 2:
542                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
543            elif M == 3:
544                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
545            else:
546                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
547        covMatrix = covData['covMatrix']
548        varyList = covData['varyList']
549        Vcov = getVCov(names,varyList,covMatrix)
550        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
551        if sig < sigVals[M-3]:
552            sig = sigVals[M-3]
553   
554    return Val,sig
555       
556   
557def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
558    # returns value(esd) string; nTZ=True for no trailing zeros
559    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
560    #get the 2 significant digits in the esd
561    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
562    #get the number of digits to represent them
563    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
564   
565    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
566    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
567    if esd > 0:
568        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
569        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
570    elif esd < 0:
571         return str(round(value,mdec(esd)-1))
572    else:
573        text = str("%f"%(value))
574        if nTZ:
575            return text.rstrip('0')
576        else:
577            return text
578           
579def adjHKLmax(SGData,Hmax):
580    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
581        Hmax[0] = ((Hmax[0]+3)/6)*6
582        Hmax[1] = ((Hmax[1]+3)/6)*6
583        Hmax[2] = ((Hmax[2]+1)/4)*4
584    else:
585        Hmax[0] = ((Hmax[0]+2)/4)*4
586        Hmax[1] = ((Hmax[1]+2)/4)*4
587        Hmax[2] = ((Hmax[2]+1)/4)*4
588
589def FourierMap(data,reflData):
590   
591    generalData = data['General']
592    if not generalData['Map']['MapType']:
593        print '**** ERROR - Fourier map not defined'
594        return
595    mapData = generalData['Map']
596    dmin = mapData['Resolution']
597    SGData = generalData['SGData']
598    cell = generalData['Cell'][1:8]       
599    A = G2lat.cell2A(cell[:6])
600    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
601    adjHKLmax(SGData,Hmax)
602    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
603#    Fhkl[0,0,0] = generalData['F000X']
604    time0 = time.time()
605    for ref in reflData:
606        if ref[4] >= dmin:
607            Fosq,Fcsq,ph = ref[8:11]
608            for i,hkl in enumerate(ref[11]):
609                hkl = np.asarray(hkl,dtype='i')
610                dp = 360.*ref[12][i]
611                a = cosd(ph+dp)
612                b = sind(ph+dp)
613                phasep = complex(a,b)
614                phasem = complex(a,-b)
615                if 'Fobs' in mapData['MapType']:
616                    F = np.sqrt(Fosq)
617                    h,k,l = hkl+Hmax
618                    Fhkl[h,k,l] = F*phasep
619                    h,k,l = -hkl+Hmax
620                    Fhkl[h,k,l] = F*phasem
621                elif 'Fcalc' in mapData['MapType']:
622                    F = np.sqrt(Fcsq)
623                    h,k,l = hkl+Hmax
624                    Fhkl[h,k,l] = F*phasep
625                    h,k,l = -hkl+Hmax
626                    Fhkl[h,k,l] = F*phasem
627                elif 'delt-F' in mapData['MapType']:
628                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
629                    h,k,l = hkl+Hmax
630                    Fhkl[h,k,l] = dF*phasep
631                    h,k,l = -hkl+Hmax
632                    Fhkl[h,k,l] = dF*phasem
633                elif '2*Fo-Fc' in mapData['MapType']:
634                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
635                    h,k,l = hkl+Hmax
636                    Fhkl[h,k,l] = F*phasep
637                    h,k,l = -hkl+Hmax
638                    Fhkl[h,k,l] = F*phasem
639                elif 'Patterson' in mapData['MapType']:
640                    h,k,l = hkl+Hmax
641                    Fhkl[h,k,l] = complex(Fosq,0.)
642                    h,k,l = -hkl+Hmax
643                    Fhkl[h,k,l] = complex(Fosq,0.)
644    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
645    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
646    mapData['rho'] = np.real(rho)
647    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
648    return mapData
649   
650# map printing for testing purposes
651def printRho(SGLaue,rho,rhoMax):                         
652    dim = len(rho.shape)
653    if dim == 2:
654        ix,jy = rho.shape
655        for j in range(jy):
656            line = ''
657            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
658                line += (jy-j)*'  '
659            for i in range(ix):
660                r = int(100*rho[i,j]/rhoMax)
661                line += '%4d'%(r)
662            print line+'\n'
663    else:
664        ix,jy,kz = rho.shape
665        for k in range(kz):
666            print 'k = ',k
667            for j in range(jy):
668                line = ''
669                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
670                    line += (jy-j)*'  '
671                for i in range(ix):
672                    r = int(100*rho[i,j,k]/rhoMax)
673                    line += '%4d'%(r)
674                print line+'\n'
675## keep this
676               
677def findOffset(SGData,A,Fhkl):   
678    if SGData['SpGrp'] == 'P 1':
679        return [0,0,0]   
680    hklShape = Fhkl.shape
681    steps = np.array(hklShape)
682    Hmax = 2*np.asarray(G2lat.getHKLmax(4.5,SGData,A),dtype='i')
683    Fmax = np.max(np.absolute(Fhkl))
684    hklHalf = np.array(hklShape)/2
685    sortHKL = np.argsort(Fhkl.flatten())
686    Fdict = {}
687    for hkl in sortHKL:
688        HKL = np.unravel_index(hkl,hklShape)
689        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
690        if F == 0.:
691            break
692        Fdict['%.6f'%(np.absolute(F))] = hkl
693    Flist = np.flipud(np.sort(Fdict.keys()))
694    F = str(1.e6)
695    i = 0
696    DH = []
697    Dphi = []
698    while i < 20 and len(DH) < 30:
699        F = Flist[i]
700        hkl = np.unravel_index(Fdict[F],hklShape)
701        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
702        Uniq = np.array(Uniq,dtype='i')
703        Phi = np.array(Phi)
704        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
705        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
706        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
707        ang0 = np.angle(Fh0,deg=True)/360.
708        for j,H in enumerate(Uniq[1:]):
709            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
710            dH = H-hkl
711            dang = ang-ang0
712            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
713                continue
714            DH.append(dH)
715            Dphi.append((dang+0.5) % 1.0)
716        i += 1
717    DH = np.array(DH)
718    print ' map offset no.of terms: %d'%(len(DH))
719    Dphi = np.array(Dphi)
720    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
721    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
722    Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)
723    chisq = np.min(Mmap)
724    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
725    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
726    return DX
727   
728def ChargeFlip(data,reflData,pgbar):
729    generalData = data['General']
730    mapData = generalData['Map']
731    flipData = generalData['Flip']
732    FFtable = {}
733    if 'None' not in flipData['Norm element']:
734        normElem = flipData['Norm element'].upper()
735        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
736        for ff in FFs:
737            if ff['Symbol'] == normElem:
738                FFtable.update(ff)
739    dmin = flipData['Resolution']
740    SGData = generalData['SGData']
741    cell = generalData['Cell'][1:8]       
742    A = G2lat.cell2A(cell[:6])
743    Vol = cell[6]
744    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
745    adjHKLmax(SGData,Hmax)
746    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
747    time0 = time.time()
748    for ref in reflData:
749        dsp = ref[4]
750        if dsp >= dmin:
751            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
752            if FFtable:
753                SQ = 0.25/dsp**2
754                ff *= G2el.ScatFac(FFtable,SQ)[0]
755            if ref[8] > 0.:
756                E = np.sqrt(ref[8])/ff
757            else:
758                E = 0.
759            ph = ref[10]
760            ph = rn.uniform(0.,360.)
761            for i,hkl in enumerate(ref[11]):
762                hkl = np.asarray(hkl,dtype='i')
763                dp = 360.*ref[12][i]
764                a = cosd(ph+dp)
765                b = sind(ph+dp)
766                phasep = complex(a,b)
767                phasem = complex(a,-b)
768                h,k,l = hkl+Hmax
769                Ehkl[h,k,l] = E*phasep
770                h,k,l = -hkl+Hmax       #Friedel pair refl.
771                Ehkl[h,k,l] = E*phasem
772#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
773    CEhkl = copy.copy(Ehkl)
774    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
775    Emask = ma.getmask(MEhkl)
776    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
777    Ncyc = 0
778    old = np.seterr(all='raise')
779    while True:       
780        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
781        CEsig = np.std(CErho)
782        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
783        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
784        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
785        phase = CFhkl/np.absolute(CFhkl)
786        CEhkl = np.absolute(Ehkl)*phase
787        Ncyc += 1
788        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
789        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
790        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
791        if Rcf < 5.:
792            break
793        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
794        if not GoOn or Ncyc > 10000:
795            break
796    np.seterr(**old)
797    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
798    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
799    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
800    roll = findOffset(SGData,A,CEhkl)
801       
802    mapData['Rcf'] = Rcf
803    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
804    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
805    mapData['rollMap'] = [0,0,0]
806    return mapData
807   
808def SearchMap(data):
809    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
810   
811    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
812   
813    def noDuplicate(xyz,peaks,Amat):
814        XYZ = np.inner(Amat,xyz)
815        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
816            print ' Peak',xyz,' <0.5A from another peak'
817            return False
818        return True
819                           
820    def fixSpecialPos(xyz,SGData,Amat):
821        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
822        X = []
823        xyzs = [equiv[0] for equiv in equivs]
824        for x in xyzs:
825            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
826                X.append(x)
827        if len(X) > 1:
828            return np.average(X,axis=0)
829        else:
830            return xyz
831       
832    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
833        Mag,x0,y0,z0,sig = parms
834        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
835       
836    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
837        Mag,x0,y0,z0,sig = parms
838        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
839        return M
840       
841    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
842        Mag,x0,y0,z0,sig = parms
843        dMdv = np.zeros(([5,]+list(rX.shape)))
844        delt = .01
845        for i in range(5):
846            parms[i] -= delt
847            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
848            parms[i] += 2.*delt
849            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
850            parms[i] -= delt
851            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
852        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
853        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
854        dMdv = np.reshape(dMdv,(5,rX.size))
855        Hess = np.inner(dMdv,dMdv)
856       
857        return Vec,Hess
858       
859    generalData = data['General']
860    phaseName = generalData['Name']
861    SGData = generalData['SGData']
862    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
863    drawingData = data['Drawing']
864    peaks = []
865    mags = []
866    dzeros = []
867    try:
868        mapData = generalData['Map']
869        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
870        rho = copy.copy(mapData['rho'])     #don't mess up original
871        mapHalf = np.array(rho.shape)/2
872        res = mapData['Resolution']
873        incre = np.array(rho.shape,dtype=np.float)
874        step = max(1.0,1./res)+1
875        steps = np.array(3*[step,])
876    except KeyError:
877        print '**** ERROR - Fourier map not defined'
878        return peaks,mags
879    rhoMask = ma.array(rho,mask=(rho<contLevel))
880    indices = (-1,0,1)
881    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
882    for roll in rolls:
883        if np.any(roll):
884            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
885    indx = np.transpose(rhoMask.nonzero())
886    peaks = indx/incre
887    mags = rhoMask[rhoMask.nonzero()]
888    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
889        rho = rollMap(rho,ind)
890        rMM = mapHalf-steps
891        rMP = mapHalf+steps+1
892        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
893        peakInt = np.sum(rhoPeak)*res**3
894        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
895        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
896        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
897            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
898        x1 = result[0]
899        if not np.any(x1 < 0):
900            mag = x1[0]
901            peak = (np.array(x1[1:4])-ind)/incre
902        peak = fixSpecialPos(peak,SGData,Amat)
903        rho = rollMap(rho,-ind)       
904    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
905    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
906   
907def sortArray(data,pos,reverse=False):
908    #data is a list of items
909    #sort by pos in list; reverse if True
910    T = []
911    for i,M in enumerate(data):
912        T.append((M[pos],i))
913    D = dict(zip(T,data))
914    T.sort()
915    if reverse:
916        T.reverse()
917    X = []
918    for key in T:
919        X.append(D[key])
920    return X
921               
922def PeaksUnique(data,Ind):
923
924    def noDuplicate(xyz,peaks,Amat):
925        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
926            return False
927        return True
928                           
929    generalData = data['General']
930    cell = generalData['Cell'][1:7]
931    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
932    A = G2lat.cell2A(cell)
933    SGData = generalData['SGData']
934    mapPeaks = data['Map Peaks']
935    Indx = {}
936    XYZ = {}
937    for ind in Ind:
938        XYZ[ind] = np.array(mapPeaks[ind][1:4])
939        Indx[ind] = True
940    for ind in Ind:
941        if Indx[ind]:
942            xyz = XYZ[ind]
943            for jnd in Ind:
944                if ind != jnd and Indx[jnd]:                       
945                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
946                    xyzs = np.array([equiv[0] for equiv in Equiv])
947                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
948    Ind = []
949    for ind in Indx:
950        if Indx[ind]:
951            Ind.append(ind)
952    return Ind
953   
954def prodQQ(QA,QB):
955    ''' Grassman quaternion product
956        QA,QB quaternions; q=r+ai+bj+ck
957    '''
958    D = np.zeros(4)
959    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
960    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
961    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
962    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
963    return D
964   
965def normQ(QA):
966    ''' get length of quaternion & normalize it
967        q=r+ai+bj+ck
968    '''
969    n = np.sqrt(np.sum(np.array(QA)**2))
970    return QA/n
971   
972def invQ(Q):
973    '''
974        get inverse of quaternion
975        q=r+ai+bj+ck; q* = r-ai-bj-ck
976    '''
977    return Q*np.array([1,-1,-1,-1])
978   
979def prodQVQ(Q,V):
980    ''' compute the quaternion vector rotation qvq-1 = v'
981        q=r+ai+bj+ck
982    '''
983    VP = np.zeros(3)
984    T2 = Q[0]*Q[1]
985    T3 = Q[0]*Q[2]
986    T4 = Q[0]*Q[3]
987    T5 = -Q[1]*Q[1]
988    T6 = Q[1]*Q[2]
989    T7 = Q[1]*Q[3]
990    T8 = -Q[2]*Q[2]
991    T9 = Q[2]*Q[3]
992    T10 = -Q[3]*Q[3]
993    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
994    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
995    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
996    return VP   
997   
998def Q2Mat(Q):
999    ''' make rotation matrix from quaternion
1000        q=r+ai+bj+ck
1001    '''
1002    aa = Q[0]**2
1003    ab = Q[0]*Q[1]
1004    ac = Q[0]*Q[2]
1005    ad = Q[0]*Q[3]
1006    bb = Q[1]**2
1007    bc = Q[1]*Q[2]
1008    bd = Q[1]*Q[3]
1009    cc = Q[2]**2
1010    cd = Q[2]*Q[3]
1011    dd = Q[3]**2
1012    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
1013        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
1014        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
1015    return np.array(M)
1016   
1017def AV2Q(A,V):
1018    ''' convert angle (radians -pi to pi) & vector to quaternion
1019        q=r+ai+bj+ck
1020    '''
1021    Q = np.zeros(4)
1022    d = np.sqrt(np.sum(np.array(V)**2))
1023    if d:
1024        V /= d
1025    else:
1026        return [1.,0.,0.,0.]    #identity
1027    p = A/2.
1028    Q[0] = np.cos(p)
1029    s = np.sin(p)
1030    Q[1:4] = V*s
1031    return Q
1032   
1033def AVdeg2Q(A,V):
1034    ''' convert angle (degrees -180 to 180) & vector to quaternion
1035        q=r+ai+bj+ck
1036    '''
1037    Q = np.zeros(4)
1038    d = np.sqrt(np.sum(np.array(V)**2))
1039    if d:
1040        V /= d
1041    else:
1042        return [1.,0.,0.,0.]    #identity
1043    p = A/2.
1044    Q[0] = cosd(p)
1045    S = sind(p)
1046    Q[1:4] = V*S
1047    return Q
1048   
1049=======
1050# -*- coding: utf-8 -*-
1051#GSASIImath - major mathematics routines
1052########### SVN repository information ###################
1053# $Date: 2012-09-26 06:34:54 +0000 (Wed, 26 Sep 2012) $
1054# $Author: vondreele $
1055# $Revision: 763 $
1056# $URL: trunk/GSASIImath.py $
1057# $Id: GSASIImath.py 763 2012-09-26 06:34:54Z vondreele $
1058########### SVN repository information ###################
1059import sys
1060import os
1061import os.path as ospath
1062import random as rn
1063import numpy as np
1064import numpy.linalg as nl
1065import numpy.ma as ma
1066import cPickle
1067import time
1068import math
1069import copy
1070import GSASIIpath
1071GSASIIpath.SetVersionNumber("$Revision: 763 $")
1072import GSASIIElem as G2el
1073import GSASIIlattice as G2lat
1074import GSASIIspc as G2spc
1075import scipy.optimize as so
1076import scipy.linalg as sl
1077import numpy.fft as fft
1078
1079sind = lambda x: np.sin(x*np.pi/180.)
1080cosd = lambda x: np.cos(x*np.pi/180.)
1081tand = lambda x: np.tan(x*np.pi/180.)
1082asind = lambda x: 180.*np.arcsin(x)/np.pi
1083acosd = lambda x: 180.*np.arccos(x)/np.pi
1084atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
1085
1086def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
1087   
1088    """
1089    Minimize the sum of squares of a set of equations.
1090
1091    ::
1092   
1093                    Nobs
1094        x = arg min(sum(func(y)**2,axis=0))
1095                    y=0
1096
1097    Parameters
1098    ----------
1099    func : callable
1100        should take at least one (possibly length N vector) argument and
1101        returns M floating point numbers.
1102    x0 : ndarray
1103        The starting estimate for the minimization of length N
1104    Hess : callable
1105        A required function or method to compute the weighted vector and Hessian for func.
1106        It must be a symmetric NxN array
1107    args : tuple
1108        Any extra arguments to func are placed in this tuple.
1109    ftol : float
1110        Relative error desired in the sum of squares.
1111    xtol : float
1112        Relative error desired in the approximate solution.
1113    maxcyc : int
1114        The maximum number of cycles of refinement to execute, if -1 refine
1115        until other limits are met (ftol, xtol)
1116
1117    Returns
1118    -------
1119    x : ndarray
1120        The solution (or the result of the last iteration for an unsuccessful
1121        call).
1122    cov_x : ndarray
1123        Uses the fjac and ipvt optional outputs to construct an
1124        estimate of the jacobian around the solution.  ``None`` if a
1125        singular matrix encountered (indicates very flat curvature in
1126        some direction).  This matrix must be multiplied by the
1127        residual standard deviation to get the covariance of the
1128        parameter estimates -- see curve_fit.
1129    infodict : dict
1130        a dictionary of optional outputs with the key s::
1131
1132            - 'fvec' : the function evaluated at the output
1133
1134
1135    Notes
1136    -----
1137
1138    """
1139               
1140    x0 = np.array(x0, ndmin=1)      #might be redundant?
1141    n = len(x0)
1142    if type(args) != type(()):
1143        args = (args,)
1144       
1145    icycle = 0
1146    One = np.ones((n,n))
1147    lam = 0.001
1148    lamMax = lam
1149    nfev = 0
1150    while icycle < maxcyc:
1151        lamMax = max(lamMax,lam)
1152        M = func(x0,*args)
1153        nfev += 1
1154        chisq0 = np.sum(M**2)
1155        Yvec,Amat = Hess(x0,*args)
1156        Adiag = np.sqrt(np.diag(Amat))
1157        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
1158        if np.any(psing):                #hard singularity in matrix
1159            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
1160        Anorm = np.outer(Adiag,Adiag)
1161        Yvec /= Adiag
1162        Amat /= Anorm       
1163        while True:
1164            Lam = np.eye(Amat.shape[0])*lam
1165            Amatlam = Amat*(One+Lam)
1166            try:
1167                Xvec = nl.solve(Amatlam,Yvec)
1168            except nl.LinAlgError:
1169                print 'ouch #1'
1170                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
1171                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
1172            Xvec /= Adiag
1173            M2 = func(x0+Xvec,*args)
1174            nfev += 1
1175            chisq1 = np.sum(M2**2)
1176            if chisq1 > chisq0:
1177                lam *= 10.
1178            else:
1179                x0 += Xvec
1180                lam /= 10.
1181                break
1182        if (chisq0-chisq1)/chisq0 < ftol:
1183            break
1184        icycle += 1
1185    M = func(x0,*args)
1186    nfev += 1
1187    Yvec,Amat = Hess(x0,*args)
1188    try:
1189        Bmat = nl.inv(Amat)
1190        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
1191    except nl.LinAlgError:
1192        print 'ouch #2 linear algebra error in LS'
1193        psing = []
1194        if maxcyc:
1195            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
1196        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
1197   
1198def getVCov(varyNames,varyList,covMatrix):
1199    vcov = np.zeros((len(varyNames),len(varyNames)))
1200    for i1,name1 in enumerate(varyNames):
1201        for i2,name2 in enumerate(varyNames):
1202            try:
1203                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
1204            except ValueError:
1205                vcov[i1][i2] = 0.0
1206    return vcov
1207
1208def getMass(generalData):
1209    mass = 0.
1210    for i,elem in enumerate(generalData['AtomTypes']):
1211        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
1212    return mass   
1213
1214def getDensity(generalData):
1215   
1216    mass = getMass(generalData)
1217    Volume = generalData['Cell'][7]
1218    density = mass/(0.6022137*Volume)
1219    return density,Volume/mass
1220   
1221def getRestDist(XYZ,Amat):
1222    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1223
1224def getRestAngle(XYZ,Amat):
1225   
1226    def calcVec(Ox,Tx,Amat):
1227        return np.inner(Amat,(Tx-Ox))
1228
1229    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1230    VecA /= np.sqrt(np.sum(VecA**2))
1231    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1232    VecB /= np.sqrt(np.sum(VecB**2))
1233    edge = VecB-VecA
1234    edge = np.sum(edge**2)
1235    angle = (2.-edge)/2.
1236    angle = max(angle,-1.)
1237    return acosd(angle)
1238   
1239def getRestPlane(XYZ,Amat):
1240    sumXYZ = np.zeros(3)
1241    for xyz in XYZ:
1242        sumXYZ += xyz
1243    sumXYZ /= len(XYZ)
1244    XYZ = np.array(XYZ)-sumXYZ
1245    XYZ = np.inner(Amat,XYZ).T
1246    Zmat = np.zeros((3,3))
1247    for i,xyz in enumerate(XYZ):
1248        Zmat += np.outer(xyz.T,xyz)
1249    Evec,Emat = nl.eig(Zmat)
1250    Evec = np.sqrt(Evec)/(len(XYZ)-3)
1251    Order = np.argsort(Evec)
1252    return Evec[Order[0]]
1253   
1254def getRestChiral(XYZ,Amat):
1255   
1256    VecA = np.empty((3,3))   
1257    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1258    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
1259    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
1260    return nl.det(VecA)
1261       
1262def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
1263   
1264    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
1265        TxT = inv*(np.inner(M,Tx)+T)+C+U
1266        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
1267       
1268    inv = Top/abs(Top)
1269    cent = abs(Top)/100
1270    op = abs(Top)%100-1
1271    M,T = SGData['SGOps'][op]
1272    C = SGData['SGCen'][cent]
1273    dx = .00001
1274    deriv = np.zeros(6)
1275    for i in [0,1,2]:
1276        Oxyz[i] += dx
1277        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1278        Oxyz[i] -= 2*dx
1279        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1280        Oxyz[i] += dx
1281        Txyz[i] += dx
1282        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1283        Txyz[i] -= 2*dx
1284        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1285        Txyz[i] += dx
1286    return deriv
1287   
1288def getAngSig(VA,VB,Amat,SGData,covData={}):
1289   
1290    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1291        TxT = inv*(np.inner(M,Tx)+T)+C
1292        TxT = G2spc.MoveToUnitCell(TxT)+U
1293        return np.inner(Amat,(TxT-Ox))
1294       
1295    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1296        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1297        VecA /= np.sqrt(np.sum(VecA**2))
1298        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1299        VecB /= np.sqrt(np.sum(VecB**2))
1300        edge = VecB-VecA
1301        edge = np.sum(edge**2)
1302        angle = (2.-edge)/2.
1303        angle = max(angle,-1.)
1304        return acosd(angle)
1305       
1306    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1307    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1308    invA = invB = 1
1309    invA = TopA/abs(TopA)
1310    invB = TopB/abs(TopB)
1311    centA = abs(TopA)/100
1312    centB = abs(TopB)/100
1313    opA = abs(TopA)%100-1
1314    opB = abs(TopB)%100-1
1315    MA,TA = SGData['SGOps'][opA]
1316    MB,TB = SGData['SGOps'][opB]
1317    CA = SGData['SGCen'][centA]
1318    CB = SGData['SGCen'][centB]
1319    if 'covMatrix' in covData:
1320        covMatrix = covData['covMatrix']
1321        varyList = covData['varyList']
1322        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1323        dx = .00001
1324        dadx = np.zeros(9)
1325        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1326        for i in [0,1,2]:
1327            OxA[i] += dx
1328            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1329            OxA[i] -= 2*dx
1330            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
1331            OxA[i] += dx
1332           
1333            TxA[i] += dx
1334            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1335            TxA[i] -= 2*dx
1336            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
1337            TxA[i] += dx
1338           
1339            TxB[i] += dx
1340            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1341            TxB[i] -= 2*dx
1342            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
1343            TxB[i] += dx
1344           
1345        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1346        if sigAng < 0.01:
1347            sigAng = 0.0
1348        return Ang,sigAng
1349    else:
1350        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1351
1352def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1353
1354    def calcDist(Atoms,SyOps,Amat):
1355        XYZ = []
1356        for i,atom in enumerate(Atoms):
1357            Inv,M,T,C,U = SyOps[i]
1358            XYZ.append(np.array(atom[1:4]))
1359            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1360            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1361        V1 = XYZ[1]-XYZ[0]
1362        return np.sqrt(np.sum(V1**2))
1363       
1364    Inv = []
1365    SyOps = []
1366    names = []
1367    for i,atom in enumerate(Oatoms):
1368        names += atom[-1]
1369        Op,unit = Atoms[i][-1]
1370        inv = Op/abs(Op)
1371        m,t = SGData['SGOps'][abs(Op)%100-1]
1372        c = SGData['SGCen'][abs(Op)/100]
1373        SyOps.append([inv,m,t,c,unit])
1374    Dist = calcDist(Oatoms,SyOps,Amat)
1375   
1376    sig = -0.001
1377    if 'covMatrix' in covData:
1378        parmNames = []
1379        dx = .00001
1380        dadx = np.zeros(6)
1381        for i in range(6):
1382            ia = i/3
1383            ix = i%3
1384            Oatoms[ia][ix+1] += dx
1385            a0 = calcDist(Oatoms,SyOps,Amat)
1386            Oatoms[ia][ix+1] -= 2*dx
1387            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1388        covMatrix = covData['covMatrix']
1389        varyList = covData['varyList']
1390        DistVcov = getVCov(names,varyList,covMatrix)
1391        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1392        if sig < 0.001:
1393            sig = -0.001
1394   
1395    return Dist,sig
1396
1397def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1398       
1399    def calcAngle(Atoms,SyOps,Amat):
1400        XYZ = []
1401        for i,atom in enumerate(Atoms):
1402            Inv,M,T,C,U = SyOps[i]
1403            XYZ.append(np.array(atom[1:4]))
1404            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1405            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1406        V1 = XYZ[1]-XYZ[0]
1407        V1 /= np.sqrt(np.sum(V1**2))
1408        V2 = XYZ[1]-XYZ[2]
1409        V2 /= np.sqrt(np.sum(V2**2))
1410        V3 = V2-V1
1411        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1412        return acosd(cang)
1413
1414    Inv = []
1415    SyOps = []
1416    names = []
1417    for i,atom in enumerate(Oatoms):
1418        names += atom[-1]
1419        Op,unit = Atoms[i][-1]
1420        inv = Op/abs(Op)
1421        m,t = SGData['SGOps'][abs(Op)%100-1]
1422        c = SGData['SGCen'][abs(Op)/100]
1423        SyOps.append([inv,m,t,c,unit])
1424    Angle = calcAngle(Oatoms,SyOps,Amat)
1425   
1426    sig = -0.01
1427    if 'covMatrix' in covData:
1428        parmNames = []
1429        dx = .00001
1430        dadx = np.zeros(9)
1431        for i in range(9):
1432            ia = i/3
1433            ix = i%3
1434            Oatoms[ia][ix+1] += dx
1435            a0 = calcAngle(Oatoms,SyOps,Amat)
1436            Oatoms[ia][ix+1] -= 2*dx
1437            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1438        covMatrix = covData['covMatrix']
1439        varyList = covData['varyList']
1440        AngVcov = getVCov(names,varyList,covMatrix)
1441        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1442        if sig < 0.01:
1443            sig = -0.01
1444   
1445    return Angle,sig
1446
1447def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1448   
1449    def calcTorsion(Atoms,SyOps,Amat):
1450       
1451        XYZ = []
1452        for i,atom in enumerate(Atoms):
1453            Inv,M,T,C,U = SyOps[i]
1454            XYZ.append(np.array(atom[1:4]))
1455            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1456            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1457        V1 = XYZ[1]-XYZ[0]
1458        V2 = XYZ[2]-XYZ[1]
1459        V3 = XYZ[3]-XYZ[2]
1460        V1 /= np.sqrt(np.sum(V1**2))
1461        V2 /= np.sqrt(np.sum(V2**2))
1462        V3 /= np.sqrt(np.sum(V3**2))
1463        M = np.array([V1,V2,V3])
1464        D = nl.det(M)
1465        Ang = 1.0
1466        P12 = np.dot(V1,V2)
1467        P13 = np.dot(V1,V3)
1468        P23 = np.dot(V2,V3)
1469        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1470        return Tors
1471           
1472    Inv = []
1473    SyOps = []
1474    names = []
1475    for i,atom in enumerate(Oatoms):
1476        names += atom[-1]
1477        Op,unit = Atoms[i][-1]
1478        inv = Op/abs(Op)
1479        m,t = SGData['SGOps'][abs(Op)%100-1]
1480        c = SGData['SGCen'][abs(Op)/100]
1481        SyOps.append([inv,m,t,c,unit])
1482    Tors = calcTorsion(Oatoms,SyOps,Amat)
1483   
1484    sig = -0.01
1485    if 'covMatrix' in covData:
1486        parmNames = []
1487        dx = .00001
1488        dadx = np.zeros(12)
1489        for i in range(12):
1490            ia = i/3
1491            ix = i%3
1492            Oatoms[ia][ix+1] += dx
1493            a0 = calcTorsion(Oatoms,SyOps,Amat)
1494            Oatoms[ia][ix+1] -= 2*dx
1495            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1496        covMatrix = covData['covMatrix']
1497        varyList = covData['varyList']
1498        TorVcov = getVCov(names,varyList,covMatrix)
1499        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1500        if sig < 0.01:
1501            sig = -0.01
1502   
1503    return Tors,sig
1504       
1505def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1506   
1507    def calcDist(Atoms,SyOps,Amat):
1508        XYZ = []
1509        for i,atom in enumerate(Atoms):
1510            Inv,M,T,C,U = SyOps[i]
1511            XYZ.append(np.array(atom[1:4]))
1512            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1513            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1514        V1 = XYZ[1]-XYZ[0]
1515        return np.sqrt(np.sum(V1**2))
1516       
1517    def calcAngle(Atoms,SyOps,Amat):
1518        XYZ = []
1519        for i,atom in enumerate(Atoms):
1520            Inv,M,T,C,U = SyOps[i]
1521            XYZ.append(np.array(atom[1:4]))
1522            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1523            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1524        V1 = XYZ[1]-XYZ[0]
1525        V1 /= np.sqrt(np.sum(V1**2))
1526        V2 = XYZ[1]-XYZ[2]
1527        V2 /= np.sqrt(np.sum(V2**2))
1528        V3 = V2-V1
1529        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1530        return acosd(cang)
1531
1532    def calcTorsion(Atoms,SyOps,Amat):
1533       
1534        XYZ = []
1535        for i,atom in enumerate(Atoms):
1536            Inv,M,T,C,U = SyOps[i]
1537            XYZ.append(np.array(atom[1:4]))
1538            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1539            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1540        V1 = XYZ[1]-XYZ[0]
1541        V2 = XYZ[2]-XYZ[1]
1542        V3 = XYZ[3]-XYZ[2]
1543        V1 /= np.sqrt(np.sum(V1**2))
1544        V2 /= np.sqrt(np.sum(V2**2))
1545        V3 /= np.sqrt(np.sum(V3**2))
1546        M = np.array([V1,V2,V3])
1547        D = nl.det(M)
1548        Ang = 1.0
1549        P12 = np.dot(V1,V2)
1550        P13 = np.dot(V1,V3)
1551        P23 = np.dot(V2,V3)
1552        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1553        return Tors
1554           
1555    Inv = []
1556    SyOps = []
1557    names = []
1558    for i,atom in enumerate(Oatoms):
1559        names += atom[-1]
1560        Op,unit = Atoms[i][-1]
1561        inv = Op/abs(Op)
1562        m,t = SGData['SGOps'][abs(Op)%100-1]
1563        c = SGData['SGCen'][abs(Op)/100]
1564        SyOps.append([inv,m,t,c,unit])
1565    M = len(Oatoms)
1566    if M == 2:
1567        Val = calcDist(Oatoms,SyOps,Amat)
1568    elif M == 3:
1569        Val = calcAngle(Oatoms,SyOps,Amat)
1570    else:
1571        Val = calcTorsion(Oatoms,SyOps,Amat)
1572   
1573    sigVals = [-0.001,-0.01,-0.01]
1574    sig = sigVals[M-3]
1575    if 'covMatrix' in covData:
1576        parmNames = []
1577        dx = .00001
1578        N = M*3
1579        dadx = np.zeros(N)
1580        for i in range(N):
1581            ia = i/3
1582            ix = i%3
1583            Oatoms[ia][ix+1] += dx
1584            if M == 2:
1585                a0 = calcDist(Oatoms,SyOps,Amat)
1586            elif M == 3:
1587                a0 = calcAngle(Oatoms,SyOps,Amat)
1588            else:
1589                a0 = calcTorsion(Oatoms,SyOps,Amat)
1590            Oatoms[ia][ix+1] -= 2*dx
1591            if M == 2:
1592                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1593            elif M == 3:
1594                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1595            else:
1596                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1597        covMatrix = covData['covMatrix']
1598        varyList = covData['varyList']
1599        Vcov = getVCov(names,varyList,covMatrix)
1600        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1601        if sig < sigVals[M-3]:
1602            sig = sigVals[M-3]
1603   
1604    return Val,sig
1605       
1606   
1607def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
1608    # returns value(esd) string; nTZ=True for no trailing zeros
1609    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
1610    #get the 2 significant digits in the esd
1611    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
1612    #get the number of digits to represent them
1613    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
1614   
1615    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
1616    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
1617    if esd > 0:
1618        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
1619        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
1620    elif esd < 0:
1621         return str(round(value,mdec(esd)-1))
1622    else:
1623        text = str("%f"%(value))
1624        if nTZ:
1625            return text.rstrip('0')
1626        else:
1627            return text
1628           
1629def adjHKLmax(SGData,Hmax):
1630    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1631        Hmax[0] = ((Hmax[0]+3)/6)*6
1632        Hmax[1] = ((Hmax[1]+3)/6)*6
1633        Hmax[2] = ((Hmax[2]+1)/4)*4
1634    else:
1635        Hmax[0] = ((Hmax[0]+2)/4)*4
1636        Hmax[1] = ((Hmax[1]+2)/4)*4
1637        Hmax[2] = ((Hmax[2]+1)/4)*4
1638
1639def FourierMap(data,reflData):
1640   
1641    generalData = data['General']
1642    if not generalData['Map']['MapType']:
1643        print '**** ERROR - Fourier map not defined'
1644        return
1645    mapData = generalData['Map']
1646    dmin = mapData['Resolution']
1647    SGData = generalData['SGData']
1648    cell = generalData['Cell'][1:8]       
1649    A = G2lat.cell2A(cell[:6])
1650    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1651    adjHKLmax(SGData,Hmax)
1652    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1653#    Fhkl[0,0,0] = generalData['F000X']
1654    time0 = time.time()
1655    for ref in reflData:
1656        if ref[4] >= dmin:
1657            Fosq,Fcsq,ph = ref[8:11]
1658            for i,hkl in enumerate(ref[11]):
1659                hkl = np.asarray(hkl,dtype='i')
1660                dp = 360.*ref[12][i]
1661                a = cosd(ph+dp)
1662                b = sind(ph+dp)
1663                phasep = complex(a,b)
1664                phasem = complex(a,-b)
1665                if 'Fobs' in mapData['MapType']:
1666                    F = np.sqrt(Fosq)
1667                    h,k,l = hkl+Hmax
1668                    Fhkl[h,k,l] = F*phasep
1669                    h,k,l = -hkl+Hmax
1670                    Fhkl[h,k,l] = F*phasem
1671                elif 'Fcalc' in mapData['MapType']:
1672                    F = np.sqrt(Fcsq)
1673                    h,k,l = hkl+Hmax
1674                    Fhkl[h,k,l] = F*phasep
1675                    h,k,l = -hkl+Hmax
1676                    Fhkl[h,k,l] = F*phasem
1677                elif 'delt-F' in mapData['MapType']:
1678                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
1679                    h,k,l = hkl+Hmax
1680                    Fhkl[h,k,l] = dF*phasep
1681                    h,k,l = -hkl+Hmax
1682                    Fhkl[h,k,l] = dF*phasem
1683                elif '2*Fo-Fc' in mapData['MapType']:
1684                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1685                    h,k,l = hkl+Hmax
1686                    Fhkl[h,k,l] = F*phasep
1687                    h,k,l = -hkl+Hmax
1688                    Fhkl[h,k,l] = F*phasem
1689                elif 'Patterson' in mapData['MapType']:
1690                    h,k,l = hkl+Hmax
1691                    Fhkl[h,k,l] = complex(Fosq,0.)
1692                    h,k,l = -hkl+Hmax
1693                    Fhkl[h,k,l] = complex(Fosq,0.)
1694    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1695    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1696    mapData['rho'] = np.real(rho)
1697    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1698    return mapData
1699   
1700# map printing for testing purposes
1701def printRho(SGLaue,rho,rhoMax):                         
1702    dim = len(rho.shape)
1703    if dim == 2:
1704        ix,jy = rho.shape
1705        for j in range(jy):
1706            line = ''
1707            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1708                line += (jy-j)*'  '
1709            for i in range(ix):
1710                r = int(100*rho[i,j]/rhoMax)
1711                line += '%4d'%(r)
1712            print line+'\n'
1713    else:
1714        ix,jy,kz = rho.shape
1715        for k in range(kz):
1716            print 'k = ',k
1717            for j in range(jy):
1718                line = ''
1719                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1720                    line += (jy-j)*'  '
1721                for i in range(ix):
1722                    r = int(100*rho[i,j,k]/rhoMax)
1723                    line += '%4d'%(r)
1724                print line+'\n'
1725## keep this
1726               
1727def findOffset(SGData,A,Fhkl):   
1728    if SGData['SpGrp'] == 'P 1':
1729        return [0,0,0]   
1730    hklShape = Fhkl.shape
1731    steps = np.array(hklShape)
1732    Hmax = 2*np.asarray(G2lat.getHKLmax(4.5,SGData,A),dtype='i')
1733    Fmax = np.max(np.absolute(Fhkl))
1734    hklHalf = np.array(hklShape)/2
1735    sortHKL = np.argsort(Fhkl.flatten())
1736    Fdict = {}
1737    for hkl in sortHKL:
1738        HKL = np.unravel_index(hkl,hklShape)
1739        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1740        if F == 0.:
1741            break
1742        Fdict['%.6f'%(np.absolute(F))] = hkl
1743    Flist = np.flipud(np.sort(Fdict.keys()))
1744    F = str(1.e6)
1745    i = 0
1746    DH = []
1747    Dphi = []
1748    while i < 20 and len(DH) < 30:
1749        F = Flist[i]
1750        hkl = np.unravel_index(Fdict[F],hklShape)
1751        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1752        Uniq = np.array(Uniq,dtype='i')
1753        Phi = np.array(Phi)
1754        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1755        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1756        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1757        ang0 = np.angle(Fh0,deg=True)/360.
1758        for j,H in enumerate(Uniq[1:]):
1759            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
1760            dH = H-hkl
1761            dang = ang-ang0
1762            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
1763                continue
1764            DH.append(dH)
1765            Dphi.append((dang+0.5) % 1.0)
1766        i += 1
1767    DH = np.array(DH)
1768    print ' map offset no.of terms: %d'%(len(DH))
1769    Dphi = np.array(Dphi)
1770    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1771    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1772    Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)
1773    chisq = np.min(Mmap)
1774    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1775    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1776    return DX
1777   
1778def ChargeFlip(data,reflData,pgbar):
1779    generalData = data['General']
1780    mapData = generalData['Map']
1781    flipData = generalData['Flip']
1782    FFtable = {}
1783    if 'None' not in flipData['Norm element']:
1784        normElem = flipData['Norm element'].upper()
1785        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1786        for ff in FFs:
1787            if ff['Symbol'] == normElem:
1788                FFtable.update(ff)
1789    dmin = flipData['Resolution']
1790    SGData = generalData['SGData']
1791    cell = generalData['Cell'][1:8]       
1792    A = G2lat.cell2A(cell[:6])
1793    Vol = cell[6]
1794    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1795    adjHKLmax(SGData,Hmax)
1796    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1797    time0 = time.time()
1798    for ref in reflData:
1799        dsp = ref[4]
1800        if dsp >= dmin:
1801            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1802            if FFtable:
1803                SQ = 0.25/dsp**2
1804                ff *= G2el.ScatFac(FFtable,SQ)[0]
1805            if ref[8] > 0.:
1806                E = np.sqrt(ref[8])/ff
1807            else:
1808                E = 0.
1809            ph = ref[10]
1810            ph = rn.uniform(0.,360.)
1811            for i,hkl in enumerate(ref[11]):
1812                hkl = np.asarray(hkl,dtype='i')
1813                dp = 360.*ref[12][i]
1814                a = cosd(ph+dp)
1815                b = sind(ph+dp)
1816                phasep = complex(a,b)
1817                phasem = complex(a,-b)
1818                h,k,l = hkl+Hmax
1819                Ehkl[h,k,l] = E*phasep
1820                h,k,l = -hkl+Hmax       #Friedel pair refl.
1821                Ehkl[h,k,l] = E*phasem
1822#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1823    CEhkl = copy.copy(Ehkl)
1824    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1825    Emask = ma.getmask(MEhkl)
1826    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1827    Ncyc = 0
1828    old = np.seterr(all='raise')
1829    while True:       
1830        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1831        CEsig = np.std(CErho)
1832        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1833        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1834        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1835        phase = CFhkl/np.absolute(CFhkl)
1836        CEhkl = np.absolute(Ehkl)*phase
1837        Ncyc += 1
1838        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1839        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1840        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1841        if Rcf < 5.:
1842            break
1843        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1844        if not GoOn or Ncyc > 10000:
1845            break
1846    np.seterr(**old)
1847    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1848    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1849    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1850    roll = findOffset(SGData,A,CEhkl)
1851       
1852    mapData['Rcf'] = Rcf
1853    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1854    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1855    mapData['rollMap'] = [0,0,0]
1856    return mapData
1857   
1858def SearchMap(data,keepDup=False,Pgbar=None):
1859   
1860    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1861   
1862    def noEquivalent(xyz,peaks,SGData):                  #be careful where this is used - it's slow
1863        equivs = G2spc.GenAtom(xyz,SGData)
1864        xyzs = [equiv[0] for equiv in equivs]
1865        for x in xyzs:
1866            if True in [np.allclose(x,peak,atol=0.02) for peak in peaks]:
1867                return False
1868        return True
1869       
1870    def noDuplicate(xyz,peaks,Amat):
1871        XYZ = np.inner(Amat,xyz)
1872        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1873            print ' Peak',xyz,' <0.5A from another peak'
1874            return False
1875        return True
1876                           
1877    def fixSpecialPos(xyz,SGData,Amat):
1878        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1879        X = []
1880        xyzs = [equiv[0] for equiv in equivs]
1881        for x in xyzs:
1882            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1883                X.append(x)
1884        if len(X) > 1:
1885            return np.average(X,axis=0)
1886        else:
1887            return xyz
1888       
1889    def findRoll(rhoMask,mapHalf):
1890        indx = np.array(ma.nonzero(rhoMask)).T
1891        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
1892        rhoMax = np.max(rhoList)
1893        return mapHalf-indx[np.argmax(rhoList)]
1894       
1895    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1896        Mag,x0,y0,z0,sig = parms
1897#        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1898#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
1899#        else:
1900#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
1901        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
1902       
1903    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1904        Mag,x0,y0,z0,sig = parms
1905        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1906        return M
1907       
1908    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1909        Mag,x0,y0,z0,sig = parms
1910        dMdv = np.zeros(([5,]+list(rX.shape)))
1911        delt = .01
1912        for i in range(5):
1913            parms[i] -= delt
1914            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1915            parms[i] += 2.*delt
1916            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1917            parms[i] -= delt
1918            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1919        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1920        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1921        dMdv = np.reshape(dMdv,(5,rX.size))
1922        Hess = np.inner(dMdv,dMdv)
1923       
1924        return Vec,Hess
1925       
1926    generalData = data['General']
1927    phaseName = generalData['Name']
1928    SGData = generalData['SGData']
1929    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1930    drawingData = data['Drawing']
1931    peaks = []
1932    mags = []
1933    dzeros = []
1934    try:
1935        mapData = generalData['Map']
1936        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1937        rho = copy.copy(mapData['rho'])     #don't mess up original
1938        mapHalf = np.array(rho.shape)/2
1939        res = mapData['Resolution']
1940        incre = np.array(rho.shape)
1941        step = max(1.0,1./res)+1
1942        steps = np.array(3*[step,])
1943    except KeyError:
1944        print '**** ERROR - Fourier map not defined'
1945        return peaks,mags
1946    while True:
1947        rhoMask = ma.array(rho,mask=(rho<contLevel))
1948        if not ma.count(rhoMask):
1949            break
1950        rMI = findRoll(rhoMask,mapHalf)
1951        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
1952        rMM = mapHalf-steps
1953        rMP = mapHalf+steps+1
1954        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1955        peakInt = np.sum(rhoPeak)*res**3
1956        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1957        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1958        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1959            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1960        x1 = result[0]
1961        if np.any(x1 < 0):
1962            break
1963        peak = (np.array(x1[1:4])-rMI)/incre
1964        dzero = np.sqrt(np.sum(np.inner(Amat,peak)**2))
1965        peak = fixSpecialPos(peak,SGData,Amat)
1966        if not len(peaks):
1967            peaks.append(peak)
1968            mags.append(x1[0])
1969            dzeros.append(dzero)
1970        else:
1971            if keepDup:
1972                if noDuplicate(peak,peaks,Amat):
1973                    peaks.append(peak)
1974                    mags.append(x1[0])
1975                    dzeros.append(dzero)
1976            elif noEquivalent(peak,peaks,SGData) and x1[0] > 0.:
1977                peaks.append(peak)
1978                mags.append(x1[0])
1979                dzeros.append(dzero)
1980            GoOn = Pgbar.Update(len(peaks),newmsg='%s%d'%('No. Peaks found =',len(peaks)))[0]
1981            if not GoOn or len(peaks) > 1000:
1982                break
1983        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(x1,rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
1984        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
1985#    if SGData['SGInv']:                 #check origin location & fix it if needed - centrosymmetric only
1986#        Ocheck = np.zeros_like(rho)
1987#        for ipeak in peaks:
1988#            for opeak in peaks:           
1989#                dx = ((opeak-ipeak)*incre)%incre
1990#                if np.any(dx):      #avoid self vector
1991#                    Ocheck[dx[0],dx[1],dx[2]] += 1
1992#        dxMax = np.unravel_index(np.argmax(Ocheck),rho.shape)
1993#        print ' Inversion at:',dxMax,' shift by ;',dxMax-mapHalf
1994#        rho = np.roll(np.roll(np.roll(rho,dxMax[2],axis=2),dxMax[1],axis=1),dxMax[0],axis=0)
1995#        for peak in peaks:
1996#            peak = (peak-(dxMax+mapHalf)/incre)%1.0
1997       
1998    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1999   
2000def sortArray(data,pos,reverse=False):
2001    #data is a list of items
2002    #sort by pos in list; reverse if True
2003    T = []
2004    for i,M in enumerate(data):
2005        T.append((M[pos],i))
2006    D = dict(zip(T,data))
2007    T.sort()
2008    if reverse:
2009        T.reverse()
2010    X = []
2011    for key in T:
2012        X.append(D[key])
2013    return X
2014               
2015def PeaksUnique(data,Ind):
2016
2017    def noDuplicate(xyz,peaks,Amat):
2018        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2019            return False
2020        return True
2021                           
2022    generalData = data['General']
2023    cell = generalData['Cell'][1:7]
2024    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2025    A = G2lat.cell2A(cell)
2026    SGData = generalData['SGData']
2027    mapPeaks = data['Map Peaks']
2028    Indx = {}
2029    XYZ = {}
2030    for ind in Ind:
2031        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2032        Indx[ind] = True
2033    for ind in Ind:
2034        if Indx[ind]:
2035            xyz = XYZ[ind]
2036            for jnd in Ind:
2037                if ind != jnd and Indx[jnd]:                       
2038                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2039                    xyzs = np.array([equiv[0] for equiv in Equiv])
2040                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2041    Ind = []
2042    for ind in Indx:
2043        if Indx[ind]:
2044            Ind.append(ind)
2045    return Ind
2046   
2047def prodQQ(QA,QB):
2048    ''' Grassman quaternion product
2049        QA,QB quaternions; q=r+ai+bj+ck
2050    '''
2051    D = np.zeros(4)
2052    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2053    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2054    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2055    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2056    return D
2057   
2058def normQ(QA):
2059    ''' get length of quaternion & normalize it
2060        q=r+ai+bj+ck
2061    '''
2062    n = np.sqrt(np.sum(np.array(QA)**2))
2063    return QA/n
2064   
2065def invQ(Q):
2066    '''
2067        get inverse of quaternion
2068        q=r+ai+bj+ck; q* = r-ai-bj-ck
2069    '''
2070    return Q*np.array([1,-1,-1,-1])
2071   
2072def prodQVQ(Q,V):
2073    ''' compute the quaternion vector rotation qvq-1 = v'
2074        q=r+ai+bj+ck
2075    '''
2076    VP = np.zeros(3)
2077    T2 = Q[0]*Q[1]
2078    T3 = Q[0]*Q[2]
2079    T4 = Q[0]*Q[3]
2080    T5 = -Q[1]*Q[1]
2081    T6 = Q[1]*Q[2]
2082    T7 = Q[1]*Q[3]
2083    T8 = -Q[2]*Q[2]
2084    T9 = Q[2]*Q[3]
2085    T10 = -Q[3]*Q[3]
2086    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2087    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2088    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2089    return VP   
2090   
2091def Q2Mat(Q):
2092    ''' make rotation matrix from quaternion
2093        q=r+ai+bj+ck
2094    '''
2095    aa = Q[0]**2
2096    ab = Q[0]*Q[1]
2097    ac = Q[0]*Q[2]
2098    ad = Q[0]*Q[3]
2099    bb = Q[1]**2
2100    bc = Q[1]*Q[2]
2101    bd = Q[1]*Q[3]
2102    cc = Q[2]**2
2103    cd = Q[2]*Q[3]
2104    dd = Q[3]**2
2105    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2106        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2107        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2108    return np.array(M)
2109   
2110def AV2Q(A,V):
2111    ''' convert angle (radians -pi to pi) & vector to quaternion
2112        q=r+ai+bj+ck
2113    '''
2114    Q = np.zeros(4)
2115    d = np.sqrt(np.sum(np.array(V)**2))
2116    if d:
2117        V /= d
2118    else:
2119        return [1.,0.,0.,0.]    #identity
2120    p = A/2.
2121    Q[0] = np.cos(p)
2122    s = np.sin(p)
2123    Q[1:4] = V*s
2124    return Q
2125   
2126def AVdeg2Q(A,V):
2127    ''' convert angle (degrees -180 to 180) & vector to quaternion
2128        q=r+ai+bj+ck
2129    '''
2130    Q = np.zeros(4)
2131    d = np.sqrt(np.sum(np.array(V)**2))
2132    if d:
2133        V /= d
2134    else:
2135        return [1.,0.,0.,0.]    #identity
2136    p = A/2.
2137    Q[0] = cosd(p)
2138    S = sind(p)
2139    Q[1:4] = V*S
2140    return Q
2141   
2142>>>>>>> .r762
Note: See TracBrowser for help on using the repository browser.