source: trunk/GSASIIstrMath.py @ 939

Last change on this file since 939 was 939, checked in by toby, 9 years ago

fix & cleanup unit tests; add/change doc strings for sphinx; add all G2 py files to sphinx

File size: 95.1 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
8# $Author: vondreele $
9# $Revision: 920 $
10# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 920 2013-05-17 17:13:15Z vondreele $
12########### SVN repository information ###################
13import time
14import math
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.optimize as so
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 920 $")
22import GSASIIElem as G2el
23import GSASIIlattice as G2lat
24import GSASIIspc as G2spc
25import GSASIIpwd as G2pwd
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28
29sind = lambda x: np.sin(x*np.pi/180.)
30cosd = lambda x: np.cos(x*np.pi/180.)
31tand = lambda x: np.tan(x*np.pi/180.)
32asind = lambda x: 180.*np.arcsin(x)/np.pi
33acosd = lambda x: 180.*np.arccos(x)/np.pi
34atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35   
36ateln2 = 8.0*math.log(2.0)
37
38################################################################################
39##### Rigid Body Models
40################################################################################
41       
42def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
43    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
44    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
45    '''
46    atxIds = ['Ax:','Ay:','Az:']
47    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
48    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
49    if not RBIds['Vector'] and not RBIds['Residue']:
50        return
51    VRBIds = RBIds['Vector']
52    RRBIds = RBIds['Residue']
53    if Update:
54        RBData = rigidbodyDict
55    else:
56        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
57    if RBIds['Vector']:                       # first update the vector magnitudes
58        VRBData = RBData['Vector']
59        for i,rbId in enumerate(VRBIds):
60            for j in range(len(VRBData[rbId]['VectMag'])):
61                name = '::RBV;'+str(j)+':'+str(i)
62                VRBData[rbId]['VectMag'][j] = parmDict[name]
63    for phase in Phases:
64        Phase = Phases[phase]
65        General = Phase['General']
66        cell = General['Cell'][1:7]
67        Amat,Bmat = G2lat.cell2AB(cell)
68        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
69        pfx = str(Phase['pId'])+'::'
70        if Update:
71            RBModels = Phase['RBModels']
72        else:
73            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
74        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
75            jrb = VRBIds.index(RBObj['RBId'])
76            rbsx = str(irb)+':'+str(jrb)
77            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
78                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
79            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
80                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
81            TLS = RBObj['ThermalMotion']
82            if 'T' in TLS[0]:
83                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
84                    TLS[1][i] = parmDict[pfx+pt+rbsx]
85            if 'L' in TLS[0]:
86                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
87                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
88            if 'S' in TLS[0]:
89                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
90                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
91            if 'U' in TLS[0]:
92                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
93            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
94            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
95            for i,x in enumerate(XYZ):
96                atId = RBObj['Ids'][i]
97                for j in [0,1,2]:
98                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
99                if UIJ[i][0] == 'A':
100                    for j in range(6):
101                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
102                elif UIJ[i][0] == 'I':
103                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
104           
105        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
106            jrb = RRBIds.index(RBObj['RBId'])
107            rbsx = str(irb)+':'+str(jrb)
108            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
109                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
110            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
111                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
112            TLS = RBObj['ThermalMotion']
113            if 'T' in TLS[0]:
114                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
115                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
116            if 'L' in TLS[0]:
117                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
118                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
119            if 'S' in TLS[0]:
120                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
121                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
122            if 'U' in TLS[0]:
123                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
124            for itors,tors in enumerate(RBObj['Torsions']):
125                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
126            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
127            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
128            for i,x in enumerate(XYZ):
129                atId = RBObj['Ids'][i]
130                for j in [0,1,2]:
131                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
132                if UIJ[i][0] == 'A':
133                    for j in range(6):
134                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
135                elif UIJ[i][0] == 'I':
136                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
137                   
138def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
139    'Needs a doc string'
140    atxIds = ['dAx:','dAy:','dAz:']
141    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
142    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
143    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
144    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
145    PIds = ['Px:','Py:','Pz:']
146    OIds = ['Oa:','Oi:','Oj:','Ok:']
147    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
148    if not RBIds['Vector'] and not RBIds['Residue']:
149        return
150    VRBIds = RBIds['Vector']
151    RRBIds = RBIds['Residue']
152    RBData = rigidbodyDict
153    for item in parmDict:
154        if 'RB' in item:
155            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
156    General = Phase['General']
157    cell = General['Cell'][1:7]
158    Amat,Bmat = G2lat.cell2AB(cell)
159    rpd = np.pi/180.
160    rpd2 = rpd**2
161    g = nl.inv(np.inner(Bmat,Bmat))
162    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
163        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
164    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
165    pfx = str(Phase['pId'])+'::'
166    RBModels =  Phase['RBModels']
167    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
168        VModel = RBData['Vector'][RBObj['RBId']]
169        Q = RBObj['Orient'][0]
170        Pos = RBObj['Orig'][0]
171        jrb = VRBIds.index(RBObj['RBId'])
172        rbsx = str(irb)+':'+str(jrb)
173        dXdv = []
174        for iv in range(len(VModel['VectMag'])):
175            dCdv = []
176            for vec in VModel['rbVect'][iv]:
177                dCdv.append(G2mth.prodQVQ(Q,vec))
178            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
179        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
180        for ia,atId in enumerate(RBObj['Ids']):
181            atNum = AtLookup[atId]
182            dx = 0.0001
183            for iv in range(len(VModel['VectMag'])):
184                for ix in [0,1,2]:
185                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
186            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
187                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
188            for iv in range(4):
189                Q[iv] -= dx
190                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
191                Q[iv] += 2.*dx
192                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
193                Q[iv] -= dx
194                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
195                for ix in [0,1,2]:
196                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
197            X = G2mth.prodQVQ(Q,Cart[ia])
198            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
199            dFdu = G2lat.U6toUij(dFdu.T)
200            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
201            dFdu = G2lat.UijtoU6(dFdu)
202            atNum = AtLookup[atId]
203            if 'T' in RBObj['ThermalMotion'][0]:
204                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
205                    dFdvDict[pfx+name+rbsx] += dFdu[i]
206            if 'L' in RBObj['ThermalMotion'][0]:
207                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
208                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
209                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
210                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
211                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
212                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
213                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
214                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
215                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
216            if 'S' in RBObj['ThermalMotion'][0]:
217                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
218                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
219                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
220                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
221                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
222                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
223                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
224                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
225            if 'U' in RBObj['ThermalMotion'][0]:
226                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
227
228
229    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
230        Q = RBObj['Orient'][0]
231        Pos = RBObj['Orig'][0]
232        jrb = RRBIds.index(RBObj['RBId'])
233        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
234        rbsx = str(irb)+':'+str(jrb)
235        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
236        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
237            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
238            orId,pvId = torData[itors][:2]
239            pivotVec = Cart[orId]-Cart[pvId]
240            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
241            QB = G2mth.AVdeg2Q(0.001,pivotVec)
242            for ir in torData[itors][3]:
243                atNum = AtLookup[RBObj['Ids'][ir]]
244                rVec = Cart[ir]-Cart[pvId]
245                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
246                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
247                for ix in [0,1,2]:
248                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
249        for ia,atId in enumerate(RBObj['Ids']):
250            atNum = AtLookup[atId]
251            dx = 0.0001
252            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
253                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
254            for iv in range(4):
255                Q[iv] -= dx
256                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
257                Q[iv] += 2.*dx
258                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
259                Q[iv] -= dx
260                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
261                for ix in [0,1,2]:
262                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
263            X = G2mth.prodQVQ(Q,Cart[ia])
264            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
265            dFdu = G2lat.U6toUij(dFdu.T)
266            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
267            dFdu = G2lat.UijtoU6(dFdu)
268            atNum = AtLookup[atId]
269            if 'T' in RBObj['ThermalMotion'][0]:
270                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
271                    dFdvDict[pfx+name+rbsx] += dFdu[i]
272            if 'L' in RBObj['ThermalMotion'][0]:
273                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
274                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
275                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
276                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
277                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
278                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
279                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
280                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
281                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
282            if 'S' in RBObj['ThermalMotion'][0]:
283                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
284                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
285                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
286                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
287                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
288                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
289                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
290                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
291            if 'U' in RBObj['ThermalMotion'][0]:
292                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
293   
294################################################################################
295##### Penalty & restraint functions
296################################################################################
297
298def penaltyFxn(HistoPhases,parmDict,varyList):
299    'Needs a doc string'
300    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
301    pNames = []
302    pVals = []
303    pWt = []
304    negWt = {}
305    for phase in Phases:
306        pId = Phases[phase]['pId']
307        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
308        General = Phases[phase]['General']
309        textureData = General['SH Texture']
310        SGData = General['SGData']
311        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
312        cell = General['Cell'][1:7]
313        Amat,Bmat = G2lat.cell2AB(cell)
314        if phase not in restraintDict:
315            continue
316        phaseRest = restraintDict[phase]
317        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
318            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
319            ['ChemComp','Sites'],['Texture','HKLs']]
320        for name,rest in names:
321            itemRest = phaseRest[name]
322            if itemRest[rest] and itemRest['Use']:
323                wt = itemRest['wtFactor']
324                if name in ['Bond','Angle','Plane','Chiral']:
325                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
326                        pNames.append(str(pId)+':'+name+':'+str(i))
327                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
328                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
329                        if name == 'Bond':
330                            calc = G2mth.getRestDist(XYZ,Amat)
331                        elif name == 'Angle':
332                            calc = G2mth.getRestAngle(XYZ,Amat)
333                        elif name == 'Plane':
334                            calc = G2mth.getRestPlane(XYZ,Amat)
335                        elif name == 'Chiral':
336                            calc = G2mth.getRestChiral(XYZ,Amat)
337                        pVals.append(obs-calc)
338                        pWt.append(wt/esd**2)
339                elif name in ['Torsion','Rama']:
340                    coeffDict = itemRest['Coeff']
341                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
342                        pNames.append(str(pId)+':'+name+':'+str(i))
343                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
344                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
345                        if name == 'Torsion':
346                            tor = G2mth.getRestTorsion(XYZ,Amat)
347                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
348                        else:
349                            phi,psi = G2mth.getRestRama(XYZ,Amat)
350                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
351                        pVals.append(obs-calc)
352                        pWt.append(wt/esd**2)
353                elif name == 'ChemComp':
354                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
355                        pNames.append(str(pId)+':'+name+':'+str(i))
356                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
357                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
358                        calc = np.sum(mul*frac*factors)
359                        pVals.append(obs-calc)
360                        pWt.append(wt/esd**2)                   
361                elif name == 'Texture':
362                    SHkeys = textureData['SH Coeff'][1].keys()
363                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
364                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
365                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
366                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
367                        PH = np.array(hkl)
368                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
369                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
370                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
371                        Z1 = -ma.masked_greater(Z,0.0)
372                        IndZ1 = np.array(ma.nonzero(Z1))
373                        for ind in IndZ1.T:
374                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
375                            pVals.append(Z1[ind[0]][ind[1]])
376                            pWt.append(wt/esd1**2)
377                        if ifesd2:
378                            Z2 = 1.-Z
379                            for ind in np.ndindex(grid,grid):
380                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
381                                pVals.append(Z1[ind[0]][ind[1]])
382                                pWt.append(wt/esd2**2)
383         
384    for item in varyList:
385        if 'PWLref' in item and parmDict[item] < 0.:
386            pId = int(item.split(':')[0])
387            if negWt[pId]:
388                pNames.append(item)
389                pVals.append(-parmDict[item])
390                pWt.append(negWt[pId])
391    pVals = np.array(pVals)
392    pWt = np.array(pWt)         #should this be np.sqrt?
393    return pNames,pVals,pWt
394   
395def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
396    'Needs a doc string'
397    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
398    pDerv = np.zeros((len(varyList),len(pVal)))
399    for phase in Phases:
400#        if phase not in restraintDict:
401#            continue
402        pId = Phases[phase]['pId']
403        General = Phases[phase]['General']
404        SGData = General['SGData']
405        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
406        cell = General['Cell'][1:7]
407        Amat,Bmat = G2lat.cell2AB(cell)
408        textureData = General['SH Texture']
409
410        SHkeys = textureData['SH Coeff'][1].keys()
411        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
412        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
413        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
414        sam = SamSym[textureData['Model']]
415        phaseRest = restraintDict.get(phase,{})
416        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
417            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
418            'ChemComp':'Sites','Texture':'HKLs'}
419        lasthkl = np.array([0,0,0])
420        for ip,pName in enumerate(pNames):
421            pnames = pName.split(':')
422            if pId == int(pnames[0]):
423                name = pnames[1]
424                if 'PWL' in pName:
425                    pDerv[varyList.index(pName)][ip] += 1.
426                    continue
427                id = int(pnames[2]) 
428                itemRest = phaseRest[name]
429                if name in ['Bond','Angle','Plane','Chiral']:
430                    indx,ops,obs,esd = itemRest[names[name]][id]
431                    dNames = []
432                    for ind in indx:
433                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
434                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
435                    if name == 'Bond':
436                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
437                    elif name == 'Angle':
438                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
439                    elif name == 'Plane':
440                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
441                    elif name == 'Chiral':
442                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
443                elif name in ['Torsion','Rama']:
444                    coffDict = itemRest['Coeff']
445                    indx,ops,cofName,esd = itemRest[names[name]][id]
446                    dNames = []
447                    for ind in indx:
448                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
449                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
450                    if name == 'Torsion':
451                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
452                    else:
453                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
454                elif name == 'ChemComp':
455                    indx,factors,obs,esd = itemRest[names[name]][id]
456                    dNames = []
457                    for ind in indx:
458                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
459                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
460                        deriv = mul*factors
461                elif 'Texture' in name:
462                    deriv = []
463                    dNames = []
464                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
465                    hkl = np.array(hkl)
466                    if np.any(lasthkl-hkl):
467                        PH = np.array(hkl)
468                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
469                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
470                        lasthkl = copy.copy(hkl)                       
471                    if 'unit' in name:
472                        pass
473                    else:
474                        gam = float(pnames[3])
475                        psi = float(pnames[4])
476                        for SHname in ODFln:
477                            l,m,n = eval(SHname[1:])
478                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
479                            dNames += [str(pId)+'::'+SHname]
480                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
481                for dName,drv in zip(dNames,deriv):
482                    try:
483                        ind = varyList.index(dName)
484                        pDerv[ind][ip] += drv
485                    except ValueError:
486                        pass
487    return pDerv
488
489################################################################################
490##### Function & derivative calculations
491################################################################################       
492                   
493def GetAtomFXU(pfx,calcControls,parmDict):
494    'Needs a doc string'
495    Natoms = calcControls['Natoms'][pfx]
496    Tdata = Natoms*[' ',]
497    Mdata = np.zeros(Natoms)
498    IAdata = Natoms*[' ',]
499    Fdata = np.zeros(Natoms)
500    FFdata = []
501    BLdata = []
502    Xdata = np.zeros((3,Natoms))
503    dXdata = np.zeros((3,Natoms))
504    Uisodata = np.zeros(Natoms)
505    Uijdata = np.zeros((6,Natoms))
506    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
507        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
508        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
509        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
510        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
511    for iatm in range(Natoms):
512        for key in keys:
513            parm = pfx+key+str(iatm)
514            if parm in parmDict:
515                keys[key][iatm] = parmDict[parm]
516    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
517   
518def getFFvalues(FFtables,SQ):
519    'Needs a doc string'
520    FFvals = {}
521    for El in FFtables:
522        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
523    return FFvals
524   
525def getBLvalues(BLtables):
526    'Needs a doc string'
527    BLvals = {}
528    for El in BLtables:
529        BLvals[El] = BLtables[El][1][1]
530    return BLvals
531       
532def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
533    ''' Compute structure factors for all h,k,l for phase
534    puts the result, F^2, in each ref[8] in refList
535    input:
536   
537    :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
538    :param np.array G:      reciprocal metric tensor
539    :param str pfx:    phase id string
540    :param dict SGData: space group info. dictionary output from SpcGroup
541    :param dict calcControls:
542    :param dict ParmDict:
543
544    '''       
545    twopi = 2.0*np.pi
546    twopisq = 2.0*np.pi**2
547    phfx = pfx.split(':')[0]+hfx
548    ast = np.sqrt(np.diag(G))
549    Mast = twopisq*np.multiply.outer(ast,ast)
550    FFtables = calcControls['FFtables']
551    BLtables = calcControls['BLtables']
552    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
553    FF = np.zeros(len(Tdata))
554    if 'N' in parmDict[hfx+'Type']:
555        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
556    else:
557        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
558        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
559    maxPos = len(SGData['SGOps'])
560    Uij = np.array(G2lat.U6toUij(Uijdata))
561    bij = Mast*Uij.T
562    for refl in refList:
563        fbs = np.array([0,0])
564        H = refl[:3]
565        SQ = 1./(2.*refl[4])**2
566        SQfactor = 4.0*SQ*twopisq
567        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
568        if not len(refl[-1]):                #no form factors
569            if 'N' in parmDict[hfx+'Type']:
570                refl[-1] = getBLvalues(BLtables)
571            else:       #'X'
572                refl[-1] = getFFvalues(FFtables,SQ)
573        for i,El in enumerate(Tdata):
574            FF[i] = refl[-1][El]           
575        Uniq = refl[11]
576        phi = refl[12]
577        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
578        sinp = np.sin(phase)
579        cosp = np.cos(phase)
580        occ = Mdata*Fdata/len(Uniq)
581        biso = -SQfactor*Uisodata
582        Tiso = np.where(biso<1.,np.exp(biso),1.0)
583        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
584        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
585        Tcorr = Tiso*Tuij
586        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
587        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
588        if not SGData['SGInv']:
589            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
590            fbs = np.sum(np.sum(fb,axis=1),axis=1)
591        fasq = fas**2
592        fbsq = fbs**2        #imaginary
593        refl[9] = np.sum(fasq)+np.sum(fbsq)
594        refl[10] = atan2d(fbs[0],fas[0])
595   
596def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
597    'Needs a doc string'
598    twopi = 2.0*np.pi
599    twopisq = 2.0*np.pi**2
600    phfx = pfx.split(':')[0]+hfx
601    ast = np.sqrt(np.diag(G))
602    Mast = twopisq*np.multiply.outer(ast,ast)
603    FFtables = calcControls['FFtables']
604    BLtables = calcControls['BLtables']
605    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
606    FF = np.zeros(len(Tdata))
607    if 'N' in parmDict[hfx+'Type']:
608        FP = 0.
609        FPP = 0.
610    else:
611        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
612        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
613    maxPos = len(SGData['SGOps'])       
614    Uij = np.array(G2lat.U6toUij(Uijdata))
615    bij = Mast*Uij.T
616    dFdvDict = {}
617    dFdfr = np.zeros((len(refList),len(Mdata)))
618    dFdx = np.zeros((len(refList),len(Mdata),3))
619    dFdui = np.zeros((len(refList),len(Mdata)))
620    dFdua = np.zeros((len(refList),len(Mdata),6))
621    dFdbab = np.zeros((len(refList),2))
622    for iref,refl in enumerate(refList):
623        H = np.array(refl[:3])
624        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
625        SQfactor = 8.0*SQ*np.pi**2
626        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
627        Bab = parmDict[phfx+'BabA']*dBabdA
628        for i,El in enumerate(Tdata):           
629            FF[i] = refl[-1][El]           
630        Uniq = refl[11]
631        phi = refl[12]
632        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
633        sinp = np.sin(phase)
634        cosp = np.cos(phase)
635        occ = Mdata*Fdata/len(Uniq)
636        biso = -SQfactor*Uisodata
637        Tiso = np.where(biso<1.,np.exp(biso),1.0)
638        HbH = -np.inner(H,np.inner(bij,H))
639        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
640        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
641        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
642        Tcorr = Tiso*Tuij
643        fot = (FF+FP-Bab)*occ*Tcorr
644        fotp = FPP*occ*Tcorr
645        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
646        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
647       
648        fas = np.sum(np.sum(fa,axis=1),axis=1)
649        fbs = np.sum(np.sum(fb,axis=1),axis=1)
650        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
651        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
652        #sum below is over Uniq
653        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
654        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
655        dfadui = np.sum(-SQfactor*fa,axis=2)
656        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
657        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
658        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
659        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
660        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
661        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
662        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
663        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
664        if not SGData['SGInv']:
665            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
666            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
667            dfbdui = np.sum(-SQfactor*fb,axis=2)
668            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
669            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
670            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
671            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
672            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
673            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
674            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
675        #loop over atoms - each dict entry is list of derivatives for all the reflections
676    for i in range(len(Mdata)):     
677        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
678        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
679        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
680        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
681        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
682        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
683        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
684        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
685        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
686        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
687        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
688        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
689        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
690    return dFdvDict
691   
692def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
693    ''' Single crystal extinction function; puts correction in ref[13] and returns
694    corrections needed for derivatives
695    '''
696    ref[13] = 1.0
697    dervCor = 1.0
698    dervDict = {}
699    if calcControls[phfx+'EType'] != 'None':
700        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
701        if 'SXC' in parmDict[hfx+'Type']:
702            AV = 7.9406e5/parmDict[pfx+'Vol']**2
703            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
704            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
705        elif 'SNT' in parmDict[hfx+'Type']:
706            AV = 1.e7/parmDict[pfx+'Vol']**2
707            PL = 1./(4.*refl[4]**2)
708            P12 = 1.0
709        elif 'SNC' in parmDict[hfx+'Type']:
710            AV = 1.e7/parmDict[pfx+'Vol']**2
711            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
712            P12 = 1.0
713           
714        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
715        if 'Primary' in calcControls[phfx+'EType']:
716            PLZ *= 1.5
717        else:
718            PLZ *= calcControls[phfx+'Tbar']
719                       
720        if 'Primary' in calcControls[phfx+'EType']:
721            PSIG = parmDict[phfx+'Ep']
722        elif 'I & II' in calcControls[phfx+'EType']:
723            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
724        elif 'Type II' in calcControls[phfx+'EType']:
725            PSIG = parmDict[phfx+'Es']
726        else:       # 'Secondary Type I'
727            PSIG = parmDict[phfx+'Eg']/PL
728           
729        AG = 0.58+0.48*cos2T+0.24*cos2T**2
730        AL = 0.025+0.285*cos2T
731        BG = 0.02-0.025*cos2T
732        BL = 0.15-0.2*(0.75-cos2T)**2
733        if cos2T < 0.:
734            BL = -0.45*cos2T
735        CG = 2.
736        CL = 2.
737        PF = PLZ*PSIG
738       
739        if 'Gaussian' in calcControls[phfx+'EApprox']:
740            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
741            extCor = np.sqrt(PF4)
742            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
743        else:
744            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
745            extCor = np.sqrt(PF4)
746            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
747
748        dervCor = (1.+PF)*PF3
749        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
750            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
751        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
752            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
753        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
754            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
755               
756        ref[13] = 1./extCor
757    return dervCor,dervDict
758       
759   
760def Dict2Values(parmdict, varylist):
761    '''Use before call to leastsq to setup list of values for the parameters
762    in parmdict, as selected by key in varylist'''
763    return [parmdict[key] for key in varylist] 
764   
765def Values2Dict(parmdict, varylist, values):
766    ''' Use after call to leastsq to update the parameter dictionary with
767    values corresponding to keys in varylist'''
768    parmdict.update(zip(varylist,values))
769   
770def GetNewCellParms(parmDict,varyList):
771    'Needs a doc string'
772    newCellDict = {}
773    Anames = ['A'+str(i) for i in range(6)]
774    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
775    for item in varyList:
776        keys = item.split(':')
777        if keys[2] in Ddict:
778            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
779            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
780            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
781    return newCellDict          # is e.g. {'0::D11':A0+D11}
782   
783def ApplyXYZshifts(parmDict,varyList):
784    '''
785    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
786   
787    :param dict parmDict: parameter dictionary
788    :param list varyList: list of variables
789    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
790
791    '''
792    newAtomDict = {}
793    for item in parmDict:
794        if 'dA' in item:
795            parm = ''.join(item.split('d'))
796            parmDict[parm] += parmDict[item]
797            newAtomDict[item] = [parm,parmDict[parm]]
798    return newAtomDict
799   
800def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
801    'Spherical harmonics texture'
802    IFCoup = 'Bragg' in calcControls[hfx+'instType']
803    odfCor = 1.0
804    H = refl[:3]
805    cell = G2lat.Gmat2cell(g)
806    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
807    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
808    phi,beta = G2lat.CrsAng(H,cell,SGData)
809    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
810    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
811    for item in SHnames:
812        L,M,N = eval(item.strip('C'))
813        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
814        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
815        Lnorm = G2lat.Lnorm(L)
816        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
817    return odfCor
818   
819def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
820    'Spherical harmonics texture derivatives'
821    FORPI = 4.0*np.pi
822    IFCoup = 'Bragg' in calcControls[hfx+'instType']
823    odfCor = 1.0
824    dFdODF = {}
825    dFdSA = [0,0,0]
826    H = refl[:3]
827    cell = G2lat.Gmat2cell(g)
828    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
829    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
830    phi,beta = G2lat.CrsAng(H,cell,SGData)
831    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
832    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
833    for item in SHnames:
834        L,M,N = eval(item.strip('C'))
835        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
836        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
837        Lnorm = G2lat.Lnorm(L)
838        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
839        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
840        for i in range(3):
841            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
842    return odfCor,dFdODF,dFdSA
843   
844def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
845    'spherical harmonics preferred orientation (cylindrical symmetry only)'
846    odfCor = 1.0
847    H = refl[:3]
848    cell = G2lat.Gmat2cell(g)
849    Sangl = [0.,0.,0.]
850    if 'Bragg' in calcControls[hfx+'instType']:
851        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
852        IFCoup = True
853    else:
854        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
855        IFCoup = False
856    phi,beta = G2lat.CrsAng(H,cell,SGData)
857    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
858    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
859    for item in SHnames:
860        L,N = eval(item.strip('C'))
861        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
862        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
863    return odfCor
864   
865def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
866    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
867    FORPI = 12.5663706143592
868    odfCor = 1.0
869    dFdODF = {}
870    H = refl[:3]
871    cell = G2lat.Gmat2cell(g)
872    Sangl = [0.,0.,0.]
873    if 'Bragg' in calcControls[hfx+'instType']:
874        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
875        IFCoup = True
876    else:
877        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
878        IFCoup = False
879    phi,beta = G2lat.CrsAng(H,cell,SGData)
880    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
881    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
882    for item in SHnames:
883        L,N = eval(item.strip('C'))
884        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
885        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
886        dFdODF[phfx+item] = Kcsl*Lnorm
887    return odfCor,dFdODF
888   
889def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
890    'Needs a doc string'
891    POcorr = 1.0
892    if calcControls[phfx+'poType'] == 'MD':
893        MD = parmDict[phfx+'MD']
894        if MD != 1.0:
895            MDAxis = calcControls[phfx+'MDAxis']
896            sumMD = 0
897            for H in refl[11]:           
898                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
899                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
900                sumMD += A**3
901            POcorr = sumMD/len(refl[11])
902    else:   #spherical harmonics
903        if calcControls[phfx+'SHord']:
904            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
905    return POcorr
906   
907def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
908    'Needs a doc string'
909    POcorr = 1.0
910    POderv = {}
911    if calcControls[phfx+'poType'] == 'MD':
912        MD = parmDict[phfx+'MD']
913        MDAxis = calcControls[phfx+'MDAxis']
914        sumMD = 0
915        sumdMD = 0
916        for H in refl[11]:           
917            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
918            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
919            sumMD += A**3
920            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
921        POcorr = sumMD/len(refl[11])
922        POderv[phfx+'MD'] = sumdMD/len(refl[11])
923    else:   #spherical harmonics
924        if calcControls[phfx+'SHord']:
925            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
926    return POcorr,POderv
927   
928def GetAbsorb(refl,hfx,calcControls,parmDict):
929    'Needs a doc string'
930    if 'Debye' in calcControls[hfx+'instType']:
931        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
932    else:
933        return 1.0
934   
935def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
936    'Needs a doc string'
937    if 'Debye' in calcControls[hfx+'instType']:
938        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
939    else:
940        return 0.0
941   
942def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
943    'Needs a doc string'
944    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
945    if 'X' in parmDict[hfx+'Type']:
946        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
947    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
948    if pfx+'SHorder' in parmDict:
949        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
950    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
951    refl[13] = Icorr       
952   
953def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
954    'Needs a doc string'
955    dIdsh = 1./parmDict[hfx+'Scale']
956    dIdsp = 1./parmDict[phfx+'Scale']
957    if 'X' in parmDict[hfx+'Type']:
958        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
959        dIdPola /= pola
960    else:       #'N'
961        dIdPola = 0.0
962    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
963    for iPO in dIdPO:
964        dIdPO[iPO] /= POcorr
965    dFdODF = {}
966    dFdSA = [0,0,0]
967    if pfx+'SHorder' in parmDict:
968        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
969        for iSH in dFdODF:
970            dFdODF[iSH] /= odfCor
971        for i in range(3):
972            dFdSA[i] /= odfCor
973    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
974    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
975       
976def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
977    'Needs a doc string'
978    costh = cosd(refl[5]/2.)
979    #crystallite size
980    if calcControls[phfx+'SizeType'] == 'isotropic':
981        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
982    elif calcControls[phfx+'SizeType'] == 'uniaxial':
983        H = np.array(refl[:3])
984        P = np.array(calcControls[phfx+'SizeAxis'])
985        cosP,sinP = G2lat.CosSinAngle(H,P,G)
986        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
987        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
988    else:           #ellipsoidal crystallites
989        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
990        H = np.array(refl[:3])
991        lenR = G2pwd.ellipseSize(H,Sij,GB)
992        Sgam = 1.8*wave/(np.pi*costh*lenR)
993    #microstrain               
994    if calcControls[phfx+'MustrainType'] == 'isotropic':
995        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
996    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
997        H = np.array(refl[:3])
998        P = np.array(calcControls[phfx+'MustrainAxis'])
999        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1000        Si = parmDict[phfx+'Mustrain;i']
1001        Sa = parmDict[phfx+'Mustrain;a']
1002        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1003    else:       #generalized - P.W. Stephens model
1004        pwrs = calcControls[phfx+'MuPwrs']
1005        sum = 0
1006        for i,pwr in enumerate(pwrs):
1007            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1008        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1009    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1010    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1011    sig /= ateln2
1012    return sig,gam
1013       
1014def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1015    'Needs a doc string'
1016    gamDict = {}
1017    sigDict = {}
1018    costh = cosd(refl[5]/2.)
1019    tanth = tand(refl[5]/2.)
1020    #crystallite size derivatives
1021    if calcControls[phfx+'SizeType'] == 'isotropic':
1022        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1023        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1024        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1025    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1026        H = np.array(refl[:3])
1027        P = np.array(calcControls[phfx+'SizeAxis'])
1028        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1029        Si = parmDict[phfx+'Size;i']
1030        Sa = parmDict[phfx+'Size;a']
1031        gami = (1.8*wave/np.pi)/(Si*Sa)
1032        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1033        Sgam = gami*sqtrm
1034        gam = Sgam/costh
1035        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
1036        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
1037        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1038        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1039        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1040        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1041    else:           #ellipsoidal crystallites
1042        const = 1.8*wave/(np.pi*costh)
1043        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1044        H = np.array(refl[:3])
1045        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1046        Sgam = 1.8*wave/(np.pi*costh*lenR)
1047        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1048            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1049            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1050    gamDict[phfx+'Size;mx'] = Sgam
1051    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1052           
1053    #microstrain derivatives               
1054    if calcControls[phfx+'MustrainType'] == 'isotropic':
1055        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1056        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1057        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1058    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1059        H = np.array(refl[:3])
1060        P = np.array(calcControls[phfx+'MustrainAxis'])
1061        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1062        Si = parmDict[phfx+'Mustrain;i']
1063        Sa = parmDict[phfx+'Mustrain;a']
1064        gami = 0.018*Si*Sa*tanth/np.pi
1065        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1066        Mgam = gami/sqtrm
1067        dsi = -gami*Si*cosP**2/sqtrm**3
1068        dsa = -gami*Sa*sinP**2/sqtrm**3
1069        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1070        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1071        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1072        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1073    else:       #generalized - P.W. Stephens model
1074        pwrs = calcControls[phfx+'MuPwrs']
1075        const = 0.018*refl[4]**2*tanth
1076        sum = 0
1077        for i,pwr in enumerate(pwrs):
1078            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1079            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1080            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1081            sigDict[phfx+'Mustrain:'+str(i)] = \
1082                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1083        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1084        for i in range(len(pwrs)):
1085            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1086    gamDict[phfx+'Mustrain;mx'] = Mgam
1087    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1088    return sigDict,gamDict
1089       
1090def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1091    'Needs a doc string'
1092    h,k,l = refl[:3]
1093    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1094    d = np.sqrt(dsq)
1095
1096    refl[4] = d
1097    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1098    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1099    if 'Bragg' in calcControls[hfx+'instType']:
1100        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1101            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1102    else:               #Debye-Scherrer - simple but maybe not right
1103        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1104    return pos
1105
1106def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1107    'Needs a doc string'
1108    dpr = 180./np.pi
1109    h,k,l = refl[:3]
1110    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1111    dst = np.sqrt(dstsq)
1112    pos = refl[5]-parmDict[hfx+'Zero']
1113    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1114    dpdw = const*dst
1115    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1116    dpdA *= const*wave/(2.0*dst)
1117    dpdZ = 1.0
1118    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1119    if 'Bragg' in calcControls[hfx+'instType']:
1120        dpdSh = -4.*const*cosd(pos/2.0)
1121        dpdTr = -const*sind(pos)*100.0
1122        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1123    else:               #Debye-Scherrer - simple but maybe not right
1124        dpdXd = -const*cosd(pos)
1125        dpdYd = -const*sind(pos)
1126        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1127           
1128def GetHStrainShift(refl,SGData,phfx,parmDict):
1129    'Needs a doc string'
1130    laue = SGData['SGLaue']
1131    uniq = SGData['SGUniq']
1132    h,k,l = refl[:3]
1133    if laue in ['m3','m3m']:
1134        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1135            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1136    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1137        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1138    elif laue in ['3R','3mR']:
1139        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1140    elif laue in ['4/m','4/mmm']:
1141        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1142    elif laue in ['mmm']:
1143        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1144    elif laue in ['2/m']:
1145        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1146        if uniq == 'a':
1147            Dij += parmDict[phfx+'D23']*k*l
1148        elif uniq == 'b':
1149            Dij += parmDict[phfx+'D13']*h*l
1150        elif uniq == 'c':
1151            Dij += parmDict[phfx+'D12']*h*k
1152    else:
1153        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1154            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1155    return -Dij*refl[4]**2*tand(refl[5]/2.0)
1156           
1157def GetHStrainShiftDerv(refl,SGData,phfx):
1158    'Needs a doc string'
1159    laue = SGData['SGLaue']
1160    uniq = SGData['SGUniq']
1161    h,k,l = refl[:3]
1162    if laue in ['m3','m3m']:
1163        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1164            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1165    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1166        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1167    elif laue in ['3R','3mR']:
1168        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1169    elif laue in ['4/m','4/mmm']:
1170        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1171    elif laue in ['mmm']:
1172        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1173    elif laue in ['2/m']:
1174        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1175        if uniq == 'a':
1176            dDijDict[phfx+'D23'] = k*l
1177        elif uniq == 'b':
1178            dDijDict[phfx+'D13'] = h*l
1179        elif uniq == 'c':
1180            dDijDict[phfx+'D12'] = h*k
1181            names.append()
1182    else:
1183        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1184            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1185    for item in dDijDict:
1186        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
1187    return dDijDict
1188   
1189def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1190    'Needs a doc string'
1191    histoList = Histograms.keys()
1192    histoList.sort()
1193    for histogram in histoList:
1194        if 'PWDR' in histogram[:4]:
1195            Histogram = Histograms[histogram]
1196            hId = Histogram['hId']
1197            hfx = ':%d:'%(hId)
1198            Limits = calcControls[hfx+'Limits']
1199            shl = max(parmDict[hfx+'SH/L'],0.0005)
1200            Ka2 = False
1201            kRatio = 0.0
1202            if hfx+'Lam1' in parmDict.keys():
1203                Ka2 = True
1204                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1205                kRatio = parmDict[hfx+'I(L2)/I(L1)']
1206            x,y,w,yc,yb,yd = Histogram['Data']
1207            xB = np.searchsorted(x,Limits[0])
1208            xF = np.searchsorted(x,Limits[1])
1209            ymb = np.array(y-yb)
1210            ymb = np.where(ymb,ymb,1.0)
1211            ycmb = np.array(yc-yb)
1212            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1213            refLists = Histogram['Reflection Lists']
1214            for phase in refLists:
1215                Phase = Phases[phase]
1216                pId = Phase['pId']
1217                phfx = '%d:%d:'%(pId,hId)
1218                refList = refLists[phase]
1219                sumFo = 0.0
1220                sumdF = 0.0
1221                sumFosq = 0.0
1222                sumdFsq = 0.0
1223                for refl in refList:
1224                    if 'C' in calcControls[hfx+'histType']:
1225                        yp = np.zeros_like(yb)
1226                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1227                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1228                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1229                        iFin2 = iFin
1230                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1231                        if Ka2:
1232                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1233                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1234                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1235                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1236                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1237                                continue
1238                            elif not iBeg2-iFin2:     #peak above high limit - done
1239                                break
1240                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
1241                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
1242                    elif 'T' in calcControls[hfx+'histType']:
1243                        print 'TOF Undefined at present'
1244                        raise Exception    #no TOF yet
1245                    Fo = np.sqrt(np.abs(refl[8]))
1246                    Fc = np.sqrt(np.abs(refl[9]))
1247                    sumFo += Fo
1248                    sumFosq += refl[8]**2
1249                    sumdF += np.abs(Fo-Fc)
1250                    sumdFsq += (refl[8]-refl[9])**2
1251                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1252                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1253                Histogram[phfx+'Nref'] = len(refList)
1254               
1255def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1256    'Needs a doc string'
1257   
1258    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1259        U = parmDict[hfx+'U']
1260        V = parmDict[hfx+'V']
1261        W = parmDict[hfx+'W']
1262        X = parmDict[hfx+'X']
1263        Y = parmDict[hfx+'Y']
1264        tanPos = tand(refl[5]/2.0)
1265        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1266        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1267        sig = max(0.001,sig)
1268        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1269        gam = max(0.001,gam)
1270        return sig,gam
1271               
1272    hId = Histogram['hId']
1273    hfx = ':%d:'%(hId)
1274    bakType = calcControls[hfx+'bakType']
1275    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1276    yc = np.zeros_like(yb)
1277       
1278    if 'C' in calcControls[hfx+'histType']:   
1279        shl = max(parmDict[hfx+'SH/L'],0.002)
1280        Ka2 = False
1281        if hfx+'Lam1' in parmDict.keys():
1282            wave = parmDict[hfx+'Lam1']
1283            Ka2 = True
1284            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1285            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1286        else:
1287            wave = parmDict[hfx+'Lam']
1288    else:
1289        print 'TOF Undefined at present'
1290        raise ValueError
1291    for phase in Histogram['Reflection Lists']:
1292        refList = Histogram['Reflection Lists'][phase]
1293        Phase = Phases[phase]
1294        pId = Phase['pId']
1295        pfx = '%d::'%(pId)
1296        phfx = '%d:%d:'%(pId,hId)
1297        hfx = ':%d:'%(hId)
1298        SGData = Phase['General']['SGData']
1299        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1300        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1301        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1302        Vst = np.sqrt(nl.det(G))    #V*
1303        if not Phase['General'].get('doPawley'):
1304            time0 = time.time()
1305            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1306#            print 'sf calc time: %.3fs'%(time.time()-time0)
1307        time0 = time.time()
1308        for refl in refList:
1309            if 'C' in calcControls[hfx+'histType']:
1310                h,k,l = refl[:3]
1311                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1312                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1313                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1314                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1315                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
1316                refl[13] *= Vst*Lorenz
1317                if Phase['General'].get('doPawley'):
1318                    try:
1319                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1320                        refl[9] = parmDict[pInd]
1321                    except KeyError:
1322#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1323                        continue
1324                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1325                iBeg = np.searchsorted(x,refl[5]-fmin)
1326                iFin = np.searchsorted(x,refl[5]+fmax)
1327                if not iBeg+iFin:       #peak below low limit - skip peak
1328                    continue
1329                elif not iBeg-iFin:     #peak above high limit - done
1330                    break
1331                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1332                if Ka2:
1333                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1334                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1335                    iBeg = np.searchsorted(x,pos2-fmin)
1336                    iFin = np.searchsorted(x,pos2+fmax)
1337                    if not iBeg+iFin:       #peak below low limit - skip peak
1338                        continue
1339                    elif not iBeg-iFin:     #peak above high limit - done
1340                        return yc,yb
1341                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1342            elif 'T' in calcControls[hfx+'histType']:
1343                print 'TOF Undefined at present'
1344                raise Exception    #no TOF yet
1345#        print 'profile calc time: %.3fs'%(time.time()-time0)
1346    return yc,yb
1347   
1348def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1349    'Needs a doc string'
1350   
1351    def cellVaryDerv(pfx,SGData,dpdA): 
1352        if SGData['SGLaue'] in ['-1',]:
1353            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1354                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1355        elif SGData['SGLaue'] in ['2/m',]:
1356            if SGData['SGUniq'] == 'a':
1357                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1358            elif SGData['SGUniq'] == 'b':
1359                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1360            else:
1361                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1362        elif SGData['SGLaue'] in ['mmm',]:
1363            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1364        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1365            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1366        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1367            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1368        elif SGData['SGLaue'] in ['3R', '3mR']:
1369            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1370        elif SGData['SGLaue'] in ['m3m','m3']:
1371            return [[pfx+'A0',dpdA[0]]]
1372           
1373    # create a list of dependent variables and set up a dictionary to hold their derivatives
1374    dependentVars = G2mv.GetDependentVars()
1375    depDerivDict = {}
1376    for j in dependentVars:
1377        depDerivDict[j] = np.zeros(shape=(len(x)))
1378    #print 'dependent vars',dependentVars
1379    lenX = len(x)               
1380    hId = Histogram['hId']
1381    hfx = ':%d:'%(hId)
1382    bakType = calcControls[hfx+'bakType']
1383    dMdv = np.zeros(shape=(len(varylist),len(x)))
1384    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1385    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1386        bBpos =varylist.index(hfx+'Back:0')
1387        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1388    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1389    for name in varylist:
1390        if 'Debye' in name:
1391            id = int(name.split(':')[-1])
1392            parm = name[:int(name.rindex(':'))]
1393            ip = names.index(parm)
1394            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1395    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1396    for name in varylist:
1397        if 'BkPk' in name:
1398            id = int(name.split(':')[-1])
1399            parm = name[:int(name.rindex(':'))]
1400            ip = names.index(parm)
1401            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1402    cw = np.diff(x)
1403    cw = np.append(cw,cw[-1])
1404    if 'C' in calcControls[hfx+'histType']:   
1405        shl = max(parmDict[hfx+'SH/L'],0.002)
1406        Ka2 = False
1407        if hfx+'Lam1' in parmDict.keys():
1408            wave = parmDict[hfx+'Lam1']
1409            Ka2 = True
1410            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1411            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1412        else:
1413            wave = parmDict[hfx+'Lam']
1414    else:
1415        print 'TOF Undefined at present'
1416        raise ValueError
1417    for phase in Histogram['Reflection Lists']:
1418        refList = Histogram['Reflection Lists'][phase]
1419        Phase = Phases[phase]
1420        SGData = Phase['General']['SGData']
1421        pId = Phase['pId']
1422        pfx = '%d::'%(pId)
1423        phfx = '%d:%d:'%(pId,hId)
1424        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1425        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1426        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1427        if not Phase['General'].get('doPawley'):
1428            time0 = time.time()
1429            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1430#            print 'sf-derv time %.3fs'%(time.time()-time0)
1431            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1432        time0 = time.time()
1433        for iref,refl in enumerate(refList):
1434            if 'C' in calcControls[hfx+'histType']:        #CW powder
1435                h,k,l = refl[:3]
1436                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1437                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1438                iBeg = np.searchsorted(x,refl[5]-fmin)
1439                iFin = np.searchsorted(x,refl[5]+fmax)
1440                if not iBeg+iFin:       #peak below low limit - skip peak
1441                    continue
1442                elif not iBeg-iFin:     #peak above high limit - done
1443                    break
1444                pos = refl[5]
1445                tanth = tand(pos/2.0)
1446                costh = cosd(pos/2.0)
1447                lenBF = iFin-iBeg
1448                dMdpk = np.zeros(shape=(6,lenBF))
1449                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1450                for i in range(5):
1451                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
1452                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1453                if Ka2:
1454                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1455                    iBeg2 = np.searchsorted(x,pos2-fmin)
1456                    iFin2 = np.searchsorted(x,pos2+fmax)
1457                    if iBeg2-iFin2:
1458                        lenBF2 = iFin2-iBeg2
1459                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1460                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
1461                        for i in range(5):
1462                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
1463                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
1464                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1465                if Phase['General'].get('doPawley'):
1466                    dMdpw = np.zeros(len(x))
1467                    try:
1468                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1469                        idx = varylist.index(pIdx)
1470                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1471                        if Ka2:
1472                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1473                        dMdv[idx] = dMdpw
1474                    except: # ValueError:
1475                        pass
1476                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1477                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1478                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1479                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1480                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1481                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1482                    hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
1483                for name in names:
1484                    item = names[name]
1485                    if name in varylist:
1486                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1487                        if Ka2:
1488                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1489                    elif name in dependentVars:
1490                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1491                        if Ka2:
1492                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1493                for iPO in dIdPO:
1494                    if iPO in varylist:
1495                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1496                        if Ka2:
1497                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1498                    elif iPO in dependentVars:
1499                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1500                        if Ka2:
1501                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1502                for i,name in enumerate(['omega','chi','phi']):
1503                    aname = pfx+'SH '+name
1504                    if aname in varylist:
1505                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1506                        if Ka2:
1507                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1508                    elif aname in dependentVars:
1509                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1510                        if Ka2:
1511                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1512                for iSH in dFdODF:
1513                    if iSH in varylist:
1514                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1515                        if Ka2:
1516                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1517                    elif iSH in dependentVars:
1518                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1519                        if Ka2:
1520                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1521                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1522                for name,dpdA in cellDervNames:
1523                    if name in varylist:
1524                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1525                        if Ka2:
1526                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1527                    elif name in dependentVars:
1528                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1529                        if Ka2:
1530                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1531                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1532                for name in dDijDict:
1533                    if name in varylist:
1534                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1535                        if Ka2:
1536                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1537                    elif name in dependentVars:
1538                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1539                        if Ka2:
1540                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1541                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1542                for name in gamDict:
1543                    if name in varylist:
1544                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1545                        if Ka2:
1546                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1547                    elif name in dependentVars:
1548                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1549                        if Ka2:
1550                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1551                for name in sigDict:
1552                    if name in varylist:
1553                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1554                        if Ka2:
1555                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1556                    elif name in dependentVars:
1557                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1558                        if Ka2:
1559                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1560                for name in ['BabA','BabU']:
1561                    if phfx+name in varylist:
1562                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1563                        if Ka2:
1564                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
1565                    elif phfx+name in dependentVars:                   
1566                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1567                        if Ka2:
1568                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
1569            elif 'T' in calcControls[hfx+'histType']:
1570                print 'TOF Undefined at present'
1571                raise Exception    #no TOF yet
1572            if not Phase['General'].get('doPawley'):
1573                #do atom derivatives -  for RB,F,X & U so far             
1574                corr = dervDict['int']/refl[9]
1575                if Ka2:
1576                    corr2 = dervDict2['int']/refl[9]
1577                for name in varylist+dependentVars:
1578                    if '::RBV;' in name:
1579                        pass
1580                    else:
1581                        try:
1582                            aname = name.split(pfx)[1][:2]
1583                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1584                        except IndexError:
1585                            continue
1586                    if name in varylist:
1587                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1588                        if Ka2:
1589                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1590                    elif name in dependentVars:
1591                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1592                        if Ka2:
1593                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1594#        print 'profile derv time: %.3fs'%(time.time()-time0)
1595    # now process derivatives in constraints
1596    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1597    return dMdv
1598
1599def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1600    'Needs a doc string'
1601    parmDict.update(zip(varylist,values))
1602    G2mv.Dict2Map(parmDict,varylist)
1603    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1604    nvar = len(varylist)
1605    dMdv = np.empty(0)
1606    histoList = Histograms.keys()
1607    histoList.sort()
1608    for histogram in histoList:
1609        if 'PWDR' in histogram[:4]:
1610            Histogram = Histograms[histogram]
1611            hId = Histogram['hId']
1612            hfx = ':%d:'%(hId)
1613            wtFactor = calcControls[hfx+'wtFactor']
1614            Limits = calcControls[hfx+'Limits']
1615            x,y,w,yc,yb,yd = Histogram['Data']
1616            W = wtFactor*w
1617            xB = np.searchsorted(x,Limits[0])
1618            xF = np.searchsorted(x,Limits[1])
1619            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1620                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1621        elif 'HKLF' in histogram[:4]:
1622            Histogram = Histograms[histogram]
1623            nobs = Histogram['Nobs']
1624            phase = Histogram['Reflection Lists']
1625            Phase = Phases[phase]
1626            hId = Histogram['hId']
1627            hfx = ':%d:'%(hId)
1628            wtFactor = calcControls[hfx+'wtFactor']
1629            pfx = '%d::'%(Phase['pId'])
1630            phfx = '%d:%d:'%(Phase['pId'],hId)
1631            SGData = Phase['General']['SGData']
1632            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1633            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1634            refList = Histogram['Data']
1635            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1636            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1637            dMdvh = np.zeros((len(varylist),len(refList)))
1638            if calcControls['F**2']:
1639                for iref,ref in enumerate(refList):
1640                    if ref[6] > 0:
1641                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1642                        w = 1.0/ref[6]
1643                        if w*ref[5] >= calcControls['minF/sig']:
1644                            for j,var in enumerate(varylist):
1645                                if var in dFdvDict:
1646                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1647                            if phfx+'Scale' in varylist:
1648                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1649                            for item in ['Ep','Es','Eg']:
1650                                if phfx+item in varylist:
1651                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1652                            for item in ['BabA','BabU']:
1653                                if phfx+item in varylist:
1654                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1655            else:
1656                for iref,ref in enumerate(refList):
1657                    if ref[5] > 0.:
1658                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1659                        Fo = np.sqrt(ref[5])
1660                        Fc = np.sqrt(ref[7])
1661                        w = 1.0/ref[6]
1662                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1663                            for j,var in enumerate(varylist):
1664                                if var in dFdvDict:
1665                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1666                            if phfx+'Scale' in varylist:
1667                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor                           
1668                            for item in ['Ep','Es','Eg']:
1669                                if phfx+item in varylist:
1670                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1671                            for item in ['BabA','BabU']:
1672                                if phfx+item in varylist:
1673                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1674        else:
1675            continue        #skip non-histogram entries
1676        if len(dMdv):
1677            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1678        else:
1679            dMdv = np.sqrt(wtFactor)*dMdvh
1680           
1681    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1682    if np.any(pVals):
1683        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1684        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1685       
1686    return dMdv
1687
1688def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1689    'Needs a doc string'
1690    parmDict.update(zip(varylist,values))
1691    G2mv.Dict2Map(parmDict,varylist)
1692    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1693    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1694    nvar = len(varylist)
1695    Hess = np.empty(0)
1696    histoList = Histograms.keys()
1697    histoList.sort()
1698    for histogram in histoList:
1699        if 'PWDR' in histogram[:4]:
1700            Histogram = Histograms[histogram]
1701            hId = Histogram['hId']
1702            hfx = ':%d:'%(hId)
1703            wtFactor = calcControls[hfx+'wtFactor']
1704            Limits = calcControls[hfx+'Limits']
1705            x,y,w,yc,yb,yd = Histogram['Data']
1706            W = wtFactor*w
1707            dy = y-yc
1708            xB = np.searchsorted(x,Limits[0])
1709            xF = np.searchsorted(x,Limits[1])
1710            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1711                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1712            Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
1713            Dy = dy[xB:xF][np.newaxis,:]
1714            dMdvh *= Wt
1715            if dlg:
1716                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1717            if len(Hess):
1718                Hess += np.inner(dMdvh,dMdvh)
1719                dMdvh *= Wt*Dy
1720                Vec += np.sum(dMdvh,axis=1)
1721            else:
1722                Hess = np.inner(dMdvh,dMdvh)
1723                dMdvh *= Wt*Dy
1724                Vec = np.sum(dMdvh,axis=1)
1725        elif 'HKLF' in histogram[:4]:
1726            Histogram = Histograms[histogram]
1727            nobs = Histogram['Nobs']
1728            phase = Histogram['Reflection Lists']
1729            Phase = Phases[phase]
1730            hId = Histogram['hId']
1731            hfx = ':%d:'%(hId)
1732            wtFactor = calcControls[hfx+'wtFactor']
1733            pfx = '%d::'%(Phase['pId'])
1734            phfx = '%d:%d:'%(Phase['pId'],hId)
1735            SGData = Phase['General']['SGData']
1736            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1737            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1738            refList = Histogram['Data']
1739            time0 = time.time()
1740            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
1741            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1742            dMdvh = np.zeros((len(varylist),len(refList)))
1743            wdf = np.zeros(len(refList))
1744            if calcControls['F**2']:
1745                for iref,ref in enumerate(refList):
1746                    if ref[6] > 0:
1747                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1748                        w =  1.0/ref[6]
1749                        if w*ref[5] >= calcControls['minF/sig']:
1750                            wdf[iref] = w*(ref[5]-ref[7])
1751                            for j,var in enumerate(varylist):
1752                                if var in dFdvDict:
1753                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1754                            if phfx+'Scale' in varylist:
1755                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1756                            for item in ['Ep','Es','Eg']:
1757                                if phfx+item in varylist:
1758                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1759                            for item in ['BabA','BabU']:
1760                                if phfx+item in varylist:
1761                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1762            else:
1763                for iref,ref in enumerate(refList):
1764                    if ref[5] > 0.:
1765                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1766                        Fo = np.sqrt(ref[5])
1767                        Fc = np.sqrt(ref[7])
1768                        w = 1.0/ref[6]
1769                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1770                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1771                            for j,var in enumerate(varylist):
1772                                if var in dFdvDict:
1773                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1774                            if phfx+'Scale' in varylist:
1775                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor                           
1776                            for item in ['Ep','Es','Eg']:
1777                                if phfx+item in varylist:
1778                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1779                            for item in ['BabA','BabU']:
1780                                if phfx+item in varylist:
1781                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1782                       
1783            if dlg:
1784                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1785            if len(Hess):
1786                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1787                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1788            else:
1789                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1790                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1791        else:
1792            continue        #skip non-histogram entries
1793    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1794    if np.any(pVals):
1795        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1796        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1797        Hess += np.inner(dpdv*pWt,dpdv)
1798    return Vec,Hess
1799
1800def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1801    'Needs a doc string'
1802    parmDict.update(zip(varylist,values))
1803    Values2Dict(parmDict, varylist, values)
1804    G2mv.Dict2Map(parmDict,varylist)
1805    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1806    M = np.empty(0)
1807    SumwYo = 0
1808    Nobs = 0
1809    ApplyRBModels(parmDict,Phases,rigidbodyDict)
1810    histoList = Histograms.keys()
1811    histoList.sort()
1812    for histogram in histoList:
1813        if 'PWDR' in histogram[:4]:
1814            Histogram = Histograms[histogram]
1815            hId = Histogram['hId']
1816            hfx = ':%d:'%(hId)
1817            wtFactor = calcControls[hfx+'wtFactor']
1818            Limits = calcControls[hfx+'Limits']
1819            x,y,w,yc,yb,yd = Histogram['Data']
1820            W = wtFactor*w
1821            yc *= 0.0                           #zero full calcd profiles
1822            yb *= 0.0
1823            yd *= 0.0
1824            xB = np.searchsorted(x,Limits[0])
1825            xF = np.searchsorted(x,Limits[1])
1826            Histogram['Nobs'] = xF-xB
1827            Nobs += Histogram['Nobs']
1828            Histogram['sumwYo'] = np.sum(W[xB:xF]*y[xB:xF]**2)
1829            SumwYo += Histogram['sumwYo']
1830            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1831                varylist,Histogram,Phases,calcControls,pawleyLookup)
1832            yc[xB:xF] += yb[xB:xF]
1833            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1834            Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
1835            wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
1836            Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
1837            if dlg:
1838                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1839            M = np.concatenate((M,wdy))
1840#end of PWDR processing
1841        elif 'HKLF' in histogram[:4]:
1842            Histogram = Histograms[histogram]
1843            phase = Histogram['Reflection Lists']
1844            Phase = Phases[phase]
1845            hId = Histogram['hId']
1846            hfx = ':%d:'%(hId)
1847            wtFactor = calcControls[hfx+'wtFactor']
1848            pfx = '%d::'%(Phase['pId'])
1849            phfx = '%d:%d:'%(Phase['pId'],hId)
1850            SGData = Phase['General']['SGData']
1851            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1852            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1853            refList = Histogram['Data']
1854            time0 = time.time()
1855            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1856#            print 'sf-calc time: %.3f'%(time.time()-time0)
1857            df = np.zeros(len(refList))
1858            sumwYo = 0
1859            sumFo = 0
1860            sumFo2 = 0
1861            sumdF = 0
1862            sumdF2 = 0
1863            nobs = 0
1864            if calcControls['F**2']:
1865                for i,ref in enumerate(refList):
1866                    if ref[6] > 0:
1867                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1868                        w = 1.0/ref[6]
1869                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1870                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1871                        ref[8] = ref[5]/parmDict[phfx+'Scale']
1872                        if w*ref[5] >= calcControls['minF/sig']:
1873                            sumFo2 += ref[5]
1874                            Fo = np.sqrt(ref[5])
1875                            sumFo += Fo
1876                            sumFo2 += ref[5]
1877                            sumdF += abs(Fo-np.sqrt(ref[7]))
1878                            sumdF2 += abs(ref[5]-ref[7])
1879                            nobs += 1
1880                            df[i] = -w*(ref[5]-ref[7])
1881                            sumwYo += (w*ref[5])**2
1882            else:
1883                for i,ref in enumerate(refList):
1884                    if ref[5] > 0.:
1885                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1886                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1887                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1888                        Fo = np.sqrt(ref[5])
1889                        Fc = np.sqrt(ref[7])
1890                        w = 2.0*Fo/ref[6]
1891                        if w*Fo >= calcControls['minF/sig']:
1892                            sumFo += Fo
1893                            sumFo2 += ref[5]
1894                            sumdF += abs(Fo-Fc)
1895                            sumdF2 += abs(ref[5]-ref[7])
1896                            nobs += 1
1897                            df[i] = -w*(Fo-Fc)
1898                            sumwYo += (w*Fo)**2
1899            Histogram['Nobs'] = nobs
1900            Histogram['sumwYo'] = sumwYo
1901            SumwYo += sumwYo
1902            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
1903            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
1904            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
1905            Histogram[phfx+'Nref'] = nobs
1906            Nobs += nobs
1907            if dlg:
1908                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1909            M = np.concatenate((M,wtFactor*df))
1910# end of HKLF processing
1911    Histograms['sumwYo'] = SumwYo
1912    Histograms['Nobs'] = Nobs
1913    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
1914    if dlg:
1915        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
1916        if not GoOn:
1917            parmDict['saved values'] = values
1918            dlg.Destroy()
1919            raise Exception         #Abort!!
1920    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1921    if np.any(pVals):
1922        pSum = np.sum(pWt*pVals**2)
1923        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
1924        Nobs += len(pVals)
1925        M = np.concatenate((M,np.sqrt(pWt)*pVals))
1926    return M
1927                       
Note: See TracBrowser for help on using the repository browser.