source: trunk/GSASIImath.py @ 2321

Last change on this file since 2321 was 2321, checked in by vondreele, 7 years ago

start on CalcDistSig? & AddAngleDialog? from tablet
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 166.0 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2016-06-15 12:10:38 +0000 (Wed, 15 Jun 2016) $
5# $Author: vondreele $
6# $Revision: 2321 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 2321 2016-06-15 12:10:38Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 2321 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIIpwd as G2pwd
34import numpy.fft as fft
35import scipy.optimize as so
36import pypowder as pwd
37
38sind = lambda x: np.sin(x*np.pi/180.)
39cosd = lambda x: np.cos(x*np.pi/180.)
40tand = lambda x: np.tan(x*np.pi/180.)
41asind = lambda x: 180.*np.arcsin(x)/np.pi
42acosd = lambda x: 180.*np.arccos(x)/np.pi
43atand = lambda x: 180.*np.arctan(x)/np.pi
44atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
45twopi = 2.0*np.pi
46twopisq = 2.0*np.pi**2
47nxs = np.newaxis
48   
49################################################################################
50##### Hessian least-squares Levenberg-Marquardt routine
51################################################################################
52
53def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
54   
55    """
56    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
57    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
58   
59    ::
60   
61                    Nobs
62        x = arg min(sum(func(y)**2,axis=0))
63                    y=0
64
65    :param function func: callable method or function
66        should take at least one (possibly length N vector) argument and
67        returns M floating point numbers.
68    :param np.ndarray x0: The starting estimate for the minimization of length N
69    :param function Hess: callable method or function
70        A required function or method to compute the weighted vector and Hessian for func.
71        It must be a symmetric NxN array
72    :param tuple args: Any extra arguments to func are placed in this tuple.
73    :param float ftol: Relative error desired in the sum of squares.
74    :param float xtol: Relative error desired in the approximate solution.
75    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
76        until other limits are met (ftol, xtol)
77    :param bool Print: True for printing results (residuals & times) by cycle
78
79    :returns: (x,cov_x,infodict) where
80
81      * x : ndarray
82        The solution (or the result of the last iteration for an unsuccessful
83        call).
84      * cov_x : ndarray
85        Uses the fjac and ipvt optional outputs to construct an
86        estimate of the jacobian around the solution.  ``None`` if a
87        singular matrix encountered (indicates very flat curvature in
88        some direction).  This matrix must be multiplied by the
89        residual standard deviation to get the covariance of the
90        parameter estimates -- see curve_fit.
91      * infodict : dict
92        a dictionary of optional outputs with the keys:
93
94         * 'fvec' : the function evaluated at the output
95         * 'num cyc':
96         * 'nfev':
97         * 'lamMax':
98         * 'psing':
99           
100    """
101               
102    ifConverged = False
103    deltaChi2 = -10.
104    x0 = np.array(x0, ndmin=1)      #might be redundant?
105    n = len(x0)
106    if type(args) != type(()):
107        args = (args,)
108       
109    icycle = 0
110    One = np.ones((n,n))
111    lam = 0.001
112    lamMax = lam
113    nfev = 0
114    if Print:
115        print ' Hessian refinement on %d variables:'%(n)
116    Lam = np.zeros((n,n))
117    while icycle < maxcyc:
118        time0 = time.time()
119        M = func(x0,*args)
120        nfev += 1
121        chisq0 = np.sum(M**2)
122        Yvec,Amat = Hess(x0,*args)
123        Adiag = np.sqrt(np.diag(Amat))
124        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
125        if np.any(psing):                #hard singularity in matrix
126            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
127        Anorm = np.outer(Adiag,Adiag)
128        Yvec /= Adiag
129        Amat /= Anorm
130        while True:
131            Lam = np.eye(Amat.shape[0])*lam
132            Amatlam = Amat*(One+Lam)
133            try:
134                Xvec = nl.solve(Amatlam,Yvec)
135            except nl.LinAlgError:
136                print 'ouch #1'
137                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
138                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
139            Xvec /= Adiag
140            M2 = func(x0+Xvec,*args)
141            nfev += 1
142            chisq1 = np.sum(M2**2)
143            if chisq1 > chisq0*(1.+ftol):
144                lam *= 10.
145                if Print:
146                    print 'matrix modification needed; lambda now %.1e'%(lam)
147            else:
148                x0 += Xvec
149                lam /= 10.
150                break
151            if lam > 10.e3:
152                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
153                break
154        lamMax = max(lamMax,lam)
155        deltaChi2 = (chisq0-chisq1)/chisq0
156        if Print:
157            print ' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Lambda: %.3g,  Delta: %.3g'%(
158                icycle,time.time()-time0,chisq1,lamMax,deltaChi2)
159        if deltaChi2 < ftol:
160            ifConverged = True
161            if Print: print "converged"
162            break
163        icycle += 1
164    M = func(x0,*args)
165    nfev += 1
166    Yvec,Amat = Hess(x0,*args)
167    Adiag = np.sqrt(np.diag(Amat))
168    Anorm = np.outer(Adiag,Adiag)
169    Lam = np.eye(Amat.shape[0])*lam
170    Amatlam = Amat/Anorm  #*(One+Lam)              #don't scale Amat to Marquardt array       
171    try:
172        Bmat = nl.inv(Amatlam)/Anorm  #*(One+Lam)      #don't rescale Bmat to Marquardt array
173        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
174    except nl.LinAlgError:
175        print 'ouch #2 linear algebra error in making v-cov matrix'
176        psing = []
177        if maxcyc:
178            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
179        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
180   
181def getVCov(varyNames,varyList,covMatrix):
182    '''obtain variance-covariance terms for a set of variables. NB: the varyList
183    and covMatrix were saved by the last least squares refinement so they must match.
184   
185    :param list varyNames: variable names to find v-cov matric for
186    :param list varyList: full list of all variables in v-cov matrix
187    :param nparray covMatrix: full variance-covariance matrix from the last
188     least squares refinement
189   
190    :returns: nparray vcov: variance-covariance matrix for the variables given
191     in varyNames
192   
193    '''
194    vcov = np.zeros((len(varyNames),len(varyNames)))
195    for i1,name1 in enumerate(varyNames):
196        for i2,name2 in enumerate(varyNames):
197            try:
198                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
199            except ValueError:
200                vcov[i1][i2] = 0.0
201#                if i1 == i2:
202#                    vcov[i1][i2] = 1e-20
203#                else:
204#                    vcov[i1][i2] = 0.0
205    return vcov
206   
207################################################################################
208##### Atom manipulations
209################################################################################
210
211def FindMolecule(ind,generalData,atomData):                    #uses numpy & masks - very fast even for proteins!
212
213    def getNeighbors(atom,radius):
214        neighList = [] 
215        Dx = UAtoms-np.array(atom[cx:cx+3])
216        dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
217        sumR = Radii+radius
218        return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0])                #get indices of bonded atoms
219
220    import numpy.ma as ma
221    indices = (-1,0,1)
222    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
223    cx,ct,cs,ci = generalData['AtomPtrs']
224    DisAglCtls = generalData['DisAglCtls']
225    SGData = generalData['SGData']
226    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
227    radii = DisAglCtls['BondRadii']
228    atomTypes = DisAglCtls['AtomTypes']
229    factor = DisAglCtls['Factors'][0]
230    unit = np.zeros(3)
231    try:
232        indH = atomTypes.index('H')
233        radii[indH] = 0.5
234    except:
235        pass           
236    nAtom = len(atomData)
237    Indx = range(nAtom)
238    UAtoms = []
239    Radii = []
240    for atom in atomData:
241        UAtoms.append(np.array(atom[cx:cx+3]))
242        Radii.append(radii[atomTypes.index(atom[ct])])
243    UAtoms = np.array(UAtoms)
244    Radii = np.array(Radii)
245    for nOp,Op in enumerate(SGData['SGOps'][1:]):
246        UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
247        Radii = np.concatenate((Radii,Radii[:nAtom]))
248        Indx += Indx[:nAtom]
249    for icen,cen in enumerate(SGData['SGCen'][1:]):
250        UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
251        Radii = np.concatenate((Radii,Radii))
252        Indx += Indx[:nAtom]
253    if SGData['SGInv']:
254        UAtoms = np.concatenate((UAtoms,-UAtoms))
255        Radii = np.concatenate((Radii,Radii))
256        Indx += Indx
257    UAtoms %= 1.
258    mAtoms = len(UAtoms)
259    for unit in Units:
260        if np.any(unit):    #skip origin cell
261            UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
262            Radii = np.concatenate((Radii,Radii[:mAtoms]))
263            Indx += Indx[:mAtoms]
264    UAtoms = np.array(UAtoms)
265    Radii = np.array(Radii)
266    newAtoms = [atomData[ind],]
267    atomData[ind] = None
268    radius = Radii[ind]
269    IndB = getNeighbors(newAtoms[-1],radius)
270    while True:
271        if not len(IndB):
272            break
273        indb = IndB.pop()
274        if atomData[Indx[indb]] == None:
275            continue
276        while True:
277            try:
278                jndb = Indb.index(indb)
279                Indb.remove(jndb)
280            except:
281                break
282        newAtom = copy.copy(atomData[Indx[indb]])
283        newAtom[cx:cx+3] = UAtoms[indb]     #NB: thermal Uij, etc. not transformed!
284        newAtoms.append(newAtom)
285        atomData[Indx[indb]] = None
286        IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
287        if len(IndB) > nAtom:
288            return 'Assemble molecule cannot be used on extended structures'
289    for atom in atomData:
290        if atom != None:
291            newAtoms.append(atom)
292    return newAtoms
293       
294def FindAtomIndexByIDs(atomData,loc,IDs,Draw=True):
295    '''finds the set of atom array indices for a list of atom IDs. Will search
296    either the Atom table or the drawAtom table.
297   
298    :param list atomData: Atom or drawAtom table containting coordinates, etc.
299    :param int loc: location of atom id in atomData record
300    :param list IDs: atom IDs to be found
301    :param bool Draw: True if drawAtom table to be searched; False if Atom table
302      is searched
303   
304    :returns: list indx: atom (or drawAtom) indices
305   
306    '''
307    indx = []
308    for i,atom in enumerate(atomData):
309        if Draw and atom[loc] in IDs:
310            indx.append(i)
311        elif atom[loc] in IDs:
312            indx.append(i)
313    return indx
314
315def FillAtomLookUp(atomData,indx):
316    '''create a dictionary of atom indexes with atom IDs as keys
317   
318    :param list atomData: Atom table to be used
319   
320    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
321   
322    '''
323    atomLookUp = {}
324    for iatm,atom in enumerate(atomData):
325        atomLookUp[atom[indx]] = iatm
326    return atomLookUp
327
328def GetAtomsById(atomData,atomLookUp,IdList):
329    '''gets a list of atoms from Atom table that match a set of atom IDs
330   
331    :param list atomData: Atom table to be used
332    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
333    :param list IdList: atom IDs to be found
334   
335    :returns: list atoms: list of atoms found
336   
337    '''
338    atoms = []
339    for id in IdList:
340        atoms.append(atomData[atomLookUp[id]])
341    return atoms
342   
343def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
344    '''gets atom parameters for atoms using atom IDs
345   
346    :param list atomData: Atom table to be used
347    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
348    :param list IdList: atom IDs to be found
349    :param int itemLoc: pointer to desired 1st item in an atom table entry
350    :param int numItems: number of items to be retrieved
351   
352    :returns: type name: description
353   
354    '''
355    Items = []
356    if not isinstance(IdList,list):
357        IdList = [IdList,]
358    for id in IdList:
359        if numItems == 1:
360            Items.append(atomData[atomLookUp[id]][itemLoc])
361        else:
362            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
363    return Items
364   
365def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
366    '''default doc string
367   
368    :param type name: description
369   
370    :returns: type name: description
371   
372    '''
373    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
374    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
375    XYZ = []
376    for ind in indx:
377        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
378        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
379        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
380    return XYZ
381   
382def TransformAtoms(Atoms,cx,cia,Trans,Vec):
383    for Atom in Atoms:
384        XYZ = Atom[cx:cx+3]
385        if 'A' in Atom[cia]:
386            U6 = Atom[cia+2:cia+8]
387   
388
389def FindNeighbors(phase,FrstName,AtNames,notName=''):
390    General = phase['General']
391    cx,ct,cs,cia = General['AtomPtrs']
392    Atoms = phase['Atoms']
393    atNames = [atom[ct-1] for atom in Atoms]
394    Cell = General['Cell'][1:7]
395    Amat,Bmat = G2lat.cell2AB(Cell)
396    atTypes = General['AtomTypes']
397    Radii = np.array(General['BondRadii'])
398    DisAglCtls = General['DisAglCtls']   
399    radiusFactor = DisAglCtls['Factors'][0]
400    AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
401    Orig = atNames.index(FrstName)
402    OId = Atoms[Orig][cia+8]
403    OType = Atoms[Orig][ct]
404    XYZ = getAtomXYZ(Atoms,cx)       
405    Neigh = []
406    Ids = []
407    Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
408    dist = np.sqrt(np.sum(Dx**2,axis=1))
409    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
410    IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))
411    for j in IndB[0]:
412        if j != Orig:
413            if AtNames[j] != notName:
414                Neigh.append([AtNames[j],dist[j],True])
415                Ids.append(Atoms[j][cia+8])
416    return Neigh,[OId,Ids]
417   
418def FindAllNeighbors(phase,FrstName,AtNames,notName=''):
419    General = phase['General']
420    cx,ct,cs,cia = General['AtomPtrs']
421    Atoms = phase['Atoms']
422    atNames = [atom[ct-1] for atom in Atoms]
423    Cell = General['Cell'][1:7]
424    Amat,Bmat = G2lat.cell2AB(Cell)
425    SGData = General['SGData']
426    indices = (-1,0,1)
427    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
428    atTypes = General['AtomTypes']
429    Radii = np.array(General['BondRadii'])
430    DisAglCtls = General['DisAglCtls']   
431    radiusFactor = DisAglCtls['Factors'][0]
432    AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
433    Orig = atNames.index(FrstName)
434    OId = Atoms[Orig][cia+8]
435    OType = Atoms[Orig][ct]
436    XYZ = getAtomXYZ(Atoms,cx)       
437    Oxyz = XYZ[Orig]
438    Neigh = []
439    Ids = []
440    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
441    sumR = np.reshape(np.tile(sumR,27),(27,-1))
442    results = []
443    for xyz in XYZ:
444        results.append(G2spc.GenAtom(xyz,SGData,False,Move=False))
445    for iA,result in enumerate(results):
446        if iA != Orig:               
447            for [Txyz,Top,Tunit] in result:
448                Dx = np.array([Txyz-Oxyz+unit for unit in Units])
449                dx = np.inner(Dx,Amat)
450                dist = np.sqrt(np.sum(dx**2,axis=1))
451                IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR[:,iA],0.))
452        #        GSASIIpath.IPyBreak()
453                for iU in IndB[0]:
454                    if AtNames[iA] != notName:
455                        unit = Units[iU]
456                        if np.any(unit):
457                            Topstr = ' +(%4d)[%2d,%2d,%2d]'%(Top,unit[0],unit[1],unit[2])
458                        else:
459                            Topstr = ' +(%4d)'%(Top)
460                        Neigh.append([AtNames[iA]+Topstr,dist[iU]])
461                        Ids.append(Atoms[iA][cia+8])
462    return Neigh,[OId,Ids]
463   
464def calcBond(A,Ax,Bx,MTCU):
465    cell = G2lat.A2cell(A)
466    Amat,Bmat = G2lat.cell2AB(cell)
467    M,T,C,U = MTCU
468    Btx = np.inner(M,Bx)+T+C+U
469    Dx = Btx-Ax
470    dist = np.sqrt(np.inner(Amat,Dx))
471    return dist
472   
473def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
474   
475    def getTransMat(RXYZ,OXYZ,TXYZ,Amat):
476        Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
477        Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
478        Mat2 = np.cross(Vec[0],Vec[1])      #UxV
479        Mat2 /= np.sqrt(np.sum(Mat2**2))
480        Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
481        return nl.inv(np.array([Vec[0],Mat2,Mat3]))       
482   
483    cx,ct,cs,cia = General['AtomPtrs']
484    Cell = General['Cell'][1:7]
485    Amat,Bmat = G2lat.cell2AB(Cell)
486    nBonds = AddHydId[-1]+len(AddHydId[1])
487    Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
488    OXYZ = np.array(Oatom[cx:cx+3])
489    if 'I' in Oatom[cia]:
490        Uiso = Oatom[cia+1]
491    else:
492        Uiso = (Oatom[cia+2]+Oatom[cia+3]+Oatom[cia+4])/3.0       #simple average
493    Uiso = max(Uiso,0.005)                      #set floor!
494    Tatoms = GetAtomsById(Atoms,AtLookUp,AddHydId[1])
495    TXYZ = np.array([tatom[cx:cx+3] for tatom in Tatoms]) #3 x xyz
496    DX = np.inner(Amat,TXYZ-OXYZ).T
497    if nBonds == 4:
498        if AddHydId[-1] == 1:
499            Vec = TXYZ-OXYZ
500            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
501            Vec = np.sum(Vec/Len,axis=0)
502            Len = np.sqrt(np.sum(Vec**2))
503            Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
504            HU = 1.1*Uiso
505            return [Hpos,],[HU,]
506        elif AddHydId[-1] == 2:
507            Vec = np.inner(Amat,TXYZ-OXYZ).T
508            Vec[0] += Vec[1]            #U - along bisector
509            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
510            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
511            Mat2 /= np.sqrt(np.sum(Mat2**2))
512            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
513            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
514            Hpos = np.array([[-0.97*cosd(54.75),0.97*sind(54.75),0.],
515                [-0.97*cosd(54.75),-0.97*sind(54.75),0.]])
516            HU = 1.2*Uiso*np.ones(2)
517            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
518            return Hpos,HU
519        else:
520            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
521            RXYZ = np.array(Ratom[cx:cx+3])
522            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
523            a = 0.96*cosd(70.5)
524            b = 0.96*sind(70.5)
525            Hpos = np.array([[a,0.,-b],[a,-b*cosd(30.),0.5*b],[a,b*cosd(30.),0.5*b]])
526            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
527            HU = 1.5*Uiso*np.ones(3)
528            return Hpos,HU         
529    elif nBonds == 3:
530        if AddHydId[-1] == 1:
531            Vec = np.sum(TXYZ-OXYZ,axis=0)               
532            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
533            Vec = -0.93*Vec/Len
534            Hpos = OXYZ+Vec
535            HU = 1.1*Uiso
536            return [Hpos,],[HU,]
537        elif AddHydId[-1] == 2:
538            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
539            RXYZ = np.array(Ratom[cx:cx+3])
540            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
541            a = 0.93*cosd(60.)
542            b = 0.93*sind(60.)
543            Hpos = [[a,b,0],[a,-b,0]]
544            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
545            HU = 1.2*Uiso*np.ones(2)
546            return Hpos,HU
547    else:   #2 bonds
548        if 'C' in Oatom[ct]:
549            Vec = TXYZ[0]-OXYZ
550            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
551            Vec = -0.93*Vec/Len
552            Hpos = OXYZ+Vec
553            HU = 1.1*Uiso
554            return [Hpos,],[HU,]
555        elif 'O' in Oatom[ct]:
556            mapData = General['Map']
557            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
558            RXYZ = np.array(Ratom[cx:cx+3])
559            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
560            a = 0.82*cosd(70.5)
561            b = 0.82*sind(70.5)
562            azm = np.arange(0.,360.,5.)
563            Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm])
564            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
565            Rhos = np.array([getRho(pos,mapData) for pos in Hpos])
566            imax = np.argmax(Rhos)
567            HU = 1.5*Uiso
568            return [Hpos[imax],],[HU,]
569    return [],[]
570       
571def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
572    '''default doc string
573   
574    :param type name: description
575   
576    :returns: type name: description
577   
578    '''
579    for atom in atomData:
580        XYZ = np.inner(Amat,atom[cx:cx+3])
581        if atom[cia] == 'A':
582            UIJ = atom[cia+2:cia+8]
583               
584def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
585    '''default doc string
586   
587    :param type name: description
588   
589    :returns: type name: description
590   
591    '''
592    TLStype,TLS = rbObj['ThermalMotion'][:2]
593    Tmat = np.zeros((3,3))
594    Lmat = np.zeros((3,3))
595    Smat = np.zeros((3,3))
596    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
597        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
598    if 'T' in TLStype:
599        Tmat = G2lat.U6toUij(TLS[:6])
600    if 'L' in TLStype:
601        Lmat = G2lat.U6toUij(TLS[6:12])
602    if 'S' in TLStype:
603        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
604    XYZ = np.inner(Amat,xyz)
605    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
606    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
607    beta = np.inner(np.inner(g,Umat),g)
608    return G2lat.UijtoU6(beta)*gvec   
609       
610def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
611    '''default doc string
612   
613    :param type name: description
614   
615    :returns: type name: description
616   
617    '''
618    cx,ct,cs,cia = atPtrs
619    TLStype,TLS = rbObj['ThermalMotion'][:2]
620    Tmat = np.zeros((3,3))
621    Lmat = np.zeros((3,3))
622    Smat = np.zeros((3,3))
623    G,g = G2lat.A2Gmat(Amat)
624    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
625    if 'T' in TLStype:
626        Tmat = G2lat.U6toUij(TLS[:6])
627    if 'L' in TLStype:
628        Lmat = G2lat.U6toUij(TLS[6:12])
629    if 'S' in TLStype:
630        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
631    for atom in atomData:
632        XYZ = np.inner(Amat,atom[cx:cx+3])
633        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
634        if 'U' in TSLtype:
635            atom[cia+1] = TLS[0]
636            atom[cia] = 'I'
637        else:
638            atom[cia] = 'A'
639            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
640            beta = np.inner(np.inner(g,Umat),g)
641            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
642
643def GetXYZDist(xyz,XYZ,Amat):
644    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
645        and are in crystal coordinates; Amat is crystal to Cart matrix
646   
647    :param type name: description
648   
649    :returns: type name: description
650   
651    '''
652    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
653
654def getAtomXYZ(atoms,cx):
655    '''default doc string
656   
657    :param type name: description
658   
659    :returns: type name: description
660   
661    '''
662    XYZ = []
663    for atom in atoms:
664        XYZ.append(atom[cx:cx+3])
665    return np.array(XYZ)
666
667def RotateRBXYZ(Bmat,Cart,oriQ):
668    '''rotate & transform cartesian coordinates to crystallographic ones
669    no translation applied. To be used for numerical derivatives
670   
671    :param type name: description
672   
673    :returns: type name: description
674   
675    '''
676    ''' returns crystal coordinates for atoms described by RBObj
677    '''
678    XYZ = np.zeros_like(Cart)
679    for i,xyz in enumerate(Cart):
680        XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
681    return XYZ
682
683def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
684    '''default doc string
685   
686    :param type name: description
687   
688    :returns: type name: description
689   
690    '''
691    ''' returns crystal coordinates for atoms described by RBObj
692    '''
693    RBRes = RBData[RBType][RBObj['RBId']]
694    if RBType == 'Vector':
695        vecs = RBRes['rbVect']
696        mags = RBRes['VectMag']
697        Cart = np.zeros_like(vecs[0])
698        for vec,mag in zip(vecs,mags):
699            Cart += vec*mag
700    elif RBType == 'Residue':
701        Cart = np.array(RBRes['rbXYZ'])
702        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
703            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
704            Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
705    XYZ = np.zeros_like(Cart)
706    for i,xyz in enumerate(Cart):
707        XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
708    return XYZ,Cart
709
710def UpdateMCSAxyz(Bmat,MCSA):
711    '''default doc string
712   
713    :param type name: description
714   
715    :returns: type name: description
716   
717    '''
718    xyz = []
719    atTypes = []
720    iatm = 0
721    for model in MCSA['Models'][1:]:        #skip the MD model
722        if model['Type'] == 'Atom':
723            xyz.append(model['Pos'][0])
724            atTypes.append(model['atType'])
725            iatm += 1
726        else:
727            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
728            Pos = np.array(model['Pos'][0])
729            Ori = np.array(model['Ori'][0])
730            Qori = AVdeg2Q(Ori[0],Ori[1:])
731            if model['Type'] == 'Vector':
732                vecs = RBRes['rbVect']
733                mags = RBRes['VectMag']
734                Cart = np.zeros_like(vecs[0])
735                for vec,mag in zip(vecs,mags):
736                    Cart += vec*mag
737            elif model['Type'] == 'Residue':
738                Cart = np.array(RBRes['rbXYZ'])
739                for itor,seq in enumerate(RBRes['rbSeq']):
740                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
741                    Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
742            if model['MolCent'][1]:
743                Cart -= model['MolCent'][0]
744            for i,x in enumerate(Cart):
745                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
746                atType = RBRes['rbTypes'][i]
747                atTypes.append(atType)
748                iatm += 1
749    return np.array(xyz),atTypes
750   
751def SetMolCent(model,RBData):
752    '''default doc string
753   
754    :param type name: description
755   
756    :returns: type name: description
757   
758    '''
759    rideList = []
760    RBRes = RBData[model['Type']][model['RBId']]
761    if model['Type'] == 'Vector':
762        vecs = RBRes['rbVect']
763        mags = RBRes['VectMag']
764        Cart = np.zeros_like(vecs[0])
765        for vec,mag in zip(vecs,mags):
766            Cart += vec*mag
767    elif model['Type'] == 'Residue':
768        Cart = np.array(RBRes['rbXYZ'])
769        for seq in RBRes['rbSeq']:
770            rideList += seq[3]
771    centList = set(range(len(Cart)))-set(rideList)
772    cent = np.zeros(3)
773    for i in centList:
774        cent += Cart[i]
775    model['MolCent'][0] = cent/len(centList)   
776   
777def UpdateRBUIJ(Bmat,Cart,RBObj):
778    '''default doc string
779   
780    :param type name: description
781   
782    :returns: type name: description
783   
784    '''
785    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
786    '''
787    TLStype,TLS = RBObj['ThermalMotion'][:2]
788    T = np.zeros(6)
789    L = np.zeros(6)
790    S = np.zeros(8)
791    if 'T' in TLStype:
792        T = TLS[:6]
793    if 'L' in TLStype:
794        L = np.array(TLS[6:12])*(np.pi/180.)**2
795    if 'S' in TLStype:
796        S = np.array(TLS[12:])*(np.pi/180.)
797    g = nl.inv(np.inner(Bmat,Bmat))
798    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
799        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
800    Uout = []
801    Q = RBObj['Orient'][0]
802    for X in Cart:
803        X = prodQVQ(Q,X)
804        if 'U' in TLStype:
805            Uout.append(['I',TLS[0],0,0,0,0,0,0])
806        elif not 'N' in TLStype:
807            U = [0,0,0,0,0,0]
808            U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
809            U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
810            U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
811            U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+  \
812                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
813            U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+  \
814                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
815            U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+  \
816                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
817            Umat = G2lat.U6toUij(U)
818            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
819            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
820        else:
821            Uout.append(['N',])
822    return Uout
823   
824def GetSHCoeff(pId,parmDict,SHkeys):
825    '''default doc string
826   
827    :param type name: description
828   
829    :returns: type name: description
830   
831    '''
832    SHCoeff = {}
833    for shkey in SHkeys:
834        shname = str(pId)+'::'+shkey
835        SHCoeff[shkey] = parmDict[shname]
836    return SHCoeff
837       
838def getMass(generalData):
839    '''Computes mass of unit cell contents
840   
841    :param dict generalData: The General dictionary in Phase
842   
843    :returns: float mass: Crystal unit cell mass in AMU.
844   
845    '''
846    mass = 0.
847    for i,elem in enumerate(generalData['AtomTypes']):
848        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
849    return max(mass,1.0)   
850
851def getDensity(generalData):
852    '''calculate crystal structure density
853   
854    :param dict generalData: The General dictionary in Phase
855   
856    :returns: float density: crystal density in gm/cm^3
857   
858    '''
859    mass = getMass(generalData)
860    Volume = generalData['Cell'][7]
861    density = mass/(0.6022137*Volume)
862    return density,Volume/mass
863   
864def getWave(Parms):
865    '''returns wavelength from Instrument parameters dictionary
866   
867    :param dict Parms: Instrument parameters;
868        must contain:
869        Lam: single wavelength
870        or
871        Lam1: Ka1 radiation wavelength
872   
873    :returns: float wave: wavelength
874   
875    '''
876    try:
877        return Parms['Lam'][1]
878    except KeyError:
879        return Parms['Lam1'][1]
880       
881def getMeanWave(Parms):
882    '''returns mean wavelength from Instrument parameters dictionary
883   
884    :param dict Parms: Instrument parameters;
885        must contain:
886        Lam: single wavelength
887        or
888        Lam1,Lam2: Ka1,Ka2 radiation wavelength
889        I(L2)/I(L1): Ka2/Ka1 ratio
890   
891    :returns: float wave: mean wavelength
892   
893    '''
894    try:
895        return Parms['Lam'][1]
896    except KeyError:
897        meanLam = (Parms['Lam1'][1]+Parms['I(L2)/I(L1)'][1]*Parms['Lam2'][1])/   \
898            (1.+Parms['I(L2)/I(L1)'][1])
899        return meanLam
900   
901       
902def El2Mass(Elements):
903    '''compute molecular weight from Elements
904   
905    :param dict Elements: elements in molecular formula;
906        each must contain
907        Num: number of atoms in formula
908        Mass: at. wt.
909   
910    :returns: float mass: molecular weight.
911   
912    '''
913    mass = 0
914    for El in Elements:
915        mass += Elements[El]['Num']*Elements[El]['Mass']
916    return mass
917       
918def Den2Vol(Elements,density):
919    '''converts density to molecular volume
920   
921    :param dict Elements: elements in molecular formula;
922        each must contain
923        Num: number of atoms in formula
924        Mass: at. wt.
925    :param float density: material density in gm/cm^3
926   
927    :returns: float volume: molecular volume in A^3
928   
929    '''
930    return El2Mass(Elements)/(density*0.6022137)
931   
932def Vol2Den(Elements,volume):
933    '''converts volume to density
934   
935    :param dict Elements: elements in molecular formula;
936        each must contain
937        Num: number of atoms in formula
938        Mass: at. wt.
939    :param float volume: molecular volume in A^3
940   
941    :returns: float density: material density in gm/cm^3
942   
943    '''
944    return El2Mass(Elements)/(volume*0.6022137)
945   
946def El2EstVol(Elements):
947    '''Estimate volume from molecular formula; assumes atom volume = 10A^3
948   
949    :param dict Elements: elements in molecular formula;
950        each must contain
951        Num: number of atoms in formula
952   
953    :returns: float volume: estimate of molecular volume in A^3
954   
955    '''
956    vol = 0
957    for El in Elements:
958        vol += 10.*Elements[El]['Num']
959    return vol
960   
961def XScattDen(Elements,vol,wave=0.):
962    '''Estimate X-ray scattering density from molecular formula & volume;
963    ignores valence, but includes anomalous effects
964   
965    :param dict Elements: elements in molecular formula;
966        each element must contain
967        Num: number of atoms in formula
968        Z: atomic number
969    :param float vol: molecular volume in A^3
970    :param float wave: optional wavelength in A
971   
972    :returns: float rho: scattering density in 10^10cm^-2;
973        if wave > 0 the includes f' contribution
974    :returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
975   
976    '''
977    rho = 0
978    mu = 0
979    if wave:
980        Xanom = XAnomAbs(Elements,wave)
981    for El in Elements:
982        f0 = Elements[El]['Z']
983        if wave:
984            f0 += Xanom[El][0]
985            mu += Xanom[El][2]*Elements[El]['Num']
986        rho += Elements[El]['Num']*f0
987    return 28.179*rho/vol,0.1*mu/vol
988   
989def wavekE(wavekE):
990    '''Convert wavelength to energy & vise versa
991   
992    :param float waveKe:wavelength in A or energy in kE
993   
994    :returns float waveKe:the other one
995   
996    '''
997    return 12.397639/wavekE
998       
999def XAnomAbs(Elements,wave):
1000    kE = wavekE(wave)
1001    Xanom = {}
1002    for El in Elements:
1003        Orbs = G2el.GetXsectionCoeff(El)
1004        Xanom[El] = G2el.FPcalc(Orbs, kE)
1005    return Xanom
1006   
1007################################################################################
1008#### Modulation math
1009################################################################################
1010
1011def makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast):
1012    '''
1013    waveTypes: array nAtoms: 'Fourier','ZigZag' or 'Block'
1014    FSSdata: array 2 x atoms x waves    (sin,cos terms)
1015    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1016    USSdata: array 2x6 x atoms X waves (sin,cos terms)
1017    Mast: array orthogonalization matrix for Uij
1018    '''
1019    ngl = 32
1020    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
1021    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
1022    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
1023    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
1024    Bf = np.array(FSSdata[1]).T    #cos frac mods...
1025    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods as betaij
1026    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods as betaij
1027    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 
1028    if nWaves[0]:
1029        tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl
1030        FmodA = Af[:,:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X ngl
1031        FmodB = Bf[:,:,nxs]*np.cos(twopi*tauF[nxs,:,:])
1032        Fmod = np.sum(1.0+FmodA+FmodB,axis=1)             #atoms X ngl; sum waves
1033    else:
1034        Fmod = 1.0           
1035    XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1036    XmodA = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1037    XmodB = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1038    for iatm in range(Ax.shape[0]):
1039        nx = 0
1040        if 'ZigZag' in waveTypes[iatm]:
1041            nx = 1
1042            Tmm = Ax[iatm][0][:2]                       
1043            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1044            XmodZ[iatm][0] += posZigZag(glTau,Tmm,XYZmax).T
1045        elif 'Block' in waveTypes[iatm]:
1046            nx = 1
1047            Tmm = Ax[iatm][0][:2]                       
1048            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1049            XmodZ[iatm][0] += posBlock(glTau,Tmm,XYZmax).T
1050        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
1051        if nx:   
1052            XmodA[iatm][:-nx] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1053            XmodB[iatm][:-nx] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
1054        else:
1055            XmodA[iatm] = Ax[iatm,:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1056            XmodB[iatm] = Bx[iatm,:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
1057    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X ngl; sum waves
1058    Xmod = np.swapaxes(Xmod,1,2)
1059    if nWaves[2]:
1060        tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
1061        UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x ngl
1062        UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
1063        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x ngl; sum waves
1064    else:
1065        Umod = 1.0
1066#    GSASIIpath.IPyBreak()
1067    return ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt
1068       
1069def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1070    '''
1071    H: array nRefBlk x ops X hklt
1072    HP: array nRefBlk x ops X hklt proj to hkl
1073    Fmod: array 2 x atoms x waves    (sin,cos terms)
1074    Xmod: array atoms X 3 X ngl
1075    Umod: array atoms x 3x3 x ngl
1076    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1077    '''
1078   
1079    if nWaves[2]:
1080        if len(HP.shape) > 2:
1081            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1082        else:
1083            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1084    else:
1085        HbH = 1.0
1086    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
1087    if len(H.shape) > 2:
1088        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1089        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1090    else:
1091        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1092        HdotXD = twopi*(HdotX+D[:,nxs,:])
1093    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1094    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1095    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1096   
1097def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1098    '''
1099    H: array nRefBlk x tw x ops X hklt
1100    HP: array nRefBlk x tw x ops X hklt proj to hkl
1101    Fmod: array 2 x atoms x waves    (sin,cos terms)
1102    Xmod: array atoms X ngl X 3
1103    Umod: array atoms x ngl x 3x3
1104    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1105    '''
1106   
1107    if nWaves[2]:
1108        if len(HP.shape) > 3:   #Blocks of reflections
1109            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1110        else:   #single reflections
1111            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1112    else:
1113        HbH = 1.0
1114    HdotX = np.inner(HP,Xmod)                   #refBlk x tw x ops x atoms X ngl
1115    if len(H.shape) > 3:
1116        D = glTau*H[:,:,:,3:,nxs]              #m*e*tau; refBlk x tw x ops X ngl
1117        HdotXD = twopi*(HdotX+D[:,:,:,nxs,:])
1118    else:
1119        D = H*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1120        HdotXD = twopi*(HdotX+D[:,nxs,:])
1121    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1122    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1123    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1124   
1125def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast):
1126    '''
1127    FSSdata: array 2 x atoms x waves    (sin,cos terms)
1128    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1129    USSdata: array 2x6 x atoms X waves (sin,cos terms)
1130    Mast: array orthogonalization matrix for Uij
1131    '''
1132    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
1133    dT = 2./ngl
1134    dX = 0.0001
1135    waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
1136    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
1137    Bf = np.array(FSSdata[1]).T    #cos frac mods...
1138    Ax = np.array(XSSdata[:3]).T   #...cos pos mods x waves x atoms
1139    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
1140    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
1141    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
1142    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 
1143    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
1144    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1145    ZtauXt = np.zeros((Ax.shape[0],2,3,ngl))               #atoms x Tminmax x 3 x ngl
1146    ZtauXx = np.zeros((Ax.shape[0],3,ngl))               #atoms x XYZmax x ngl
1147    for iatm in range(Ax.shape[0]):
1148        nx = 0
1149        if 'ZigZag' in waveTypes[iatm]:
1150            nx = 1
1151            Tmm = Ax[iatm][0][:2]                       
1152            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
1153            ZtauXt[iatm],ZtauXx[iatm] = posZigZagDerv(glTau,Tmm,XYZmax)
1154        elif 'Block' in waveTypes[iatm]:
1155            nx = 1
1156            Tmm = Ax[iatm][0][:2]                       
1157            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
1158            ZtauXt[iatm],ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax)
1159        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
1160        if nx:   
1161            StauX[iatm][:-nx] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1162            CtauX[iatm][:-nx] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1163        else:
1164            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1165            CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1166#    GSASIIpath.IPyBreak()
1167    if nWaves[0]:
1168        tauF = np.arange(1.,nWaves[0]+1-nf)[:,nxs]*glTau  #Fwaves x ngl
1169        StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
1170        CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
1171    else:
1172        StauF = 1.0
1173        CtauF = 1.0
1174    if nWaves[2]:
1175        tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
1176        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodA/dAu
1177        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodB/dBu
1178        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x ngl
1179        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
1180#derivs need to be ops x atoms x waves x 6uij; ops x atoms x waves x ngl x 6uij before sum
1181        StauU = np.rollaxis(np.rollaxis(np.swapaxes(StauU,2,4),-1),-1)
1182        CtauU = np.rollaxis(np.rollaxis(np.swapaxes(CtauU,2,4),-1),-1)
1183    else:
1184        StauU = 1.0
1185        CtauU = 1.0
1186        UmodA = 0.
1187        UmodB = 0.
1188    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
1189   
1190def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1191    '''
1192    H: array ops X hklt proj to hkl
1193    HP: array ops X hklt proj to hkl
1194    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
1195    '''
1196   
1197    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
1198    dGdMfC = np.zeros(Mf)
1199    dGdMfS = np.zeros(Mf)
1200    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
1201    dGdMxC = np.zeros(Mx)
1202    dGdMxS = np.zeros(Mx)
1203    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
1204    dGdMuC = np.zeros(Mu)
1205    dGdMuS = np.zeros(Mu)
1206   
1207    D = twopi*H[:,3][:,nxs]*glTau[nxs,:]              #m*e*tau; ops X ngl
1208    HdotX = twopi*np.inner(HP,Xmod)        #ops x atoms X ngl
1209    HdotXD = HdotX+D[:,nxs,:]
1210    if nWaves[2]:
1211        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1212        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1213        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1214        HbH = np.exp(-np.sum(HuH,axis=-2)) # ops x atoms x ngl; sum waves - OK vs Modulation version
1215        part1 = -np.exp(-HuH)*Fmod    #ops x atoms x waves x ngl
1216        dUdAu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
1217        dUdBu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
1218        dGdMuCa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1219        dGdMuCb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1220        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
1221        dGdMuSa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1222        dGdMuSb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1223        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
1224    else:
1225        HbH = np.ones_like(HdotX)
1226    dHdXA = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
1227    dHdXB = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:]    #ditto - cos waves
1228# ops x atoms x waves x 2xyz - real part - good
1229    dGdMxCa = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1230    dGdMxCb = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1231    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1232# ops x atoms x waves x 2xyz - imag part - good
1233    dGdMxSa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1234    dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1235    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1236# ZigZag/Block waves - problems here?
1237    dHdXZt = -twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
1238    dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
1239    dGdMzCt = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1240    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
1241    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
1242    dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1243    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
1244    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
1245#    GSASIIpath.IPyBreak()
1246    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
1247   
1248def ModulationDerv2(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1249    '''
1250    H: array refBlk x ops X hklt proj to hkl
1251    HP: array refBlk x ops X hklt proj to hkl
1252    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
1253    '''
1254   
1255    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
1256    dGdMfC = np.zeros(Mf)
1257    dGdMfS = np.zeros(Mf)
1258    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
1259    dGdMxC = np.zeros(Mx)
1260    dGdMxS = np.zeros(Mx)
1261    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
1262    dGdMuC = np.zeros(Mu)
1263    dGdMuS = np.zeros(Mu)
1264   
1265    D = twopi*H[:,:,3,nxs]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1266    HdotX = twopi*np.inner(HP,Xmod)        #refBlk x ops x atoms X ngl
1267    HdotXD = HdotX+D[:,:,nxs,:]
1268    if nWaves[2]:
1269        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1270        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
1271        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
1272        HbH = np.exp(-np.sum(HuH,axis=-2)) #refBlk x ops x atoms x ngl; sum waves - OK vs Modulation version
1273        part1 = -np.exp(-HuH)*Fmod    #refBlk x ops x atoms x waves x ngl
1274        dUdAu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
1275        dUdBu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
1276        dGdMuCa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1277        dGdMuCb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1278        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
1279        dGdMuSa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1280        dGdMuSb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1281        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
1282    else:
1283        HbH = np.ones_like(HdotX)
1284    dHdXA = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
1285    dHdXB = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,nxs,:,:,:,:]    #ditto - cos waves
1286# ops x atoms x waves x 2xyz - real part - good
1287    dGdMxCa = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1288    dGdMxCb = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1289    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1290# ops x atoms x waves x 2xyz - imag part - good
1291    dGdMxSa = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
1292    dGdMxSb = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
1293    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1294# ZigZag/Block waves - problems here?
1295    dHdXZt = -twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
1296    dHdXZx = twopi*HP[:,:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
1297    dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1298    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1299    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
1300    dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1301    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1302    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
1303#    GSASIIpath.IPyBreak()
1304    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
1305   
1306def posFourier(tau,psin,pcos):
1307    A = np.array([ps[:,nxs]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
1308    B = np.array([pc[:,nxs]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
1309    return np.sum(A,axis=0)+np.sum(B,axis=0)
1310   
1311def posZigZag(T,Tmm,Xmax):
1312    DT = Tmm[1]-Tmm[0]
1313    Su = 2.*Xmax/DT
1314    Sd = 2.*Xmax/(1.-DT)
1315    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax+Su*((t-Tmm[0])%1.),Xmax-Sd*((t-Tmm[1])%1.)) for t in T])
1316    return A
1317   
1318def posZigZagDerv(T,Tmm,Xmax):
1319    DT = Tmm[1]-Tmm[0]
1320    Su = 2.*Xmax/DT
1321    Sd = 2.*Xmax/(1.-DT)
1322    dAdT = np.zeros((2,3,len(T)))
1323    dAdT[0] = np.array([np.where(Tmm[0] < t <= Tmm[1],Su*(t-Tmm[0]-1)/DT,-Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1324    dAdT[1] = np.array([np.where(Tmm[0] < t <= Tmm[1],-Su*(t-Tmm[0])/DT,Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1325    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.+2.*(t-Tmm[0])/DT,1.-2.*(t-Tmm[1])%1./DT) for t in T])
1326    return dAdT,dAdX
1327
1328def posBlock(T,Tmm,Xmax):
1329    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T])
1330    return A
1331   
1332def posBlockDerv(T,Tmm,Xmax):
1333    dAdT = np.zeros((2,3,len(T)))
1334    ind = np.searchsorted(T,Tmm)
1335    dAdT[0,:,ind[0]] = -Xmax/len(T)
1336    dAdT[1,:,ind[1]] = Xmax/len(T)
1337    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T])  #OK
1338    return dAdT,dAdX
1339   
1340def fracCrenel(tau,Toff,Twid):
1341    Tau = (tau-Toff)%1.
1342    A = np.where(Tau<Twid,1.,0.)
1343    return A
1344   
1345def fracFourier(tau,fsin,fcos):
1346    A = np.array([fs[:,nxs]*np.sin(2.*np.pi*(i+1)*tau) for i,fs in enumerate(fsin)])
1347    B = np.array([fc[:,nxs]*np.cos(2.*np.pi*(i+1)*tau) for i,fc in enumerate(fcos)])
1348    return np.sum(A,axis=0)+np.sum(B,axis=0)
1349   
1350def ApplyModulation(data,tau):
1351    '''Applies modulation to drawing atom positions & Uijs for given tau
1352    '''
1353    generalData = data['General']
1354    SGData = generalData['SGData']
1355    SSGData = generalData['SSGData']
1356    cx,ct,cs,cia = generalData['AtomPtrs']
1357    drawingData = data['Drawing']
1358    dcx,dct,dcs,dci = drawingData['atomPtrs']
1359    atoms = data['Atoms']
1360    drawAtoms = drawingData['Atoms']
1361    Fade = np.zeros(len(drawAtoms))
1362    for atom in atoms:   
1363        atxyz = G2spc.MoveToUnitCell(np.array(atom[cx:cx+3]))[0]
1364        atuij = np.array(atom[cia+2:cia+8])
1365        waveType = atom[-1]['SS1']['waveType']
1366        Sfrac = atom[-1]['SS1']['Sfrac']
1367        Spos = atom[-1]['SS1']['Spos']
1368        Sadp = atom[-1]['SS1']['Sadp']
1369        indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
1370        for ind in indx:
1371            drawatom = drawAtoms[ind]
1372            opr = drawatom[dcs-1]
1373            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
1374            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
1375            tauT *= icent       #invert wave on -1
1376            wave = np.zeros(3)
1377            uwave = np.zeros(6)
1378            #how do I handle Sfrac? - fade the atoms?
1379            if len(Sfrac):
1380                scof = []
1381                ccof = []
1382                for i,sfrac in enumerate(Sfrac):
1383                    if not i and 'Crenel' in waveType:
1384                        Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
1385                    else:
1386                        scof.append(sfrac[0][0])
1387                        ccof.append(sfrac[0][1])
1388                    if len(scof):
1389                        Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))                           
1390            if len(Spos):
1391                scof = []
1392                ccof = []
1393                for i,spos in enumerate(Spos):
1394                    if waveType in ['ZigZag','Block'] and not i:
1395                        Tminmax = spos[0][:2]
1396                        XYZmax = np.array(spos[0][2:])
1397                        if waveType == 'Block':
1398                            wave = np.array(posBlock([tauT,],Tminmax,XYZmax))[0]
1399                        elif waveType == 'ZigZag':
1400                            wave = np.array(posZigZag([tauT,],Tminmax,XYZmax))[0]
1401                    else:
1402                        scof.append(spos[0][:3])
1403                        ccof.append(spos[0][3:])
1404                if len(scof):
1405                    wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1406            if len(Sadp):
1407                scof = []
1408                ccof = []
1409                for i,sadp in enumerate(Sadp):
1410                    scof.append(sadp[0][:6])
1411                    ccof.append(sadp[0][6:])
1412                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1413            if atom[cia] == 'A':                   
1414                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1415                drawatom[dcx:dcx+3] = X
1416                drawatom[dci-6:dci] = U
1417            else:
1418                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1419                drawatom[dcx:dcx+3] = X
1420    return drawAtoms,Fade
1421   
1422# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1423# integrate from a to b using n evaluations of the function f(x) 
1424# usage: from gauleg import gaulegf         
1425#        x,w = gaulegf( a, b, n)                               
1426#        area = 0.0                                           
1427#        for i in range(1,n+1):          #  yes, 1..n                   
1428#          area += w[i]*f(x[i])                                   
1429
1430import math
1431def gaulegf(a, b, n):
1432  x = range(n+1) # x[0] unused
1433  w = range(n+1) # w[0] unused
1434  eps = 3.0E-14
1435  m = (n+1)/2
1436  xm = 0.5*(b+a)
1437  xl = 0.5*(b-a)
1438  for i in range(1,m+1):
1439    z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1440    while True:
1441      p1 = 1.0
1442      p2 = 0.0
1443      for j in range(1,n+1):
1444        p3 = p2
1445        p2 = p1
1446        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
1447
1448      pp = n*(z*p1-p2)/(z*z-1.0)
1449      z1 = z
1450      z = z1 - p1/pp
1451      if abs(z-z1) <= eps:
1452            break
1453
1454    x[i] = xm - xl*z
1455    x[n+1-i] = xm + xl*z
1456    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
1457    w[n+1-i] = w[i]
1458  return np.array(x), np.array(w)
1459# end gaulegf
1460   
1461   
1462def BessJn(nmax,x):
1463    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1464    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1465   
1466    :param  integer nmax: maximul order for Jn(x)
1467    :param  float x: argument for Jn(x)
1468   
1469    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1470   
1471    '''
1472    import scipy.special as sp
1473    bessJn = np.zeros(2*nmax+1)
1474    bessJn[nmax] = sp.j0(x)
1475    bessJn[nmax+1] = sp.j1(x)
1476    bessJn[nmax-1] = -bessJn[nmax+1]
1477    for i in range(2,nmax+1):
1478        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1479        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
1480    return bessJn
1481   
1482def BessIn(nmax,x):
1483    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1484    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1485   
1486    :param  integer nmax: maximul order for In(x)
1487    :param  float x: argument for In(x)
1488   
1489    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1490   
1491    '''
1492    import scipy.special as sp
1493    bessIn = np.zeros(2*nmax+1)
1494    bessIn[nmax] = sp.i0(x)
1495    bessIn[nmax+1] = sp.i1(x)
1496    bessIn[nmax-1] = bessIn[nmax+1]
1497    for i in range(2,nmax+1):
1498        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
1499        bessIn[nmax-i] = bessIn[i+nmax]
1500    return bessIn
1501       
1502   
1503################################################################################
1504##### distance, angle, planes, torsion stuff
1505################################################################################
1506
1507def CalcDist(distance_dict, distance_atoms, parmDict):
1508    if not len(parmDict):
1509        return 0.
1510    pId = distance_dict['pId']
1511    pfx = '%d::'%(pId)
1512    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1513    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1514    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1515    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1516    inv = 1
1517    symNo = distance_dict['symNo']
1518    if symNo < 0:
1519        inv = -1
1520        symNo *= -1
1521    cen = symNo/100
1522    op = symNo%100-1
1523    M,T = distance_dict['SGData']['SGOps'][op]
1524    D = T*inv+distance_dict['SGData']['SGCen'][cen]
1525    D += distance_dict['cellNo']
1526    Txyz = np.inner(M*inv,Txyz)+D
1527    dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1528#    GSASIIpath.IPyBreak()
1529    return dist   
1530   
1531def CalcDistSig(distance_dict, distance_atoms, parmDict,covData={}):
1532    '''default doc string
1533   
1534    :param type name: description
1535   
1536    :returns: type name: description
1537   
1538    '''
1539    def calcDist(Atoms,SyOps,Amat):
1540        XYZ = []
1541        for i,atom in enumerate(Atoms):
1542            Inv,M,T,C,U = SyOps[i]
1543            XYZ.append(np.array(atom[1:4]))
1544            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1545            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1546        V1 = XYZ[1]-XYZ[0]
1547        return np.sqrt(np.sum(V1**2))
1548       
1549    if not len(parmDict):
1550        return 0.
1551    pId = distance_dict['pId']
1552    pfx = '%d::'%(pId)
1553    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1554    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1555    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1556    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1557    inv = 1
1558    symNo = distance_dict['symNo']
1559    if symNo < 0:
1560        inv = -1
1561        symNo *= -1
1562    cen = symNo/100
1563    op = symNo%100-1
1564    M,T = distance_dict['SGData']['SGOps'][op]
1565    D = T*inv+distance_dict['SGData']['SGCen'][cen]
1566    D += distance_dict['cellNo']
1567    Txyz = np.inner(M*inv,Txyz)+D
1568    dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1569
1570    sig = -0.001
1571    if 'covMatrix' in covData:
1572        parmNames = []
1573        dx = .00001
1574        dadx = np.zeros(6)
1575        for i in range(6):
1576            ia = i/3
1577            ix = i%3
1578            Oxyz[ix+1] += dx
1579            a0 = calcDist(Oatoms,SyOps,Amat)
1580            Oatoms[ia][ix+1] -= 2*dx
1581            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1582        covMatrix = covData['covMatrix']
1583        varyList = covData['varyList']
1584        DistVcov = getVCov(names,varyList,covMatrix)
1585        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1586        if sig < 0.001:
1587            sig = -0.001
1588   
1589    return Dist,sig
1590
1591def CalcAngle(distance_dict, distance_atoms, parmDict):
1592    if not len(parmDict):
1593        return 0.
1594
1595    return 0.   #angle
1596
1597def getSyXYZ(XYZ,ops,SGData):
1598    '''default doc string
1599   
1600    :param type name: description
1601   
1602    :returns: type name: description
1603   
1604    '''
1605    XYZout = np.zeros_like(XYZ)
1606    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1607        if op == '1':
1608            XYZout[i] = xyz
1609        else:
1610            oprs = op.split('+')
1611            unit = [0,0,0]
1612            if len(oprs)>1:
1613                unit = np.array(list(eval(oprs[1])))
1614            syop =int(oprs[0])
1615            inv = syop/abs(syop)
1616            syop *= inv
1617            cent = syop/100
1618            syop %= 100
1619            syop -= 1
1620            M,T = SGData['SGOps'][syop]
1621            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1622    return XYZout
1623   
1624def getRestDist(XYZ,Amat):
1625    '''default doc string
1626   
1627    :param type name: description
1628   
1629    :returns: type name: description
1630   
1631    '''
1632    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1633   
1634def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1635    '''default doc string
1636   
1637    :param type name: description
1638   
1639    :returns: type name: description
1640   
1641    '''
1642    deriv = np.zeros((len(XYZ),3))
1643    dx = 0.00001
1644    for j,xyz in enumerate(XYZ):
1645        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1646            XYZ[j] -= x
1647            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1648            XYZ[j] += 2*x
1649            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1650            XYZ[j] -= x
1651            deriv[j][i] = (d1-d2)/(2*dx)
1652    return deriv.flatten()
1653
1654def getRestAngle(XYZ,Amat):
1655    '''default doc string
1656   
1657    :param type name: description
1658   
1659    :returns: type name: description
1660   
1661    '''
1662   
1663    def calcVec(Ox,Tx,Amat):
1664        return np.inner(Amat,(Tx-Ox))
1665
1666    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1667    VecA /= np.sqrt(np.sum(VecA**2))
1668    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1669    VecB /= np.sqrt(np.sum(VecB**2))
1670    edge = VecB-VecA
1671    edge = np.sum(edge**2)
1672    angle = (2.-edge)/2.
1673    angle = max(angle,-1.)
1674    return acosd(angle)
1675   
1676def getRestPlane(XYZ,Amat):
1677    '''default doc string
1678   
1679    :param type name: description
1680   
1681    :returns: type name: description
1682   
1683    '''
1684    sumXYZ = np.zeros(3)
1685    for xyz in XYZ:
1686        sumXYZ += xyz
1687    sumXYZ /= len(XYZ)
1688    XYZ = np.array(XYZ)-sumXYZ
1689    XYZ = np.inner(Amat,XYZ).T
1690    Zmat = np.zeros((3,3))
1691    for i,xyz in enumerate(XYZ):
1692        Zmat += np.outer(xyz.T,xyz)
1693    Evec,Emat = nl.eig(Zmat)
1694    Evec = np.sqrt(Evec)/(len(XYZ)-3)
1695    Order = np.argsort(Evec)
1696    return Evec[Order[0]]
1697   
1698def getRestChiral(XYZ,Amat):   
1699    '''default doc string
1700   
1701    :param type name: description
1702   
1703    :returns: type name: description
1704   
1705    '''
1706    VecA = np.empty((3,3))   
1707    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1708    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
1709    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
1710    return nl.det(VecA)
1711   
1712def getRestTorsion(XYZ,Amat):
1713    '''default doc string
1714   
1715    :param type name: description
1716   
1717    :returns: type name: description
1718   
1719    '''
1720    VecA = np.empty((3,3))
1721    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1722    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
1723    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
1724    D = nl.det(VecA)
1725    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
1726    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
1727    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
1728    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
1729    Ang = 1.0
1730    if abs(P12) < 1.0 and abs(P23) < 1.0:
1731        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
1732    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
1733    return TOR
1734   
1735def calcTorsionEnergy(TOR,Coeff=[]):
1736    '''default doc string
1737   
1738    :param type name: description
1739   
1740    :returns: type name: description
1741   
1742    '''
1743    sum = 0.
1744    if len(Coeff):
1745        cof = np.reshape(Coeff,(3,3)).T
1746        delt = TOR-cof[1]
1747        delt = np.where(delt<-180.,delt+360.,delt)
1748        delt = np.where(delt>180.,delt-360.,delt)
1749        term = -cof[2]*delt**2
1750        val = cof[0]*np.exp(term/1000.0)
1751        pMax = cof[0][np.argmin(val)]
1752        Eval = np.sum(val)
1753        sum = Eval-pMax
1754    return sum,Eval
1755
1756def getTorsionDeriv(XYZ,Amat,Coeff):
1757    '''default doc string
1758   
1759    :param type name: description
1760   
1761    :returns: type name: description
1762   
1763    '''
1764    deriv = np.zeros((len(XYZ),3))
1765    dx = 0.00001
1766    for j,xyz in enumerate(XYZ):
1767        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1768            XYZ[j] -= x
1769            tor = getRestTorsion(XYZ,Amat)
1770            p1,d1 = calcTorsionEnergy(tor,Coeff)
1771            XYZ[j] += 2*x
1772            tor = getRestTorsion(XYZ,Amat)
1773            p2,d2 = calcTorsionEnergy(tor,Coeff)           
1774            XYZ[j] -= x
1775            deriv[j][i] = (p2-p1)/(2*dx)
1776    return deriv.flatten()
1777
1778def getRestRama(XYZ,Amat):
1779    '''Computes a pair of torsion angles in a 5 atom string
1780   
1781    :param nparray XYZ: crystallographic coordinates of 5 atoms
1782    :param nparray Amat: crystal to cartesian transformation matrix
1783   
1784    :returns: list (phi,psi) two torsion angles in degrees
1785   
1786    '''
1787    phi = getRestTorsion(XYZ[:5],Amat)
1788    psi = getRestTorsion(XYZ[1:],Amat)
1789    return phi,psi
1790   
1791def calcRamaEnergy(phi,psi,Coeff=[]):
1792    '''Computes pseudo potential energy from a pair of torsion angles and a
1793    numerical description of the potential energy surface. Used to create
1794    penalty function in LS refinement:     
1795    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
1796
1797    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
1798   
1799    :param float phi: first torsion angle (:math:`\\phi`)
1800    :param float psi: second torsion angle (:math:`\\psi`)
1801    :param list Coeff: pseudo potential coefficients
1802   
1803    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
1804      sum is used for penalty function.
1805   
1806    '''
1807    sum = 0.
1808    Eval = 0.
1809    if len(Coeff):
1810        cof = Coeff.T
1811        dPhi = phi-cof[1]
1812        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
1813        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
1814        dPsi = psi-cof[2]
1815        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
1816        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
1817        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
1818        val = cof[0]*np.exp(val/1000.)
1819        pMax = cof[0][np.argmin(val)]
1820        Eval = np.sum(val)
1821        sum = Eval-pMax
1822    return sum,Eval
1823
1824def getRamaDeriv(XYZ,Amat,Coeff):
1825    '''Computes numerical derivatives of torsion angle pair pseudo potential
1826    with respect of crystallographic atom coordinates of the 5 atom sequence
1827   
1828    :param nparray XYZ: crystallographic coordinates of 5 atoms
1829    :param nparray Amat: crystal to cartesian transformation matrix
1830    :param list Coeff: pseudo potential coefficients
1831   
1832    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
1833     crystallographic xyz coordinates.
1834   
1835    '''
1836    deriv = np.zeros((len(XYZ),3))
1837    dx = 0.00001
1838    for j,xyz in enumerate(XYZ):
1839        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1840            XYZ[j] -= x
1841            phi,psi = getRestRama(XYZ,Amat)
1842            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
1843            XYZ[j] += 2*x
1844            phi,psi = getRestRama(XYZ,Amat)
1845            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
1846            XYZ[j] -= x
1847            deriv[j][i] = (p2-p1)/(2*dx)
1848    return deriv.flatten()
1849
1850def getRestPolefig(ODFln,SamSym,Grid):
1851    '''default doc string
1852   
1853    :param type name: description
1854   
1855    :returns: type name: description
1856   
1857    '''
1858    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
1859    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
1860    R = np.where(R <= 1.,2.*atand(R),0.0)
1861    Z = np.zeros_like(R)
1862    Z = G2lat.polfcal(ODFln,SamSym,R,P)
1863    Z = np.reshape(Z,(Grid,Grid))
1864    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
1865
1866def getRestPolefigDerv(HKL,Grid,SHCoeff):
1867    '''default doc string
1868   
1869    :param type name: description
1870   
1871    :returns: type name: description
1872   
1873    '''
1874    pass
1875       
1876def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
1877    '''default doc string
1878   
1879    :param type name: description
1880   
1881    :returns: type name: description
1882   
1883    '''
1884    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
1885        TxT = inv*(np.inner(M,Tx)+T)+C+U
1886        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
1887       
1888    inv = Top/abs(Top)
1889    cent = abs(Top)/100
1890    op = abs(Top)%100-1
1891    M,T = SGData['SGOps'][op]
1892    C = SGData['SGCen'][cent]
1893    dx = .00001
1894    deriv = np.zeros(6)
1895    for i in [0,1,2]:
1896        Oxyz[i] -= dx
1897        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1898        Oxyz[i] += 2*dx
1899        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1900        Oxyz[i] -= dx
1901        Txyz[i] -= dx
1902        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1903        Txyz[i] += 2*dx
1904        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1905        Txyz[i] -= dx
1906    return deriv
1907   
1908def getAngSig(VA,VB,Amat,SGData,covData={}):
1909    '''default doc string
1910   
1911    :param type name: description
1912   
1913    :returns: type name: description
1914   
1915    '''
1916    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1917        TxT = inv*(np.inner(M,Tx)+T)+C+U
1918        return np.inner(Amat,(TxT-Ox))
1919       
1920    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1921        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1922        VecA /= np.sqrt(np.sum(VecA**2))
1923        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1924        VecB /= np.sqrt(np.sum(VecB**2))
1925        edge = VecB-VecA
1926        edge = np.sum(edge**2)
1927        angle = (2.-edge)/2.
1928        angle = max(angle,-1.)
1929        return acosd(angle)
1930       
1931    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1932    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1933    invA = invB = 1
1934    invA = TopA/abs(TopA)
1935    invB = TopB/abs(TopB)
1936    centA = abs(TopA)/100
1937    centB = abs(TopB)/100
1938    opA = abs(TopA)%100-1
1939    opB = abs(TopB)%100-1
1940    MA,TA = SGData['SGOps'][opA]
1941    MB,TB = SGData['SGOps'][opB]
1942    CA = SGData['SGCen'][centA]
1943    CB = SGData['SGCen'][centB]
1944    if 'covMatrix' in covData:
1945        covMatrix = covData['covMatrix']
1946        varyList = covData['varyList']
1947        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1948        dx = .00001
1949        dadx = np.zeros(9)
1950        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1951        for i in [0,1,2]:
1952            OxA[i] -= dx
1953            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1954            OxA[i] += 2*dx
1955            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1956            OxA[i] -= dx
1957           
1958            TxA[i] -= dx
1959            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1960            TxA[i] += 2*dx
1961            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1962            TxA[i] -= dx
1963           
1964            TxB[i] -= dx
1965            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1966            TxB[i] += 2*dx
1967            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1968            TxB[i] -= dx
1969           
1970        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1971        if sigAng < 0.01:
1972            sigAng = 0.0
1973        return Ang,sigAng
1974    else:
1975        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1976
1977def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1978    '''default doc string
1979   
1980    :param type name: description
1981   
1982    :returns: type name: description
1983   
1984    '''
1985    def calcDist(Atoms,SyOps,Amat):
1986        XYZ = []
1987        for i,atom in enumerate(Atoms):
1988            Inv,M,T,C,U = SyOps[i]
1989            XYZ.append(np.array(atom[1:4]))
1990            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1991            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1992        V1 = XYZ[1]-XYZ[0]
1993        return np.sqrt(np.sum(V1**2))
1994       
1995    Inv = []
1996    SyOps = []
1997    names = []
1998    for i,atom in enumerate(Oatoms):
1999        names += atom[-1]
2000        Op,unit = Atoms[i][-1]
2001        inv = Op/abs(Op)
2002        m,t = SGData['SGOps'][abs(Op)%100-1]
2003        c = SGData['SGCen'][abs(Op)/100]
2004        SyOps.append([inv,m,t,c,unit])
2005    Dist = calcDist(Oatoms,SyOps,Amat)
2006   
2007    sig = -0.001
2008    if 'covMatrix' in covData:
2009        parmNames = []
2010        dx = .00001
2011        dadx = np.zeros(6)
2012        for i in range(6):
2013            ia = i/3
2014            ix = i%3
2015            Oatoms[ia][ix+1] += dx
2016            a0 = calcDist(Oatoms,SyOps,Amat)
2017            Oatoms[ia][ix+1] -= 2*dx
2018            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2019        covMatrix = covData['covMatrix']
2020        varyList = covData['varyList']
2021        DistVcov = getVCov(names,varyList,covMatrix)
2022        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
2023        if sig < 0.001:
2024            sig = -0.001
2025   
2026    return Dist,sig
2027
2028def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
2029    '''default doc string
2030   
2031    :param type name: description
2032   
2033    :returns: type name: description
2034   
2035    '''
2036
2037    def calcAngle(Atoms,SyOps,Amat):
2038        XYZ = []
2039        for i,atom in enumerate(Atoms):
2040            Inv,M,T,C,U = SyOps[i]
2041            XYZ.append(np.array(atom[1:4]))
2042            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2043            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2044        V1 = XYZ[1]-XYZ[0]
2045        V1 /= np.sqrt(np.sum(V1**2))
2046        V2 = XYZ[1]-XYZ[2]
2047        V2 /= np.sqrt(np.sum(V2**2))
2048        V3 = V2-V1
2049        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2050        return acosd(cang)
2051
2052    Inv = []
2053    SyOps = []
2054    names = []
2055    for i,atom in enumerate(Oatoms):
2056        names += atom[-1]
2057        Op,unit = Atoms[i][-1]
2058        inv = Op/abs(Op)
2059        m,t = SGData['SGOps'][abs(Op)%100-1]
2060        c = SGData['SGCen'][abs(Op)/100]
2061        SyOps.append([inv,m,t,c,unit])
2062    Angle = calcAngle(Oatoms,SyOps,Amat)
2063   
2064    sig = -0.01
2065    if 'covMatrix' in covData:
2066        parmNames = []
2067        dx = .00001
2068        dadx = np.zeros(9)
2069        for i in range(9):
2070            ia = i/3
2071            ix = i%3
2072            Oatoms[ia][ix+1] += dx
2073            a0 = calcAngle(Oatoms,SyOps,Amat)
2074            Oatoms[ia][ix+1] -= 2*dx
2075            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2076        covMatrix = covData['covMatrix']
2077        varyList = covData['varyList']
2078        AngVcov = getVCov(names,varyList,covMatrix)
2079        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2080        if sig < 0.01:
2081            sig = -0.01
2082   
2083    return Angle,sig
2084
2085def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
2086    '''default doc string
2087   
2088    :param type name: description
2089   
2090    :returns: type name: description
2091   
2092    '''
2093
2094    def calcTorsion(Atoms,SyOps,Amat):
2095       
2096        XYZ = []
2097        for i,atom in enumerate(Atoms):
2098            Inv,M,T,C,U = SyOps[i]
2099            XYZ.append(np.array(atom[1:4]))
2100            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2101            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2102        V1 = XYZ[1]-XYZ[0]
2103        V2 = XYZ[2]-XYZ[1]
2104        V3 = XYZ[3]-XYZ[2]
2105        V1 /= np.sqrt(np.sum(V1**2))
2106        V2 /= np.sqrt(np.sum(V2**2))
2107        V3 /= np.sqrt(np.sum(V3**2))
2108        M = np.array([V1,V2,V3])
2109        D = nl.det(M)
2110        Ang = 1.0
2111        P12 = np.dot(V1,V2)
2112        P13 = np.dot(V1,V3)
2113        P23 = np.dot(V2,V3)
2114        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2115        return Tors
2116           
2117    Inv = []
2118    SyOps = []
2119    names = []
2120    for i,atom in enumerate(Oatoms):
2121        names += atom[-1]
2122        Op,unit = Atoms[i][-1]
2123        inv = Op/abs(Op)
2124        m,t = SGData['SGOps'][abs(Op)%100-1]
2125        c = SGData['SGCen'][abs(Op)/100]
2126        SyOps.append([inv,m,t,c,unit])
2127    Tors = calcTorsion(Oatoms,SyOps,Amat)
2128   
2129    sig = -0.01
2130    if 'covMatrix' in covData:
2131        parmNames = []
2132        dx = .00001
2133        dadx = np.zeros(12)
2134        for i in range(12):
2135            ia = i/3
2136            ix = i%3
2137            Oatoms[ia][ix+1] -= dx
2138            a0 = calcTorsion(Oatoms,SyOps,Amat)
2139            Oatoms[ia][ix+1] += 2*dx
2140            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2141            Oatoms[ia][ix+1] -= dx           
2142        covMatrix = covData['covMatrix']
2143        varyList = covData['varyList']
2144        TorVcov = getVCov(names,varyList,covMatrix)
2145        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
2146        if sig < 0.01:
2147            sig = -0.01
2148   
2149    return Tors,sig
2150       
2151def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
2152    '''default doc string
2153   
2154    :param type name: description
2155   
2156    :returns: type name: description
2157   
2158    '''
2159
2160    def calcDist(Atoms,SyOps,Amat):
2161        XYZ = []
2162        for i,atom in enumerate(Atoms):
2163            Inv,M,T,C,U = SyOps[i]
2164            XYZ.append(np.array(atom[1:4]))
2165            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2166            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2167        V1 = XYZ[1]-XYZ[0]
2168        return np.sqrt(np.sum(V1**2))
2169       
2170    def calcAngle(Atoms,SyOps,Amat):
2171        XYZ = []
2172        for i,atom in enumerate(Atoms):
2173            Inv,M,T,C,U = SyOps[i]
2174            XYZ.append(np.array(atom[1:4]))
2175            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2176            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2177        V1 = XYZ[1]-XYZ[0]
2178        V1 /= np.sqrt(np.sum(V1**2))
2179        V2 = XYZ[1]-XYZ[2]
2180        V2 /= np.sqrt(np.sum(V2**2))
2181        V3 = V2-V1
2182        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2183        return acosd(cang)
2184
2185    def calcTorsion(Atoms,SyOps,Amat):
2186       
2187        XYZ = []
2188        for i,atom in enumerate(Atoms):
2189            Inv,M,T,C,U = SyOps[i]
2190            XYZ.append(np.array(atom[1:4]))
2191            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2192            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2193        V1 = XYZ[1]-XYZ[0]
2194        V2 = XYZ[2]-XYZ[1]
2195        V3 = XYZ[3]-XYZ[2]
2196        V1 /= np.sqrt(np.sum(V1**2))
2197        V2 /= np.sqrt(np.sum(V2**2))
2198        V3 /= np.sqrt(np.sum(V3**2))
2199        M = np.array([V1,V2,V3])
2200        D = nl.det(M)
2201        Ang = 1.0
2202        P12 = np.dot(V1,V2)
2203        P13 = np.dot(V1,V3)
2204        P23 = np.dot(V2,V3)
2205        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2206        return Tors
2207           
2208    Inv = []
2209    SyOps = []
2210    names = []
2211    for i,atom in enumerate(Oatoms):
2212        names += atom[-1]
2213        Op,unit = Atoms[i][-1]
2214        inv = Op/abs(Op)
2215        m,t = SGData['SGOps'][abs(Op)%100-1]
2216        c = SGData['SGCen'][abs(Op)/100]
2217        SyOps.append([inv,m,t,c,unit])
2218    M = len(Oatoms)
2219    if M == 2:
2220        Val = calcDist(Oatoms,SyOps,Amat)
2221    elif M == 3:
2222        Val = calcAngle(Oatoms,SyOps,Amat)
2223    else:
2224        Val = calcTorsion(Oatoms,SyOps,Amat)
2225   
2226    sigVals = [-0.001,-0.01,-0.01]
2227    sig = sigVals[M-3]
2228    if 'covMatrix' in covData:
2229        parmNames = []
2230        dx = .00001
2231        N = M*3
2232        dadx = np.zeros(N)
2233        for i in range(N):
2234            ia = i/3
2235            ix = i%3
2236            Oatoms[ia][ix+1] += dx
2237            if M == 2:
2238                a0 = calcDist(Oatoms,SyOps,Amat)
2239            elif M == 3:
2240                a0 = calcAngle(Oatoms,SyOps,Amat)
2241            else:
2242                a0 = calcTorsion(Oatoms,SyOps,Amat)
2243            Oatoms[ia][ix+1] -= 2*dx
2244            if M == 2:
2245                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2246            elif M == 3:
2247                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2248            else:
2249                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2250        covMatrix = covData['covMatrix']
2251        varyList = covData['varyList']
2252        Vcov = getVCov(names,varyList,covMatrix)
2253        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
2254        if sig < sigVals[M-3]:
2255            sig = sigVals[M-3]
2256   
2257    return Val,sig
2258       
2259def ValEsd(value,esd=0,nTZ=False):
2260    '''Format a floating point number with a given level of precision or
2261    with in crystallographic format with a "esd", as value(esd). If esd is
2262    negative the number is formatted with the level of significant figures
2263    appropriate if abs(esd) were the esd, but the esd is not included.
2264    if the esd is zero, approximately 6 significant figures are printed.
2265    nTZ=True causes "extra" zeros to be removed after the decimal place.
2266    for example:
2267
2268      * "1.235(3)" for value=1.2346 & esd=0.003
2269      * "1.235(3)e4" for value=12346. & esd=30
2270      * "1.235(3)e6" for value=0.12346e7 & esd=3000
2271      * "1.235" for value=1.2346 & esd=-0.003
2272      * "1.240" for value=1.2395 & esd=-0.003
2273      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
2274      * "1.23460" for value=1.2346 & esd=0.0
2275
2276    :param float value: number to be formatted
2277    :param float esd: uncertainty or if esd < 0, specifies level of
2278      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
2279    :param bool nTZ: True to remove trailing zeros (default is False)
2280    :returns: value(esd) or value as a string
2281
2282    '''
2283    # Note: this routine is Python 3 compatible -- I think
2284    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
2285    if math.isnan(value): # invalid value, bail out
2286        return '?'
2287    if math.isnan(esd): # invalid esd, treat as zero
2288        esd = 0
2289        esdoff = 5
2290    elif esd != 0:
2291        # transform the esd to a one or two digit integer
2292        l = math.log10(abs(esd)) % 1.
2293        if l < math.log10(cutoff): l+= 1.       
2294        intesd = int(round(10**l)) # esd as integer
2295        # determine the number of digits offset for the esd
2296        esdoff = int(round(math.log10(intesd*1./abs(esd))))
2297    else:
2298        esdoff = 5
2299    valoff = 0
2300    if abs(value) < abs(esdoff): # value is effectively zero
2301        pass
2302    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
2303        # where the digit offset is to the left of the decimal place or where too many
2304        # digits are needed
2305        if abs(value) > 1:
2306            valoff = int(math.log10(abs(value)))
2307        elif abs(value) > 0:
2308            valoff = int(math.log10(abs(value))-0.9999999)
2309        else:
2310            valoff = 0
2311    if esd != 0:
2312        if valoff+esdoff < 0:
2313            valoff = esdoff = 0
2314        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
2315    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
2316        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
2317    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
2318        if abs(value) > 0:           
2319            extra = -math.log10(abs(value))
2320        else:
2321            extra = 0
2322        if extra > 0: extra += 1
2323        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
2324    if esd > 0:
2325        out += ("({:d})").format(intesd)  # add the esd
2326    elif nTZ and '.' in out:
2327        out = out.rstrip('0')  # strip zeros to right of decimal
2328        out = out.rstrip('.')  # and decimal place when not needed
2329    if valoff != 0:
2330        out += ("e{:d}").format(valoff) # add an exponent, when needed
2331    return out
2332   
2333################################################################################
2334##### Texture fitting stuff
2335################################################################################
2336
2337def FitTexture(General,Gangls,refData,keyList,pgbar):
2338    import pytexture as ptx
2339    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2340   
2341    def printSpHarm(textureData,SHtextureSig):
2342        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
2343        names = ['omega','chi','phi']
2344        namstr = '  names :'
2345        ptstr =  '  values:'
2346        sigstr = '  esds  :'
2347        for name in names:
2348            namstr += '%12s'%('Sample '+name)
2349            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2350            if 'Sample '+name in SHtextureSig:
2351                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2352            else:
2353                sigstr += 12*' '
2354        print namstr
2355        print ptstr
2356        print sigstr
2357        print '\n Texture coefficients:'
2358        SHcoeff = textureData['SH Coeff'][1]
2359        SHkeys = SHcoeff.keys()
2360        nCoeff = len(SHcoeff)
2361        nBlock = nCoeff/10+1
2362        iBeg = 0
2363        iFin = min(iBeg+10,nCoeff)
2364        for block in range(nBlock):
2365            namstr = '  names :'
2366            ptstr =  '  values:'
2367            sigstr = '  esds  :'
2368            for name in SHkeys[iBeg:iFin]:
2369                if 'C' in name:
2370                    namstr += '%12s'%(name)
2371                    ptstr += '%12.3f'%(SHcoeff[name])
2372                    if name in SHtextureSig:
2373                        sigstr += '%12.3f'%(SHtextureSig[name])
2374                    else:
2375                        sigstr += 12*' '
2376            print namstr
2377            print ptstr
2378            print sigstr
2379            iBeg += 10
2380            iFin = min(iBeg+10,nCoeff)
2381           
2382    def Dict2Values(parmdict, varylist):
2383        '''Use before call to leastsq to setup list of values for the parameters
2384        in parmdict, as selected by key in varylist'''
2385        return [parmdict[key] for key in varylist] 
2386       
2387    def Values2Dict(parmdict, varylist, values):
2388        ''' Use after call to leastsq to update the parameter dictionary with
2389        values corresponding to keys in varylist'''
2390        parmdict.update(zip(varylist,values))
2391       
2392    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2393        parmDict.update(zip(varyList,values))
2394        Mat = np.empty(0)
2395        sumObs = 0
2396        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2397        for hist in Gangls.keys():
2398            Refs = refData[hist]
2399            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2400            wt = 1./np.sqrt(np.max(Refs[:,4],.25))
2401#            wt = 1./np.max(Refs[:,4],.25)
2402            sumObs += np.sum(wt*Refs[:,5])
2403            Refs[:,6] = 1.
2404            H = Refs[:,:3]
2405            phi,beta = G2lat.CrsAng(H,cell,SGData)
2406            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2407            for item in parmDict:
2408                if 'C' in item:
2409                    L,M,N = eval(item.strip('C'))
2410                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2411                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2412                    Lnorm = G2lat.Lnorm(L)
2413                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2414            mat = wt*(Refs[:,5]-Refs[:,6])
2415            Mat = np.concatenate((Mat,mat))
2416        sumD = np.sum(np.abs(Mat))
2417        R = min(100.,100.*sumD/sumObs)
2418        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2419        print ' Residual: %.3f%%'%(R)
2420        return Mat
2421       
2422    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2423        Mat = np.empty(0)
2424        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2425        for hist in Gangls.keys():
2426            mat = np.zeros((len(varyList),len(refData[hist])))
2427            Refs = refData[hist]
2428            H = Refs[:,:3]
2429            wt = 1./np.sqrt(np.max(Refs[:,4],.25))
2430#            wt = 1./np.max(Refs[:,4],.25)
2431            phi,beta = G2lat.CrsAng(H,cell,SGData)
2432            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2433            for j,item in enumerate(varyList):
2434                if 'C' in item:
2435                    L,M,N = eval(item.strip('C'))
2436                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2437                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2438                    Lnorm = G2lat.Lnorm(L)
2439                    mat[j] = -wt*Lnorm*Kcl*Ksl
2440                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2441                        try:
2442                            l = varyList.index(itema)
2443                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2444                        except ValueError:
2445                            pass
2446            if len(Mat):
2447                Mat = np.concatenate((Mat,mat.T))
2448            else:
2449                Mat = mat.T
2450        print 'deriv'
2451        return Mat
2452
2453    print ' Fit texture for '+General['Name']
2454    SGData = General['SGData']
2455    cell = General['Cell'][1:7]
2456    Texture = General['SH Texture']
2457    if not Texture['Order']:
2458        return 'No spherical harmonics coefficients'
2459    varyList = []
2460    parmDict = copy.copy(Texture['SH Coeff'][1])
2461    for item in ['Sample omega','Sample chi','Sample phi']:
2462        parmDict[item] = Texture[item][1]
2463        if Texture[item][0]:
2464            varyList.append(item)
2465    if Texture['SH Coeff'][0]:
2466        varyList += Texture['SH Coeff'][1].keys()
2467    while True:
2468        begin = time.time()
2469        values =  np.array(Dict2Values(parmDict, varyList))
2470        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2471            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2472        ncyc = int(result[2]['nfev']/2)
2473        if ncyc:
2474            runtime = time.time()-begin   
2475            chisq = np.sum(result[2]['fvec']**2)
2476            Values2Dict(parmDict, varyList, result[0])
2477            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2478            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
2479            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2480            try:
2481                sig = np.sqrt(np.diag(result[1])*GOF)
2482                if np.any(np.isnan(sig)):
2483                    print '*** Least squares aborted - some invalid esds possible ***'
2484                break                   #refinement succeeded - finish up!
2485            except ValueError:          #result[1] is None on singular matrix
2486                print '**** Refinement failed - singular matrix ****'
2487                return None
2488        else:
2489            break
2490   
2491    if ncyc:
2492        for parm in parmDict:
2493            if 'C' in parm:
2494                Texture['SH Coeff'][1][parm] = parmDict[parm]
2495            else:
2496                Texture[parm][1] = parmDict[parm] 
2497        sigDict = dict(zip(varyList,sig))
2498        printSpHarm(Texture,sigDict)
2499       
2500    return None
2501   
2502################################################################################
2503##### Fourier & charge flip stuff
2504################################################################################
2505
2506def adjHKLmax(SGData,Hmax):
2507    '''default doc string
2508   
2509    :param type name: description
2510   
2511    :returns: type name: description
2512   
2513    '''
2514    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2515        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
2516        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
2517        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2518    else:
2519        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
2520        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
2521        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2522
2523def OmitMap(data,reflDict,pgbar=None):
2524    '''default doc string
2525   
2526    :param type name: description
2527   
2528    :returns: type name: description
2529   
2530    '''
2531    generalData = data['General']
2532    if not generalData['Map']['MapType']:
2533        print '**** ERROR - Fourier map not defined'
2534        return
2535    mapData = generalData['Map']
2536    dmin = mapData['Resolution']
2537    SGData = generalData['SGData']
2538    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2539    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2540    cell = generalData['Cell'][1:8]       
2541    A = G2lat.cell2A(cell[:6])
2542    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2543    adjHKLmax(SGData,Hmax)
2544    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2545    time0 = time.time()
2546    for iref,ref in enumerate(reflDict['RefList']):
2547        if ref[4] >= dmin:
2548            Fosq,Fcsq,ph = ref[8:11]
2549            Uniq = np.inner(ref[:3],SGMT)
2550            Phi = np.inner(ref[:3],SGT)
2551            for i,hkl in enumerate(Uniq):        #uses uniq
2552                hkl = np.asarray(hkl,dtype='i')
2553                dp = 360.*Phi[i]                #and phi
2554                a = cosd(ph+dp)
2555                b = sind(ph+dp)
2556                phasep = complex(a,b)
2557                phasem = complex(a,-b)
2558                Fo = np.sqrt(Fosq)
2559                if '2Fo-Fc' in mapData['MapType']:
2560                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
2561                else:
2562                    F = np.sqrt(Fosq)
2563                h,k,l = hkl+Hmax
2564                Fhkl[h,k,l] = F*phasep
2565                h,k,l = -hkl+Hmax
2566                Fhkl[h,k,l] = F*phasem
2567    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2568    M = np.mgrid[0:4,0:4,0:4]
2569    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
2570    iBeg = blkIds*rho0.shape/4
2571    iFin = (blkIds+1)*rho0.shape/4
2572    rho_omit = np.zeros_like(rho0)
2573    nBlk = 0
2574    for iB,iF in zip(iBeg,iFin):
2575        rho1 = np.copy(rho0)
2576        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
2577        Fnew = fft.ifftshift(fft.ifftn(rho1))
2578        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
2579        phase = Fnew/np.absolute(Fnew)
2580        OFhkl = np.absolute(Fhkl)*phase
2581        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
2582        rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]])
2583        nBlk += 1
2584        pgbar.Update(nBlk)
2585    mapData['rho'] = np.real(rho_omit)/cell[6]
2586    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2587    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2588    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2589    return mapData
2590   
2591def FourierMap(data,reflDict):
2592    '''default doc string
2593   
2594    :param type name: description
2595   
2596    :returns: type name: description
2597   
2598    '''
2599    generalData = data['General']
2600    mapData = generalData['Map']
2601    dmin = mapData['Resolution']
2602    SGData = generalData['SGData']
2603    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2604    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2605    cell = generalData['Cell'][1:8]       
2606    A = G2lat.cell2A(cell[:6])
2607    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2608    adjHKLmax(SGData,Hmax)
2609    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2610#    Fhkl[0,0,0] = generalData['F000X']
2611    time0 = time.time()
2612    for iref,ref in enumerate(reflDict['RefList']):
2613        if ref[4] > dmin:
2614            Fosq,Fcsq,ph = ref[8:11]
2615            Uniq = np.inner(ref[:3],SGMT)
2616            Phi = np.inner(ref[:3],SGT)
2617            for i,hkl in enumerate(Uniq):        #uses uniq
2618                hkl = np.asarray(hkl,dtype='i')
2619                dp = 360.*Phi[i]                #and phi
2620                a = cosd(ph+dp)
2621                b = sind(ph+dp)
2622                phasep = complex(a,b)
2623                phasem = complex(a,-b)
2624                if 'Fobs' in mapData['MapType']:
2625                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
2626                    h,k,l = hkl+Hmax
2627                    Fhkl[h,k,l] = F*phasep
2628                    h,k,l = -hkl+Hmax
2629                    Fhkl[h,k,l] = F*phasem
2630                elif 'Fcalc' in mapData['MapType']:
2631                    F = np.sqrt(Fcsq)
2632                    h,k,l = hkl+Hmax
2633                    Fhkl[h,k,l] = F*phasep
2634                    h,k,l = -hkl+Hmax
2635                    Fhkl[h,k,l] = F*phasem
2636                elif 'delt-F' in mapData['MapType']:
2637                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2638                    h,k,l = hkl+Hmax
2639                    Fhkl[h,k,l] = dF*phasep
2640                    h,k,l = -hkl+Hmax
2641                    Fhkl[h,k,l] = dF*phasem
2642                elif '2*Fo-Fc' in mapData['MapType']:
2643                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2644                    h,k,l = hkl+Hmax
2645                    Fhkl[h,k,l] = F*phasep
2646                    h,k,l = -hkl+Hmax
2647                    Fhkl[h,k,l] = F*phasem
2648                elif 'Patterson' in mapData['MapType']:
2649                    h,k,l = hkl+Hmax
2650                    Fhkl[h,k,l] = complex(Fosq,0.)
2651                    h,k,l = -hkl+Hmax
2652                    Fhkl[h,k,l] = complex(Fosq,0.)
2653    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2654    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2655    mapData['Type'] = reflDict['Type']
2656    mapData['rho'] = np.real(rho)
2657    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2658    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2659   
2660def Fourier4DMap(data,reflDict):
2661    '''default doc string
2662   
2663    :param type name: description
2664   
2665    :returns: type name: description
2666   
2667    '''
2668    generalData = data['General']
2669    map4DData = generalData['4DmapData']
2670    mapData = generalData['Map']
2671    dmin = mapData['Resolution']
2672    SGData = generalData['SGData']
2673    SSGData = generalData['SSGData']
2674    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2675    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2676    cell = generalData['Cell'][1:8]       
2677    A = G2lat.cell2A(cell[:6])
2678    maxM = 4
2679    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
2680    adjHKLmax(SGData,Hmax)
2681    Hmax = np.asarray(Hmax,dtype='i')+1
2682    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2683    time0 = time.time()
2684    for iref,ref in enumerate(reflDict['RefList']):
2685        if ref[5] > dmin:
2686            Fosq,Fcsq,ph = ref[9:12]
2687            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
2688            Uniq = np.inner(ref[:4],SSGMT)
2689            Phi = np.inner(ref[:4],SSGT)
2690            for i,hkl in enumerate(Uniq):        #uses uniq
2691                hkl = np.asarray(hkl,dtype='i')
2692                dp = 360.*Phi[i]                #and phi
2693                a = cosd(ph+dp)
2694                b = sind(ph+dp)
2695                phasep = complex(a,b)
2696                phasem = complex(a,-b)
2697                if 'Fobs' in mapData['MapType']:
2698                    F = np.sqrt(Fosq)
2699                    h,k,l,m = hkl+Hmax
2700                    Fhkl[h,k,l,m] = F*phasep
2701                    h,k,l,m = -hkl+Hmax
2702                    Fhkl[h,k,l,m] = F*phasem
2703                elif 'Fcalc' in mapData['MapType']:
2704                    F = np.sqrt(Fcsq)
2705                    h,k,l,m = hkl+Hmax
2706                    Fhkl[h,k,l,m] = F*phasep
2707                    h,k,l,m = -hkl+Hmax
2708                    Fhkl[h,k,l,m] = F*phasem                   
2709                elif 'delt-F' in mapData['MapType']:
2710                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
2711                    h,k,l,m = hkl+Hmax
2712                    Fhkl[h,k,l,m] = dF*phasep
2713                    h,k,l,m = -hkl+Hmax
2714                    Fhkl[h,k,l,m] = dF*phasem
2715    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
2716    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
2717    map4DData['rho'] = np.real(SSrho)
2718    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2719    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2720    map4DData['Type'] = reflDict['Type']
2721    mapData['Type'] = reflDict['Type']
2722    mapData['rho'] = np.real(rho)
2723    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2724    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2725    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2726
2727# map printing for testing purposes
2728def printRho(SGLaue,rho,rhoMax):                         
2729    '''default doc string
2730   
2731    :param type name: description
2732   
2733    :returns: type name: description
2734   
2735    '''
2736    dim = len(rho.shape)
2737    if dim == 2:
2738        ix,jy = rho.shape
2739        for j in range(jy):
2740            line = ''
2741            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2742                line += (jy-j)*'  '
2743            for i in range(ix):
2744                r = int(100*rho[i,j]/rhoMax)
2745                line += '%4d'%(r)
2746            print line+'\n'
2747    else:
2748        ix,jy,kz = rho.shape
2749        for k in range(kz):
2750            print 'k = ',k
2751            for j in range(jy):
2752                line = ''
2753                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2754                    line += (jy-j)*'  '
2755                for i in range(ix):
2756                    r = int(100*rho[i,j,k]/rhoMax)
2757                    line += '%4d'%(r)
2758                print line+'\n'
2759## keep this
2760               
2761def findOffset(SGData,A,Fhkl):   
2762    '''default doc string
2763   
2764    :param type name: description
2765   
2766    :returns: type name: description
2767   
2768    '''
2769    if SGData['SpGrp'] == 'P 1':
2770        return [0,0,0]   
2771    hklShape = Fhkl.shape
2772    hklHalf = np.array(hklShape)/2
2773    sortHKL = np.argsort(Fhkl.flatten())
2774    Fdict = {}
2775    for hkl in sortHKL:
2776        HKL = np.unravel_index(hkl,hklShape)
2777        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
2778        if F == 0.:
2779            break
2780        Fdict['%.6f'%(np.absolute(F))] = hkl
2781    Flist = np.flipud(np.sort(Fdict.keys()))
2782    F = str(1.e6)
2783    i = 0
2784    DH = []
2785    Dphi = []
2786    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2787    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2788    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2789    for F in Flist:
2790        hkl = np.unravel_index(Fdict[F],hklShape)
2791        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
2792            continue
2793        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
2794        Uniq = np.array(Uniq,dtype='i')
2795        Phi = np.array(Phi)
2796        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
2797        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2798        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
2799        ang0 = np.angle(Fh0,deg=True)/360.
2800        for H,phi in zip(Uniq,Phi)[1:]:
2801            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
2802            dH = H-hkl
2803            dang = ang-ang0
2804            DH.append(dH)
2805            Dphi.append((dang+.5) % 1.0)
2806        if i > 20 or len(DH) > 30:
2807            break
2808        i += 1
2809    DH = np.array(DH)
2810    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2811    Dphi = np.array(Dphi)
2812    steps = np.array(hklShape)
2813    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
2814    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
2815    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
2816    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2817    hist,bins = np.histogram(Mmap,bins=1000)
2818#    for i,item in enumerate(hist[:10]):
2819#        print item,bins[i]
2820    chisq = np.min(Mmap)
2821    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2822    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
2823#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2824    return DX
2825   
2826def ChargeFlip(data,reflDict,pgbar):
2827    '''default doc string
2828   
2829    :param type name: description
2830   
2831    :returns: type name: description
2832   
2833    '''
2834    generalData = data['General']
2835    mapData = generalData['Map']
2836    flipData = generalData['Flip']
2837    FFtable = {}
2838    if 'None' not in flipData['Norm element']:
2839        normElem = flipData['Norm element'].upper()
2840        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2841        for ff in FFs:
2842            if ff['Symbol'] == normElem:
2843                FFtable.update(ff)
2844    dmin = flipData['Resolution']
2845    SGData = generalData['SGData']
2846    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2847    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2848    cell = generalData['Cell'][1:8]       
2849    A = G2lat.cell2A(cell[:6])
2850    Vol = cell[6]
2851    im = 0
2852    if generalData['Type'] in ['modulated','magnetic',]:
2853        im = 1
2854    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2855    adjHKLmax(SGData,Hmax)
2856    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2857    time0 = time.time()
2858    for iref,ref in enumerate(reflDict['RefList']):
2859        dsp = ref[4+im]
2860        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
2861            continue
2862        if dsp > dmin:
2863            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2864            if FFtable:
2865                SQ = 0.25/dsp**2
2866                ff *= G2el.ScatFac(FFtable,SQ)[0]
2867            if ref[8+im] > 0.:         #use only +ve Fobs**2
2868                E = np.sqrt(ref[8+im])/ff
2869            else:
2870                E = 0.
2871            ph = ref[10]
2872            ph = rn.uniform(0.,360.)
2873            Uniq = np.inner(ref[:3],SGMT)
2874            Phi = np.inner(ref[:3],SGT)
2875            for i,hkl in enumerate(Uniq):        #uses uniq
2876                hkl = np.asarray(hkl,dtype='i')
2877                dp = 360.*Phi[i]                #and phi
2878                a = cosd(ph+dp)
2879                b = sind(ph+dp)
2880                phasep = complex(a,b)
2881                phasem = complex(a,-b)
2882                h,k,l = hkl+Hmax
2883                Ehkl[h,k,l] = E*phasep
2884                h,k,l = -hkl+Hmax       #Friedel pair refl.
2885                Ehkl[h,k,l] = E*phasem
2886#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2887    testHKL = np.array(flipData['testHKL'])+Hmax
2888    CEhkl = copy.copy(Ehkl)
2889    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2890    Emask = ma.getmask(MEhkl)
2891    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2892    Ncyc = 0
2893    old = np.seterr(all='raise')
2894    twophases = []
2895    while True:       
2896        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2897        CEsig = np.std(CErho)
2898        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2899        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2900        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2901        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2902        phase = CFhkl/np.absolute(CFhkl)
2903        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
2904        CEhkl = np.absolute(Ehkl)*phase
2905        Ncyc += 1
2906        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2907        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2908        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2909        if Rcf < 5.:
2910            break
2911        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2912        if not GoOn or Ncyc > 10000:
2913            break
2914    np.seterr(**old)
2915    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2916    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
2917    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2918    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2919       
2920    mapData['Rcf'] = Rcf
2921    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2922    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2923    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2924    mapData['Type'] = reflDict['Type']
2925    return mapData,twophases
2926   
2927def findSSOffset(SGData,SSGData,A,Fhklm):   
2928    '''default doc string
2929   
2930    :param type name: description
2931   
2932    :returns: type name: description
2933   
2934    '''
2935    if SGData['SpGrp'] == 'P 1':
2936        return [0,0,0,0]   
2937    hklmShape = Fhklm.shape
2938    hklmHalf = np.array(hklmShape)/2
2939    sortHKLM = np.argsort(Fhklm.flatten())
2940    Fdict = {}
2941    for hklm in sortHKLM:
2942        HKLM = np.unravel_index(hklm,hklmShape)
2943        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
2944        if F == 0.:
2945            break
2946        Fdict['%.6f'%(np.absolute(F))] = hklm
2947    Flist = np.flipud(np.sort(Fdict.keys()))
2948    F = str(1.e6)
2949    i = 0
2950    DH = []
2951    Dphi = []
2952    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2953    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2954    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2955    for F in Flist:
2956        hklm = np.unravel_index(Fdict[F],hklmShape)
2957        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
2958            continue
2959        Uniq = np.inner(hklm-hklmHalf,SSGMT)
2960        Phi = np.inner(hklm-hklmHalf,SSGT)
2961        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
2962        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2963        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
2964        ang0 = np.angle(Fh0,deg=True)/360.
2965        for H,phi in zip(Uniq,Phi)[1:]:
2966            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
2967            dH = H-hklm
2968            dang = ang-ang0
2969            DH.append(dH)
2970            Dphi.append((dang+.5) % 1.0)
2971        if i > 20 or len(DH) > 30:
2972            break
2973        i += 1
2974    DH = np.array(DH)
2975    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2976    Dphi = np.array(Dphi)
2977    steps = np.array(hklmShape)
2978    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
2979    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
2980    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
2981    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2982    hist,bins = np.histogram(Mmap,bins=1000)
2983#    for i,item in enumerate(hist[:10]):
2984#        print item,bins[i]
2985    chisq = np.min(Mmap)
2986    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2987    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
2988#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2989    return DX
2990   
2991def SSChargeFlip(data,reflDict,pgbar):
2992    '''default doc string
2993   
2994    :param type name: description
2995   
2996    :returns: type name: description
2997   
2998    '''
2999    generalData = data['General']
3000    mapData = generalData['Map']
3001    map4DData = {}
3002    flipData = generalData['Flip']
3003    FFtable = {}
3004    if 'None' not in flipData['Norm element']:
3005        normElem = flipData['Norm element'].upper()
3006        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3007        for ff in FFs:
3008            if ff['Symbol'] == normElem:
3009                FFtable.update(ff)
3010    dmin = flipData['Resolution']
3011    SGData = generalData['SGData']
3012    SSGData = generalData['SSGData']
3013    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3014    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3015    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3016    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3017    cell = generalData['Cell'][1:8]       
3018    A = G2lat.cell2A(cell[:6])
3019    Vol = cell[6]
3020    maxM = 4
3021    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3022    adjHKLmax(SGData,Hmax)
3023    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3024    time0 = time.time()
3025    for iref,ref in enumerate(reflDict['RefList']):
3026        dsp = ref[5]
3027        if dsp > dmin:
3028            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3029            if FFtable:
3030                SQ = 0.25/dsp**2
3031                ff *= G2el.ScatFac(FFtable,SQ)[0]
3032            if ref[9] > 0.:         #use only +ve Fobs**2
3033                E = np.sqrt(ref[9])/ff
3034            else:
3035                E = 0.
3036            ph = ref[11]
3037            ph = rn.uniform(0.,360.)
3038            Uniq = np.inner(ref[:4],SSGMT)
3039            Phi = np.inner(ref[:4],SSGT)
3040            for i,hklm in enumerate(Uniq):        #uses uniq
3041                hklm = np.asarray(hklm,dtype='i')
3042                dp = 360.*Phi[i]                #and phi
3043                a = cosd(ph+dp)
3044                b = sind(ph+dp)
3045                phasep = complex(a,b)
3046                phasem = complex(a,-b)
3047                h,k,l,m = hklm+Hmax
3048                Ehkl[h,k,l,m] = E*phasep
3049                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3050                Ehkl[h,k,l,m] = E*phasem
3051#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3052    CEhkl = copy.copy(Ehkl)
3053    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3054    Emask = ma.getmask(MEhkl)
3055    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3056    Ncyc = 0
3057    old = np.seterr(all='raise')
3058    while True:       
3059        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3060        CEsig = np.std(CErho)
3061        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3062        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3063        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3064        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3065        phase = CFhkl/np.absolute(CFhkl)
3066        CEhkl = np.absolute(Ehkl)*phase
3067        Ncyc += 1
3068        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3069        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3070        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3071        if Rcf < 5.:
3072            break
3073        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3074        if not GoOn or Ncyc > 10000:
3075            break
3076    np.seterr(**old)
3077    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3078    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3079    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3080    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3081    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3082
3083    mapData['Rcf'] = Rcf
3084    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3085    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3086    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3087    mapData['Type'] = reflDict['Type']
3088
3089    map4DData['Rcf'] = Rcf
3090    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
3091    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3092    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3093    map4DData['Type'] = reflDict['Type']
3094    return mapData,map4DData
3095   
3096def getRho(xyz,mapData):
3097    ''' get scattering density at a point by 8-point interpolation
3098    param xyz: coordinate to be probed
3099    param: mapData: dict of map data
3100   
3101    :returns: density at xyz
3102    '''
3103    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3104    if not len(mapData):
3105        return 0.0
3106    rho = copy.copy(mapData['rho'])     #don't mess up original
3107    if not len(rho):
3108        return 0.0
3109    mapShape = np.array(rho.shape)
3110    mapStep = 1./mapShape
3111    X = np.array(xyz)%1.    #get into unit cell
3112    I = np.array(X*mapShape,dtype='int')
3113    D = X-I*mapStep         #position inside map cell
3114    D12 = D[0]*D[1]
3115    D13 = D[0]*D[2]
3116    D23 = D[1]*D[2]
3117    D123 = np.prod(D)
3118    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3119    R = Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \
3120        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3121        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3122        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3123    return R
3124       
3125def SearchMap(generalData,drawingData,Neg=False):
3126    '''Does a search of a density map for peaks meeting the criterion of peak
3127    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3128    mapData is data['General']['mapData']; the map is also in mapData.
3129
3130    :param generalData: the phase data structure; includes the map
3131    :param drawingData: the drawing data structure
3132    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3133
3134    :returns: (peaks,mags,dzeros) where
3135
3136        * peaks : ndarray
3137          x,y,z positions of the peaks found in the map
3138        * mags : ndarray
3139          the magnitudes of the peaks
3140        * dzeros : ndarray
3141          the distance of the peaks from  the unit cell origin
3142        * dcent : ndarray
3143          the distance of the peaks from  the unit cell center
3144
3145    '''       
3146    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3147   
3148    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3149   
3150#    def noDuplicate(xyz,peaks,Amat):
3151#        XYZ = np.inner(Amat,xyz)
3152#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3153#            print ' Peak',xyz,' <0.5A from another peak'
3154#            return False
3155#        return True
3156#                           
3157    def fixSpecialPos(xyz,SGData,Amat):
3158        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3159        X = []
3160        xyzs = [equiv[0] for equiv in equivs]
3161        for x in xyzs:
3162            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3163                X.append(x)
3164        if len(X) > 1:
3165            return np.average(X,axis=0)
3166        else:
3167            return xyz
3168       
3169    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3170        Mag,x0,y0,z0,sig = parms
3171        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3172#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3173        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3174       
3175    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3176        Mag,x0,y0,z0,sig = parms
3177        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3178        return M
3179       
3180    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3181        Mag,x0,y0,z0,sig = parms
3182        dMdv = np.zeros(([5,]+list(rX.shape)))
3183        delt = .01
3184        for i in range(5):
3185            parms[i] -= delt
3186            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3187            parms[i] += 2.*delt
3188            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3189            parms[i] -= delt
3190            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3191        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3192        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3193        dMdv = np.reshape(dMdv,(5,rX.size))
3194        Hess = np.inner(dMdv,dMdv)
3195       
3196        return Vec,Hess
3197       
3198    phaseName = generalData['Name']
3199    SGData = generalData['SGData']
3200    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3201    peaks = []
3202    mags = []
3203    dzeros = []
3204    dcent = []
3205    try:
3206        mapData = generalData['Map']
3207        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3208        if Neg:
3209            rho = -copy.copy(mapData['rho'])     #flip +/-
3210        else:
3211            rho = copy.copy(mapData['rho'])     #don't mess up original
3212        mapHalf = np.array(rho.shape)/2
3213        res = mapData['Resolution']
3214        incre = np.array(rho.shape,dtype=np.float)
3215        step = max(1.0,1./res)+1
3216        steps = np.array((3*[step,]),dtype='int32')
3217    except KeyError:
3218        print '**** ERROR - Fourier map not defined'
3219        return peaks,mags
3220    rhoMask = ma.array(rho,mask=(rho<contLevel))
3221    indices = (-1,0,1)
3222    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3223    for roll in rolls:
3224        if np.any(roll):
3225            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3226    indx = np.transpose(rhoMask.nonzero())
3227    peaks = indx/incre
3228    mags = rhoMask[rhoMask.nonzero()]
3229    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3230        rho = rollMap(rho,ind)
3231        rMM = mapHalf-steps
3232        rMP = mapHalf+steps+1
3233        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3234        peakInt = np.sum(rhoPeak)*res**3
3235        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3236        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3237        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3238            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3239        x1 = result[0]
3240        if not np.any(x1 < 0):
3241            mag = x1[0]
3242            peak = (np.array(x1[1:4])-ind)/incre
3243        peak = fixSpecialPos(peak,SGData,Amat)
3244        rho = rollMap(rho,-ind)
3245    cent = np.ones(3)*.5     
3246    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3247    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3248    if Neg:     #want negative magnitudes for negative peaks
3249        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3250    else:
3251        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3252   
3253def sortArray(data,pos,reverse=False):
3254    '''data is a list of items
3255    sort by pos in list; reverse if True
3256    '''
3257    T = []
3258    for i,M in enumerate(data):
3259        try:
3260            T.append((M[pos],i))
3261        except IndexError:
3262            return data
3263    D = dict(zip(T,data))
3264    T.sort()
3265    if reverse:
3266        T.reverse()
3267    X = []
3268    for key in T:
3269        X.append(D[key])
3270    return X
3271
3272def PeaksEquiv(data,Ind):
3273    '''Find the equivalent map peaks for those selected. Works on the
3274    contents of data['Map Peaks'].
3275
3276    :param data: the phase data structure
3277    :param list Ind: list of selected peak indices
3278    :returns: augmented list of peaks including those related by symmetry to the
3279      ones in Ind
3280
3281    '''       
3282    def Duplicate(xyz,peaks,Amat):
3283        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3284            return True
3285        return False
3286                           
3287    generalData = data['General']
3288    cell = generalData['Cell'][1:7]
3289    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3290    A = G2lat.cell2A(cell)
3291    SGData = generalData['SGData']
3292    mapPeaks = data['Map Peaks']
3293    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3294    Indx = {}
3295    for ind in Ind:
3296        xyz = np.array(mapPeaks[ind][1:4])
3297        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3298        for jnd,xyz in enumerate(XYZ):       
3299            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3300    Ind = []
3301    for ind in Indx:
3302        if Indx[ind]:
3303            Ind.append(ind)
3304    return Ind
3305               
3306def PeaksUnique(data,Ind):
3307    '''Finds the symmetry unique set of peaks from those selected. Works on the
3308    contents of data['Map Peaks'].
3309
3310    :param data: the phase data structure
3311    :param list Ind: list of selected peak indices
3312    :returns: the list of symmetry unique peaks from among those given in Ind
3313
3314    '''       
3315#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3316
3317    def noDuplicate(xyz,peaks,Amat):
3318        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3319            return False
3320        return True
3321                           
3322    generalData = data['General']
3323    cell = generalData['Cell'][1:7]
3324    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3325    A = G2lat.cell2A(cell)
3326    SGData = generalData['SGData']
3327    mapPeaks = data['Map Peaks']
3328    Indx = {}
3329    XYZ = {}
3330    for ind in Ind:
3331        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3332        Indx[ind] = True
3333    for ind in Ind:
3334        if Indx[ind]:
3335            xyz = XYZ[ind]
3336            for jnd in Ind:
3337                if ind != jnd and Indx[jnd]:                       
3338                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3339                    xyzs = np.array([equiv[0] for equiv in Equiv])
3340                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3341    Ind = []
3342    for ind in Indx:
3343        if Indx[ind]:
3344            Ind.append(ind)
3345    return Ind
3346   
3347################################################################################
3348##### single peak fitting profile fxn stuff
3349################################################################################
3350
3351def getCWsig(ins,pos):
3352    '''get CW peak profile sigma^2
3353   
3354    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3355      as values only
3356    :param float pos: 2-theta of peak
3357    :returns: float getCWsig: peak sigma^2
3358   
3359    '''
3360    tp = tand(pos/2.0)
3361    return ins['U']*tp**2+ins['V']*tp+ins['W']
3362   
3363def getCWsigDeriv(pos):
3364    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3365   
3366    :param float pos: 2-theta of peak
3367   
3368    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3369   
3370    '''
3371    tp = tand(pos/2.0)
3372    return tp**2,tp,1.0
3373   
3374def getCWgam(ins,pos):
3375    '''get CW peak profile gamma
3376   
3377    :param dict ins: instrument parameters with at least 'X' & 'Y'
3378      as values only
3379    :param float pos: 2-theta of peak
3380    :returns: float getCWgam: peak gamma
3381   
3382    '''
3383    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
3384   
3385def getCWgamDeriv(pos):
3386    '''get derivatives of CW peak profile gamma wrt X & Y
3387   
3388    :param float pos: 2-theta of peak
3389   
3390    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3391   
3392    '''
3393    return 1./cosd(pos/2.0),tand(pos/2.0)
3394   
3395def getTOFsig(ins,dsp):
3396    '''get TOF peak profile sigma^2
3397   
3398    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3399      as values only
3400    :param float dsp: d-spacing of peak
3401   
3402    :returns: float getTOFsig: peak sigma^2
3403   
3404    '''
3405    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
3406   
3407def getTOFsigDeriv(dsp):
3408    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3409   
3410    :param float dsp: d-spacing of peak
3411   
3412    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3413   
3414    '''
3415    return 1.0,dsp**2,dsp**4,1./dsp**2
3416   
3417def getTOFgamma(ins,dsp):
3418    '''get TOF peak profile gamma
3419   
3420    :param dict ins: instrument parameters with at least 'X' & 'Y'
3421      as values only
3422    :param float dsp: d-spacing of peak
3423   
3424    :returns: float getTOFgamma: peak gamma
3425   
3426    '''
3427    return ins['X']*dsp+ins['Y']*dsp**2
3428   
3429def getTOFgammaDeriv(dsp):
3430    '''get derivatives of TOF peak profile gamma wrt X & Y
3431   
3432    :param float dsp: d-spacing of peak
3433   
3434    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3435   
3436    '''
3437    return dsp,dsp**2
3438   
3439def getTOFbeta(ins,dsp):
3440    '''get TOF peak profile beta
3441   
3442    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3443      as values only
3444    :param float dsp: d-spacing of peak
3445   
3446    :returns: float getTOFbeta: peak beat
3447   
3448    '''
3449    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3450   
3451def getTOFbetaDeriv(dsp):
3452    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3453   
3454    :param float dsp: d-spacing of peak
3455   
3456    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3457   
3458    '''
3459    return 1.0,1./dsp**4,1./dsp**2
3460   
3461def getTOFalpha(ins,dsp):
3462    '''get TOF peak profile alpha
3463   
3464    :param dict ins: instrument parameters with at least 'alpha'
3465      as values only
3466    :param float dsp: d-spacing of peak
3467   
3468    :returns: flaot getTOFalpha: peak alpha
3469   
3470    '''
3471    return ins['alpha']/dsp
3472   
3473def getTOFalphaDeriv(dsp):
3474    '''get derivatives of TOF peak profile beta wrt alpha
3475   
3476    :param float dsp: d-spacing of peak
3477   
3478    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3479   
3480    '''
3481    return 1./dsp
3482   
3483def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3484    '''set starting peak parameters for single peak fits from plot selection or auto selection
3485   
3486    :param dict Parms: instrument parameters dictionary
3487    :param dict Parms2: table lookup for TOF profile coefficients
3488    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3489    :param float mag: peak top magnitude from pick
3490    :param bool ifQ: True if pos in Q
3491    :param bool useFit: True if use fitted CW Parms values (not defaults)
3492   
3493    :returns: list XY: peak list entry:
3494        for CW: [pos,0,mag,1,sig,0,gam,0]
3495        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3496        NB: mag refinement set by default, all others off
3497   
3498    '''
3499    ind = 0
3500    if useFit:
3501        ind = 1
3502    ins = {}
3503    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
3504        for x in ['U','V','W','X','Y']:
3505            ins[x] = Parms[x][ind]
3506        if ifQ:                              #qplot - convert back to 2-theta
3507            pos = 2.0*asind(pos*wave/(4*math.pi))
3508        sig = getCWsig(ins,pos)
3509        gam = getCWgam(ins,pos)           
3510        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
3511    else:
3512        if ifQ:
3513            dsp = 2.*np.pi/pos
3514            pos = Parms['difC']*dsp
3515        else:
3516            dsp = pos/Parms['difC'][1]
3517        if 'Pdabc' in Parms2:
3518            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
3519                ins[x] = Parms[x][ind]
3520            Pdabc = Parms2['Pdabc'].T
3521            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
3522            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
3523        else:
3524            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
3525                ins[x] = Parms[x][ind]
3526            alp = getTOFalpha(ins,dsp)
3527            bet = getTOFbeta(ins,dsp)
3528        sig = getTOFsig(ins,dsp)
3529        gam = getTOFgamma(ins,dsp)
3530        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3531    return XY
3532   
3533################################################################################
3534##### MC/SA stuff
3535################################################################################
3536
3537#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
3538# Original Author: Travis Oliphant 2002
3539# Bug-fixes in 2006 by Tim Leslie
3540
3541
3542import numpy
3543from numpy import asarray, tan, exp, ones, squeeze, sign, \
3544     all, log, sqrt, pi, shape, array, minimum, where
3545from numpy import random
3546
3547#__all__ = ['anneal']
3548
3549_double_min = numpy.finfo(float).min
3550_double_max = numpy.finfo(float).max
3551class base_schedule(object):
3552    def __init__(self):
3553        self.dwell = 20
3554        self.learn_rate = 0.5
3555        self.lower = -10
3556        self.upper = 10
3557        self.Ninit = 50
3558        self.accepted = 0
3559        self.tests = 0
3560        self.feval = 0
3561        self.k = 0
3562        self.T = None
3563
3564    def init(self, **options):
3565        self.__dict__.update(options)
3566        self.lower = asarray(self.lower)
3567        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3568        self.upper = asarray(self.upper)
3569        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3570        self.k = 0
3571        self.accepted = 0
3572        self.feval = 0
3573        self.tests = 0
3574
3575    def getstart_temp(self, best_state):
3576        """ Find a matching starting temperature and starting parameters vector
3577        i.e. find x0 such that func(x0) = T0.
3578
3579        Parameters
3580        ----------
3581        best_state : _state
3582            A _state object to store the function value and x0 found.
3583
3584        returns
3585        -------
3586        x0 : array
3587            The starting parameters vector.
3588        """
3589
3590        assert(not self.dims is None)
3591        lrange = self.lower
3592        urange = self.upper
3593        fmax = _double_min
3594        fmin = _double_max
3595        for _ in range(self.Ninit):
3596            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3597            fval = self.func(x0, *self.args)
3598            self.feval += 1
3599            if fval > fmax:
3600                fmax = fval
3601            if fval < fmin:
3602                fmin = fval
3603                best_state.cost = fval
3604                best_state.x = array(x0)
3605
3606        self.T0 = (fmax-fmin)*1.5
3607        return best_state.x
3608       
3609    def set_range(self,x0,frac):
3610        delrange = frac*(self.upper-self.lower)
3611        self.upper = x0+delrange
3612        self.lower = x0-delrange
3613
3614    def accept_test(self, dE):
3615        T = self.T
3616        self.tests += 1
3617        if dE < 0:
3618            self.accepted += 1
3619            return 1
3620        p = exp(-dE*1.0/self.boltzmann/T)
3621        if (p > random.uniform(0.0, 1.0)):
3622            self.accepted += 1
3623            return 1
3624        return 0
3625
3626    def update_guess(self, x0):
3627        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3628
3629    def update_temp(self, x0):
3630        pass
3631
3632
3633#  A schedule due to Lester Ingber modified to use bounds - OK
3634class fast_sa(base_schedule):
3635    def init(self, **options):
3636        self.__dict__.update(options)
3637        if self.m is None:
3638            self.m = 1.0
3639        if self.n is None:
3640            self.n = 1.0
3641        self.c = self.m * exp(-self.n * self.quench)
3642
3643    def update_guess(self, x0):
3644        x0 = asarray(x0)
3645        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
3646        T = self.T
3647        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
3648        xnew = xc*(self.upper - self.lower)+self.lower
3649        return xnew
3650#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
3651#        xc = y*(self.upper - self.lower)
3652#        xnew = x0 + xc
3653#        return xnew
3654
3655    def update_temp(self):
3656        self.T = self.T0*exp(-self.c * self.k**(self.quench))
3657        self.k += 1
3658        return
3659
3660class cauchy_sa(base_schedule):     #modified to use bounds - not good
3661    def update_guess(self, x0):
3662        x0 = asarray(x0)
3663        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
3664        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
3665        xnew = xc*(self.upper - self.lower)+self.lower
3666        return xnew
3667#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
3668#        xc = self.learn_rate * self.T * tan(numbers)
3669#        xnew = x0 + xc
3670#        return xnew
3671
3672    def update_temp(self):
3673        self.T = self.T0/(1+self.k)
3674        self.k += 1
3675        return
3676
3677class boltzmann_sa(base_schedule):
3678#    def update_guess(self, x0):
3679#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
3680#        x0 = asarray(x0)
3681#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
3682#
3683#        xnew = x0 + xc*std*self.learn_rate
3684#        return xnew
3685
3686    def update_temp(self):
3687        self.k += 1
3688        self.T = self.T0 / log(self.k+1.0)
3689        return
3690
3691class log_sa(base_schedule):        #OK
3692
3693    def init(self,**options):
3694        self.__dict__.update(options)
3695       
3696    def update_guess(self,x0):     #same as default
3697        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3698       
3699    def update_temp(self):
3700        self.k += 1
3701        self.T = self.T0*self.slope**self.k
3702       
3703class _state(object):
3704    def __init__(self):
3705        self.x = None
3706        self.cost = None
3707
3708# TODO:
3709#     allow for general annealing temperature profile
3710#     in that case use update given by alpha and omega and
3711#     variation of all previous updates and temperature?
3712
3713# Simulated annealing
3714
3715def anneal(func, x0, args=(), schedule='fast', full_output=0,
3716           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
3717           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
3718           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
3719           ranRange=0.10,autoRan=False,dlg=None):
3720    """Minimize a function using simulated annealing.
3721
3722    Schedule is a schedule class implementing the annealing schedule.
3723    Available ones are 'fast', 'cauchy', 'boltzmann'
3724
3725    :param callable func: f(x, \*args)
3726        Function to be optimized.
3727    :param ndarray x0:
3728        Initial guess.
3729    :param tuple args:
3730        Extra parameters to `func`.
3731    :param base_schedule schedule:
3732        Annealing schedule to use (a class).
3733    :param bool full_output:
3734        Whether to return optional outputs.
3735    :param float T0:
3736        Initial Temperature (estimated as 1.2 times the largest
3737        cost-function deviation over random points in the range).
3738    :param float Tf:
3739        Final goal temperature.
3740    :param int maxeval:
3741        Maximum function evaluations.
3742    :param int maxaccept:
3743        Maximum changes to accept.
3744    :param int maxiter:
3745        Maximum cooling iterations.
3746    :param float learn_rate:
3747        Scale constant for adjusting guesses.
3748    :param float boltzmann:
3749        Boltzmann constant in acceptance test
3750        (increase for less stringent test at each temperature).
3751    :param float feps:
3752        Stopping relative error tolerance for the function value in
3753        last four coolings.
3754    :param float quench,m,n:
3755        Parameters to alter fast_sa schedule.
3756    :param float/ndarray lower,upper:
3757        Lower and upper bounds on `x`.
3758    :param int dwell:
3759        The number of times to search the space at each temperature.
3760    :param float slope:
3761        Parameter for log schedule
3762    :param bool ranStart=False:
3763        True for set 10% of ranges about x
3764
3765    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
3766
3767     * xmin (ndarray): Point giving smallest value found.
3768     * Jmin (float): Minimum value of function found.
3769     * T (float): Final temperature.
3770     * feval (int): Number of function evaluations.
3771     * iters (int): Number of cooling iterations.
3772     * accept (int): Number of tests accepted.
3773     * retval (int): Flag indicating stopping condition:
3774
3775              *  0: Points no longer changing
3776              *  1: Cooled to final temperature
3777              *  2: Maximum function evaluations
3778              *  3: Maximum cooling iterations reached
3779              *  4: Maximum accepted query locations reached
3780              *  5: Final point not the minimum amongst encountered points
3781
3782    *Notes*:
3783    Simulated annealing is a random algorithm which uses no derivative
3784    information from the function being optimized. In practice it has
3785    been more useful in discrete optimization than continuous
3786    optimization, as there are usually better algorithms for continuous
3787    optimization problems.
3788
3789    Some experimentation by trying the difference temperature
3790    schedules and altering their parameters is likely required to
3791    obtain good performance.
3792
3793    The randomness in the algorithm comes from random sampling in numpy.
3794    To obtain the same results you can call numpy.random.seed with the
3795    same seed immediately before calling scipy.optimize.anneal.
3796
3797    We give a brief description of how the three temperature schedules
3798    generate new points and vary their temperature. Temperatures are
3799    only updated with iterations in the outer loop. The inner loop is
3800    over xrange(dwell), and new points are generated for
3801    every iteration in the inner loop. (Though whether the proposed
3802    new points are accepted is probabilistic.)
3803
3804    For readability, let d denote the dimension of the inputs to func.
3805    Also, let x_old denote the previous state, and k denote the
3806    iteration number of the outer loop. All other variables not
3807    defined below are input variables to scipy.optimize.anneal itself.
3808
3809    In the 'fast' schedule the updates are ::
3810
3811        u ~ Uniform(0, 1, size=d)
3812        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
3813        xc = y * (upper - lower)
3814        x_new = x_old + xc
3815
3816        c = n * exp(-n * quench)
3817        T_new = T0 * exp(-c * k**quench)
3818
3819
3820    In the 'cauchy' schedule the updates are ::
3821
3822        u ~ Uniform(-pi/2, pi/2, size=d)
3823        xc = learn_rate * T * tan(u)
3824        x_new = x_old + xc
3825
3826        T_new = T0 / (1+k)
3827
3828    In the 'boltzmann' schedule the updates are ::
3829
3830        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
3831        y ~ Normal(0, std, size=d)
3832        x_new = x_old + learn_rate * y
3833
3834        T_new = T0 / log(1+k)
3835
3836    """
3837    x0 = asarray(x0)
3838    lower = asarray(lower)
3839    upper = asarray(upper)
3840
3841    schedule = eval(schedule+'_sa()')
3842    #   initialize the schedule
3843    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
3844                  learn_rate=learn_rate, lower=lower, upper=upper,
3845                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
3846
3847    current_state, last_state, best_state = _state(), _state(), _state()
3848    if ranStart:
3849        schedule.set_range(x0,ranRange)
3850    if T0 is None:
3851        x0 = schedule.getstart_temp(best_state)
3852    else:
3853        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
3854        best_state.x = None
3855        best_state.cost = numpy.Inf
3856
3857    last_state.x = asarray(x0).copy()
3858    fval = func(x0,*args)
3859    schedule.feval += 1
3860    last_state.cost = fval
3861    if last_state.cost < best_state.cost:
3862        best_state.cost = fval
3863        best_state.x = asarray(x0).copy()
3864    schedule.T = schedule.T0
3865    fqueue = [100, 300, 500, 700]
3866    iters = 1
3867    keepGoing = True
3868    bestn = 0
3869    while keepGoing:
3870        retval = 0
3871        for n in xrange(dwell):
3872            current_state.x = schedule.update_guess(last_state.x)
3873            current_state.cost = func(current_state.x,*args)
3874            schedule.feval += 1
3875
3876            dE = current_state.cost - last_state.cost
3877            if schedule.accept_test(dE):
3878                last_state.x = current_state.x.copy()
3879                last_state.cost = current_state.cost
3880                if last_state.cost < best_state.cost:
3881                    best_state.x = last_state.x.copy()
3882                    best_state.cost = last_state.cost
3883                    bestn = n
3884                    if best_state.cost < 1.0 and autoRan:
3885                        schedule.set_range(x0,best_state.cost/2.)                       
3886        if dlg:
3887            GoOn = dlg.Update(min(100.,best_state.cost*100),
3888                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
3889                    'Best trial:',bestn,  \
3890                    'MC/SA Residual:',best_state.cost*100,'%', \
3891                    ))[0]
3892            if not GoOn:
3893                best_state.x = last_state.x.copy()
3894                best_state.cost = last_state.cost
3895                retval = 5
3896        schedule.update_temp()
3897        iters += 1
3898        # Stopping conditions
3899        # 0) last saved values of f from each cooling step
3900        #     are all very similar (effectively cooled)
3901        # 1) Tf is set and we are below it
3902        # 2) maxeval is set and we are past it
3903        # 3) maxiter is set and we are past it
3904        # 4) maxaccept is set and we are past it
3905        # 5) user canceled run via progress bar
3906
3907        fqueue.append(squeeze(last_state.cost))
3908        fqueue.pop(0)
3909        af = asarray(fqueue)*1.0
3910        if retval == 5:
3911            print ' User terminated run; incomplete MC/SA'
3912            keepGoing = False
3913            break
3914        if all(abs((af-af[0])/af[0]) < feps):
3915            retval = 0
3916            if abs(af[-1]-best_state.cost) > feps*10:
3917                retval = 5
3918#                print "Warning: Cooled to %f at %s but this is not" \
3919#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
3920#                      + " the smallest point found."
3921            break
3922        if (Tf is not None) and (schedule.T < Tf):
3923            retval = 1
3924            break
3925        if (maxeval is not None) and (schedule.feval > maxeval):
3926            retval = 2
3927            break
3928        if (iters > maxiter):
3929            print "Warning: Maximum number of iterations exceeded."
3930            retval = 3
3931            break
3932        if (maxaccept is not None) and (schedule.accepted > maxaccept):
3933            retval = 4
3934            break
3935
3936    if full_output:
3937        return best_state.x, best_state.cost, schedule.T, \
3938               schedule.feval, iters, schedule.accepted, retval
3939    else:
3940        return best_state.x, retval
3941
3942def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
3943    outlist = []
3944    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
3945    for n in range(iCyc):
3946        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
3947        outlist.append(result[0])
3948        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
3949    out_q.put(outlist)
3950
3951def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
3952    import multiprocessing as mp
3953   
3954    nprocs = mp.cpu_count()
3955    out_q = mp.Queue()
3956    procs = []
3957    iCyc = np.zeros(nprocs)
3958    for i in range(nCyc):
3959        iCyc[i%nprocs] += 1
3960    for i in range(nprocs):
3961        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
3962        procs.append(p)
3963        p.start()
3964    resultlist = []
3965    for i in range(nprocs):
3966        resultlist += out_q.get()
3967    for p in procs:
3968        p.join()
3969    return resultlist
3970
3971def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
3972    '''default doc string
3973   
3974    :param type name: description
3975   
3976    :returns: type name: description
3977    '''
3978   
3979    global tsum
3980    tsum = 0.
3981   
3982    def getMDparms(item,pfx,parmDict,varyList):
3983        parmDict[pfx+'MDaxis'] = item['axis']
3984        parmDict[pfx+'MDval'] = item['Coef'][0]
3985        if item['Coef'][1]:
3986            varyList += [pfx+'MDval',]
3987            limits = item['Coef'][2]
3988            lower.append(limits[0])
3989            upper.append(limits[1])
3990                       
3991    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
3992        parmDict[pfx+'Atype'] = item['atType']
3993        aTypes |= set([item['atType'],]) 
3994        pstr = ['Ax','Ay','Az']
3995        XYZ = [0,0,0]
3996        for i in range(3):
3997            name = pfx+pstr[i]
3998            parmDict[name] = item['Pos'][0][i]
3999            XYZ[i] = parmDict[name]
4000            if item['Pos'][1][i]:
4001                varyList += [name,]
4002                limits = item['Pos'][2][i]
4003                lower.append(limits[0])
4004                upper.append(limits[1])
4005        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4006           
4007    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4008        parmDict[mfx+'MolCent'] = item['MolCent']
4009        parmDict[mfx+'RBId'] = item['RBId']
4010        pstr = ['Px','Py','Pz']
4011        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4012        for i in range(3):
4013            name = mfx+pstr[i]
4014            parmDict[name] = item['Pos'][0][i]
4015            if item['Pos'][1][i]:
4016                varyList += [name,]
4017                limits = item['Pos'][2][i]
4018                lower.append(limits[0])
4019                upper.append(limits[1])
4020        AV = item['Ori'][0]
4021        A = AV[0]
4022        V = AV[1:]
4023        for i in range(4):
4024            name = mfx+ostr[i]
4025            if i:
4026                parmDict[name] = V[i-1]
4027            else:
4028                parmDict[name] = A
4029            if item['Ovar'] == 'AV':
4030                varyList += [name,]
4031                limits = item['Ori'][2][i]
4032                lower.append(limits[0])
4033                upper.append(limits[1])
4034            elif item['Ovar'] == 'A' and not i:
4035                varyList += [name,]
4036                limits = item['Ori'][2][i]
4037                lower.append(limits[0])
4038                upper.append(limits[1])
4039        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
4040            for i in range(len(item['Tor'][0])):
4041                name = mfx+'Tor'+str(i)
4042                parmDict[name] = item['Tor'][0][i]
4043                if item['Tor'][1][i]:
4044                    varyList += [name,]
4045                    limits = item['Tor'][2][i]
4046                    lower.append(limits[0])
4047                    upper.append(limits[1])
4048        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
4049        aTypes |= set(atypes)
4050        atNo += len(atypes)
4051        return atNo
4052               
4053    def GetAtomM(Xdata,SGData):
4054        Mdata = []
4055        for xyz in Xdata:
4056            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
4057        return np.array(Mdata)
4058       
4059    def GetAtomT(RBdata,parmDict):
4060        'Needs a doc string'
4061        atNo = parmDict['atNo']
4062        nfixAt = parmDict['nfixAt']
4063        Tdata = atNo*[' ',]
4064        for iatm in range(nfixAt):
4065            parm = ':'+str(iatm)+':Atype'
4066            if parm in parmDict:
4067                Tdata[iatm] = aTypes.index(parmDict[parm])
4068        iatm = nfixAt
4069        for iObj in range(parmDict['nObj']):
4070            pfx = str(iObj)+':'
4071            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4072                if parmDict[pfx+'Type'] == 'Vector':
4073                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4074                    nAtm = len(RBRes['rbVect'][0])
4075                else:       #Residue
4076                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4077                    nAtm = len(RBRes['rbXYZ'])
4078                for i in range(nAtm):
4079                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
4080                    iatm += 1
4081            elif parmDict[pfx+'Type'] == 'Atom':
4082                atNo = parmDict[pfx+'atNo']
4083                parm = pfx+'Atype'              #remove extra ':'
4084                if parm in parmDict:
4085                    Tdata[atNo] = aTypes.index(parmDict[parm])
4086                iatm += 1
4087            else:
4088                continue        #skips March Dollase
4089        return Tdata
4090       
4091    def GetAtomX(RBdata,parmDict):
4092        'Needs a doc string'
4093        Bmat = parmDict['Bmat']
4094        atNo = parmDict['atNo']
4095        nfixAt = parmDict['nfixAt']
4096        Xdata = np.zeros((3,atNo))
4097        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
4098        for iatm in range(nfixAt):
4099            for key in keys:
4100                parm = ':'+str(iatm)+key
4101                if parm in parmDict:
4102                    keys[key][iatm] = parmDict[parm]
4103        iatm = nfixAt
4104        for iObj in range(parmDict['nObj']):
4105            pfx = str(iObj)+':'
4106            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4107                if parmDict[pfx+'Type'] == 'Vector':
4108                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4109                    vecs = RBRes['rbVect']
4110                    mags = RBRes['VectMag']
4111                    Cart = np.zeros_like(vecs[0])
4112                    for vec,mag in zip(vecs,mags):
4113                        Cart += vec*mag
4114                elif parmDict[pfx+'Type'] == 'Residue':
4115                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4116                    Cart = np.array(RBRes['rbXYZ'])
4117                    for itor,seq in enumerate(RBRes['rbSeq']):
4118                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
4119                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
4120                if parmDict[pfx+'MolCent'][1]:
4121                    Cart -= parmDict[pfx+'MolCent'][0]
4122                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4123                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
4124                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
4125                iatm += len(Cart)
4126            elif parmDict[pfx+'Type'] == 'Atom':
4127                atNo = parmDict[pfx+'atNo']
4128                for key in keys:
4129                    parm = pfx+key[1:]              #remove extra ':'
4130                    if parm in parmDict:
4131                        keys[key][atNo] = parmDict[parm]
4132                iatm += 1
4133            else:
4134                continue        #skips March Dollase
4135        return Xdata.T
4136       
4137    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
4138        allX = np.inner(Xdata,SGM)+SGT
4139        allT = np.repeat(Tdata,allX.shape[1])
4140        allM = np.repeat(Mdata,allX.shape[1])
4141        allX = np.reshape(allX,(-1,3))
4142        return allT,allM,allX
4143
4144    def getAllX(Xdata,SGM,SGT):
4145        allX = np.inner(Xdata,SGM)+SGT
4146        allX = np.reshape(allX,(-1,3))
4147        return allX
4148       
4149    def normQuaternions(RBdata,parmDict,varyList,values):
4150        for iObj in range(parmDict['nObj']):
4151            pfx = str(iObj)+':'
4152            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4153                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4154                A,V = Q2AVdeg(Qori)
4155                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
4156                    if i:
4157                        parmDict[pfx+name] = V[i-1]
4158                    else:
4159                        parmDict[pfx+name] = A
4160       
4161    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
4162        ''' Compute structure factors for all h,k,l for phase
4163        input:
4164            refList: [ref] where each ref = h,k,l,m,d,...
4165            rcov:   array[nref,nref] covariance terms between Fo^2 values
4166            ifInv:  bool True if centrosymmetric
4167            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
4168            RBdata: [dict] rigid body dictionary
4169            varyList: [list] names of varied parameters in MC/SA (not used here)           
4170            ParmDict: [dict] problem parameters
4171        puts result F^2 in each ref[5] in refList
4172        returns:
4173            delt-F*rcov*delt-F/sum(Fo^2)^2
4174        '''       
4175        global tsum
4176        t0 = time.time()
4177        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
4178        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
4179        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
4180        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
4181        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
4182        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
4183        refList[5] = Aterm
4184        if not ifInv:
4185            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
4186        if len(cosTable):        #apply MD correction
4187            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
4188        sumFcsq = np.sum(refList[5])
4189        scale = parmDict['sumFosq']/sumFcsq
4190        refList[5] *= scale
4191        refList[6] = refList[4]-refList[5]
4192        M = np.inner(refList[6],np.inner(rcov,refList[6]))
4193        tsum += (time.time()-t0)
4194        return M/np.sum(refList[4]**2)
4195
4196    sq8ln2 = np.sqrt(8*np.log(2))
4197    sq2pi = np.sqrt(2*np.pi)
4198    sq4pi = np.sqrt(4*np.pi)
4199    generalData = data['General']
4200    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
4201    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
4202    SGData = generalData['SGData']
4203    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
4204    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
4205    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
4206    fixAtoms = data['Atoms']                       #if any
4207    cx,ct,cs = generalData['AtomPtrs'][:3]
4208    aTypes = set([])
4209    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
4210    varyList = []
4211    atNo = 0
4212    for atm in fixAtoms:
4213        pfx = ':'+str(atNo)+':'
4214        parmDict[pfx+'Atype'] = atm[ct]
4215        aTypes |= set([atm[ct],])
4216        pstr = ['Ax','Ay','Az']
4217        parmDict[pfx+'Amul'] = atm[cs+1]
4218        for i in range(3):
4219            name = pfx+pstr[i]
4220            parmDict[name] = atm[cx+i]
4221        atNo += 1
4222    parmDict['nfixAt'] = len(fixAtoms)       
4223    MCSA = generalData['MCSA controls']
4224    reflName = MCSA['Data source']
4225    phaseName = generalData['Name']
4226    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
4227    upper = []
4228    lower = []
4229    MDvec = np.zeros(3)
4230    for i,item in enumerate(MCSAObjs):
4231        mfx = str(i)+':'
4232        parmDict[mfx+'Type'] = item['Type']
4233        if item['Type'] == 'MD':
4234            getMDparms(item,mfx,parmDict,varyList)
4235            MDvec = np.array(item['axis'])
4236        elif item['Type'] == 'Atom':
4237            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
4238            parmDict[mfx+'atNo'] = atNo
4239            atNo += 1
4240        elif item['Type'] in ['Residue','Vector']:
4241            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
4242    parmDict['atNo'] = atNo                 #total no. of atoms
4243    parmDict['nObj'] = len(MCSAObjs)
4244    aTypes = list(aTypes)
4245    Tdata = GetAtomT(RBdata,parmDict)
4246    Xdata = GetAtomX(RBdata,parmDict)
4247    Mdata = GetAtomM(Xdata,SGData)
4248    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
4249    FFtables = G2el.GetFFtable(aTypes)
4250    refs = []
4251    allFF = []
4252    cosTable = []
4253    sumFosq = 0
4254    if 'PWDR' in reflName:
4255        for ref in reflData:
4256            h,k,l,m,d,pos,sig,gam,f = ref[:9]
4257            if d >= MCSA['dmin']:
4258                sig = np.sqrt(sig)      #var -> sig in centideg
4259                sig = G2pwd.getgamFW(gam,sig)/sq8ln2        #FWHM -> sig in centideg
4260                SQ = 0.25/d**2
4261                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
4262                refs.append([h,k,l,m,f*m,pos,sig])
4263                sumFosq += f*m
4264                Heqv = np.inner(np.array([h,k,l]),SGMT)
4265                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
4266        nRef = len(refs)
4267        cosTable = np.array(cosTable)**2
4268        rcov = np.zeros((nRef,nRef))
4269        for iref,refI in enumerate(refs):
4270            rcov[iref][iref] = 1./(sq4pi*refI[6])
4271            for jref,refJ in enumerate(refs[:iref]):
4272                t1 = refI[6]**2+refJ[6]**2
4273                t2 = (refJ[5]-refI[5])**2/(2.*t1)
4274                if t2 > 10.:
4275                    rcov[iref][jref] = 0.
4276                else:
4277                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
4278        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
4279        Rdiag = np.sqrt(np.diag(rcov))
4280        Rnorm = np.outer(Rdiag,Rdiag)
4281        rcov /= Rnorm
4282    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
4283        vNames = []
4284        pfx = str(data['pId'])+'::PWLref:'
4285        for iref,refI in enumerate(reflData):           #Pawley reflection set
4286            h,k,l,m,d,v,f,s = refI
4287            if d >= MCSA['dmin'] and v:       #skip unrefined ones
4288                vNames.append(pfx+str(iref))
4289                SQ = 0.25/d**2
4290                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
4291                refs.append([h,k,l,m,f*m,iref,0.])
4292                sumFosq += f*m
4293                Heqv = np.inner(np.array([h,k,l]),SGMT)
4294                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
4295        cosTable = np.array(cosTable)**2
4296        nRef = len(refs)
4297#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
4298        if covData['freshCOV'] or  MCSA['newDmin']:
4299            vList = covData['varyList']
4300            covMatrix = covData['covMatrix']
4301            rcov = getVCov(vNames,vList,covMatrix)
4302            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
4303            Rdiag = np.sqrt(np.diag(rcov))
4304            Rdiag = np.where(Rdiag,Rdiag,1.0)
4305            Rnorm = np.outer(Rdiag,Rdiag)
4306            rcov /= Rnorm
4307            MCSA['rcov'] = rcov
4308            covData['freshCOV'] = False
4309            MCSA['newDmin'] = False
4310        else:
4311            rcov = MCSA['rcov']
4312    elif 'HKLF' in reflName:
4313        for ref in reflData:
4314            [h,k,l,m,d],f = ref[:5],ref[6]
4315            if d >= MCSA['dmin']:
4316                SQ = 0.25/d**2
4317                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
4318                refs.append([h,k,l,m,f*m,0.,0.])
4319                sumFosq += f*m
4320        nRef = len(refs)
4321        rcov = np.identity(len(refs))
4322    allFF = np.array(allFF).T
4323    refs = np.array(refs).T
4324    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
4325    print ' Number of parameters varied: %d'%(len(varyList))
4326    parmDict['sumFosq'] = sumFosq
4327    x0 = [parmDict[val] for val in varyList]
4328    ifInv = SGData['SGInv']
4329    # consider replacing anneal with scipy.optimize.basinhopping
4330    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
4331        schedule=MCSA['Algorithm'], full_output=True,
4332        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
4333        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
4334        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
4335        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
4336        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
4337    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
4338#    for ref in refs.T:
4339#        print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
4340#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
4341    Result = [False,False,results[1],results[2],]+list(results[0])
4342    Result.append(varyList)
4343    return Result,tsum
4344
4345       
4346################################################################################
4347##### Quaternion stuff
4348################################################################################
4349
4350def prodQQ(QA,QB):
4351    ''' Grassman quaternion product
4352        QA,QB quaternions; q=r+ai+bj+ck
4353    '''
4354    D = np.zeros(4)
4355    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
4356    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
4357    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
4358    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
4359   
4360#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
4361#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
4362   
4363    return D
4364   
4365def normQ(QA):
4366    ''' get length of quaternion & normalize it
4367        q=r+ai+bj+ck
4368    '''
4369    n = np.sqrt(np.sum(np.array(QA)**2))
4370    return QA/n
4371   
4372def invQ(Q):
4373    '''
4374        get inverse of quaternion
4375        q=r+ai+bj+ck; q* = r-ai-bj-ck
4376    '''
4377    return Q*np.array([1,-1,-1,-1])
4378   
4379def prodQVQ(Q,V):
4380    """
4381    compute the quaternion vector rotation qvq-1 = v'
4382    q=r+ai+bj+ck
4383    """
4384    T2 = Q[0]*Q[1]
4385    T3 = Q[0]*Q[2]
4386    T4 = Q[0]*Q[3]
4387    T5 = -Q[1]*Q[1]
4388    T6 = Q[1]*Q[2]
4389    T7 = Q[1]*Q[3]
4390    T8 = -Q[2]*Q[2]
4391    T9 = Q[2]*Q[3]
4392    T10 = -Q[3]*Q[3]
4393    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
4394    VP = 2.*np.inner(V,M)
4395    return VP+V
4396   
4397def Q2Mat(Q):
4398    ''' make rotation matrix from quaternion
4399        q=r+ai+bj+ck
4400    '''
4401    QN = normQ(Q)
4402    aa = QN[0]**2
4403    ab = QN[0]*QN[1]
4404    ac = QN[0]*QN[2]
4405    ad = QN[0]*QN[3]
4406    bb = QN[1]**2
4407    bc = QN[1]*QN[2]
4408    bd = QN[1]*QN[3]
4409    cc = QN[2]**2
4410    cd = QN[2]*QN[3]
4411    dd = QN[3]**2
4412    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
4413        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
4414        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
4415    return np.array(M)
4416   
4417def AV2Q(A,V):
4418    ''' convert angle (radians) & vector to quaternion
4419        q=r+ai+bj+ck
4420    '''
4421    Q = np.zeros(4)
4422    d = nl.norm(np.array(V))
4423    if d:
4424        V /= d
4425        if not A:       #==0.
4426            A = 2.*np.pi
4427        p = A/2.
4428        Q[0] = np.cos(p)
4429        Q[1:4] = V*np.sin(p)
4430    else:
4431        Q[3] = 1.
4432    return Q
4433   
4434def AVdeg2Q(A,V):
4435    ''' convert angle (degrees) & vector to quaternion
4436        q=r+ai+bj+ck
4437    '''
4438    Q = np.zeros(4)
4439    d = nl.norm(np.array(V))
4440    if not A:       #== 0.!
4441        A = 360.
4442    if d:
4443        V /= d
4444        p = A/2.
4445        Q[0] = cosd(p)
4446        Q[1:4] = V*sind(p)
4447    else:
4448        Q[3] = 1.
4449    return Q
4450   
4451def Q2AVdeg(Q):
4452    ''' convert quaternion to angle (degrees 0-360) & normalized vector
4453        q=r+ai+bj+ck
4454    '''
4455    A = 2.*acosd(Q[0])
4456    V = np.array(Q[1:])
4457    V /= sind(A/2.)
4458    return A,V
4459   
4460def Q2AV(Q):
4461    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
4462        q=r+ai+bj+ck
4463    '''
4464    A = 2.*np.arccos(Q[0])
4465    V = np.array(Q[1:])
4466    V /= np.sin(A/2.)
4467    return A,V
4468   
4469def randomQ(r0,r1,r2,r3):
4470    ''' create random quaternion from 4 random numbers in range (-1,1)
4471    '''
4472    sum = 0
4473    Q = np.array(4)
4474    Q[0] = r0
4475    sum += Q[0]**2
4476    Q[1] = np.sqrt(1.-sum)*r1
4477    sum += Q[1]**2
4478    Q[2] = np.sqrt(1.-sum)*r2
4479    sum += Q[2]**2
4480    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
4481    return Q
4482   
4483def randomAVdeg(r0,r1,r2,r3):
4484    ''' create random angle (deg),vector from 4 random number in range (-1,1)
4485    '''
4486    return Q2AVdeg(randomQ(r0,r1,r2,r3))
4487   
4488def makeQuat(A,B,C):
4489    ''' Make quaternion from rotation of A vector to B vector about C axis
4490
4491        :param np.array A,B,C: Cartesian 3-vectors
4492        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
4493    '''
4494
4495    V1 = np.cross(A,C)
4496    V2 = np.cross(B,C)
4497    if nl.norm(V1)*nl.norm(V2):
4498        V1 /= nl.norm(V1)
4499        V2 /= nl.norm(V2)
4500        V3 = np.cross(V1,V2)
4501    else:
4502        V3 = np.zeros(3)
4503    Q = np.array([0.,0.,0.,1.])
4504    D = 0.
4505    if nl.norm(V3):
4506        V3 /= nl.norm(V3)
4507        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
4508        D = np.arccos(D1)/2.0
4509        V1 = C-V3
4510        V2 = C+V3
4511        DM = nl.norm(V1)
4512        DP = nl.norm(V2)
4513        S = np.sin(D)
4514        Q[0] = np.cos(D)
4515        Q[1:] = V3*S
4516        D *= 2.
4517        if DM > DP:
4518            D *= -1.
4519    return Q,D
4520   
4521def annealtests():
4522    from numpy import cos
4523    # minimum expected at ~-0.195
4524    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
4525    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
4526    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
4527    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
4528
4529    # minimum expected at ~[-0.195, -0.1]
4530    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
4531    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
4532    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
4533    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
4534
4535
4536if __name__ == '__main__':
4537    annealtests()
Note: See TracBrowser for help on using the repository browser.