source: trunk/GSASIIstrMath.py @ 1884

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

remove some unused imports
add merohedral/pseudomerohedral Twin Laws to G2ddataGUI and G2strIO (not in G2strmath yet).
allow ReImport? atoms to fill otherwise empty Atom List
clarify HKL importers as Shelx HKL 4 & HKL 5 files.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 141.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-06-09 21:02:06 +0000 (Tue, 09 Jun 2015) $
8# $Author: vondreele $
9# $Revision: 1884 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1884 2015-06-09 21:02:06Z vondreele $
12########### SVN repository information ###################
13import time
14import copy
15import numpy as np
16import numpy.ma as ma
17import numpy.linalg as nl
18import scipy.optimize as so
19import scipy.stats as st
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 1884 $")
22import GSASIIElem as G2el
23import GSASIIlattice as G2lat
24import GSASIIspc as G2spc
25import GSASIIpwd as G2pwd
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28import GSASIIobj as G2obj
29
30sind = lambda x: np.sin(x*np.pi/180.)
31cosd = lambda x: np.cos(x*np.pi/180.)
32tand = lambda x: np.tan(x*np.pi/180.)
33asind = lambda x: 180.*np.arcsin(x)/np.pi
34acosd = lambda x: 180.*np.arccos(x)/np.pi
35atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
36   
37ateln2 = 8.0*np.log(2.0)
38twopi = 2.0*np.pi
39twopisq = 2.0*np.pi**2
40
41################################################################################
42##### Rigid Body Models
43################################################################################
44       
45def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
46    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
47    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
48    '''
49    atxIds = ['Ax:','Ay:','Az:']
50    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
51    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
52    if not RBIds['Vector'] and not RBIds['Residue']:
53        return
54    VRBIds = RBIds['Vector']
55    RRBIds = RBIds['Residue']
56    if Update:
57        RBData = rigidbodyDict
58    else:
59        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
60    if RBIds['Vector']:                       # first update the vector magnitudes
61        VRBData = RBData['Vector']
62        for i,rbId in enumerate(VRBIds):
63            if VRBData[rbId]['useCount']:
64                for j in range(len(VRBData[rbId]['VectMag'])):
65                    name = '::RBV;'+str(j)+':'+str(i)
66                    VRBData[rbId]['VectMag'][j] = parmDict[name]
67    for phase in Phases:
68        Phase = Phases[phase]
69        General = Phase['General']
70        cx,ct,cs,cia = General['AtomPtrs']
71        cell = General['Cell'][1:7]
72        Amat,Bmat = G2lat.cell2AB(cell)
73        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
74        pfx = str(Phase['pId'])+'::'
75        if Update:
76            RBModels = Phase['RBModels']
77        else:
78            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
79        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
80            jrb = VRBIds.index(RBObj['RBId'])
81            rbsx = str(irb)+':'+str(jrb)
82            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
83                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
84            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
85                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
86            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
87            TLS = RBObj['ThermalMotion']
88            if 'T' in TLS[0]:
89                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
90                    TLS[1][i] = parmDict[pfx+pt+rbsx]
91            if 'L' in TLS[0]:
92                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
93                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
94            if 'S' in TLS[0]:
95                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
96                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
97            if 'U' in TLS[0]:
98                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
99            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
100            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
101            for i,x in enumerate(XYZ):
102                atId = RBObj['Ids'][i]
103                for j in [0,1,2]:
104                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
105                if UIJ[i][0] == 'A':
106                    for j in range(6):
107                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
108                elif UIJ[i][0] == 'I':
109                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
110           
111        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
112            jrb = RRBIds.index(RBObj['RBId'])
113            rbsx = str(irb)+':'+str(jrb)
114            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
115                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
116            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
117                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
118            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
119            TLS = RBObj['ThermalMotion']
120            if 'T' in TLS[0]:
121                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
122                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
123            if 'L' in TLS[0]:
124                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
125                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
126            if 'S' in TLS[0]:
127                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
128                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
129            if 'U' in TLS[0]:
130                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
131            for itors,tors in enumerate(RBObj['Torsions']):
132                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
133            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
134            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
135            for i,x in enumerate(XYZ):
136                atId = RBObj['Ids'][i]
137                for j in [0,1,2]:
138                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
139                if UIJ[i][0] == 'A':
140                    for j in range(6):
141                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
142                elif UIJ[i][0] == 'I':
143                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
144                   
145def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
146    'Needs a doc string'
147    atxIds = ['dAx:','dAy:','dAz:']
148    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
149    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
150    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
151    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
152    PIds = ['Px:','Py:','Pz:']
153    OIds = ['Oa:','Oi:','Oj:','Ok:']
154    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
155    if not RBIds['Vector'] and not RBIds['Residue']:
156        return
157    VRBIds = RBIds['Vector']
158    RRBIds = RBIds['Residue']
159    RBData = rigidbodyDict
160    for item in parmDict:
161        if 'RB' in item:
162            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
163    General = Phase['General']
164    cx,ct,cs,cia = General['AtomPtrs']
165    cell = General['Cell'][1:7]
166    Amat,Bmat = G2lat.cell2AB(cell)
167    rpd = np.pi/180.
168    rpd2 = rpd**2
169    g = nl.inv(np.inner(Bmat,Bmat))
170    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
171        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
172    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
173    pfx = str(Phase['pId'])+'::'
174    RBModels =  Phase['RBModels']
175    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
176        VModel = RBData['Vector'][RBObj['RBId']]
177        Q = RBObj['Orient'][0]
178        Pos = RBObj['Orig'][0]
179        jrb = VRBIds.index(RBObj['RBId'])
180        rbsx = str(irb)+':'+str(jrb)
181        dXdv = []
182        for iv in range(len(VModel['VectMag'])):
183            dCdv = []
184            for vec in VModel['rbVect'][iv]:
185                dCdv.append(G2mth.prodQVQ(Q,vec))
186            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
187        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
188        for ia,atId in enumerate(RBObj['Ids']):
189            atNum = AtLookup[atId]
190            dx = 0.00001
191            for iv in range(len(VModel['VectMag'])):
192                for ix in [0,1,2]:
193                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
194            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
195                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
196            for iv in range(4):
197                Q[iv] -= dx
198                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
199                Q[iv] += 2.*dx
200                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201                Q[iv] -= dx
202                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
203                for ix in [0,1,2]:
204                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
205            X = G2mth.prodQVQ(Q,Cart[ia])
206            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
207            dFdu = G2lat.U6toUij(dFdu.T)
208            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
209            dFdu = G2lat.UijtoU6(dFdu)
210            atNum = AtLookup[atId]
211            if 'T' in RBObj['ThermalMotion'][0]:
212                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
213                    dFdvDict[pfx+name+rbsx] += dFdu[i]
214            if 'L' in RBObj['ThermalMotion'][0]:
215                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
216                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
217                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
218                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
219                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
220                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
221                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
222                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
223                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
224            if 'S' in RBObj['ThermalMotion'][0]:
225                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
226                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
227                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
228                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
229                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
230                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
231                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
232                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
233            if 'U' in RBObj['ThermalMotion'][0]:
234                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
235
236
237    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
238        Q = RBObj['Orient'][0]
239        Pos = RBObj['Orig'][0]
240        jrb = RRBIds.index(RBObj['RBId'])
241        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
242        rbsx = str(irb)+':'+str(jrb)
243        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
244        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
245            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
246            orId,pvId = torData[itors][:2]
247            pivotVec = Cart[orId]-Cart[pvId]
248            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
249            QB = G2mth.AVdeg2Q(0.001,pivotVec)
250            for ir in torData[itors][3]:
251                atNum = AtLookup[RBObj['Ids'][ir]]
252                rVec = Cart[ir]-Cart[pvId]
253                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
254                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
255                for ix in [0,1,2]:
256                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
257        for ia,atId in enumerate(RBObj['Ids']):
258            atNum = AtLookup[atId]
259            dx = 0.00001
260            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
261                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
262            for iv in range(4):
263                Q[iv] -= dx
264                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
265                Q[iv] += 2.*dx
266                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
267                Q[iv] -= dx
268                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
269                for ix in [0,1,2]:
270                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
271            X = G2mth.prodQVQ(Q,Cart[ia])
272            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
273            dFdu = G2lat.U6toUij(dFdu.T)
274            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
275            dFdu = G2lat.UijtoU6(dFdu)
276            atNum = AtLookup[atId]
277            if 'T' in RBObj['ThermalMotion'][0]:
278                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
279                    dFdvDict[pfx+name+rbsx] += dFdu[i]
280            if 'L' in RBObj['ThermalMotion'][0]:
281                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
282                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
283                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
284                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
285                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
286                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
287                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
288                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
289                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
290            if 'S' in RBObj['ThermalMotion'][0]:
291                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
292                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
293                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
294                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
295                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
296                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
297                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
298                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
299            if 'U' in RBObj['ThermalMotion'][0]:
300                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
301   
302################################################################################
303##### Penalty & restraint functions
304################################################################################
305
306def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
307    'Needs a doc string'
308    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
309    pNames = []
310    pVals = []
311    pWt = []
312    negWt = {}
313    pWsum = {}
314    for phase in Phases:
315        pId = Phases[phase]['pId']
316        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
317        General = Phases[phase]['General']
318        cx,ct,cs,cia = General['AtomPtrs']
319        textureData = General['SH Texture']
320        SGData = General['SGData']
321        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
322        cell = General['Cell'][1:7]
323        Amat,Bmat = G2lat.cell2AB(cell)
324        if phase not in restraintDict:
325            continue
326        phaseRest = restraintDict[phase]
327        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
328            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
329            ['ChemComp','Sites'],['Texture','HKLs'],]
330        for name,rest in names:
331            pWsum[name] = 0.
332            itemRest = phaseRest[name]
333            if itemRest[rest] and itemRest['Use']:
334                wt = itemRest['wtFactor']
335                if name in ['Bond','Angle','Plane','Chiral']:
336                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
337                        pNames.append(str(pId)+':'+name+':'+str(i))
338                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
339                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
340                        if name == 'Bond':
341                            calc = G2mth.getRestDist(XYZ,Amat)
342                        elif name == 'Angle':
343                            calc = G2mth.getRestAngle(XYZ,Amat)
344                        elif name == 'Plane':
345                            calc = G2mth.getRestPlane(XYZ,Amat)
346                        elif name == 'Chiral':
347                            calc = G2mth.getRestChiral(XYZ,Amat)
348                        pVals.append(obs-calc)
349                        pWt.append(wt/esd**2)
350                        pWsum[name] += wt*((obs-calc)/esd)**2
351                elif name in ['Torsion','Rama']:
352                    coeffDict = itemRest['Coeff']
353                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
354                        pNames.append(str(pId)+':'+name+':'+str(i))
355                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
356                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
357                        if name == 'Torsion':
358                            tor = G2mth.getRestTorsion(XYZ,Amat)
359                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
360                        else:
361                            phi,psi = G2mth.getRestRama(XYZ,Amat)
362                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
363                        pVals.append(restr)
364                        pWt.append(wt/esd**2)
365                        pWsum[name] += wt*(restr/esd)**2
366                elif name == 'ChemComp':
367                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
368                        pNames.append(str(pId)+':'+name+':'+str(i))
369                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
370                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
371                        calc = np.sum(mul*frac*factors)
372                        pVals.append(obs-calc)
373                        pWt.append(wt/esd**2)                   
374                        pWsum[name] += wt*((obs-calc)/esd)**2
375                elif name == 'Texture':
376                    SHkeys = textureData['SH Coeff'][1].keys()
377                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
378                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
379                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
380                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
381                        PH = np.array(hkl)
382                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
383                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
384                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
385                        Z1 = -ma.masked_greater(Z,0.0)
386                        IndZ1 = np.array(ma.nonzero(Z1))
387                        for ind in IndZ1.T:
388                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
389                            pVals.append(Z1[ind[0]][ind[1]])
390                            pWt.append(wt/esd1**2)
391                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
392                        if ifesd2:
393                            Z2 = 1.-Z
394                            for ind in np.ndindex(grid,grid):
395                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
396                                pVals.append(Z1[ind[0]][ind[1]])
397                                pWt.append(wt/esd2**2)
398                                pWsum[name] += wt*(Z2/esd2)**2
399       
400    for phase in Phases:
401        name = 'SH-Pref.Ori.'
402        pId = Phases[phase]['pId']
403        General = Phases[phase]['General']
404        SGData = General['SGData']
405        cell = General['Cell'][1:7]
406        pWsum[name] = 0.0
407        for hist in Phases[phase]['Histograms']:
408            if hist in Histograms and 'PWDR' in hist:
409                hId = Histograms[hist]['hId']
410                phfx = '%d:%d:'%(pId,hId)
411                if calcControls[phfx+'poType'] == 'SH':
412                    toler = calcControls[phfx+'SHtoler']
413                    wt = 1./toler**2
414                    HKLs = np.array(calcControls[phfx+'SHhkl'])
415                    SHnames = calcControls[phfx+'SHnames']
416                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
417                    for i,PH in enumerate(HKLs):
418                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
419                        SH3Coef = {}
420                        for item in SHcof:
421                            L,N = eval(item.strip('C'))
422                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
423                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
424                        X = np.linspace(0,90.0,26)
425                        Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
426                        IndY = ma.nonzero(Y)
427                        for ind in IndY[0]:
428                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
429                            pVals.append(Y[ind])
430                            pWt.append(wt)
431                            pWsum[name] += wt*(Y[ind])**2
432    pWsum['PWLref'] = 0.
433    for item in varyList:
434        if 'PWLref' in item and parmDict[item] < 0.:
435            pId = int(item.split(':')[0])
436            if negWt[pId]:
437                pNames.append(item)
438                pVals.append(-parmDict[item])
439                pWt.append(negWt[pId])
440                pWsum['PWLref'] += negWt[pId]*(-parmDict[item])**2
441    pVals = np.array(pVals)
442    pWt = np.array(pWt)         #should this be np.sqrt?
443    return pNames,pVals,pWt,pWsum
444   
445def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
446    'Needs a doc string'
447    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
448    pDerv = np.zeros((len(varyList),len(pVal)))
449    for phase in Phases:
450#        if phase not in restraintDict:
451#            continue
452        pId = Phases[phase]['pId']
453        General = Phases[phase]['General']
454        cx,ct,cs,cia = General['AtomPtrs']
455        SGData = General['SGData']
456        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
457        cell = General['Cell'][1:7]
458        Amat,Bmat = G2lat.cell2AB(cell)
459        textureData = General['SH Texture']
460
461        SHkeys = textureData['SH Coeff'][1].keys()
462        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
463        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
464        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
465        sam = SamSym[textureData['Model']]
466        phaseRest = restraintDict.get(phase,{})
467        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
468            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
469            'ChemComp':'Sites','Texture':'HKLs'}
470        lasthkl = np.array([0,0,0])
471        for ip,pName in enumerate(pNames):
472            pnames = pName.split(':')
473            if pId == int(pnames[0]):
474                name = pnames[1]
475                if 'PWL' in pName:
476                    pDerv[varyList.index(pName)][ip] += 1.
477                    continue
478                elif 'SH-' in pName:
479                    continue
480                id = int(pnames[2]) 
481                itemRest = phaseRest[name]
482                if name in ['Bond','Angle','Plane','Chiral']:
483                    indx,ops,obs,esd = itemRest[names[name]][id]
484                    dNames = []
485                    for ind in indx:
486                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
487                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
488                    if name == 'Bond':
489                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
490                    elif name == 'Angle':
491                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
492                    elif name == 'Plane':
493                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
494                    elif name == 'Chiral':
495                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
496                elif name in ['Torsion','Rama']:
497                    coffDict = itemRest['Coeff']
498                    indx,ops,cofName,esd = itemRest[names[name]][id]
499                    dNames = []
500                    for ind in indx:
501                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
502                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
503                    if name == 'Torsion':
504                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
505                    else:
506                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
507                elif name == 'ChemComp':
508                    indx,factors,obs,esd = itemRest[names[name]][id]
509                    dNames = []
510                    for ind in indx:
511                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
512                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
513                        deriv = mul*factors
514                elif 'Texture' in name:
515                    deriv = []
516                    dNames = []
517                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
518                    hkl = np.array(hkl)
519                    if np.any(lasthkl-hkl):
520                        PH = np.array(hkl)
521                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
522                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
523                        lasthkl = copy.copy(hkl)                       
524                    if 'unit' in name:
525                        pass
526                    else:
527                        gam = float(pnames[3])
528                        psi = float(pnames[4])
529                        for SHname in ODFln:
530                            l,m,n = eval(SHname[1:])
531                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
532                            dNames += [str(pId)+'::'+SHname]
533                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
534                for dName,drv in zip(dNames,deriv):
535                    try:
536                        ind = varyList.index(dName)
537                        pDerv[ind][ip] += drv
538                    except ValueError:
539                        pass
540       
541        lasthkl = np.array([0,0,0])
542        for ip,pName in enumerate(pNames):
543            deriv = []
544            dNames = []
545            pnames = pName.split(':')
546            if 'SH-' in pName and pId == int(pnames[0]):
547                hId = int(pnames[1])
548                phfx = '%d:%d:'%(pId,hId)
549                psi = float(pnames[4])
550                HKLs = calcControls[phfx+'SHhkl']
551                SHnames = calcControls[phfx+'SHnames']
552                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
553                hkl = np.array(HKLs[int(pnames[3])])     
554                if np.any(lasthkl-hkl):
555                    PH = np.array(hkl)
556                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
557                    SH3Coef = {}
558                    for item in SHcof:
559                        L,N = eval(item.strip('C'))
560                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
561                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
562                    lasthkl = copy.copy(hkl)                       
563                for SHname in SHnames:
564                    l,n = eval(SHname[1:])
565                    SH3name = 'C%d,0,%d'%(l,n)
566                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
567                    dNames += [phfx+SHname]
568                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
569            for dName,drv in zip(dNames,deriv):
570                try:
571                    ind = varyList.index(dName)
572                    pDerv[ind][ip] += drv
573                except ValueError:
574                    pass
575    return pDerv
576
577################################################################################
578##### Function & derivative calculations
579################################################################################       
580                   
581def GetAtomFXU(pfx,calcControls,parmDict):
582    'Needs a doc string'
583    Natoms = calcControls['Natoms'][pfx]
584    Tdata = Natoms*[' ',]
585    Mdata = np.zeros(Natoms)
586    IAdata = Natoms*[' ',]
587    Fdata = np.zeros(Natoms)
588    FFdata = []
589    BLdata = []
590    Xdata = np.zeros((3,Natoms))
591    dXdata = np.zeros((3,Natoms))
592    Uisodata = np.zeros(Natoms)
593    Uijdata = np.zeros((6,Natoms))
594    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
595        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
596        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
597        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
598        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
599    for iatm in range(Natoms):
600        for key in keys:
601            parm = pfx+key+str(iatm)
602            if parm in parmDict:
603                keys[key][iatm] = parmDict[parm]
604    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.?
605    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
606   
607def GetAtomSSFXU(pfx,calcControls,parmDict):
608    'Needs a doc string'
609    Natoms = calcControls['Natoms'][pfx]
610    maxSSwave = calcControls['maxSSwave'][pfx]
611    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
612        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
613    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
614    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
615    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
616    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
617    waveTypes = []
618    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
619        'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
620        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
621        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
622        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
623        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
624    for iatm in range(Natoms):
625        waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
626        for key in keys:
627            for m in range(Nwave[key[0]]):
628                parm = pfx+key+str(iatm)+':%d'%(m)
629                if parm in parmDict:
630                    keys[key][m][iatm] = parmDict[parm]
631    return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()   
632   
633def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
634    ''' Not Used: here only for comparison the StructureFactor2 - faster version
635    Compute structure factors for all h,k,l for phase
636    puts the result, F^2, in each ref[8] in refList
637    input:
638   
639    :param dict refDict: where
640        'RefList' list where each ref = h,k,l,m,d,...
641        'FF' dict of form factors - filed in below
642    :param np.array G:      reciprocal metric tensor
643    :param str pfx:    phase id string
644    :param dict SGData: space group info. dictionary output from SpcGroup
645    :param dict calcControls:
646    :param dict ParmDict:
647
648    '''       
649    phfx = pfx.split(':')[0]+hfx
650    ast = np.sqrt(np.diag(G))
651    Mast = twopisq*np.multiply.outer(ast,ast)
652    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
653    SGT = np.array([ops[1] for ops in SGData['SGOps']])
654    FFtables = calcControls['FFtables']
655    BLtables = calcControls['BLtables']
656    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
657    FF = np.zeros(len(Tdata))
658    if 'NC' in calcControls[hfx+'histType']:
659        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
660    else:
661        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
662        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
663    Uij = np.array(G2lat.U6toUij(Uijdata))
664    bij = Mast*Uij.T
665    if not len(refDict['FF']):
666        if 'N' in calcControls[hfx+'histType']:
667            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
668        else:
669            dat = G2el.getFFvalues(FFtables,0.)       
670        refDict['FF']['El'] = dat.keys()
671        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
672    for iref,refl in enumerate(refDict['RefList']):
673        if 'NT' in calcControls[hfx+'histType']:
674            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14])
675        fbs = np.array([0,0])
676        H = refl[:3]
677        SQ = 1./(2.*refl[4])**2
678        SQfactor = 4.0*SQ*twopisq
679        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
680        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
681            if 'N' in calcControls[hfx+'histType']:
682                dat = G2el.getBLvalues(BLtables)
683                refDict['FF']['FF'][iref] = dat.values()
684            else:       #'X'
685                dat = G2el.getFFvalues(FFtables,SQ)
686                refDict['FF']['FF'][iref] = dat.values()
687        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
688        FF = refDict['FF']['FF'][iref][Tindx]
689        Uniq = np.inner(H,SGMT)
690        Phi = np.inner(H,SGT)
691        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
692        sinp = np.sin(phase)
693        cosp = np.cos(phase)
694        biso = -SQfactor*Uisodata
695        Tiso = np.where(biso<1.,np.exp(biso),1.0)
696        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
697        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
698        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
699        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
700        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
701        if not SGData['SGInv']:
702            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
703            fbs = np.sum(np.sum(fb,axis=1),axis=1)
704        fasq = fas**2
705        fbsq = fbs**2        #imaginary
706        refl[9] = np.sum(fasq)+np.sum(fbsq)
707        refl[10] = atan2d(fbs[0],fas[0])
708   
709def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
710    '''
711    Compute super structure factors for all h,k,l,m for phase
712    puts the result, F^2, in each ref[8+im] in refList
713    input:
714   
715    :param dict refDict: where
716        'RefList' list where each ref = h,k,l,m,d,...
717        'FF' dict of form factors - filed in below
718    :param np.array G:      reciprocal metric tensor
719    :param str pfx:    phase id string
720    :param dict SGData: space group info. dictionary output from SpcGroup
721    :param dict calcControls:
722    :param dict ParmDict:
723
724    '''       
725    phfx = pfx.split(':')[0]+hfx
726    ast = np.sqrt(np.diag(G))
727    Mast = twopisq*np.multiply.outer(ast,ast)
728    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
729    SGT = np.array([ops[1] for ops in SGData['SGOps']])
730    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
731    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
732    FFtables = calcControls['FFtables']
733    BLtables = calcControls['BLtables']
734    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
735    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
736    FF = np.zeros(len(Tdata))
737    if 'NC' in calcControls[hfx+'histType']:
738        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
739    else:
740        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
741        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
742    Uij = np.array(G2lat.U6toUij(Uijdata))
743    bij = Mast*Uij.T
744    if not len(refDict['FF']):
745        if 'N' in calcControls[hfx+'histType']:
746            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
747        else:
748            dat = G2el.getFFvalues(FFtables,0.)       
749        refDict['FF']['El'] = dat.keys()
750        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
751    for iref,refl in enumerate(refDict['RefList']):
752        if 'NT' in calcControls[hfx+'histType']:
753            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
754        fbs = np.array([0,0])
755        H = refl[:4]
756        SQ = 1./(2.*refl[4+im])**2
757        SQfactor = 4.0*SQ*twopisq
758        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
759        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
760            if 'N' in calcControls[hfx+'histType']:
761                dat = G2el.getBLvalues(BLtables)
762                refDict['FF']['FF'][iref] = dat.values()
763            else:       #'X'
764                dat = G2el.getFFvalues(FFtables,SQ)
765                refDict['FF']['FF'][iref] = dat.values()
766        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
767        FF = refDict['FF']['FF'][iref][Tindx]
768        Uniq = np.inner(H[:3],SGMT)
769        SSUniq = np.inner(H,SSGMT)
770        Phi = np.inner(H[:3],SGT)
771        SSPhi = np.inner(H,SSGT)
772        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)       
773        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
774        sinp = np.sin(phase)
775        cosp = np.cos(phase)
776        biso = -SQfactor*Uisodata
777        Tiso = np.where(biso<1.,np.exp(biso),1.0)
778        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
779        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
780        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
781        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
782        fb = np.zeros_like(fa)
783        if not SGData['SGInv']:
784            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
785        fa = fa*GfpuA-fb*GfpuB
786        fb = fb*GfpuA+fa*GfpuB
787        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
788        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
789           
790        fasq = fas**2
791        fbsq = fbs**2        #imaginary
792        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
793        refl[10+im] = atan2d(fbs[0],fas[0])
794   
795def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
796    ''' Compute structure factors for all h,k,l for phase
797    puts the result, F^2, in each ref[8] in refList
798    input:
799   
800    :param dict refDict: where
801        'RefList' list where each ref = h,k,l,m,d,...
802        'FF' dict of form factors - filed in below
803    :param np.array G:      reciprocal metric tensor
804    :param str pfx:    phase id string
805    :param dict SGData: space group info. dictionary output from SpcGroup
806    :param dict calcControls:
807    :param dict ParmDict:
808
809    '''       
810    phfx = pfx.split(':')[0]+hfx
811    ast = np.sqrt(np.diag(G))
812    Mast = twopisq*np.multiply.outer(ast,ast)
813    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
814    SGT = np.array([ops[1] for ops in SGData['SGOps']])
815    FFtables = calcControls['FFtables']
816    BLtables = calcControls['BLtables']
817    Flack = 1.0
818    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
819        Flack = 1.-2.*parmDict[phfx+'Flack']
820    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
821    FF = np.zeros(len(Tdata))
822    if 'NC' in calcControls[hfx+'histType']:
823        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
824    elif 'X' in calcControls[hfx+'histType']:
825        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
826        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
827    Uij = np.array(G2lat.U6toUij(Uijdata))
828    bij = Mast*Uij.T
829    blkSize = 100       #no. of reflections in a block
830    nRef = refDict['RefList'].shape[0]
831    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
832        if 'N' in calcControls[hfx+'histType']:
833            dat = G2el.getBLvalues(BLtables)
834            refDict['FF']['El'] = dat.keys()
835            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
836        else:       #'X'
837            dat = G2el.getFFvalues(FFtables,0.)
838            refDict['FF']['El'] = dat.keys()
839            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
840            for iref,ref in enumerate(refDict['RefList']):
841                SQ = 1./(2.*ref[4])**2
842                dat = G2el.getFFvalues(FFtables,SQ)
843                refDict['FF']['FF'][iref] *= dat.values()
844#reflection processing begins here - big arrays!
845    iBeg = 0           
846    while iBeg < nRef:
847        iFin = min(iBeg+blkSize,nRef)
848        refl = refDict['RefList'][iBeg:iFin]
849        H = refl.T[:3]
850        SQ = 1./(2.*refl.T[4])**2
851        SQfactor = 4.0*SQ*twopisq
852        if 'T' in calcControls[hfx+'histType']:
853            if 'P' in calcControls[hfx+'histType']:
854                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
855            else:
856                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
857            FP = np.repeat(FP.T,len(SGT),axis=0)
858            FPP = np.repeat(FPP.T,len(SGT),axis=0)
859        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
860        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
861        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
862        Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
863        Phi = np.inner(H.T,SGT).flatten()
864        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
865        sinp = np.sin(phase)
866        cosp = np.cos(phase)
867        biso = -SQfactor*Uisodata[:,np.newaxis]
868        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
869        HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
870        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
871        Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
872        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
873        fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))   #real A,-b
874        fas = np.sum(np.sum(fa,axis=2),axis=2)        #real sum over atoms & unique hkl
875        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
876        fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))   #imag -B,+a       
877        fbs = np.sum(np.sum(fb,axis=2),axis=2)  #imag sum over atoms & uniq hkl
878        if SGData['SGInv']: #centrosymmetric; B=0
879            fbs[0] *= 0.
880        if 'P' in calcControls[hfx+'histType']:
881            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
882        else:
883            refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
884        refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"?
885#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
886        iBeg += blkSize
887   
888def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
889    'Needs a doc string'
890    phfx = pfx.split(':')[0]+hfx
891    ast = np.sqrt(np.diag(G))
892    Mast = twopisq*np.multiply.outer(ast,ast)
893    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
894    SGT = np.array([ops[1] for ops in SGData['SGOps']])
895    FFtables = calcControls['FFtables']
896    BLtables = calcControls['BLtables']
897    nRef = len(refDict['RefList'])
898    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
899    mSize = len(Mdata)
900    FF = np.zeros(len(Tdata))
901    if 'NC' in calcControls[hfx+'histType']:
902        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
903    elif 'X' in calcControls[hfx+'histType']:
904        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
905        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
906    Uij = np.array(G2lat.U6toUij(Uijdata))
907    bij = Mast*Uij.T
908    dFdvDict = {}
909    dFdfr = np.zeros((nRef,mSize))
910    dFdx = np.zeros((nRef,mSize,3))
911    dFdui = np.zeros((nRef,mSize))
912    dFdua = np.zeros((nRef,mSize,6))
913    dFdbab = np.zeros((nRef,2))
914    dFdfl = np.zeros(nRef)
915    Flack = 1.0
916    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
917        Flack = 1.-2.*parmDict[phfx+'Flack']
918    for iref,refl in enumerate(refDict['RefList']):
919        if 'T' in calcControls[hfx+'histType']:
920            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
921        H = np.array(refl[:3])
922        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
923        SQfactor = 8.0*SQ*np.pi**2
924        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
925        Bab = parmDict[phfx+'BabA']*dBabdA
926        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
927        FF = refDict['FF']['FF'][iref].T[Tindx]
928        Uniq = np.inner(H,SGMT)
929        Phi = np.inner(H,SGT)
930        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
931        sinp = np.sin(phase)
932        cosp = np.cos(phase)
933        occ = Mdata*Fdata/len(Uniq)
934        biso = -SQfactor*Uisodata
935        Tiso = np.where(biso<1.,np.exp(biso),1.0)
936        HbH = -np.inner(H,np.inner(bij,H))
937        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
938        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
939        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
940        Tcorr = Tiso*Tuij
941        fot = (FF+FP-Bab)*occ*Tcorr
942        fotp = FPP*occ*Tcorr
943        fa = np.array([fot[:,np.newaxis]*cosp,-Flack*fotp[:,np.newaxis]*sinp])       #non positions
944        fb = np.array([fot[:,np.newaxis]*sinp,Flack*fotp[:,np.newaxis]*cosp])
945       
946        fas = np.sum(np.sum(fa,axis=1),axis=1)      #real sum over atoms & unique hkl
947        fbs = np.sum(np.sum(fb,axis=1),axis=1)      #imag sum over atoms & uniq hkl
948        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*cosp])   #positions
949        fbx = np.array([fot[:,np.newaxis]*cosp,-fotp[:,np.newaxis]*sinp])
950        #sum below is over Uniq
951        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
952        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
953        dfadui = np.sum(-SQfactor*fa,axis=2)
954        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
955        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
956        if not SGData['SGInv']:
957            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
958            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
959            dfbdui = np.sum(-SQfactor*fb,axis=2)
960            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
961            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
962            dfadfl = np.sum(-fotp[:,np.newaxis]*sinp)
963            dfbdfl = np.sum(fotp[:,np.newaxis]*cosp)
964        else:
965            dfbdfr = np.zeros_like(dfadfr)
966            dfbdx = np.zeros_like(dfadx)
967            dfbdui = np.zeros_like(dfadui)
968            dfbdua = np.zeros_like(dfadua)
969            dfbdba = np.zeros_like(dfadba)
970            dfadfl = 0.0
971            dfbdfl = 0.0
972        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
973        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
974            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
975                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
976            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
977                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
978            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
979                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
980            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
981                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
982        else:
983            SA = fas[0]-fbs[1]
984            SB = fbs[0]+fas[1]
985            dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
986                2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
987            dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])
988            dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])
989            dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])
990            dFdfl[iref] = -4.*SA*dfadfl-4.*SB*dfbdfl
991        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
992            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
993           
994        #loop over atoms - each dict entry is list of derivatives for all the reflections
995    for i in range(len(Mdata)):
996        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
997        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
998        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
999        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1000        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1001        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1002        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1003        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1004        dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
1005        dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
1006        dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
1007    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1008    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1009    dFdvDict[phfx+'Flack'] = dFdfl.T
1010    return dFdvDict
1011   
1012def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1013    'Needs a doc string'
1014    phfx = pfx.split(':')[0]+hfx
1015    ast = np.sqrt(np.diag(G))
1016    Mast = twopisq*np.multiply.outer(ast,ast)
1017    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1018    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1019    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1020    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1021    FFtables = calcControls['FFtables']
1022    BLtables = calcControls['BLtables']
1023    nRef = len(refDict['RefList'])
1024    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1025    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1026    mSize = len(Mdata)
1027    FF = np.zeros(len(Tdata))
1028    if 'NC' in calcControls[hfx+'histType']:
1029        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1030    elif 'X' in calcControls[hfx+'histType']:
1031        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1032        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1033    Uij = np.array(G2lat.U6toUij(Uijdata))
1034    bij = Mast*Uij.T
1035    dFdvDict = {}
1036    dFdfr = np.zeros((nRef,mSize))
1037    dFdx = np.zeros((nRef,mSize,3))
1038    dFdui = np.zeros((nRef,mSize))
1039    dFdua = np.zeros((nRef,mSize,6))
1040    dFdbab = np.zeros((nRef,2))
1041    for iref,refl in enumerate(refDict['RefList']):
1042        if 'T' in calcControls[hfx+'histType']:
1043            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1044        H = np.array(refl[:4])
1045        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1046        SQfactor = 8.0*SQ*np.pi**2
1047        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1048        Bab = parmDict[phfx+'BabA']*dBabdA
1049        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1050        FF = refDict['FF']['FF'][iref].T[Tindx]
1051        Uniq = np.inner(H[:3],SGMT)
1052        SSUniq = np.inner(H,SSGMT)
1053        Phi = np.inner(H[:3],SGT)
1054        SSPhi = np.inner(H,SSGT)
1055        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1056        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1057        #need ModulationDerv here dGAdXsin, etc 
1058        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
1059        sinp = np.sin(phase)
1060        cosp = np.cos(phase)
1061        occ = Mdata*Fdata/len(Uniq)
1062        biso = -SQfactor*Uisodata
1063        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1064        HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
1065        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1066        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1067        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1068        Tcorr = Tiso*Tuij
1069        fot = (FF+FP-Bab)*occ*Tcorr
1070        fotp = FPP*occ*Tcorr
1071        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1072        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1073        GfpuA = np.swapaxes(GfpuA,1,2)
1074        GfpuB = np.swapaxes(GfpuB,1,2)
1075        fa = fa*GfpuA-fb*GfpuB
1076        fb = fb*GfpuA+fa*GfpuB
1077       
1078        fas = np.sum(np.sum(fa,axis=1),axis=1)
1079        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1080        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1081        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1082        fax = fax*GfpuA-fbx*GfpuB
1083        fbx = fbx*GfpuA+fax*GfpuB
1084        #sum below is over Uniq
1085        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1086        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1087        dfadui = np.sum(-SQfactor*fa,axis=2)
1088        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1089        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
1090        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
1091        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1092        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1093        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1094        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1095        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1096        #need dFdXsin, etc for modulations...
1097        if not SGData['SGInv']:
1098            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1099            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
1100            dfbdui = np.sum(-SQfactor*fb,axis=2)
1101            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1102            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
1103            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1104            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1105            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1106            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1107            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1108        #need dFdXsin, etc for modulations...
1109        #loop over atoms - each dict entry is list of derivatives for all the reflections
1110    for i in range(len(Mdata)):     
1111        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1112        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1113        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1114        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1115        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1116        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1117        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1118        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1119        dFdvDict[pfx+'AU12:'+str(i)] = .5*dFdua.T[3][i]
1120        dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
1121        dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
1122        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1123    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1124    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1125    return dFdvDict
1126   
1127def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1128    ''' Single crystal extinction function; returns extinction & derivative
1129    '''
1130    extCor = 1.0
1131    dervDict = {}
1132    dervCor = 1.0
1133    if calcControls[phfx+'EType'] != 'None':
1134        SQ = 1/(4.*ref[4+im]**2)
1135        if 'C' in parmDict[hfx+'Type']:           
1136            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1137        else:   #'T'
1138            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1139        if 'SXC' in parmDict[hfx+'Type']:
1140            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1141            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1142            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1143            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1144        elif 'SNT' in parmDict[hfx+'Type']:
1145            AV = 1.e7/parmDict[pfx+'Vol']**2
1146            PL = SQ
1147            PLZ = AV*ref[9+im]*ref[12+im]**2
1148        elif 'SNC' in parmDict[hfx+'Type']:
1149            AV = 1.e7/parmDict[pfx+'Vol']**2
1150            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1151            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1152           
1153        if 'Primary' in calcControls[phfx+'EType']:
1154            PLZ *= 1.5
1155        else:
1156            if 'C' in parmDict[hfx+'Type']:
1157                PLZ *= calcControls[phfx+'Tbar']
1158            else: #'T'
1159                PLZ *= ref[13+im]      #t-bar
1160        if 'Primary' in calcControls[phfx+'EType']:
1161            PLZ *= 1.5
1162            PSIG = parmDict[phfx+'Ep']
1163        elif 'I & II' in calcControls[phfx+'EType']:
1164            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1165        elif 'Type II' in calcControls[phfx+'EType']:
1166            PSIG = parmDict[phfx+'Es']
1167        else:       # 'Secondary Type I'
1168            PSIG = parmDict[phfx+'Eg']/PL
1169           
1170        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1171        AL = 0.025+0.285*cos2T
1172        BG = 0.02-0.025*cos2T
1173        BL = 0.15-0.2*(0.75-cos2T)**2
1174        if cos2T < 0.:
1175            BL = -0.45*cos2T
1176        CG = 2.
1177        CL = 2.
1178        PF = PLZ*PSIG
1179       
1180        if 'Gaussian' in calcControls[phfx+'EApprox']:
1181            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1182            extCor = np.sqrt(PF4)
1183            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1184        else:
1185            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1186            extCor = np.sqrt(PF4)
1187            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1188
1189        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
1190        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1191            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1192        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1193            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1194        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1195            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1196               
1197    return 1./extCor,dervDict,dervCor
1198   
1199def Dict2Values(parmdict, varylist):
1200    '''Use before call to leastsq to setup list of values for the parameters
1201    in parmdict, as selected by key in varylist'''
1202    return [parmdict[key] for key in varylist] 
1203   
1204def Values2Dict(parmdict, varylist, values):
1205    ''' Use after call to leastsq to update the parameter dictionary with
1206    values corresponding to keys in varylist'''
1207    parmdict.update(zip(varylist,values))
1208   
1209def GetNewCellParms(parmDict,varyList):
1210    'Needs a doc string'
1211    newCellDict = {}
1212    Anames = ['A'+str(i) for i in range(6)]
1213    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1214    for item in varyList:
1215        keys = item.split(':')
1216        if keys[2] in Ddict:
1217            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1218            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1219            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1220    return newCellDict          # is e.g. {'0::D11':A0-D11}
1221   
1222def ApplyXYZshifts(parmDict,varyList):
1223    '''
1224    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1225   
1226    :param dict parmDict: parameter dictionary
1227    :param list varyList: list of variables (not used!)
1228    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1229
1230    '''
1231    newAtomDict = {}
1232    for item in parmDict:
1233        if 'dA' in item:
1234            parm = ''.join(item.split('d'))
1235            parmDict[parm] += parmDict[item]
1236            newAtomDict[item] = [parm,parmDict[parm]]
1237    return newAtomDict
1238   
1239def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1240    'Spherical harmonics texture'
1241    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1242    if 'T' in calcControls[hfx+'histType']:
1243        tth = parmDict[hfx+'2-theta']
1244    else:
1245        tth = refl[5+im]
1246    odfCor = 1.0
1247    H = refl[:3]
1248    cell = G2lat.Gmat2cell(g)
1249    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1250    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1251    phi,beta = G2lat.CrsAng(H,cell,SGData)
1252    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1253    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1254    for item in SHnames:
1255        L,M,N = eval(item.strip('C'))
1256        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1257        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1258        Lnorm = G2lat.Lnorm(L)
1259        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1260    return odfCor
1261   
1262def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1263    'Spherical harmonics texture derivatives'
1264    if 'T' in calcControls[hfx+'histType']:
1265        tth = parmDict[hfx+'2-theta']
1266    else:
1267        tth = refl[5+im]
1268    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1269    odfCor = 1.0
1270    dFdODF = {}
1271    dFdSA = [0,0,0]
1272    H = refl[:3]
1273    cell = G2lat.Gmat2cell(g)
1274    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1275    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1276    phi,beta = G2lat.CrsAng(H,cell,SGData)
1277    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1278    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1279    for item in SHnames:
1280        L,M,N = eval(item.strip('C'))
1281        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1282        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1283        Lnorm = G2lat.Lnorm(L)
1284        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1285        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1286        for i in range(3):
1287            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1288    return odfCor,dFdODF,dFdSA
1289   
1290def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1291    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1292    if 'T' in calcControls[hfx+'histType']:
1293        tth = parmDict[hfx+'2-theta']
1294    else:
1295        tth = refl[5+im]
1296    odfCor = 1.0
1297    H = refl[:3]
1298    cell = G2lat.Gmat2cell(g)
1299    Sangls = [0.,0.,0.]
1300    if 'Bragg' in calcControls[hfx+'instType']:
1301        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1302        IFCoup = True
1303    else:
1304        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1305        IFCoup = False
1306    phi,beta = G2lat.CrsAng(H,cell,SGData)
1307    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1308    SHnames = calcControls[phfx+'SHnames']
1309    for item in SHnames:
1310        L,N = eval(item.strip('C'))
1311        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1312        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1313        Lnorm = G2lat.Lnorm(L)
1314        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1315    return np.squeeze(odfCor)
1316   
1317def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1318    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1319    if 'T' in calcControls[hfx+'histType']:
1320        tth = parmDict[hfx+'2-theta']
1321    else:
1322        tth = refl[5+im]
1323    odfCor = 1.0
1324    dFdODF = {}
1325    H = refl[:3]
1326    cell = G2lat.Gmat2cell(g)
1327    Sangls = [0.,0.,0.]
1328    if 'Bragg' in calcControls[hfx+'instType']:
1329        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1330        IFCoup = True
1331    else:
1332        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1333        IFCoup = False
1334    phi,beta = G2lat.CrsAng(H,cell,SGData)
1335    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1336    SHnames = calcControls[phfx+'SHnames']
1337    for item in SHnames:
1338        L,N = eval(item.strip('C'))
1339        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1340        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1341        Lnorm = G2lat.Lnorm(L)
1342        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1343        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
1344    return odfCor,dFdODF
1345   
1346def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1347    'March-Dollase preferred orientation correction'
1348    POcorr = 1.0
1349    MD = parmDict[phfx+'MD']
1350    if MD != 1.0:
1351        MDAxis = calcControls[phfx+'MDAxis']
1352        sumMD = 0
1353        for H in uniq:           
1354            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1355            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1356            sumMD += A**3
1357        POcorr = sumMD/len(uniq)
1358    return POcorr
1359   
1360def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1361    'Needs a doc string'
1362    POcorr = 1.0
1363    POderv = {}
1364    if calcControls[phfx+'poType'] == 'MD':
1365        MD = parmDict[phfx+'MD']
1366        MDAxis = calcControls[phfx+'MDAxis']
1367        sumMD = 0
1368        sumdMD = 0
1369        for H in uniq:           
1370            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1371            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1372            sumMD += A**3
1373            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1374        POcorr = sumMD/len(uniq)
1375        POderv[phfx+'MD'] = sumdMD/len(uniq)
1376    else:   #spherical harmonics
1377        if calcControls[phfx+'SHord']:
1378            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1379    return POcorr,POderv
1380   
1381def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1382    'Needs a doc string'
1383    if 'Debye' in calcControls[hfx+'instType']:
1384        if 'T' in calcControls[hfx+'histType']:
1385            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1386        else:
1387            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1388    else:
1389        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1390   
1391def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1392    'Needs a doc string'
1393    if 'Debye' in calcControls[hfx+'instType']:
1394        if 'T' in calcControls[hfx+'histType']:
1395            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1396        else:
1397            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1398    else:
1399        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1400       
1401def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1402    'Needs a doc string'
1403    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1404    pi2 = np.sqrt(2./np.pi)
1405    if 'T' in calcControls[hfx+'histType']:
1406        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1407        wave = refl[14+im]
1408    else:   #'C'W
1409        sth2 = sind(refl[5+im]/2.)**2
1410        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1411    c2th = 1.-2.0*sth2
1412    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1413    if 'X' in calcControls[hfx+'histType']:
1414        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1415    xfac = flv2*parmDict[phfx+'Extinction']
1416    exb = 1.0
1417    if xfac > -1.:
1418        exb = 1./np.sqrt(1.+xfac)
1419    exl = 1.0
1420    if 0 < xfac <= 1.:
1421        xn = np.array([xfac**(i+1) for i in range(6)])
1422        exl += np.sum(xn*coef)
1423    elif xfac > 1.:
1424        xfac2 = 1./np.sqrt(xfac)
1425        exl = pi2*(1.-0.125/xfac)*xfac2
1426    return exb*sth2+exl*(1.-sth2)
1427   
1428def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1429    'Needs a doc string'
1430    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1431    pi2 = np.sqrt(2./np.pi)
1432    if 'T' in calcControls[hfx+'histType']:
1433        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1434        wave = refl[14+im]
1435    else:   #'C'W
1436        sth2 = sind(refl[5+im]/2.)**2
1437        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1438    c2th = 1.-2.0*sth2
1439    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1440    if 'X' in calcControls[hfx+'histType']:
1441        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1442    xfac = flv2*parmDict[phfx+'Extinction']
1443    dbde = -500.*flv2
1444    if xfac > -1.:
1445        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
1446    dlde = 0.
1447    if 0 < xfac <= 1.:
1448        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1449        dlde = np.sum(xn*coef)
1450    elif xfac > 1.:
1451        xfac2 = 1./np.sqrt(xfac)
1452        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1453       
1454    return dbde*sth2+dlde*(1.-sth2)
1455   
1456def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1457    'Needs a doc string'    #need powder extinction!
1458    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1459    if 'X' in parmDict[hfx+'Type']:
1460        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1461    POcorr = 1.0
1462    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1463        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1464    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1465        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1466    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1467        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1468    Icorr *= POcorr
1469    AbsCorr = 1.0
1470    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1471    Icorr *= AbsCorr
1472    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1473    Icorr *= ExtCorr
1474    return Icorr,POcorr,AbsCorr,ExtCorr
1475   
1476def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1477    'Needs a doc string'    #need powder extinction derivs!
1478    dIdsh = 1./parmDict[hfx+'Scale']
1479    dIdsp = 1./parmDict[phfx+'Scale']
1480    if 'X' in parmDict[hfx+'Type']:
1481        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1482        dIdPola /= pola
1483    else:       #'N'
1484        dIdPola = 0.0
1485    dFdODF = {}
1486    dFdSA = [0,0,0]
1487    dIdPO = {}
1488    if pfx+'SHorder' in parmDict:
1489        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1490        for iSH in dFdODF:
1491            dFdODF[iSH] /= odfCor
1492        for i in range(3):
1493            dFdSA[i] /= odfCor
1494    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1495        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1496        for iPO in dIdPO:
1497            dIdPO[iPO] /= POcorr
1498    if 'T' in parmDict[hfx+'Type']:
1499        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1500        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1501    else:
1502        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1503        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1504    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1505       
1506def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1507    'Needs a doc string'
1508    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1509        costh = cosd(refl[5+im]/2.)
1510        #crystallite size
1511        if calcControls[phfx+'SizeType'] == 'isotropic':
1512            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1513        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1514            H = np.array(refl[:3])
1515            P = np.array(calcControls[phfx+'SizeAxis'])
1516            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1517            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1518            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1519        else:           #ellipsoidal crystallites
1520            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1521            H = np.array(refl[:3])
1522            lenR = G2pwd.ellipseSize(H,Sij,GB)
1523            Sgam = 1.8*wave/(np.pi*costh*lenR)
1524        #microstrain               
1525        if calcControls[phfx+'MustrainType'] == 'isotropic':
1526            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1527        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1528            H = np.array(refl[:3])
1529            P = np.array(calcControls[phfx+'MustrainAxis'])
1530            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1531            Si = parmDict[phfx+'Mustrain;i']
1532            Sa = parmDict[phfx+'Mustrain;a']
1533            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1534        else:       #generalized - P.W. Stephens model
1535            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1536            Sum = 0
1537            for i,strm in enumerate(Strms):
1538                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1539            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1540    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1541        #crystallite size
1542        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1543            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1544        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1545            H = np.array(refl[:3])
1546            P = np.array(calcControls[phfx+'SizeAxis'])
1547            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1548            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1549            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1550        else:           #ellipsoidal crystallites   #OK
1551            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1552            H = np.array(refl[:3])
1553            lenR = G2pwd.ellipseSize(H,Sij,GB)
1554            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1555        #microstrain               
1556        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1557            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1558        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1559            H = np.array(refl[:3])
1560            P = np.array(calcControls[phfx+'MustrainAxis'])
1561            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1562            Si = parmDict[phfx+'Mustrain;i']
1563            Sa = parmDict[phfx+'Mustrain;a']
1564            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1565        else:       #generalized - P.W. Stephens model  OK
1566            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1567            Sum = 0
1568            for i,strm in enumerate(Strms):
1569                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1570            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1571           
1572    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1573    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1574    sig /= ateln2
1575    return sig,gam
1576       
1577def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1578    'Needs a doc string'
1579    gamDict = {}
1580    sigDict = {}
1581    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1582        costh = cosd(refl[5+im]/2.)
1583        tanth = tand(refl[5+im]/2.)
1584        #crystallite size derivatives
1585        if calcControls[phfx+'SizeType'] == 'isotropic':
1586            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1587            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1588            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1589        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1590            H = np.array(refl[:3])
1591            P = np.array(calcControls[phfx+'SizeAxis'])
1592            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1593            Si = parmDict[phfx+'Size;i']
1594            Sa = parmDict[phfx+'Size;a']
1595            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1596            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1597            Sgam = gami*sqtrm
1598            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1599            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1600            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1601            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1602            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1603            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1604        else:           #ellipsoidal crystallites
1605            const = 1.8*wave/(np.pi*costh)
1606            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1607            H = np.array(refl[:3])
1608            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1609            Sgam = const/lenR
1610            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1611                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1612                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1613        gamDict[phfx+'Size;mx'] = Sgam
1614        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1615               
1616        #microstrain derivatives               
1617        if calcControls[phfx+'MustrainType'] == 'isotropic':
1618            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1619            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1620            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1621        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1622            H = np.array(refl[:3])
1623            P = np.array(calcControls[phfx+'MustrainAxis'])
1624            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1625            Si = parmDict[phfx+'Mustrain;i']
1626            Sa = parmDict[phfx+'Mustrain;a']
1627            gami = 0.018*Si*Sa*tanth/np.pi
1628            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1629            Mgam = gami/sqtrm
1630            dsi = -gami*Si*cosP**2/sqtrm**3
1631            dsa = -gami*Sa*sinP**2/sqtrm**3
1632            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1633            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1634            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1635            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1636        else:       #generalized - P.W. Stephens model
1637            const = 0.018*refl[4+im]**2*tanth/np.pi
1638            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1639            Sum = 0
1640            for i,strm in enumerate(Strms):
1641                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1642                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1643                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1644            Mgam = const*np.sqrt(Sum)
1645            for i in range(len(Strms)):
1646                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1647                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1648        gamDict[phfx+'Mustrain;mx'] = Mgam
1649        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1650    else:   #'T'OF - All checked & OK
1651        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1652            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1653            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1654            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1655        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1656            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1657            H = np.array(refl[:3])
1658            P = np.array(calcControls[phfx+'SizeAxis'])
1659            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1660            Si = parmDict[phfx+'Size;i']
1661            Sa = parmDict[phfx+'Size;a']
1662            gami = const/(Si*Sa)
1663            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1664            Sgam = gami*sqtrm
1665            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1666            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1667            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1668            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1669            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1670            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1671        else:           #OK  ellipsoidal crystallites
1672            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1673            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1674            H = np.array(refl[:3])
1675            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1676            Sgam = const/lenR
1677            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1678                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1679                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1680        gamDict[phfx+'Size;mx'] = Sgam  #OK
1681        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1682               
1683        #microstrain derivatives               
1684        if calcControls[phfx+'MustrainType'] == 'isotropic':
1685            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1686            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1687            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1688        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1689            H = np.array(refl[:3])
1690            P = np.array(calcControls[phfx+'MustrainAxis'])
1691            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1692            Si = parmDict[phfx+'Mustrain;i']
1693            Sa = parmDict[phfx+'Mustrain;a']
1694            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1695            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1696            Mgam = gami/sqtrm
1697            dsi = -gami*Si*cosP**2/sqtrm**3
1698            dsa = -gami*Sa*sinP**2/sqtrm**3
1699            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1700            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1701            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1702            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1703        else:       #generalized - P.W. Stephens model OK
1704            pwrs = calcControls[phfx+'MuPwrs']
1705            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1706            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1707            Sum = 0
1708            for i,strm in enumerate(Strms):
1709                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1710                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1711                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1712            Mgam = const*np.sqrt(Sum)
1713            for i in range(len(Strms)):
1714                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1715                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1716        gamDict[phfx+'Mustrain;mx'] = Mgam
1717        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1718       
1719    return sigDict,gamDict
1720       
1721def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1722    'Needs a doc string'
1723    if im:
1724        h,k,l,m = refl[:4]
1725        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1726        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1727    else:
1728        h,k,l = refl[:3]
1729        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1730    refl[4+im] = d
1731    if 'C' in calcControls[hfx+'histType']:
1732        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1733        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1734        if 'Bragg' in calcControls[hfx+'instType']:
1735            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1736                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1737        else:               #Debye-Scherrer - simple but maybe not right
1738            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1739    elif 'T' in calcControls[hfx+'histType']:
1740        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1741        #do I need sample position effects - maybe?
1742    return pos
1743
1744def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1745    'Needs a doc string'
1746    dpr = 180./np.pi
1747    if im:
1748        h,k,l,m = refl[:4]
1749        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1750        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1751    else:
1752        m = 0
1753        h,k,l = refl[:3]       
1754        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1755    dst = np.sqrt(dstsq)
1756    dsp = 1./dst
1757    if 'C' in calcControls[hfx+'histType']:
1758        pos = refl[5+im]-parmDict[hfx+'Zero']
1759        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1760        dpdw = const*dst
1761        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1762        dpdZ = 1.0
1763        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1764            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1765        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1766        if 'Bragg' in calcControls[hfx+'instType']:
1767            dpdSh = -4.*shft*cosd(pos/2.0)
1768            dpdTr = -shft*sind(pos)*100.0
1769            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1770        else:               #Debye-Scherrer - simple but maybe not right
1771            dpdXd = -shft*cosd(pos)
1772            dpdYd = -shft*sind(pos)
1773            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1774    elif 'T' in calcControls[hfx+'histType']:
1775        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1776        dpdZ = 1.0
1777        dpdDC = dsp
1778        dpdDA = dsp**2
1779        dpdDB = 1./dsp
1780        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1781            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1782        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1783           
1784def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1785    'Needs a doc string'
1786    laue = SGData['SGLaue']
1787    uniq = SGData['SGUniq']
1788    h,k,l = refl[:3]
1789    if laue in ['m3','m3m']:
1790        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1791            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1792    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1793        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1794    elif laue in ['3R','3mR']:
1795        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1796    elif laue in ['4/m','4/mmm']:
1797        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1798    elif laue in ['mmm']:
1799        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1800    elif laue in ['2/m']:
1801        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1802        if uniq == 'a':
1803            Dij += parmDict[phfx+'D23']*k*l
1804        elif uniq == 'b':
1805            Dij += parmDict[phfx+'D13']*h*l
1806        elif uniq == 'c':
1807            Dij += parmDict[phfx+'D12']*h*k
1808    else:
1809        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1810            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1811    if 'C' in calcControls[hfx+'histType']:
1812        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1813    else:
1814        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1815           
1816def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1817    'Needs a doc string'
1818    laue = SGData['SGLaue']
1819    uniq = SGData['SGUniq']
1820    h,k,l = refl[:3]
1821    if laue in ['m3','m3m']:
1822        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1823            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1824    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1825        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1826    elif laue in ['3R','3mR']:
1827        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1828    elif laue in ['4/m','4/mmm']:
1829        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1830    elif laue in ['mmm']:
1831        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1832    elif laue in ['2/m']:
1833        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1834        if uniq == 'a':
1835            dDijDict[phfx+'D23'] = k*l
1836        elif uniq == 'b':
1837            dDijDict[phfx+'D13'] = h*l
1838        elif uniq == 'c':
1839            dDijDict[phfx+'D12'] = h*k
1840    else:
1841        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1842            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1843    if 'C' in calcControls[hfx+'histType']:
1844        for item in dDijDict:
1845            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1846    else:
1847        for item in dDijDict:
1848            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1849    return dDijDict
1850   
1851def GetDij(phfx,SGData,parmDict):
1852    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1853    return G2spc.HStrainVals(HSvals,SGData)
1854               
1855def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1856    'Needs a doc string'
1857    histoList = Histograms.keys()
1858    histoList.sort()
1859    for histogram in histoList:
1860        if 'PWDR' in histogram[:4]:
1861            Histogram = Histograms[histogram]
1862            hId = Histogram['hId']
1863            hfx = ':%d:'%(hId)
1864            Limits = calcControls[hfx+'Limits']
1865            if 'C' in calcControls[hfx+'histType']:
1866                shl = max(parmDict[hfx+'SH/L'],0.0005)
1867                Ka2 = False
1868                kRatio = 0.0
1869                if hfx+'Lam1' in parmDict.keys():
1870                    Ka2 = True
1871                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1872                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1873            x,y,w,yc,yb,yd = Histogram['Data']
1874            xB = np.searchsorted(x,Limits[0])
1875            xF = np.searchsorted(x,Limits[1])
1876            ymb = np.array(y-yb)
1877            ymb = np.where(ymb,ymb,1.0)
1878            ycmb = np.array(yc-yb)
1879            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1880            refLists = Histogram['Reflection Lists']
1881            for phase in refLists:
1882                if phase not in Phases:     #skips deleted or renamed phases silently!
1883                    continue
1884                Phase = Phases[phase]
1885                im = 0
1886                if Phase['General']['Type'] in ['modulated','magnetic']:
1887                    im = 1
1888                pId = Phase['pId']
1889                phfx = '%d:%d:'%(pId,hId)
1890                refDict = refLists[phase]
1891                sumFo = 0.0
1892                sumdF = 0.0
1893                sumFosq = 0.0
1894                sumdFsq = 0.0
1895                sumInt = 0.0
1896                for refl in refDict['RefList']:
1897                    if 'C' in calcControls[hfx+'histType']:
1898                        yp = np.zeros_like(yb)
1899                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1900                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1901                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1902                        iFin2 = iFin
1903                        if not iBeg+iFin:       #peak below low limit - skip peak
1904                            continue
1905                        elif not iBeg-iFin:     #peak above high limit - done
1906                            break
1907                        elif iBeg < iFin:
1908                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
1909                            sumInt += refl[11+im]*refl[9+im]
1910                            if Ka2:
1911                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1912                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1913                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1914                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1915                                if iFin2 > iBeg2: 
1916                                    yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
1917                                    sumInt += refl[11+im]*refl[9+im]*kRatio
1918                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1919                               
1920                    elif 'T' in calcControls[hfx+'histType']:
1921                        yp = np.zeros_like(yb)
1922                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1923                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1924                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1925                        if iBeg < iFin:
1926                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
1927                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1928                            sumInt += refl[11+im]*refl[9+im]
1929                    Fo = np.sqrt(np.abs(refl[8+im]))
1930                    Fc = np.sqrt(np.abs(refl[9]+im))
1931                    sumFo += Fo
1932                    sumFosq += refl[8+im]**2
1933                    sumdF += np.abs(Fo-Fc)
1934                    sumdFsq += (refl[8+im]-refl[9+im])**2
1935                if sumFo:
1936                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1937                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1938                else:
1939                    Histogram['Residuals'][phfx+'Rf'] = 100.
1940                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1941                Histogram['Residuals'][phfx+'sumInt'] = sumInt
1942                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1943                Histogram['Residuals']['hId'] = hId
1944        elif 'HKLF' in histogram[:4]:
1945            Histogram = Histograms[histogram]
1946            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1947               
1948def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1949    'Needs a doc string'
1950   
1951    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1952        U = parmDict[hfx+'U']
1953        V = parmDict[hfx+'V']
1954        W = parmDict[hfx+'W']
1955        X = parmDict[hfx+'X']
1956        Y = parmDict[hfx+'Y']
1957        tanPos = tand(refl[5+im]/2.0)
1958        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1959        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1960        sig = max(0.001,sig)
1961        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1962        gam = max(0.001,gam)
1963        return sig,gam
1964               
1965    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1966        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1967            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1968        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1969        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1970        sig += Ssig
1971        gam += Sgam
1972        return sig,gam
1973       
1974    def GetReflAlpBet(refl,im,hfx,parmDict):
1975        alp = parmDict[hfx+'alpha']/refl[4+im]
1976        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1977        return alp,bet
1978       
1979    hId = Histogram['hId']
1980    hfx = ':%d:'%(hId)
1981    bakType = calcControls[hfx+'bakType']
1982    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1983    yc = np.zeros_like(yb)
1984    cw = np.diff(x)
1985    cw = np.append(cw,cw[-1])
1986       
1987    if 'C' in calcControls[hfx+'histType']:   
1988        shl = max(parmDict[hfx+'SH/L'],0.002)
1989        Ka2 = False
1990        if hfx+'Lam1' in parmDict.keys():
1991            wave = parmDict[hfx+'Lam1']
1992            Ka2 = True
1993            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1994            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1995        else:
1996            wave = parmDict[hfx+'Lam']
1997    for phase in Histogram['Reflection Lists']:
1998        refDict = Histogram['Reflection Lists'][phase]
1999        if phase not in Phases:     #skips deleted or renamed phases silently!
2000            continue
2001        Phase = Phases[phase]
2002        pId = Phase['pId']
2003        pfx = '%d::'%(pId)
2004        phfx = '%d:%d:'%(pId,hId)
2005        hfx = ':%d:'%(hId)
2006        SGData = Phase['General']['SGData']
2007        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2008        im = 0
2009        if Phase['General']['Type'] in ['modulated','magnetic']:
2010            SSGData = Phase['General']['SSGData']
2011            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2012            im = 1  #offset in SS reflection list
2013            #??
2014        Dij = GetDij(phfx,SGData,parmDict)
2015        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2016        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2017        if np.any(np.diag(G)<0.):
2018            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2019        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2020        Vst = np.sqrt(nl.det(G))    #V*
2021        if not Phase['General'].get('doPawley'):
2022            time0 = time.time()
2023            if im:
2024                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2025            else:
2026                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2027#            print 'sf calc time: %.3fs'%(time.time()-time0)
2028        time0 = time.time()
2029        badPeak = False
2030        for iref,refl in enumerate(refDict['RefList']):
2031            if 'C' in calcControls[hfx+'histType']:
2032                if im:
2033                    h,k,l,m = refl[:4]
2034                else:
2035                    h,k,l = refl[:3]
2036                Uniq = np.inner(refl[:3],SGMT)
2037                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2038                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2039                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2040                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2041                refl[11+im] *= Vst*Lorenz
2042                 
2043                if Phase['General'].get('doPawley'):
2044                    try:
2045                        if im:
2046                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2047                        else:
2048                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2049                        refl[9+im] = parmDict[pInd]
2050                    except KeyError:
2051#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2052                        continue
2053                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2054                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2055                iFin = np.searchsorted(x,refl[5+im]+fmax)
2056                if not iBeg+iFin:       #peak below low limit - skip peak
2057                    continue
2058                elif not iBeg-iFin:     #peak above high limit - done
2059                    break
2060                elif iBeg > iFin:   #bad peak coeff - skip
2061                    badPeak = True
2062                    continue
2063                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
2064                if Ka2:
2065                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2066                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2067                    iBeg = np.searchsorted(x,pos2-fmin)
2068                    iFin = np.searchsorted(x,pos2+fmax)
2069                    if not iBeg+iFin:       #peak below low limit - skip peak
2070                        continue
2071                    elif not iBeg-iFin:     #peak above high limit - done
2072                        return yc,yb
2073                    elif iBeg > iFin:   #bad peak coeff - skip
2074                        continue
2075                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
2076            elif 'T' in calcControls[hfx+'histType']:
2077                h,k,l = refl[:3]
2078                Uniq = np.inner(refl[:3],SGMT)
2079                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2080                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2081#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2082                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2083                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2084                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2085                refl[11+im] *= Vst*Lorenz
2086                if Phase['General'].get('doPawley'):
2087                    try:
2088                        if im:
2089                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2090                        else:
2091                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2092                        refl[9+im] = parmDict[pInd]
2093                    except KeyError:
2094#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2095                        continue
2096                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2097                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2098                iFin = np.searchsorted(x,refl[5+im]+fmax)
2099                if not iBeg+iFin:       #peak below low limit - skip peak
2100                    continue
2101                elif not iBeg-iFin:     #peak above high limit - done
2102                    break
2103                elif iBeg > iFin:   #bad peak coeff - skip
2104                    badPeak = True
2105                    continue
2106                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
2107#        print 'profile calc time: %.3fs'%(time.time()-time0)
2108    if badPeak:
2109        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2110    return yc,yb
2111   
2112def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2113    'Needs a doc string'
2114   
2115    def cellVaryDerv(pfx,SGData,dpdA): 
2116        if SGData['SGLaue'] in ['-1',]:
2117            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2118                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2119        elif SGData['SGLaue'] in ['2/m',]:
2120            if SGData['SGUniq'] == 'a':
2121                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2122            elif SGData['SGUniq'] == 'b':
2123                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2124            else:
2125                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2126        elif SGData['SGLaue'] in ['mmm',]:
2127            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2128        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2129            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2130        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2131            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2132        elif SGData['SGLaue'] in ['3R', '3mR']:
2133            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2134        elif SGData['SGLaue'] in ['m3m','m3']:
2135            return [[pfx+'A0',dpdA[0]]]
2136           
2137    # create a list of dependent variables and set up a dictionary to hold their derivatives
2138    dependentVars = G2mv.GetDependentVars()
2139    depDerivDict = {}
2140    for j in dependentVars:
2141        depDerivDict[j] = np.zeros(shape=(len(x)))
2142    #print 'dependent vars',dependentVars
2143    lenX = len(x)               
2144    hId = Histogram['hId']
2145    hfx = ':%d:'%(hId)
2146    bakType = calcControls[hfx+'bakType']
2147    dMdv = np.zeros(shape=(len(varylist),len(x)))
2148    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2149    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2150        bBpos =varylist.index(hfx+'Back;0')
2151        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2152    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2153    for name in varylist:
2154        if 'Debye' in name:
2155            id = int(name.split(';')[-1])
2156            parm = name[:int(name.rindex(';'))]
2157            ip = names.index(parm)
2158            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2159    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2160    for name in varylist:
2161        if 'BkPk' in name:
2162            parm,id = name.split(';')
2163            id = int(id)
2164            if parm in names:
2165                ip = names.index(parm)
2166                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2167    cw = np.diff(x)
2168    cw = np.append(cw,cw[-1])
2169    Ka2 = False #also for TOF!
2170    if 'C' in calcControls[hfx+'histType']:   
2171        shl = max(parmDict[hfx+'SH/L'],0.002)
2172        if hfx+'Lam1' in parmDict.keys():
2173            wave = parmDict[hfx+'Lam1']
2174            Ka2 = True
2175            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2176            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2177        else:
2178            wave = parmDict[hfx+'Lam']
2179    for phase in Histogram['Reflection Lists']:
2180        refDict = Histogram['Reflection Lists'][phase]
2181        if phase not in Phases:     #skips deleted or renamed phases silently!
2182            continue
2183        Phase = Phases[phase]
2184        SGData = Phase['General']['SGData']
2185        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2186        im = 0
2187        if Phase['General']['Type'] in ['modulated','magnetic']:
2188            SSGData = Phase['General']['SSGData']
2189            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2190            im = 1  #offset in SS reflection list
2191            #??
2192        pId = Phase['pId']
2193        pfx = '%d::'%(pId)
2194        phfx = '%d:%d:'%(pId,hId)
2195        Dij = GetDij(phfx,SGData,parmDict)
2196        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2197        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2198        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2199        if not Phase['General'].get('doPawley'):
2200            time0 = time.time()
2201            if im:
2202                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2203            else:
2204                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2205#            print 'sf-derv time %.3fs'%(time.time()-time0)
2206            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2207        time0 = time.time()
2208        for iref,refl in enumerate(refDict['RefList']):
2209            if im:
2210                h,k,l,m = refl[:4]
2211            else:
2212                h,k,l = refl[:3]
2213            Uniq = np.inner(refl[:3],SGMT)
2214            if 'T' in calcControls[hfx+'histType']:
2215                wave = refl[14+im]
2216            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2217            if 'C' in calcControls[hfx+'histType']:        #CW powder
2218                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2219            else: #'T'OF
2220                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2221            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2222            iFin = np.searchsorted(x,refl[5+im]+fmax)
2223            if not iBeg+iFin:       #peak below low limit - skip peak
2224                continue
2225            elif not iBeg-iFin:     #peak above high limit - done
2226                break
2227            pos = refl[5+im]
2228            if 'C' in calcControls[hfx+'histType']:
2229                tanth = tand(pos/2.0)
2230                costh = cosd(pos/2.0)
2231                lenBF = iFin-iBeg
2232                dMdpk = np.zeros(shape=(6,lenBF))
2233                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2234                for i in range(5):
2235                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2236                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2237                if Ka2:
2238                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2239                    iBeg2 = np.searchsorted(x,pos2-fmin)
2240                    iFin2 = np.searchsorted(x,pos2+fmax)
2241                    if iBeg2-iFin2:
2242                        lenBF2 = iFin2-iBeg2
2243                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2244                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2245                        for i in range(5):
2246                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2247                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2248                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2249            else:   #'T'OF
2250                lenBF = iFin-iBeg
2251                if lenBF < 0:   #bad peak coeff
2252                    break
2253                dMdpk = np.zeros(shape=(6,lenBF))
2254                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2255                for i in range(6):
2256                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2257                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2258            if Phase['General'].get('doPawley'):
2259                dMdpw = np.zeros(len(x))
2260                try:
2261                    if im:
2262                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2263                    else:
2264                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2265                    idx = varylist.index(pIdx)
2266                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2267                    if Ka2: #not for TOF either
2268                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2269                    dMdv[idx] = dMdpw
2270                except: # ValueError:
2271                    pass
2272            if 'C' in calcControls[hfx+'histType']:
2273                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2274                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2275                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2276                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2277                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2278                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2279                    hfx+'DisplaceY':[dpdY,'pos'],}
2280                if 'Bragg' in calcControls[hfx+'instType']:
2281                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2282                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2283                else:
2284                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2285            else:   #'T'OF
2286                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2287                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2288                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2289                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2290                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2291                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2292                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2293                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2294            for name in names:
2295                item = names[name]
2296                if name in varylist:
2297                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2298                    if Ka2 and iFin2-iBeg2:
2299                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2300                elif name in dependentVars:
2301                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2302                    if Ka2 and iFin2-iBeg2:
2303                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2304            for iPO in dIdPO:
2305                if iPO in varylist:
2306                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2307                    if Ka2 and iFin2-iBeg2:
2308                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2309                elif iPO in dependentVars:
2310                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2311                    if Ka2 and iFin2-iBeg2:
2312                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2313            for i,name in enumerate(['omega','chi','phi']):
2314                aname = pfx+'SH '+name
2315                if aname in varylist:
2316                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2317                    if Ka2 and iFin2-iBeg2:
2318                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2319                elif aname in dependentVars:
2320                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2321                    if Ka2 and iFin2-iBeg2:
2322                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2323            for iSH in dFdODF:
2324                if iSH in varylist:
2325                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2326                    if Ka2 and iFin2-iBeg2:
2327                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2328                elif iSH in dependentVars:
2329                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2330                    if Ka2 and iFin2-iBeg2:
2331                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2332            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2333            for name,dpdA in cellDervNames:
2334                if name in varylist:
2335                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2336                    if Ka2 and iFin2-iBeg2:
2337                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2338                elif name in dependentVars:
2339                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2340                    if Ka2 and iFin2-iBeg2:
2341                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2342            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2343            for name in dDijDict:
2344                if name in varylist:
2345                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2346                    if Ka2 and iFin2-iBeg2:
2347                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2348                elif name in dependentVars:
2349                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2350                    if Ka2 and iFin2-iBeg2:
2351                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2352            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2353                if name in varylist:
2354                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2355                    if Ka2 and iFin2-iBeg2:
2356                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2357                elif name in dependentVars:
2358                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2359                    if Ka2 and iFin2-iBeg2:
2360                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2361            if 'C' in calcControls[hfx+'histType']:
2362                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2363            else:   #'T'OF
2364                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2365            for name in gamDict:
2366                if name in varylist:
2367                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2368                    if Ka2 and iFin2-iBeg2:
2369                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2370                elif name in dependentVars:
2371                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2372                    if Ka2 and iFin2-iBeg2:
2373                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2374            for name in sigDict:
2375                if name in varylist:
2376                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2377                    if Ka2 and iFin2-iBeg2:
2378                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2379                elif name in dependentVars:
2380                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2381                    if Ka2 and iFin2-iBeg2:
2382                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2383            for name in ['BabA','BabU']:
2384                if refl[9+im]:
2385                    if phfx+name in varylist:
2386                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2387                        if Ka2 and iFin2-iBeg2:
2388                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2389                    elif phfx+name in dependentVars:                   
2390                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2391                        if Ka2 and iFin2-iBeg2:
2392                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2393            if not Phase['General'].get('doPawley'):
2394                #do atom derivatives -  for RB,F,X & U so far             
2395                corr = dervDict['int']/refl[9+im]
2396                if Ka2 and iFin2-iBeg2:
2397                    corr2 = dervDict2['int']/refl[9+im]
2398                for name in varylist+dependentVars:
2399                    if '::RBV;' in name:
2400                        pass
2401                    else:
2402                        try:
2403                            aname = name.split(pfx)[1][:2]
2404                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2405                        except IndexError:
2406                            continue
2407                    if name in varylist:
2408                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2409                        if Ka2 and iFin2-iBeg2:
2410                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2411                    elif name in dependentVars:
2412                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2413                        if Ka2 and iFin2-iBeg2:
2414                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2415    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2416    # now process derivatives in constraints
2417    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2418    return dMdv
2419   
2420def UserRejectHKL(ref,im,userReject):
2421    if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2422        return False
2423    elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2424        return False
2425    elif ref[11+im] < userReject['MinExt']:
2426        return False
2427    elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2428        return False
2429    return True
2430   
2431def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2432    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2433    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2434    are returned to either dervRefine or HessRefine.
2435
2436    :returns:
2437    '''
2438    nobs = Histogram['Residuals']['Nobs']
2439    hId = Histogram['hId']
2440    hfx = ':%d:'%(hId)
2441    pfx = '%d::'%(Phase['pId'])
2442    phfx = '%d:%d:'%(Phase['pId'],hId)
2443    SGData = Phase['General']['SGData']
2444    im = 0
2445    if Phase['General']['Type'] in ['modulated','magnetic']:
2446        SSGData = Phase['General']['SSGData']
2447        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2448        im = 1  #offset in SS reflection list
2449    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2450    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2451    refDict = Histogram['Data']
2452    if im:
2453        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2454    else:
2455        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2456    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2457    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2458    dependentVars = G2mv.GetDependentVars()
2459    depDerivDict = {}
2460    for j in dependentVars:
2461        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2462    wdf = np.zeros(len(refDict['RefList']))
2463    if calcControls['F**2']:
2464        for iref,ref in enumerate(refDict['RefList']):
2465            if ref[6+im] > 0:
2466                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2467                w = 1.0/ref[6+im]
2468                if ref[3+im] > 0:
2469                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2470                    for j,var in enumerate(varylist):
2471                        if var in dFdvDict:
2472                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2473                    for var in dependentVars:
2474                        if var in dFdvDict:
2475                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2476                    if phfx+'Scale' in varylist:
2477                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2478                    elif phfx+'Scale' in dependentVars:
2479                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK
2480                    for item in ['Ep','Es','Eg']:
2481                        if phfx+item in varylist and phfx+item in dervDict:
2482                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2483                        elif phfx+item in dependentVars and phfx+item in dervDict:
2484                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2485                    for item in ['BabA','BabU']:
2486                        if phfx+item in varylist:
2487                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2488                        elif phfx+item in dependentVars:
2489                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2490    else:   #F refinement
2491        for iref,ref in enumerate(refDict['RefList']):
2492            if ref[5+im] > 0.:
2493                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2494                Fo = np.sqrt(ref[5+im])
2495                Fc = np.sqrt(ref[7+im])
2496                w = 1.0/ref[6+im]
2497                if ref[3+im] > 0:
2498                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
2499                    for j,var in enumerate(varylist):
2500                        if var in dFdvDict:
2501                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2502                    for var in dependentVars:
2503                        if var in dFdvDict:
2504                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2505                    if phfx+'Scale' in varylist:
2506                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2507                    elif phfx+'Scale' in dependentVars:
2508                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK                   
2509                    for item in ['Ep','Es','Eg']:   #OK!
2510                        if phfx+item in varylist and phfx+item in dervDict:
2511                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] 
2512                        elif phfx+item in dependentVars and phfx+item in dervDict:
2513                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2514                    for item in ['BabA','BabU']:
2515                        if phfx+item in varylist:
2516                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2517                        elif phfx+item in dependentVars:
2518                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2519    return dMdvh,depDerivDict,wdf
2520   
2521
2522def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2523    '''Loop over histograms and compute derivatives of the fitting
2524    model (M) with respect to all parameters.  Results are returned in
2525    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2526    n is the number of parameters and m is the number of data
2527    points. This can exceed memory when m gets large. This routine is
2528    used when refinement derivatives are selected as "analtytic
2529    Jacobian" in Controls.
2530
2531    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2532    '''
2533    parmDict.update(zip(varylist,values))
2534    G2mv.Dict2Map(parmDict,varylist)
2535    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2536    nvar = len(varylist)
2537    dMdv = np.empty(0)
2538    histoList = Histograms.keys()
2539    histoList.sort()
2540    for histogram in histoList:
2541        if 'PWDR' in histogram[:4]:
2542            Histogram = Histograms[histogram]
2543            hId = Histogram['hId']
2544            hfx = ':%d:'%(hId)
2545            wtFactor = calcControls[hfx+'wtFactor']
2546            Limits = calcControls[hfx+'Limits']
2547            x,y,w,yc,yb,yd = Histogram['Data']
2548            xB = np.searchsorted(x,Limits[0])
2549            xF = np.searchsorted(x,Limits[1])
2550            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2551                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2552        elif 'HKLF' in histogram[:4]:
2553            Histogram = Histograms[histogram]
2554            phase = Histogram['Reflection Lists']
2555            Phase = Phases[phase]
2556            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2557            hfx = ':%d:'%(Histogram['hId'])
2558            wtFactor = calcControls[hfx+'wtFactor']
2559            # now process derivatives in constraints
2560            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2561        else:
2562            continue        #skip non-histogram entries
2563        if len(dMdv):
2564            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2565        else:
2566            dMdv = np.sqrt(wtFactor)*dMdvh
2567           
2568    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2569    if np.any(pVals):
2570        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2571        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2572       
2573    return dMdv
2574
2575def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2576    '''Loop over histograms and compute derivatives of the fitting
2577    model (M) with respect to all parameters.  For each histogram, the
2578    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2579    number of parameters and m is the number of data points *in the
2580    histogram*. The (n by n) Hessian is computed from each Jacobian
2581    and it is returned.  This routine is used when refinement
2582    derivatives are selected as "analtytic Hessian" in Controls.
2583
2584    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2585    '''
2586    parmDict.update(zip(varylist,values))
2587    G2mv.Dict2Map(parmDict,varylist)
2588    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2589    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2590    nvar = len(varylist)
2591    Hess = np.empty(0)
2592    histoList = Histograms.keys()
2593    histoList.sort()
2594    for histogram in histoList:
2595        if 'PWDR' in histogram[:4]:
2596            Histogram = Histograms[histogram]
2597            hId = Histogram['hId']
2598            hfx = ':%d:'%(hId)
2599            wtFactor = calcControls[hfx+'wtFactor']
2600            Limits = calcControls[hfx+'Limits']
2601            x,y,w,yc,yb,yd = Histogram['Data']
2602            W = wtFactor*w
2603            dy = y-yc
2604            xB = np.searchsorted(x,Limits[0])
2605            xF = np.searchsorted(x,Limits[1])
2606            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2607                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2608            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2609            Dy = dy[xB:xF][np.newaxis,:]
2610            dMdvh *= Wt
2611            if dlg:
2612                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2613            if len(Hess):
2614                Hess += np.inner(dMdvh,dMdvh)
2615                dMdvh *= Wt*Dy
2616                Vec += np.sum(dMdvh,axis=1)
2617            else:
2618                Hess = np.inner(dMdvh,dMdvh)
2619                dMdvh *= Wt*Dy
2620                Vec = np.sum(dMdvh,axis=1)
2621        elif 'HKLF' in histogram[:4]:
2622            Histogram = Histograms[histogram]
2623            phase = Histogram['Reflection Lists']
2624            Phase = Phases[phase]
2625            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2626            hId = Histogram['hId']
2627            hfx = ':%d:'%(Histogram['hId'])
2628            wtFactor = calcControls[hfx+'wtFactor']
2629            # now process derivatives in constraints
2630            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2631#            print 'matrix build time: %.3f'%(time.time()-time0)
2632
2633            if dlg:
2634                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2635            if len(Hess):
2636                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2637                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2638            else:
2639                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2640                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2641        else:
2642            continue        #skip non-histogram entries
2643    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2644    if np.any(pVals):
2645        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2646        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2647        Hess += np.inner(dpdv*pWt,dpdv)
2648    return Vec,Hess
2649
2650def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2651    '''Computes the point-by-point discrepancies between every data point in every histogram
2652    and the observed value
2653   
2654    :returns: an np array of differences between observed and computed diffraction values.
2655    '''
2656    Values2Dict(parmDict, varylist, values)
2657    G2mv.Dict2Map(parmDict,varylist)
2658    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2659    M = np.empty(0)
2660    SumwYo = 0
2661    Nobs = 0
2662    Nrej = 0
2663    Next = 0
2664    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2665    histoList = Histograms.keys()
2666    histoList.sort()
2667    for histogram in histoList:
2668        if 'PWDR' in histogram[:4]:
2669            Histogram = Histograms[histogram]
2670            hId = Histogram['hId']
2671            hfx = ':%d:'%(hId)
2672            wtFactor = calcControls[hfx+'wtFactor']
2673            Limits = calcControls[hfx+'Limits']
2674            x,y,w,yc,yb,yd = Histogram['Data']
2675            yc *= 0.0                           #zero full calcd profiles
2676            yb *= 0.0
2677            yd *= 0.0
2678            xB = np.searchsorted(x,Limits[0])
2679            xF = np.searchsorted(x,Limits[1])
2680            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2681                varylist,Histogram,Phases,calcControls,pawleyLookup)
2682            yc[xB:xF] += yb[xB:xF]
2683            if not np.any(y):                   #fill dummy data
2684                rv = st.poisson(yc[xB:xF])
2685                y[xB:xF] = rv.rvs()
2686                Z = np.ones_like(yc[xB:xF])
2687                Z[1::2] *= -1
2688                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2689                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2690            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2691            W = wtFactor*w
2692            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2693            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2694            Nobs += Histogram['Residuals']['Nobs']
2695            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2696            SumwYo += Histogram['Residuals']['sumwYo']
2697            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2698            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2699            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2700            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2701            sumYB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yc[xB:xF])*ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF],0.))
2702            sumwYB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(ma.abs(y[xB:xF]-yc[xB:xF])*ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF])**2,0.))
2703            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2704            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2705            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2706            if dlg:
2707                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2708            M = np.concatenate((M,wdy))
2709#end of PWDR processing
2710        elif 'HKLF' in histogram[:4]:
2711            Histogram = Histograms[histogram]
2712            Histogram['Residuals'] = {}
2713            phase = Histogram['Reflection Lists']
2714            Phase = Phases[phase]
2715            hId = Histogram['hId']
2716            hfx = ':%d:'%(hId)
2717            wtFactor = calcControls[hfx+'wtFactor']
2718            pfx = '%d::'%(Phase['pId'])
2719            phfx = '%d:%d:'%(Phase['pId'],hId)
2720            SGData = Phase['General']['SGData']
2721            im = 0
2722            if Phase['General']['Type'] in ['modulated','magnetic']:
2723                SSGData = Phase['General']['SSGData']
2724                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2725                im = 1  #offset in SS reflection list
2726            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2727            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2728            refDict = Histogram['Data']
2729            time0 = time.time()
2730            if im:
2731                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2732            else:
2733                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2734#            print 'sf-calc time: %.3f'%(time.time()-time0)
2735            df = np.zeros(len(refDict['RefList']))
2736            sumwYo = 0
2737            sumFo = 0
2738            sumFo2 = 0
2739            sumdF = 0
2740            sumdF2 = 0
2741            nobs = 0
2742            nrej = 0
2743            next = 0
2744            if calcControls['F**2']:
2745                for i,ref in enumerate(refDict['RefList']):
2746                    if ref[6+im] > 0:
2747                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2748                        w = 1.0/ref[6+im]   # 1/sig(F^2)
2749                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
2750                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2751                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2752                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2753                            Fo = np.sqrt(ref[5+im])
2754                            sumFo += Fo
2755                            sumFo2 += ref[5+im]
2756                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2757                            sumdF2 += abs(ref[5+im]-ref[7+im])
2758                            nobs += 1
2759                            df[i] = -w*(ref[5+im]-ref[7+im])
2760                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
2761                        else:
2762                            if ref[3+im]:
2763                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2764                                nrej += 1
2765                            else:   #sp.gp.extinct
2766                                next += 1
2767            else:
2768                for i,ref in enumerate(refDict['RefList']):
2769                    if ref[5+im] > 0.:
2770                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2771                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
2772                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2773                        Fo = np.sqrt(ref[5+im])
2774                        Fc = np.sqrt(ref[7+im])
2775                        w = 2.0*Fo/ref[6+im]    # 1/sig(F)?
2776                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2777                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2778                            sumFo += Fo
2779                            sumFo2 += ref[5+im]
2780                            sumdF += abs(Fo-Fc)
2781                            sumdF2 += abs(ref[5+im]-ref[7+im])
2782                            nobs += 1
2783                            df[i] = -w*(Fo-Fc)
2784                            sumwYo += (w*Fo)**2
2785                        else:
2786                            if ref[3+im]:
2787                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2788                                nrej += 1
2789                            else:   #sp.gp.extinct
2790                                next += 1
2791            Histogram['Residuals']['Nobs'] = nobs
2792            Histogram['Residuals']['sumwYo'] = sumwYo
2793            SumwYo += sumwYo
2794            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/sumwYo)*100.)
2795            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2796            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2797            Histogram['Residuals'][phfx+'Nref'] = nobs
2798            Histogram['Residuals'][phfx+'Nrej'] = nrej
2799            Histogram['Residuals'][phfx+'Next'] = next
2800            Nobs += nobs
2801            Nrej += nrej
2802            Next += next
2803            if dlg:
2804                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2805            M = np.concatenate((M,wtFactor*df))
2806# end of HKLF processing
2807    Histograms['sumwYo'] = SumwYo
2808    Histograms['Nobs'] = Nobs
2809    Histograms['Nrej'] = Nrej
2810    Histograms['Next'] = Next
2811    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2812    if dlg:
2813        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2814        if not GoOn:
2815            parmDict['saved values'] = values
2816            dlg.Destroy()
2817            raise G2obj.G2Exception('User abort')         #Abort!!
2818    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2819    if len(pVals):
2820        pSum = np.sum(pWt*pVals**2)
2821        for name in pWsum:
2822            if pWsum[name]:
2823                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2824        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2825        Nobs += len(pVals)
2826        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2827    return M
2828                       
Note: See TracBrowser for help on using the repository browser.