source: trunk/GSASIImath.py @ 2043

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

some SS deriv work
trap G2frame.PatternId? = None (or 0) in PLotPatterns
put missing import numpy in G2img_1TIF.py

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