source: trunk/GSASIImath.py @ 815

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

another hot key
more restraint development
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 39.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2012-12-14 22:28:59 +0000 (Fri, 14 Dec 2012) $
5# $Author: vondreele $
6# $Revision: 815 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 815 2012-12-14 22:28:59Z 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: 815 $")
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
33atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
34
35def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
36   
37    """
38    Minimize the sum of squares of a set of equations.
39
40    ::
41   
42                    Nobs
43        x = arg min(sum(func(y)**2,axis=0))
44                    y=0
45
46    Parameters
47    ----------
48    func : callable
49        should take at least one (possibly length N vector) argument and
50        returns M floating point numbers.
51    x0 : ndarray
52        The starting estimate for the minimization of length N
53    Hess : callable
54        A required function or method to compute the weighted vector and Hessian for func.
55        It must be a symmetric NxN array
56    args : tuple
57        Any extra arguments to func are placed in this tuple.
58    ftol : float
59        Relative error desired in the sum of squares.
60    xtol : float
61        Relative error desired in the approximate solution.
62    maxcyc : int
63        The maximum number of cycles of refinement to execute, if -1 refine
64        until other limits are met (ftol, xtol)
65
66    Returns
67    -------
68    x : ndarray
69        The solution (or the result of the last iteration for an unsuccessful
70        call).
71    cov_x : ndarray
72        Uses the fjac and ipvt optional outputs to construct an
73        estimate of the jacobian around the solution.  ``None`` if a
74        singular matrix encountered (indicates very flat curvature in
75        some direction).  This matrix must be multiplied by the
76        residual standard deviation to get the covariance of the
77        parameter estimates -- see curve_fit.
78    infodict : dict
79        a dictionary of optional outputs with the key s::
80
81            - 'fvec' : the function evaluated at the output
82
83
84    Notes
85    -----
86
87    """
88               
89    x0 = np.array(x0, ndmin=1)      #might be redundant?
90    n = len(x0)
91    if type(args) != type(()):
92        args = (args,)
93       
94    icycle = 0
95    One = np.ones((n,n))
96    lam = 0.001
97    lamMax = lam
98    nfev = 0
99    while icycle < maxcyc:
100        lamMax = max(lamMax,lam)
101        M = func(x0,*args)
102        nfev += 1
103        chisq0 = np.sum(M**2)
104        Yvec,Amat = Hess(x0,*args)
105        Adiag = np.sqrt(np.diag(Amat))
106        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
107        if np.any(psing):                #hard singularity in matrix
108            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
109        Anorm = np.outer(Adiag,Adiag)
110        Yvec /= Adiag
111        Amat /= Anorm       
112        while True:
113            Lam = np.eye(Amat.shape[0])*lam
114            Amatlam = Amat*(One+Lam)
115            try:
116                Xvec = nl.solve(Amatlam,Yvec)
117            except nl.LinAlgError:
118                print 'ouch #1'
119                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
120                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
121            Xvec /= Adiag
122            M2 = func(x0+Xvec,*args)
123            nfev += 1
124            chisq1 = np.sum(M2**2)
125            if chisq1 > chisq0:
126                lam *= 10.
127            else:
128                x0 += Xvec
129                lam /= 10.
130                break
131        if (chisq0-chisq1)/chisq0 < ftol:
132            break
133        icycle += 1
134    M = func(x0,*args)
135    nfev += 1
136    Yvec,Amat = Hess(x0,*args)
137    try:
138        Bmat = nl.inv(Amat)
139        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
140    except nl.LinAlgError:
141        print 'ouch #2 linear algebra error in LS'
142        psing = []
143        if maxcyc:
144            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
145        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
146   
147def getVCov(varyNames,varyList,covMatrix):
148    vcov = np.zeros((len(varyNames),len(varyNames)))
149    for i1,name1 in enumerate(varyNames):
150        for i2,name2 in enumerate(varyNames):
151            try:
152                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
153            except ValueError:
154                vcov[i1][i2] = 0.0
155    return vcov
156
157def FindAtomIndexByIDs(atomData,IDs,Draw=True):
158    indx = []
159    for i,atom in enumerate(atomData):
160        if Draw and atom[-3] in IDs:
161            indx.append(i)
162        elif atom[-1] in IDs:
163            indx.append(i)
164    return indx
165
166def FillAtomLookUp(atomData):
167    atomLookUp = {}
168    for iatm,atom in enumerate(atomData):
169        atomLookUp[atom[-1]] = iatm
170    return atomLookUp
171
172def GetAtomsById(atomData,atomLookUp,IdList):
173    atoms = []
174    for id in IdList:
175        atoms.append(atomData[atomLookUp[id]])
176    return atoms
177   
178def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
179    Items = []
180    if not isinstance(IdList,list):
181        IdList = [IdList,]
182    for id in IdList:
183        if numItems == 1:
184            Items.append(atomData[atomLookUp[id]][itemLoc])
185        else:
186            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
187    return Items
188       
189def getMass(generalData):
190    mass = 0.
191    for i,elem in enumerate(generalData['AtomTypes']):
192        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
193    return mass   
194
195def getDensity(generalData):
196   
197    mass = getMass(generalData)
198    Volume = generalData['Cell'][7]
199    density = mass/(0.6022137*Volume)
200    return density,Volume/mass
201   
202def getSyXYZ(XYZ,ops,SGData):
203    XYZout = np.zeros_like(XYZ)
204    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
205        if op == '1':
206            XYZout[i] = xyz
207        else:
208            oprs = op.split('+')
209            unit = [0,0,0]
210            if oprs[1]:
211                unit = np.array(list(eval(oprs[1])))
212            syop =int(oprs[0])
213            inv = syop/abs(syop)
214            syop *= inv
215            cent = syop/100
216            syop %= 100
217            syop -= 1
218            M,T = SGData['SGOps'][syop]
219            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
220    return XYZout
221   
222def getRestDist(XYZ,Amat):
223    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
224   
225def getRestDeriv(Func,XYZ,Amat):
226    deriv = np.array((3,len(XYZ)))
227    dx = 0.00001
228    for j,xyz in enumerate(XYZ):
229        for i,x in enumerate(xyz):
230            x += dx
231            d1 = Func(XYZ,Amat)
232            x -= 2*dx
233            d2 = Func(XYZ,Amat)
234            x += dx
235            deriv[i][j] = (d1-d2)/(2*dx)
236    return deriv
237
238def getRestAngle(XYZ,Amat):
239   
240    def calcVec(Ox,Tx,Amat):
241        return np.inner(Amat,(Tx-Ox))
242
243    VecA = calcVec(XYZ[1],XYZ[0],Amat)
244    VecA /= np.sqrt(np.sum(VecA**2))
245    VecB = calcVec(XYZ[1],XYZ[2],Amat)
246    VecB /= np.sqrt(np.sum(VecB**2))
247    edge = VecB-VecA
248    edge = np.sum(edge**2)
249    angle = (2.-edge)/2.
250    angle = max(angle,-1.)
251    return acosd(angle)
252   
253def getRestPlane(XYZ,Amat):
254    sumXYZ = np.zeros(3)
255    for xyz in XYZ:
256        sumXYZ += xyz
257    sumXYZ /= len(XYZ)
258    XYZ = np.array(XYZ)-sumXYZ
259    XYZ = np.inner(Amat,XYZ).T
260    Zmat = np.zeros((3,3))
261    for i,xyz in enumerate(XYZ):
262        Zmat += np.outer(xyz.T,xyz)
263    Evec,Emat = nl.eig(Zmat)
264    Evec = np.sqrt(Evec)/(len(XYZ)-3)
265    Order = np.argsort(Evec)
266    return Evec[Order[0]]
267   
268def getRestChiral(XYZ,Amat):   
269    VecA = np.empty((3,3))   
270    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
271    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
272    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
273    return nl.det(VecA)
274   
275def getRestTorsion(XYZ,Amat,Coeff=[]):
276    VecA = np.empty((3,3))
277    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
278    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
279    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
280    D = nl.det(VecA)
281    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
282    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
283    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
284    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
285    Ang = 1.0
286    if abs(P12) < 1.0 and abs(P23) < 1.0:
287        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
288    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
289    sum = 0.
290    if len(Coeff):
291        cof = np.reshape(Coeff,(3,3)).T
292        delt = TOR-cof[1]
293        delt = np.where(delt<-180.,delt+360.,delt)
294        delt = np.where(delt>180.,delt-360.,delt)
295#        pMax = np.min(cof[0])
296        pMax = cof[0][np.argmin(delt)]
297        term = -cof[2]*delt**2
298        sum = np.sum(cof[0]*np.exp(term/1000.0))-pMax
299    return TOR,sum
300   
301def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
302   
303    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
304        TxT = inv*(np.inner(M,Tx)+T)+C+U
305        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
306       
307    inv = Top/abs(Top)
308    cent = abs(Top)/100
309    op = abs(Top)%100-1
310    M,T = SGData['SGOps'][op]
311    C = SGData['SGCen'][cent]
312    dx = .00001
313    deriv = np.zeros(6)
314    for i in [0,1,2]:
315        Oxyz[i] += dx
316        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
317        Oxyz[i] -= 2*dx
318        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
319        Oxyz[i] += dx
320        Txyz[i] += dx
321        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
322        Txyz[i] -= 2*dx
323        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
324        Txyz[i] += dx
325    return deriv
326   
327def getAngSig(VA,VB,Amat,SGData,covData={}):
328   
329    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
330        TxT = inv*(np.inner(M,Tx)+T)+C
331        TxT = G2spc.MoveToUnitCell(TxT)+U
332        return np.inner(Amat,(TxT-Ox))
333       
334    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
335        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
336        VecA /= np.sqrt(np.sum(VecA**2))
337        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
338        VecB /= np.sqrt(np.sum(VecB**2))
339        edge = VecB-VecA
340        edge = np.sum(edge**2)
341        angle = (2.-edge)/2.
342        angle = max(angle,-1.)
343        return acosd(angle)
344       
345    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
346    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
347    invA = invB = 1
348    invA = TopA/abs(TopA)
349    invB = TopB/abs(TopB)
350    centA = abs(TopA)/100
351    centB = abs(TopB)/100
352    opA = abs(TopA)%100-1
353    opB = abs(TopB)%100-1
354    MA,TA = SGData['SGOps'][opA]
355    MB,TB = SGData['SGOps'][opB]
356    CA = SGData['SGCen'][centA]
357    CB = SGData['SGCen'][centB]
358    if 'covMatrix' in covData:
359        covMatrix = covData['covMatrix']
360        varyList = covData['varyList']
361        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
362        dx = .00001
363        dadx = np.zeros(9)
364        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
365        for i in [0,1,2]:
366            OxA[i] += dx
367            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
368            OxA[i] -= 2*dx
369            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
370            OxA[i] += dx
371           
372            TxA[i] += dx
373            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
374            TxA[i] -= 2*dx
375            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
376            TxA[i] += dx
377           
378            TxB[i] += dx
379            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
380            TxB[i] -= 2*dx
381            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
382            TxB[i] += dx
383           
384        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
385        if sigAng < 0.01:
386            sigAng = 0.0
387        return Ang,sigAng
388    else:
389        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
390
391def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
392
393    def calcDist(Atoms,SyOps,Amat):
394        XYZ = []
395        for i,atom in enumerate(Atoms):
396            Inv,M,T,C,U = SyOps[i]
397            XYZ.append(np.array(atom[1:4]))
398            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
399            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
400        V1 = XYZ[1]-XYZ[0]
401        return np.sqrt(np.sum(V1**2))
402       
403    Inv = []
404    SyOps = []
405    names = []
406    for i,atom in enumerate(Oatoms):
407        names += atom[-1]
408        Op,unit = Atoms[i][-1]
409        inv = Op/abs(Op)
410        m,t = SGData['SGOps'][abs(Op)%100-1]
411        c = SGData['SGCen'][abs(Op)/100]
412        SyOps.append([inv,m,t,c,unit])
413    Dist = calcDist(Oatoms,SyOps,Amat)
414   
415    sig = -0.001
416    if 'covMatrix' in covData:
417        parmNames = []
418        dx = .00001
419        dadx = np.zeros(6)
420        for i in range(6):
421            ia = i/3
422            ix = i%3
423            Oatoms[ia][ix+1] += dx
424            a0 = calcDist(Oatoms,SyOps,Amat)
425            Oatoms[ia][ix+1] -= 2*dx
426            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
427        covMatrix = covData['covMatrix']
428        varyList = covData['varyList']
429        DistVcov = getVCov(names,varyList,covMatrix)
430        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
431        if sig < 0.001:
432            sig = -0.001
433   
434    return Dist,sig
435
436def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
437       
438    def calcAngle(Atoms,SyOps,Amat):
439        XYZ = []
440        for i,atom in enumerate(Atoms):
441            Inv,M,T,C,U = SyOps[i]
442            XYZ.append(np.array(atom[1:4]))
443            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
444            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
445        V1 = XYZ[1]-XYZ[0]
446        V1 /= np.sqrt(np.sum(V1**2))
447        V2 = XYZ[1]-XYZ[2]
448        V2 /= np.sqrt(np.sum(V2**2))
449        V3 = V2-V1
450        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
451        return acosd(cang)
452
453    Inv = []
454    SyOps = []
455    names = []
456    for i,atom in enumerate(Oatoms):
457        names += atom[-1]
458        Op,unit = Atoms[i][-1]
459        inv = Op/abs(Op)
460        m,t = SGData['SGOps'][abs(Op)%100-1]
461        c = SGData['SGCen'][abs(Op)/100]
462        SyOps.append([inv,m,t,c,unit])
463    Angle = calcAngle(Oatoms,SyOps,Amat)
464   
465    sig = -0.01
466    if 'covMatrix' in covData:
467        parmNames = []
468        dx = .00001
469        dadx = np.zeros(9)
470        for i in range(9):
471            ia = i/3
472            ix = i%3
473            Oatoms[ia][ix+1] += dx
474            a0 = calcAngle(Oatoms,SyOps,Amat)
475            Oatoms[ia][ix+1] -= 2*dx
476            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
477        covMatrix = covData['covMatrix']
478        varyList = covData['varyList']
479        AngVcov = getVCov(names,varyList,covMatrix)
480        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
481        if sig < 0.01:
482            sig = -0.01
483   
484    return Angle,sig
485
486def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
487   
488    def calcTorsion(Atoms,SyOps,Amat):
489       
490        XYZ = []
491        for i,atom in enumerate(Atoms):
492            Inv,M,T,C,U = SyOps[i]
493            XYZ.append(np.array(atom[1:4]))
494            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
495            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
496        V1 = XYZ[1]-XYZ[0]
497        V2 = XYZ[2]-XYZ[1]
498        V3 = XYZ[3]-XYZ[2]
499        V1 /= np.sqrt(np.sum(V1**2))
500        V2 /= np.sqrt(np.sum(V2**2))
501        V3 /= np.sqrt(np.sum(V3**2))
502        M = np.array([V1,V2,V3])
503        D = nl.det(M)
504        Ang = 1.0
505        P12 = np.dot(V1,V2)
506        P13 = np.dot(V1,V3)
507        P23 = np.dot(V2,V3)
508        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
509        return Tors
510           
511    Inv = []
512    SyOps = []
513    names = []
514    for i,atom in enumerate(Oatoms):
515        names += atom[-1]
516        Op,unit = Atoms[i][-1]
517        inv = Op/abs(Op)
518        m,t = SGData['SGOps'][abs(Op)%100-1]
519        c = SGData['SGCen'][abs(Op)/100]
520        SyOps.append([inv,m,t,c,unit])
521    Tors = calcTorsion(Oatoms,SyOps,Amat)
522   
523    sig = -0.01
524    if 'covMatrix' in covData:
525        parmNames = []
526        dx = .00001
527        dadx = np.zeros(12)
528        for i in range(12):
529            ia = i/3
530            ix = i%3
531            Oatoms[ia][ix+1] += dx
532            a0 = calcTorsion(Oatoms,SyOps,Amat)
533            Oatoms[ia][ix+1] -= 2*dx
534            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
535        covMatrix = covData['covMatrix']
536        varyList = covData['varyList']
537        TorVcov = getVCov(names,varyList,covMatrix)
538        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
539        if sig < 0.01:
540            sig = -0.01
541   
542    return Tors,sig
543       
544def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
545   
546    def calcDist(Atoms,SyOps,Amat):
547        XYZ = []
548        for i,atom in enumerate(Atoms):
549            Inv,M,T,C,U = SyOps[i]
550            XYZ.append(np.array(atom[1:4]))
551            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
552            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
553        V1 = XYZ[1]-XYZ[0]
554        return np.sqrt(np.sum(V1**2))
555       
556    def calcAngle(Atoms,SyOps,Amat):
557        XYZ = []
558        for i,atom in enumerate(Atoms):
559            Inv,M,T,C,U = SyOps[i]
560            XYZ.append(np.array(atom[1:4]))
561            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
562            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
563        V1 = XYZ[1]-XYZ[0]
564        V1 /= np.sqrt(np.sum(V1**2))
565        V2 = XYZ[1]-XYZ[2]
566        V2 /= np.sqrt(np.sum(V2**2))
567        V3 = V2-V1
568        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
569        return acosd(cang)
570
571    def calcTorsion(Atoms,SyOps,Amat):
572       
573        XYZ = []
574        for i,atom in enumerate(Atoms):
575            Inv,M,T,C,U = SyOps[i]
576            XYZ.append(np.array(atom[1:4]))
577            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
578            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
579        V1 = XYZ[1]-XYZ[0]
580        V2 = XYZ[2]-XYZ[1]
581        V3 = XYZ[3]-XYZ[2]
582        V1 /= np.sqrt(np.sum(V1**2))
583        V2 /= np.sqrt(np.sum(V2**2))
584        V3 /= np.sqrt(np.sum(V3**2))
585        M = np.array([V1,V2,V3])
586        D = nl.det(M)
587        Ang = 1.0
588        P12 = np.dot(V1,V2)
589        P13 = np.dot(V1,V3)
590        P23 = np.dot(V2,V3)
591        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
592        return Tors
593           
594    Inv = []
595    SyOps = []
596    names = []
597    for i,atom in enumerate(Oatoms):
598        names += atom[-1]
599        Op,unit = Atoms[i][-1]
600        inv = Op/abs(Op)
601        m,t = SGData['SGOps'][abs(Op)%100-1]
602        c = SGData['SGCen'][abs(Op)/100]
603        SyOps.append([inv,m,t,c,unit])
604    M = len(Oatoms)
605    if M == 2:
606        Val = calcDist(Oatoms,SyOps,Amat)
607    elif M == 3:
608        Val = calcAngle(Oatoms,SyOps,Amat)
609    else:
610        Val = calcTorsion(Oatoms,SyOps,Amat)
611   
612    sigVals = [-0.001,-0.01,-0.01]
613    sig = sigVals[M-3]
614    if 'covMatrix' in covData:
615        parmNames = []
616        dx = .00001
617        N = M*3
618        dadx = np.zeros(N)
619        for i in range(N):
620            ia = i/3
621            ix = i%3
622            Oatoms[ia][ix+1] += dx
623            if M == 2:
624                a0 = calcDist(Oatoms,SyOps,Amat)
625            elif M == 3:
626                a0 = calcAngle(Oatoms,SyOps,Amat)
627            else:
628                a0 = calcTorsion(Oatoms,SyOps,Amat)
629            Oatoms[ia][ix+1] -= 2*dx
630            if M == 2:
631                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
632            elif M == 3:
633                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
634            else:
635                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
636        covMatrix = covData['covMatrix']
637        varyList = covData['varyList']
638        Vcov = getVCov(names,varyList,covMatrix)
639        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
640        if sig < sigVals[M-3]:
641            sig = sigVals[M-3]
642   
643    return Val,sig
644       
645   
646def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
647    # returns value(esd) string; nTZ=True for no trailing zeros
648    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
649    #get the 2 significant digits in the esd
650    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
651    #get the number of digits to represent them
652    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
653   
654    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
655    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
656    if esd > 0:
657        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
658        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
659    elif esd < 0:
660         return str(round(value,mdec(esd)-1))
661    else:
662        text = str("%f"%(value))
663        if nTZ:
664            return text.rstrip('0')
665        else:
666            return text
667           
668def adjHKLmax(SGData,Hmax):
669    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
670        Hmax[0] = ((Hmax[0]+3)/6)*6
671        Hmax[1] = ((Hmax[1]+3)/6)*6
672        Hmax[2] = ((Hmax[2]+1)/4)*4
673    else:
674        Hmax[0] = ((Hmax[0]+2)/4)*4
675        Hmax[1] = ((Hmax[1]+2)/4)*4
676        Hmax[2] = ((Hmax[2]+1)/4)*4
677
678def FourierMap(data,reflData):
679   
680    generalData = data['General']
681    if not generalData['Map']['MapType']:
682        print '**** ERROR - Fourier map not defined'
683        return
684    mapData = generalData['Map']
685    dmin = mapData['Resolution']
686    SGData = generalData['SGData']
687    cell = generalData['Cell'][1:8]       
688    A = G2lat.cell2A(cell[:6])
689    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
690    adjHKLmax(SGData,Hmax)
691    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
692#    Fhkl[0,0,0] = generalData['F000X']
693    time0 = time.time()
694    for ref in reflData:
695        if ref[4] >= dmin:
696            Fosq,Fcsq,ph = ref[8:11]
697            for i,hkl in enumerate(ref[11]):
698                hkl = np.asarray(hkl,dtype='i')
699                dp = 360.*ref[12][i]
700                a = cosd(ph+dp)
701                b = sind(ph+dp)
702                phasep = complex(a,b)
703                phasem = complex(a,-b)
704                if 'Fobs' in mapData['MapType']:
705                    F = np.sqrt(Fosq)
706                    h,k,l = hkl+Hmax
707                    Fhkl[h,k,l] = F*phasep
708                    h,k,l = -hkl+Hmax
709                    Fhkl[h,k,l] = F*phasem
710                elif 'Fcalc' in mapData['MapType']:
711                    F = np.sqrt(Fcsq)
712                    h,k,l = hkl+Hmax
713                    Fhkl[h,k,l] = F*phasep
714                    h,k,l = -hkl+Hmax
715                    Fhkl[h,k,l] = F*phasem
716                elif 'delt-F' in mapData['MapType']:
717                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
718                    h,k,l = hkl+Hmax
719                    Fhkl[h,k,l] = dF*phasep
720                    h,k,l = -hkl+Hmax
721                    Fhkl[h,k,l] = dF*phasem
722                elif '2*Fo-Fc' in mapData['MapType']:
723                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
724                    h,k,l = hkl+Hmax
725                    Fhkl[h,k,l] = F*phasep
726                    h,k,l = -hkl+Hmax
727                    Fhkl[h,k,l] = F*phasem
728                elif 'Patterson' in mapData['MapType']:
729                    h,k,l = hkl+Hmax
730                    Fhkl[h,k,l] = complex(Fosq,0.)
731                    h,k,l = -hkl+Hmax
732                    Fhkl[h,k,l] = complex(Fosq,0.)
733    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
734    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
735    mapData['rho'] = np.real(rho)
736    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
737    return mapData
738   
739# map printing for testing purposes
740def printRho(SGLaue,rho,rhoMax):                         
741    dim = len(rho.shape)
742    if dim == 2:
743        ix,jy = rho.shape
744        for j in range(jy):
745            line = ''
746            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
747                line += (jy-j)*'  '
748            for i in range(ix):
749                r = int(100*rho[i,j]/rhoMax)
750                line += '%4d'%(r)
751            print line+'\n'
752    else:
753        ix,jy,kz = rho.shape
754        for k in range(kz):
755            print 'k = ',k
756            for j in range(jy):
757                line = ''
758                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
759                    line += (jy-j)*'  '
760                for i in range(ix):
761                    r = int(100*rho[i,j,k]/rhoMax)
762                    line += '%4d'%(r)
763                print line+'\n'
764## keep this
765               
766def findOffset(SGData,A,Fhkl):   
767    if SGData['SpGrp'] == 'P 1':
768        return [0,0,0]   
769    hklShape = Fhkl.shape
770    hklHalf = np.array(hklShape)/2
771    sortHKL = np.argsort(Fhkl.flatten())
772    Fdict = {}
773    for hkl in sortHKL:
774        HKL = np.unravel_index(hkl,hklShape)
775        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
776        if F == 0.:
777            break
778        Fdict['%.6f'%(np.absolute(F))] = hkl
779    Flist = np.flipud(np.sort(Fdict.keys()))
780    F = str(1.e6)
781    i = 0
782    DH = []
783    Dphi = []
784    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
785    for F in Flist:
786        hkl = np.unravel_index(Fdict[F],hklShape)
787        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
788        Uniq = np.array(Uniq,dtype='i')
789        Phi = np.array(Phi)
790        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
791        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
792        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
793        ang0 = np.angle(Fh0,deg=True)/360.
794        for H,phi in zip(Uniq,Phi)[1:]:
795            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
796            dH = H-hkl
797            dang = ang-ang0
798            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
799                continue
800            DH.append(dH)
801            Dphi.append((dang+.5) % 1.0)
802        if i > 20 or len(DH) > 30:
803            break
804        i += 1
805    DH = np.array(DH)
806    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
807    Dphi = np.array(Dphi)
808    steps = np.array(hklShape)
809    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
810    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
811    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
812    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
813    hist,bins = np.histogram(Mmap,bins=1000)
814#    for i,item in enumerate(hist[:10]):
815#        print item,bins[i]
816    chisq = np.min(Mmap)
817    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
818    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
819#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
820    return DX
821   
822def ChargeFlip(data,reflData,pgbar):
823    generalData = data['General']
824    mapData = generalData['Map']
825    flipData = generalData['Flip']
826    FFtable = {}
827    if 'None' not in flipData['Norm element']:
828        normElem = flipData['Norm element'].upper()
829        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
830        for ff in FFs:
831            if ff['Symbol'] == normElem:
832                FFtable.update(ff)
833    dmin = flipData['Resolution']
834    SGData = generalData['SGData']
835    cell = generalData['Cell'][1:8]       
836    A = G2lat.cell2A(cell[:6])
837    Vol = cell[6]
838    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
839    adjHKLmax(SGData,Hmax)
840    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
841    time0 = time.time()
842    for ref in reflData:
843        dsp = ref[4]
844        if dsp >= dmin:
845            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
846            if FFtable:
847                SQ = 0.25/dsp**2
848                ff *= G2el.ScatFac(FFtable,SQ)[0]
849            if ref[8] > 0.:         #use only +ve Fobs**2
850                E = np.sqrt(ref[8])/ff
851            else:
852                E = 0.
853            ph = ref[10]
854            ph = rn.uniform(0.,360.)
855            for i,hkl in enumerate(ref[11]):
856                hkl = np.asarray(hkl,dtype='i')
857                dp = 360.*ref[12][i]
858                a = cosd(ph+dp)
859                b = sind(ph+dp)
860                phasep = complex(a,b)
861                phasem = complex(a,-b)
862                h,k,l = hkl+Hmax
863                Ehkl[h,k,l] = E*phasep
864                h,k,l = -hkl+Hmax       #Friedel pair refl.
865                Ehkl[h,k,l] = E*phasem
866#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
867    CEhkl = copy.copy(Ehkl)
868    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
869    Emask = ma.getmask(MEhkl)
870    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
871    Ncyc = 0
872    old = np.seterr(all='raise')
873    while True:       
874        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
875        CEsig = np.std(CErho)
876        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
877        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
878        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
879        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
880        phase = CFhkl/np.absolute(CFhkl)
881        CEhkl = np.absolute(Ehkl)*phase
882        Ncyc += 1
883        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
884        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
885        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
886        if Rcf < 5.:
887            break
888        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
889        if not GoOn or Ncyc > 10000:
890            break
891    np.seterr(**old)
892    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
893    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
894    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
895    roll = findOffset(SGData,A,CEhkl)
896       
897    mapData['Rcf'] = Rcf
898    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
899    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
900    mapData['rollMap'] = [0,0,0]
901    return mapData
902   
903def SearchMap(data):
904    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
905   
906    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
907   
908    def noDuplicate(xyz,peaks,Amat):
909        XYZ = np.inner(Amat,xyz)
910        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
911            print ' Peak',xyz,' <0.5A from another peak'
912            return False
913        return True
914                           
915    def fixSpecialPos(xyz,SGData,Amat):
916        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
917        X = []
918        xyzs = [equiv[0] for equiv in equivs]
919        for x in xyzs:
920            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
921                X.append(x)
922        if len(X) > 1:
923            return np.average(X,axis=0)
924        else:
925            return xyz
926       
927    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
928        Mag,x0,y0,z0,sig = parms
929        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
930#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
931        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
932       
933    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
934        Mag,x0,y0,z0,sig = parms
935        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
936        return M
937       
938    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
939        Mag,x0,y0,z0,sig = parms
940        dMdv = np.zeros(([5,]+list(rX.shape)))
941        delt = .01
942        for i in range(5):
943            parms[i] -= delt
944            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
945            parms[i] += 2.*delt
946            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
947            parms[i] -= delt
948            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
949        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
950        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
951        dMdv = np.reshape(dMdv,(5,rX.size))
952        Hess = np.inner(dMdv,dMdv)
953       
954        return Vec,Hess
955       
956    generalData = data['General']
957    phaseName = generalData['Name']
958    SGData = generalData['SGData']
959    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
960    drawingData = data['Drawing']
961    peaks = []
962    mags = []
963    dzeros = []
964    try:
965        mapData = generalData['Map']
966        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
967        rho = copy.copy(mapData['rho'])     #don't mess up original
968        mapHalf = np.array(rho.shape)/2
969        res = mapData['Resolution']
970        incre = np.array(rho.shape,dtype=np.float)
971        step = max(1.0,1./res)+1
972        steps = np.array(3*[step,])
973    except KeyError:
974        print '**** ERROR - Fourier map not defined'
975        return peaks,mags
976    rhoMask = ma.array(rho,mask=(rho<contLevel))
977    indices = (-1,0,1)
978    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
979    for roll in rolls:
980        if np.any(roll):
981            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
982    indx = np.transpose(rhoMask.nonzero())
983    peaks = indx/incre
984    mags = rhoMask[rhoMask.nonzero()]
985    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
986        rho = rollMap(rho,ind)
987        rMM = mapHalf-steps
988        rMP = mapHalf+steps+1
989        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
990        peakInt = np.sum(rhoPeak)*res**3
991        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
992        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
993        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
994            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
995        x1 = result[0]
996        if not np.any(x1 < 0):
997            mag = x1[0]
998            peak = (np.array(x1[1:4])-ind)/incre
999        peak = fixSpecialPos(peak,SGData,Amat)
1000        rho = rollMap(rho,-ind)       
1001    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1002    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1003   
1004def sortArray(data,pos,reverse=False):
1005    #data is a list of items
1006    #sort by pos in list; reverse if True
1007    T = []
1008    for i,M in enumerate(data):
1009        T.append((M[pos],i))
1010    D = dict(zip(T,data))
1011    T.sort()
1012    if reverse:
1013        T.reverse()
1014    X = []
1015    for key in T:
1016        X.append(D[key])
1017    return X
1018
1019def PeaksEquiv(data,Ind):
1020
1021    def Duplicate(xyz,peaks,Amat):
1022        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1023            return True
1024        return False
1025                           
1026    generalData = data['General']
1027    cell = generalData['Cell'][1:7]
1028    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1029    A = G2lat.cell2A(cell)
1030    SGData = generalData['SGData']
1031    mapPeaks = data['Map Peaks']
1032    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1033    Indx = {}
1034    for ind in Ind:
1035        xyz = np.array(mapPeaks[ind][1:4])
1036        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1037#        for x in xyzs: print x
1038        for jnd,xyz in enumerate(XYZ):       
1039            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1040    Ind = []
1041    for ind in Indx:
1042        if Indx[ind]:
1043            Ind.append(ind)
1044    return Ind
1045               
1046def PeaksUnique(data,Ind):
1047#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1048
1049    def noDuplicate(xyz,peaks,Amat):
1050        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1051            return False
1052        return True
1053                           
1054    generalData = data['General']
1055    cell = generalData['Cell'][1:7]
1056    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1057    A = G2lat.cell2A(cell)
1058    SGData = generalData['SGData']
1059    mapPeaks = data['Map Peaks']
1060    Indx = {}
1061    XYZ = {}
1062    for ind in Ind:
1063        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1064        Indx[ind] = True
1065    for ind in Ind:
1066        if Indx[ind]:
1067            xyz = XYZ[ind]
1068            for jnd in Ind:
1069                if ind != jnd and Indx[jnd]:                       
1070                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1071                    xyzs = np.array([equiv[0] for equiv in Equiv])
1072                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1073    Ind = []
1074    for ind in Indx:
1075        if Indx[ind]:
1076            Ind.append(ind)
1077    return Ind
1078   
1079def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1080    ind = 0
1081    if useFit:
1082        ind = 1
1083    ins = {}
1084    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1085        for x in ['U','V','W','X','Y']:
1086            ins[x] = Parms[x][ind]
1087        if ifQ:                              #qplot - convert back to 2-theta
1088            pos = 2.0*asind(pos*wave/(4*math.pi))
1089        sig = ins['U']*tand(pos/2.0)**2+ins['V']*tand(pos/2.0)+ins['W']
1090        gam = ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)           
1091        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1092    else:
1093        if ifQ:
1094            dsp = 2.*np.pi/pos
1095            pos = Parms['difC']*dsp
1096        else:
1097            dsp = pos/Parms['difC'][1]
1098        if 'Pdabc' in Parms2:
1099            for x in ['sig-0','sig-1','X','Y']:
1100                ins[x] = Parms[x][ind]
1101            Pdabc = Parms2['Pdabc'].T
1102            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1103            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1104        else:
1105            for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
1106                ins[x] = Parms[x][ind]
1107            alp = ins['alpha']/dsp
1108            bet = ins['beta-0']+ins['beta-1']/dsp**4
1109        sig = ins['sig-0']+ins['sig-1']*dsp**2
1110        gam = ins['X']*dsp+ins['Y']*dsp**2
1111        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
1112    return XY
1113       
1114def getWave(Parms):
1115    try:
1116        return Parms['Lam'][1]
1117    except KeyError:
1118        return Parms['Lam1'][1]
1119   
1120def prodQQ(QA,QB):
1121    ''' Grassman quaternion product
1122        QA,QB quaternions; q=r+ai+bj+ck
1123    '''
1124    D = np.zeros(4)
1125    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
1126    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
1127    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
1128    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
1129    return D
1130   
1131def normQ(QA):
1132    ''' get length of quaternion & normalize it
1133        q=r+ai+bj+ck
1134    '''
1135    n = np.sqrt(np.sum(np.array(QA)**2))
1136    return QA/n
1137   
1138def invQ(Q):
1139    '''
1140        get inverse of quaternion
1141        q=r+ai+bj+ck; q* = r-ai-bj-ck
1142    '''
1143    return Q*np.array([1,-1,-1,-1])
1144   
1145def prodQVQ(Q,V):
1146    ''' compute the quaternion vector rotation qvq-1 = v'
1147        q=r+ai+bj+ck
1148    '''
1149    VP = np.zeros(3)
1150    T2 = Q[0]*Q[1]
1151    T3 = Q[0]*Q[2]
1152    T4 = Q[0]*Q[3]
1153    T5 = -Q[1]*Q[1]
1154    T6 = Q[1]*Q[2]
1155    T7 = Q[1]*Q[3]
1156    T8 = -Q[2]*Q[2]
1157    T9 = Q[2]*Q[3]
1158    T10 = -Q[3]*Q[3]
1159    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
1160    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
1161    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
1162    return VP   
1163   
1164def Q2Mat(Q):
1165    ''' make rotation matrix from quaternion
1166        q=r+ai+bj+ck
1167    '''
1168    aa = Q[0]**2
1169    ab = Q[0]*Q[1]
1170    ac = Q[0]*Q[2]
1171    ad = Q[0]*Q[3]
1172    bb = Q[1]**2
1173    bc = Q[1]*Q[2]
1174    bd = Q[1]*Q[3]
1175    cc = Q[2]**2
1176    cd = Q[2]*Q[3]
1177    dd = Q[3]**2
1178    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
1179        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
1180        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
1181    return np.array(M)
1182   
1183def AV2Q(A,V):
1184    ''' convert angle (radians -pi to pi) & vector to quaternion
1185        q=r+ai+bj+ck
1186    '''
1187    Q = np.zeros(4)
1188    d = np.sqrt(np.sum(np.array(V)**2))
1189    if d:
1190        V /= d
1191    else:
1192        return [1.,0.,0.,0.]    #identity
1193    p = A/2.
1194    Q[0] = np.cos(p)
1195    s = np.sin(p)
1196    Q[1:4] = V*s
1197    return Q
1198   
1199def AVdeg2Q(A,V):
1200    ''' convert angle (degrees -180 to 180) & vector to quaternion
1201        q=r+ai+bj+ck
1202    '''
1203    Q = np.zeros(4)
1204    d = np.sqrt(np.sum(np.array(V)**2))
1205    if d:
1206        V /= d
1207    else:
1208        return [1.,0.,0.,0.]    #identity
1209    p = A/2.
1210    Q[0] = cosd(p)
1211    S = sind(p)
1212    Q[1:4] = V*S
1213    return Q
1214   
Note: See TracBrowser for help on using the repository browser.