source: trunk/GSASIImath.py @ 2006

Last change on this file since 2006 was 2006, checked in by vondreele, 6 years ago

minor work on SS structure factors

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