source: trunk/GSASIIstrMath.py @ 1921

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

twin derivatives now OK.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 143.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-07-06 17:48:05 +0000 (Mon, 06 Jul 2015) $
8# $Author: vondreele $
9# $Revision: 1921 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1921 2015-07-06 17:48:05Z 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: 1921 $")
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 SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
634    '''
635    Compute super structure factors for all h,k,l,m for phase
636    puts the result, F^2, in each ref[8+im] 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    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
655    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
656    FFtables = calcControls['FFtables']
657    BLtables = calcControls['BLtables']
658    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
659    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
660    FF = np.zeros(len(Tdata))
661    if 'NC' in calcControls[hfx+'histType']:
662        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
663    else:
664        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
665        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
666    Uij = np.array(G2lat.U6toUij(Uijdata))
667    bij = Mast*Uij.T
668    if not len(refDict['FF']):
669        if 'N' in calcControls[hfx+'histType']:
670            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
671        else:
672            dat = G2el.getFFvalues(FFtables,0.)       
673        refDict['FF']['El'] = dat.keys()
674        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
675    for iref,refl in enumerate(refDict['RefList']):
676        if 'NT' in calcControls[hfx+'histType']:
677            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
678        fbs = np.array([0,0])
679        H = refl[:4]
680        SQ = 1./(2.*refl[4+im])**2
681        SQfactor = 4.0*SQ*twopisq
682        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
683        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
684            if 'N' in calcControls[hfx+'histType']:
685                dat = G2el.getBLvalues(BLtables)
686                refDict['FF']['FF'][iref] = dat.values()
687            else:       #'X'
688                dat = G2el.getFFvalues(FFtables,SQ)
689                refDict['FF']['FF'][iref] = dat.values()
690        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
691        FF = refDict['FF']['FF'][iref][Tindx]
692        Uniq = np.inner(H[:3],SGMT)
693        SSUniq = np.inner(H,SSGMT)
694        Phi = np.inner(H[:3],SGT)
695        SSPhi = np.inner(H,SSGT)
696        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)       
697        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
698        sinp = np.sin(phase)
699        cosp = np.cos(phase)
700        biso = -SQfactor*Uisodata
701        Tiso = np.where(biso<1.,np.exp(biso),1.0)
702        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
703        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
704        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
705        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
706        fb = np.zeros_like(fa)
707        if not SGData['SGInv']:
708            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
709        fa = fa*GfpuA-fb*GfpuB
710        fb = fb*GfpuA+fa*GfpuB
711        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
712        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
713           
714        fasq = fas**2
715        fbsq = fbs**2        #imaginary
716        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
717        refl[10+im] = atan2d(fbs[0],fas[0])
718   
719def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
720    ''' Compute structure factors for all h,k,l for phase
721    puts the result, F^2, in each ref[8] in refList
722    input:
723   
724    :param dict refDict: where
725        'RefList' list where each ref = h,k,l,m,d,...
726        'FF' dict of form factors - filed in below
727    :param np.array G:      reciprocal metric tensor
728    :param str pfx:    phase id string
729    :param dict SGData: space group info. dictionary output from SpcGroup
730    :param dict calcControls:
731    :param dict ParmDict:
732
733    '''       
734    phfx = pfx.split(':')[0]+hfx
735    ast = np.sqrt(np.diag(G))
736    Mast = twopisq*np.multiply.outer(ast,ast)
737    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
738    SGT = np.array([ops[1] for ops in SGData['SGOps']])
739    FFtables = calcControls['FFtables']
740    BLtables = calcControls['BLtables']
741    Flack = 1.0
742    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
743        Flack = 1.-2.*parmDict[phfx+'Flack']
744    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
745    TwDict = refDict.get('TwDict',{})           
746    if 'S' in calcControls[hfx+'histType']:
747        TwinLaw = calcControls[phfx+'TwinLaw']
748        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
749        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
750    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
751    FF = np.zeros(len(Tdata))
752    if 'NC' in calcControls[hfx+'histType']:
753        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
754    elif 'X' in calcControls[hfx+'histType']:
755        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
756        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
757    Uij = np.array(G2lat.U6toUij(Uijdata))
758    bij = Mast*Uij.T
759    blkSize = 100       #no. of reflections in a block
760    nRef = refDict['RefList'].shape[0]
761    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
762        if 'N' in calcControls[hfx+'histType']:
763            dat = G2el.getBLvalues(BLtables)
764            refDict['FF']['El'] = dat.keys()
765            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
766        else:       #'X'
767            dat = G2el.getFFvalues(FFtables,0.)
768            refDict['FF']['El'] = dat.keys()
769            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
770            for iref,ref in enumerate(refDict['RefList']):
771                SQ = 1./(2.*ref[4])**2
772                dat = G2el.getFFvalues(FFtables,SQ)
773                refDict['FF']['FF'][iref] *= dat.values()
774#reflection processing begins here - big arrays!
775    iBeg = 0
776    while iBeg < nRef:
777        iFin = min(iBeg+blkSize,nRef)
778        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
779        H = refl.T[:3]                          #array(blkSize,3)
780        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
781        TwMask = np.any(H,axis=-1)
782        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
783            for ir in range(blkSize):
784                iref = ir+iBeg
785                if iref in TwDict:
786                    for i in TwDict[iref]:
787                        H[ir][i] = np.array(TwDict[iref][i])*TwinInv[i]
788            TwMask = np.any(H,axis=-1)
789        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
790        SQfactor = 4.0*SQ*twopisq               #ditto prev.
791        if 'T' in calcControls[hfx+'histType']:
792            if 'P' in calcControls[hfx+'histType']:
793                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
794            else:
795                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
796            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
797            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
798        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
799        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
800        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
801        Uniq = np.inner(H,SGMT)
802        Phi = np.inner(H,SGT)
803        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
804        sinp = np.sin(phase)
805        cosp = np.cos(phase)
806        biso = -SQfactor*Uisodata[:,np.newaxis]
807        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
808        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
809        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
810        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
811        if 'T' in calcControls[hfx+'histType']:
812            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
813            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
814        else:
815            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
816            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
817        fas = np.sum(np.sum(fa,axis=-1),axis=-1)       #real sum over atoms & unique hkl
818        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag sum over atoms & uniq hkl
819        if SGData['SGInv']: #centrosymmetric; B=0
820            fbs[0] *= 0.
821        if 'P' in calcControls[hfx+'histType']:
822            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
823            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
824        else:
825            if len(TwinLaw) > 1:
826                refl.T[9] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
827                refl.T[7] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
828                    np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
829                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
830            else:
831                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
832                refl.T[7] = np.copy(refl.T[9])               
833                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
834#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
835        iBeg += blkSize
836   
837def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
838    'Needs a doc string'
839    phfx = pfx.split(':')[0]+hfx
840    ast = np.sqrt(np.diag(G))
841    Mast = twopisq*np.multiply.outer(ast,ast)
842    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
843    SGT = np.array([ops[1] for ops in SGData['SGOps']])
844    FFtables = calcControls['FFtables']
845    BLtables = calcControls['BLtables']
846    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
847    TwDict = refDict.get('TwDict',{})           
848    if 'S' in calcControls[hfx+'histType']:
849        TwinLaw = calcControls[phfx+'TwinLaw']
850        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
851        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
852    nTwin = len(TwinLaw)       
853    nRef = len(refDict['RefList'])
854    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
855    mSize = len(Mdata)
856    FF = np.zeros(len(Tdata))
857    if 'NC' in calcControls[hfx+'histType']:
858        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
859    elif 'X' in calcControls[hfx+'histType']:
860        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
861        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
862    Uij = np.array(G2lat.U6toUij(Uijdata))
863    bij = Mast*Uij.T
864    dFdvDict = {}
865    dFdfr = np.squeeze(np.zeros((nRef,nTwin,mSize)))
866    dFdx = np.squeeze(np.zeros((nRef,nTwin,mSize,3)))
867    dFdui = np.squeeze(np.zeros((nRef,nTwin,mSize)))
868    dFdua = np.squeeze(np.zeros((nRef,nTwin,mSize,6)))
869    dFdbab = np.squeeze(np.zeros((nRef,nTwin,2)))
870    dFdfl = np.squeeze(np.zeros((nRef,nTwin)))
871    dFdtw = np.zeros((nRef,nTwin))
872    Flack = 1.0
873    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
874        Flack = 1.-2.*parmDict[phfx+'Flack']
875    for iref,refl in enumerate(refDict['RefList']):
876        if 'T' in calcControls[hfx+'histType']:
877            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
878        H = np.array(refl[:3])
879        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
880        TwMask = np.any(H,axis=-1)
881        if TwinLaw.shape[0] > 1 and TwDict:
882            if iref in TwDict:
883                for i in TwDict[iref]:
884                    H[i] = np.array(TwDict[iref][i])*TwinInv[i]
885            TwMask = np.any(H,axis=-1)
886        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
887        SQfactor = 8.0*SQ*np.pi**2
888        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
889        Bab = parmDict[phfx+'BabA']*dBabdA
890        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
891        FF = refDict['FF']['FF'][iref].T[Tindx].T
892        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
893        Phi = np.inner(H,SGT)
894        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
895        sinp = np.sin(phase)
896        cosp = np.cos(phase)
897        occ = Mdata*Fdata/len(SGT)
898        biso = -SQfactor*Uisodata[:,np.newaxis]
899        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
900        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
901        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
902        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
903        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
904        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
905        fot = (FF+FP-Bab)*Tcorr
906        fotp = FPP*Tcorr       
907        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
908        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
909        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
910        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
911        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
912        fbx = np.array([fot*cosp,-fotp*sinp])
913        #sum below is over Uniq
914        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
915        dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
916        dfadui = np.sum(-SQfactor*fa,axis=-2)
917        if nTwin > 1:
918            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
919            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
920            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
921        else:
922            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
923            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
924            # array(2,nAtom,3) & array(2,nAtom,6)
925        if not SGData['SGInv']:
926            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
927            dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
928            dfbdui = np.sum(-SQfactor*fb,axis=-2)
929            if len(TwinLaw) > 1:
930                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
931                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
932            else:
933                dfadfl = np.sum(-FPP*Tcorr*sinp)
934                dfbdfl = np.sum(FPP*Tcorr*cosp)
935                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
936                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
937        else:
938            dfbdfr = np.zeros_like(dfadfr)
939            dfbdx = np.zeros_like(dfadx)
940            dfbdui = np.zeros_like(dfadui)
941            dfbdua = np.zeros_like(dfadua)
942            dfbdba = np.zeros_like(dfadba)
943            dfadfl = 0.0
944            dfbdfl = 0.0
945        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
946        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
947            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
948                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
949            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
950                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
951            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
952                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
953            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
954                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
955        else:
956            SA = fas[0]-fbs[1]
957            SB = fbs[0]+fas[1]
958            if nTwin > 1:
959                dFdfr[iref] = [2.*TwMask[it]*SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
960                    2.*TwMask[it]*SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
961                dFdx[iref] = [2.*TwMask[it]*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*TwMask[it]*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
962                dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
963                dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
964                dFdtw[iref] = TwMask*np.sum(fas,axis=0)**2+TwMask*np.sum(fbs,axis=0)**2
965            else:   #these are good for no twin single crystals
966                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
967                    2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
968                dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])    #array(nRef,nAtom,3)
969                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])   #array(nRef,nAtom)
970                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])    #array(nRef,nAtom,6)
971                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
972                dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
973                    2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
974        #loop over atoms - each dict entry is list of derivatives for all the reflections
975    if nTwin > 1:
976        for i in range(len(Mdata)):     #these all OK?
977            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
978            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
979            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
980            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
981            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
982            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
983            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
984            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
985            dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)
986            dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)
987            dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)
988    else:
989        for i in range(len(Mdata)):
990            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
991            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
992            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
993            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
994            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
995            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
996            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
997            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
998            dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
999            dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
1000            dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
1001        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1002    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1003    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1004    if nTwin > 1:
1005        for i in range(nTwin):
1006            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1007    return dFdvDict
1008   
1009def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1010    'Needs a doc string'
1011    phfx = pfx.split(':')[0]+hfx
1012    ast = np.sqrt(np.diag(G))
1013    Mast = twopisq*np.multiply.outer(ast,ast)
1014    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1015    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1016    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1017    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1018    FFtables = calcControls['FFtables']
1019    BLtables = calcControls['BLtables']
1020    nRef = len(refDict['RefList'])
1021    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1022    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1023    mSize = len(Mdata)
1024    FF = np.zeros(len(Tdata))
1025    if 'NC' in calcControls[hfx+'histType']:
1026        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1027    elif 'X' in calcControls[hfx+'histType']:
1028        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1029        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1030    Uij = np.array(G2lat.U6toUij(Uijdata))
1031    bij = Mast*Uij.T
1032    dFdvDict = {}
1033    dFdfr = np.zeros((nRef,mSize))
1034    dFdx = np.zeros((nRef,mSize,3))
1035    dFdui = np.zeros((nRef,mSize))
1036    dFdua = np.zeros((nRef,mSize,6))
1037    dFdbab = np.zeros((nRef,2))
1038    for iref,refl in enumerate(refDict['RefList']):
1039        if 'T' in calcControls[hfx+'histType']:
1040            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1041        H = np.array(refl[:4])
1042        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1043        SQfactor = 8.0*SQ*np.pi**2
1044        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1045        Bab = parmDict[phfx+'BabA']*dBabdA
1046        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1047        FF = refDict['FF']['FF'][iref].T[Tindx]
1048        Uniq = np.inner(H[:3],SGMT)
1049        SSUniq = np.inner(H,SSGMT)
1050        Phi = np.inner(H[:3],SGT)
1051        SSPhi = np.inner(H,SSGT)
1052        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1053        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1054        #need ModulationDerv here dGAdXsin, etc 
1055        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
1056        sinp = np.sin(phase)
1057        cosp = np.cos(phase)
1058        occ = Mdata*Fdata/len(Uniq)
1059        biso = -SQfactor*Uisodata
1060        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1061        HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
1062        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1063        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1064        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1065        Tcorr = Tiso*Tuij
1066        fot = (FF+FP-Bab)*occ*Tcorr
1067        fotp = FPP*occ*Tcorr
1068        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1069        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1070        GfpuA = np.swapaxes(GfpuA,1,2)
1071        GfpuB = np.swapaxes(GfpuB,1,2)
1072        fa = fa*GfpuA-fb*GfpuB
1073        fb = fb*GfpuA+fa*GfpuB
1074       
1075        fas = np.sum(np.sum(fa,axis=1),axis=1)
1076        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1077        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1078        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1079        fax = fax*GfpuA-fbx*GfpuB
1080        fbx = fbx*GfpuA+fax*GfpuB
1081        #sum below is over Uniq
1082        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1083        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1084        dfadui = np.sum(-SQfactor*fa,axis=2)
1085        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1086        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
1087        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
1088        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1089        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1090        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1091        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1092        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1093        #need dFdXsin, etc for modulations...
1094        if not SGData['SGInv']:
1095            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1096            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
1097            dfbdui = np.sum(-SQfactor*fb,axis=2)
1098            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1099            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
1100            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1101            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1102            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1103            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1104            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1105        #need dFdXsin, etc for modulations...
1106        #loop over atoms - each dict entry is list of derivatives for all the reflections
1107    for i in range(len(Mdata)):     
1108        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1109        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1110        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1111        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1112        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1113        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1114        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1115        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1116        dFdvDict[pfx+'AU12:'+str(i)] = .5*dFdua.T[3][i]
1117        dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
1118        dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
1119        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1120    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1121    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1122    return dFdvDict
1123   
1124def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1125    ''' Single crystal extinction function; returns extinction & derivative
1126    '''
1127    extCor = 1.0
1128    dervDict = {}
1129    dervCor = 1.0
1130    if calcControls[phfx+'EType'] != 'None':
1131        SQ = 1/(4.*ref[4+im]**2)
1132        if 'C' in parmDict[hfx+'Type']:           
1133            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1134        else:   #'T'
1135            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1136        if 'SXC' in parmDict[hfx+'Type']:
1137            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1138            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1139            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1140            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1141        elif 'SNT' in parmDict[hfx+'Type']:
1142            AV = 1.e7/parmDict[pfx+'Vol']**2
1143            PL = SQ
1144            PLZ = AV*ref[9+im]*ref[12+im]**2
1145        elif 'SNC' in parmDict[hfx+'Type']:
1146            AV = 1.e7/parmDict[pfx+'Vol']**2
1147            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1148            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1149           
1150        if 'Primary' in calcControls[phfx+'EType']:
1151            PLZ *= 1.5
1152        else:
1153            if 'C' in parmDict[hfx+'Type']:
1154                PLZ *= calcControls[phfx+'Tbar']
1155            else: #'T'
1156                PLZ *= ref[13+im]      #t-bar
1157        if 'Primary' in calcControls[phfx+'EType']:
1158            PLZ *= 1.5
1159            PSIG = parmDict[phfx+'Ep']
1160        elif 'I & II' in calcControls[phfx+'EType']:
1161            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1162        elif 'Type II' in calcControls[phfx+'EType']:
1163            PSIG = parmDict[phfx+'Es']
1164        else:       # 'Secondary Type I'
1165            PSIG = parmDict[phfx+'Eg']/PL
1166           
1167        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1168        AL = 0.025+0.285*cos2T
1169        BG = 0.02-0.025*cos2T
1170        BL = 0.15-0.2*(0.75-cos2T)**2
1171        if cos2T < 0.:
1172            BL = -0.45*cos2T
1173        CG = 2.
1174        CL = 2.
1175        PF = PLZ*PSIG
1176       
1177        if 'Gaussian' in calcControls[phfx+'EApprox']:
1178            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1179            extCor = np.sqrt(PF4)
1180            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1181        else:
1182            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1183            extCor = np.sqrt(PF4)
1184            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1185
1186        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
1187        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1188            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1189        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1190            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1191        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1192            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1193               
1194    return 1./extCor,dervDict,dervCor
1195   
1196def Dict2Values(parmdict, varylist):
1197    '''Use before call to leastsq to setup list of values for the parameters
1198    in parmdict, as selected by key in varylist'''
1199    return [parmdict[key] for key in varylist] 
1200   
1201def Values2Dict(parmdict, varylist, values):
1202    ''' Use after call to leastsq to update the parameter dictionary with
1203    values corresponding to keys in varylist'''
1204    parmdict.update(zip(varylist,values))
1205   
1206def GetNewCellParms(parmDict,varyList):
1207    'Needs a doc string'
1208    newCellDict = {}
1209    Anames = ['A'+str(i) for i in range(6)]
1210    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1211    for item in varyList:
1212        keys = item.split(':')
1213        if keys[2] in Ddict:
1214            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1215            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1216            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1217    return newCellDict          # is e.g. {'0::D11':A0-D11}
1218   
1219def ApplyXYZshifts(parmDict,varyList):
1220    '''
1221    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1222   
1223    :param dict parmDict: parameter dictionary
1224    :param list varyList: list of variables (not used!)
1225    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1226
1227    '''
1228    newAtomDict = {}
1229    for item in parmDict:
1230        if 'dA' in item:
1231            parm = ''.join(item.split('d'))
1232            parmDict[parm] += parmDict[item]
1233            newAtomDict[item] = [parm,parmDict[parm]]
1234    return newAtomDict
1235   
1236def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1237    'Spherical harmonics texture'
1238    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1239    if 'T' in calcControls[hfx+'histType']:
1240        tth = parmDict[hfx+'2-theta']
1241    else:
1242        tth = refl[5+im]
1243    odfCor = 1.0
1244    H = refl[:3]
1245    cell = G2lat.Gmat2cell(g)
1246    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1247    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1248    phi,beta = G2lat.CrsAng(H,cell,SGData)
1249    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1250    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1251    for item in SHnames:
1252        L,M,N = eval(item.strip('C'))
1253        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1254        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1255        Lnorm = G2lat.Lnorm(L)
1256        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1257    return odfCor
1258   
1259def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1260    'Spherical harmonics texture derivatives'
1261    if 'T' in calcControls[hfx+'histType']:
1262        tth = parmDict[hfx+'2-theta']
1263    else:
1264        tth = refl[5+im]
1265    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1266    odfCor = 1.0
1267    dFdODF = {}
1268    dFdSA = [0,0,0]
1269    H = refl[:3]
1270    cell = G2lat.Gmat2cell(g)
1271    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1272    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1273    phi,beta = G2lat.CrsAng(H,cell,SGData)
1274    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1275    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1276    for item in SHnames:
1277        L,M,N = eval(item.strip('C'))
1278        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1279        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1280        Lnorm = G2lat.Lnorm(L)
1281        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1282        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1283        for i in range(3):
1284            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1285    return odfCor,dFdODF,dFdSA
1286   
1287def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1288    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1289    if 'T' in calcControls[hfx+'histType']:
1290        tth = parmDict[hfx+'2-theta']
1291    else:
1292        tth = refl[5+im]
1293    odfCor = 1.0
1294    H = refl[:3]
1295    cell = G2lat.Gmat2cell(g)
1296    Sangls = [0.,0.,0.]
1297    if 'Bragg' in calcControls[hfx+'instType']:
1298        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1299        IFCoup = True
1300    else:
1301        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1302        IFCoup = False
1303    phi,beta = G2lat.CrsAng(H,cell,SGData)
1304    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1305    SHnames = calcControls[phfx+'SHnames']
1306    for item in SHnames:
1307        L,N = eval(item.strip('C'))
1308        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1309        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1310        Lnorm = G2lat.Lnorm(L)
1311        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1312    return np.squeeze(odfCor)
1313   
1314def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1315    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1316    if 'T' in calcControls[hfx+'histType']:
1317        tth = parmDict[hfx+'2-theta']
1318    else:
1319        tth = refl[5+im]
1320    odfCor = 1.0
1321    dFdODF = {}
1322    H = refl[:3]
1323    cell = G2lat.Gmat2cell(g)
1324    Sangls = [0.,0.,0.]
1325    if 'Bragg' in calcControls[hfx+'instType']:
1326        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1327        IFCoup = True
1328    else:
1329        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1330        IFCoup = False
1331    phi,beta = G2lat.CrsAng(H,cell,SGData)
1332    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1333    SHnames = calcControls[phfx+'SHnames']
1334    for item in SHnames:
1335        L,N = eval(item.strip('C'))
1336        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1337        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1338        Lnorm = G2lat.Lnorm(L)
1339        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1340        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
1341    return odfCor,dFdODF
1342   
1343def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1344    'March-Dollase preferred orientation correction'
1345    POcorr = 1.0
1346    MD = parmDict[phfx+'MD']
1347    if MD != 1.0:
1348        MDAxis = calcControls[phfx+'MDAxis']
1349        sumMD = 0
1350        for H in uniq:           
1351            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1352            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1353            sumMD += A**3
1354        POcorr = sumMD/len(uniq)
1355    return POcorr
1356   
1357def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1358    'Needs a doc string'
1359    POcorr = 1.0
1360    POderv = {}
1361    if calcControls[phfx+'poType'] == 'MD':
1362        MD = parmDict[phfx+'MD']
1363        MDAxis = calcControls[phfx+'MDAxis']
1364        sumMD = 0
1365        sumdMD = 0
1366        for H in uniq:           
1367            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1368            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1369            sumMD += A**3
1370            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1371        POcorr = sumMD/len(uniq)
1372        POderv[phfx+'MD'] = sumdMD/len(uniq)
1373    else:   #spherical harmonics
1374        if calcControls[phfx+'SHord']:
1375            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1376    return POcorr,POderv
1377   
1378def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1379    'Needs a doc string'
1380    if 'Debye' in calcControls[hfx+'instType']:
1381        if 'T' in calcControls[hfx+'histType']:
1382            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1383        else:
1384            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1385    else:
1386        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1387   
1388def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1389    'Needs a doc string'
1390    if 'Debye' in calcControls[hfx+'instType']:
1391        if 'T' in calcControls[hfx+'histType']:
1392            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1393        else:
1394            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1395    else:
1396        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1397       
1398def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1399    'Needs a doc string'
1400    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1401    pi2 = np.sqrt(2./np.pi)
1402    if 'T' in calcControls[hfx+'histType']:
1403        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1404        wave = refl[14+im]
1405    else:   #'C'W
1406        sth2 = sind(refl[5+im]/2.)**2
1407        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1408    c2th = 1.-2.0*sth2
1409    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1410    if 'X' in calcControls[hfx+'histType']:
1411        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1412    xfac = flv2*parmDict[phfx+'Extinction']
1413    exb = 1.0
1414    if xfac > -1.:
1415        exb = 1./np.sqrt(1.+xfac)
1416    exl = 1.0
1417    if 0 < xfac <= 1.:
1418        xn = np.array([xfac**(i+1) for i in range(6)])
1419        exl += np.sum(xn*coef)
1420    elif xfac > 1.:
1421        xfac2 = 1./np.sqrt(xfac)
1422        exl = pi2*(1.-0.125/xfac)*xfac2
1423    return exb*sth2+exl*(1.-sth2)
1424   
1425def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1426    'Needs a doc string'
1427    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1428    pi2 = np.sqrt(2./np.pi)
1429    if 'T' in calcControls[hfx+'histType']:
1430        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1431        wave = refl[14+im]
1432    else:   #'C'W
1433        sth2 = sind(refl[5+im]/2.)**2
1434        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1435    c2th = 1.-2.0*sth2
1436    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1437    if 'X' in calcControls[hfx+'histType']:
1438        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1439    xfac = flv2*parmDict[phfx+'Extinction']
1440    dbde = -500.*flv2
1441    if xfac > -1.:
1442        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
1443    dlde = 0.
1444    if 0 < xfac <= 1.:
1445        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1446        dlde = np.sum(xn*coef)
1447    elif xfac > 1.:
1448        xfac2 = 1./np.sqrt(xfac)
1449        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1450       
1451    return dbde*sth2+dlde*(1.-sth2)
1452   
1453def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1454    'Needs a doc string'    #need powder extinction!
1455    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1456    if 'X' in parmDict[hfx+'Type']:
1457        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1458    POcorr = 1.0
1459    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1460        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1461    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1462        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1463    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1464        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1465    Icorr *= POcorr
1466    AbsCorr = 1.0
1467    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1468    Icorr *= AbsCorr
1469    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1470    Icorr *= ExtCorr
1471    return Icorr,POcorr,AbsCorr,ExtCorr
1472   
1473def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1474    'Needs a doc string'    #need powder extinction derivs!
1475    dIdsh = 1./parmDict[hfx+'Scale']
1476    dIdsp = 1./parmDict[phfx+'Scale']
1477    if 'X' in parmDict[hfx+'Type']:
1478        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1479        dIdPola /= pola
1480    else:       #'N'
1481        dIdPola = 0.0
1482    dFdODF = {}
1483    dFdSA = [0,0,0]
1484    dIdPO = {}
1485    if pfx+'SHorder' in parmDict:
1486        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1487        for iSH in dFdODF:
1488            dFdODF[iSH] /= odfCor
1489        for i in range(3):
1490            dFdSA[i] /= odfCor
1491    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1492        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1493        for iPO in dIdPO:
1494            dIdPO[iPO] /= POcorr
1495    if 'T' in parmDict[hfx+'Type']:
1496        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1497        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1498    else:
1499        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1500        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1501    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1502       
1503def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1504    'Needs a doc string'
1505    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1506        costh = cosd(refl[5+im]/2.)
1507        #crystallite size
1508        if calcControls[phfx+'SizeType'] == 'isotropic':
1509            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1510        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1511            H = np.array(refl[:3])
1512            P = np.array(calcControls[phfx+'SizeAxis'])
1513            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1514            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1515            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1516        else:           #ellipsoidal crystallites
1517            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1518            H = np.array(refl[:3])
1519            lenR = G2pwd.ellipseSize(H,Sij,GB)
1520            Sgam = 1.8*wave/(np.pi*costh*lenR)
1521        #microstrain               
1522        if calcControls[phfx+'MustrainType'] == 'isotropic':
1523            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1524        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1525            H = np.array(refl[:3])
1526            P = np.array(calcControls[phfx+'MustrainAxis'])
1527            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1528            Si = parmDict[phfx+'Mustrain;i']
1529            Sa = parmDict[phfx+'Mustrain;a']
1530            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1531        else:       #generalized - P.W. Stephens model
1532            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1533            Sum = 0
1534            for i,strm in enumerate(Strms):
1535                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1536            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1537    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1538        #crystallite size
1539        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1540            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1541        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1542            H = np.array(refl[:3])
1543            P = np.array(calcControls[phfx+'SizeAxis'])
1544            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1545            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1546            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1547        else:           #ellipsoidal crystallites   #OK
1548            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1549            H = np.array(refl[:3])
1550            lenR = G2pwd.ellipseSize(H,Sij,GB)
1551            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1552        #microstrain               
1553        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1554            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1555        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1556            H = np.array(refl[:3])
1557            P = np.array(calcControls[phfx+'MustrainAxis'])
1558            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1559            Si = parmDict[phfx+'Mustrain;i']
1560            Sa = parmDict[phfx+'Mustrain;a']
1561            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1562        else:       #generalized - P.W. Stephens model  OK
1563            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1564            Sum = 0
1565            for i,strm in enumerate(Strms):
1566                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1567            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1568           
1569    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1570    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1571    sig /= ateln2
1572    return sig,gam
1573       
1574def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1575    'Needs a doc string'
1576    gamDict = {}
1577    sigDict = {}
1578    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1579        costh = cosd(refl[5+im]/2.)
1580        tanth = tand(refl[5+im]/2.)
1581        #crystallite size derivatives
1582        if calcControls[phfx+'SizeType'] == 'isotropic':
1583            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1584            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1585            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1586        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1587            H = np.array(refl[:3])
1588            P = np.array(calcControls[phfx+'SizeAxis'])
1589            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1590            Si = parmDict[phfx+'Size;i']
1591            Sa = parmDict[phfx+'Size;a']
1592            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1593            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1594            Sgam = gami*sqtrm
1595            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1596            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1597            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1598            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1599            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1600            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1601        else:           #ellipsoidal crystallites
1602            const = 1.8*wave/(np.pi*costh)
1603            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1604            H = np.array(refl[:3])
1605            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1606            Sgam = const/lenR
1607            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1608                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1609                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1610        gamDict[phfx+'Size;mx'] = Sgam
1611        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1612               
1613        #microstrain derivatives               
1614        if calcControls[phfx+'MustrainType'] == 'isotropic':
1615            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1616            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1617            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1618        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1619            H = np.array(refl[:3])
1620            P = np.array(calcControls[phfx+'MustrainAxis'])
1621            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1622            Si = parmDict[phfx+'Mustrain;i']
1623            Sa = parmDict[phfx+'Mustrain;a']
1624            gami = 0.018*Si*Sa*tanth/np.pi
1625            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1626            Mgam = gami/sqtrm
1627            dsi = -gami*Si*cosP**2/sqtrm**3
1628            dsa = -gami*Sa*sinP**2/sqtrm**3
1629            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1630            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1631            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1632            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1633        else:       #generalized - P.W. Stephens model
1634            const = 0.018*refl[4+im]**2*tanth/np.pi
1635            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1636            Sum = 0
1637            for i,strm in enumerate(Strms):
1638                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1639                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1640                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1641            Mgam = const*np.sqrt(Sum)
1642            for i in range(len(Strms)):
1643                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1644                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1645        gamDict[phfx+'Mustrain;mx'] = Mgam
1646        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1647    else:   #'T'OF - All checked & OK
1648        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1649            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1650            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1651            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1652        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1653            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1654            H = np.array(refl[:3])
1655            P = np.array(calcControls[phfx+'SizeAxis'])
1656            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1657            Si = parmDict[phfx+'Size;i']
1658            Sa = parmDict[phfx+'Size;a']
1659            gami = const/(Si*Sa)
1660            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1661            Sgam = gami*sqtrm
1662            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1663            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1664            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1665            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1666            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1667            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1668        else:           #OK  ellipsoidal crystallites
1669            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1670            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1671            H = np.array(refl[:3])
1672            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1673            Sgam = const/lenR
1674            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1675                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1676                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1677        gamDict[phfx+'Size;mx'] = Sgam  #OK
1678        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1679               
1680        #microstrain derivatives               
1681        if calcControls[phfx+'MustrainType'] == 'isotropic':
1682            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1683            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1684            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1685        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1686            H = np.array(refl[:3])
1687            P = np.array(calcControls[phfx+'MustrainAxis'])
1688            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1689            Si = parmDict[phfx+'Mustrain;i']
1690            Sa = parmDict[phfx+'Mustrain;a']
1691            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1692            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1693            Mgam = gami/sqtrm
1694            dsi = -gami*Si*cosP**2/sqtrm**3
1695            dsa = -gami*Sa*sinP**2/sqtrm**3
1696            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1697            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1698            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1699            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1700        else:       #generalized - P.W. Stephens model OK
1701            pwrs = calcControls[phfx+'MuPwrs']
1702            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1703            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1704            Sum = 0
1705            for i,strm in enumerate(Strms):
1706                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1707                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1708                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1709            Mgam = const*np.sqrt(Sum)
1710            for i in range(len(Strms)):
1711                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1712                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1713        gamDict[phfx+'Mustrain;mx'] = Mgam
1714        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1715       
1716    return sigDict,gamDict
1717       
1718def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1719    'Needs a doc string'
1720    if im:
1721        h,k,l,m = refl[:4]
1722        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1723        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1724    else:
1725        h,k,l = refl[:3]
1726        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1727    refl[4+im] = d
1728    if 'C' in calcControls[hfx+'histType']:
1729        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1730        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1731        if 'Bragg' in calcControls[hfx+'instType']:
1732            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1733                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1734        else:               #Debye-Scherrer - simple but maybe not right
1735            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1736    elif 'T' in calcControls[hfx+'histType']:
1737        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1738        #do I need sample position effects - maybe?
1739    return pos
1740
1741def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1742    'Needs a doc string'
1743    dpr = 180./np.pi
1744    if im:
1745        h,k,l,m = refl[:4]
1746        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1747        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1748    else:
1749        m = 0
1750        h,k,l = refl[:3]       
1751        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1752    dst = np.sqrt(dstsq)
1753    dsp = 1./dst
1754    if 'C' in calcControls[hfx+'histType']:
1755        pos = refl[5+im]-parmDict[hfx+'Zero']
1756        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1757        dpdw = const*dst
1758        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1759        dpdZ = 1.0
1760        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1761            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1762        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1763        if 'Bragg' in calcControls[hfx+'instType']:
1764            dpdSh = -4.*shft*cosd(pos/2.0)
1765            dpdTr = -shft*sind(pos)*100.0
1766            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1767        else:               #Debye-Scherrer - simple but maybe not right
1768            dpdXd = -shft*cosd(pos)
1769            dpdYd = -shft*sind(pos)
1770            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1771    elif 'T' in calcControls[hfx+'histType']:
1772        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1773        dpdZ = 1.0
1774        dpdDC = dsp
1775        dpdDA = dsp**2
1776        dpdDB = 1./dsp
1777        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1778            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1779        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1780           
1781def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1782    'Needs a doc string'
1783    laue = SGData['SGLaue']
1784    uniq = SGData['SGUniq']
1785    h,k,l = refl[:3]
1786    if laue in ['m3','m3m']:
1787        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1788            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1789    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1790        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1791    elif laue in ['3R','3mR']:
1792        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1793    elif laue in ['4/m','4/mmm']:
1794        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1795    elif laue in ['mmm']:
1796        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1797    elif laue in ['2/m']:
1798        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1799        if uniq == 'a':
1800            Dij += parmDict[phfx+'D23']*k*l
1801        elif uniq == 'b':
1802            Dij += parmDict[phfx+'D13']*h*l
1803        elif uniq == 'c':
1804            Dij += parmDict[phfx+'D12']*h*k
1805    else:
1806        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1807            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1808    if 'C' in calcControls[hfx+'histType']:
1809        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1810    else:
1811        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1812           
1813def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1814    'Needs a doc string'
1815    laue = SGData['SGLaue']
1816    uniq = SGData['SGUniq']
1817    h,k,l = refl[:3]
1818    if laue in ['m3','m3m']:
1819        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1820            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1821    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1822        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1823    elif laue in ['3R','3mR']:
1824        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1825    elif laue in ['4/m','4/mmm']:
1826        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1827    elif laue in ['mmm']:
1828        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1829    elif laue in ['2/m']:
1830        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1831        if uniq == 'a':
1832            dDijDict[phfx+'D23'] = k*l
1833        elif uniq == 'b':
1834            dDijDict[phfx+'D13'] = h*l
1835        elif uniq == 'c':
1836            dDijDict[phfx+'D12'] = h*k
1837    else:
1838        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1839            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1840    if 'C' in calcControls[hfx+'histType']:
1841        for item in dDijDict:
1842            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1843    else:
1844        for item in dDijDict:
1845            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1846    return dDijDict
1847   
1848def GetDij(phfx,SGData,parmDict):
1849    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1850    return G2spc.HStrainVals(HSvals,SGData)
1851               
1852def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1853    'Needs a doc string'
1854    histoList = Histograms.keys()
1855    histoList.sort()
1856    for histogram in histoList:
1857        if 'PWDR' in histogram[:4]:
1858            Histogram = Histograms[histogram]
1859            hId = Histogram['hId']
1860            hfx = ':%d:'%(hId)
1861            Limits = calcControls[hfx+'Limits']
1862            if 'C' in calcControls[hfx+'histType']:
1863                shl = max(parmDict[hfx+'SH/L'],0.0005)
1864                Ka2 = False
1865                kRatio = 0.0
1866                if hfx+'Lam1' in parmDict.keys():
1867                    Ka2 = True
1868                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1869                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1870            x,y,w,yc,yb,yd = Histogram['Data']
1871            xB = np.searchsorted(x,Limits[0])
1872            xF = np.searchsorted(x,Limits[1])
1873            ymb = np.array(y-yb)
1874            ymb = np.where(ymb,ymb,1.0)
1875            ycmb = np.array(yc-yb)
1876            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1877            refLists = Histogram['Reflection Lists']
1878            for phase in refLists:
1879                if phase not in Phases:     #skips deleted or renamed phases silently!
1880                    continue
1881                Phase = Phases[phase]
1882                im = 0
1883                if Phase['General']['Type'] in ['modulated','magnetic']:
1884                    im = 1
1885                pId = Phase['pId']
1886                phfx = '%d:%d:'%(pId,hId)
1887                refDict = refLists[phase]
1888                sumFo = 0.0
1889                sumdF = 0.0
1890                sumFosq = 0.0
1891                sumdFsq = 0.0
1892                sumInt = 0.0
1893                for refl in refDict['RefList']:
1894                    if 'C' in calcControls[hfx+'histType']:
1895                        yp = np.zeros_like(yb)
1896                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1897                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1898                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1899                        iFin2 = iFin
1900                        if not iBeg+iFin:       #peak below low limit - skip peak
1901                            continue
1902                        elif not iBeg-iFin:     #peak above high limit - done
1903                            break
1904                        elif iBeg < iFin:
1905                            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
1906                            sumInt += refl[11+im]*refl[9+im]
1907                            if Ka2:
1908                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1909                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1910                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1911                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1912                                if iFin2 > iBeg2: 
1913                                    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
1914                                    sumInt += refl[11+im]*refl[9+im]*kRatio
1915                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1916                               
1917                    elif 'T' in calcControls[hfx+'histType']:
1918                        yp = np.zeros_like(yb)
1919                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1920                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1921                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1922                        if iBeg < iFin:
1923                            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
1924                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1925                            sumInt += refl[11+im]*refl[9+im]
1926                    Fo = np.sqrt(np.abs(refl[8+im]))
1927                    Fc = np.sqrt(np.abs(refl[9]+im))
1928                    sumFo += Fo
1929                    sumFosq += refl[8+im]**2
1930                    sumdF += np.abs(Fo-Fc)
1931                    sumdFsq += (refl[8+im]-refl[9+im])**2
1932                if sumFo:
1933                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1934                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1935                else:
1936                    Histogram['Residuals'][phfx+'Rf'] = 100.
1937                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1938                Histogram['Residuals'][phfx+'sumInt'] = sumInt
1939                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1940                Histogram['Residuals']['hId'] = hId
1941        elif 'HKLF' in histogram[:4]:
1942            Histogram = Histograms[histogram]
1943            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1944               
1945def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1946    'Needs a doc string'
1947   
1948    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1949        U = parmDict[hfx+'U']
1950        V = parmDict[hfx+'V']
1951        W = parmDict[hfx+'W']
1952        X = parmDict[hfx+'X']
1953        Y = parmDict[hfx+'Y']
1954        tanPos = tand(refl[5+im]/2.0)
1955        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1956        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1957        sig = max(0.001,sig)
1958        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1959        gam = max(0.001,gam)
1960        return sig,gam
1961               
1962    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1963        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1964            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1965        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1966        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1967        sig += Ssig
1968        gam += Sgam
1969        return sig,gam
1970       
1971    def GetReflAlpBet(refl,im,hfx,parmDict):
1972        alp = parmDict[hfx+'alpha']/refl[4+im]
1973        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1974        return alp,bet
1975       
1976    hId = Histogram['hId']
1977    hfx = ':%d:'%(hId)
1978    bakType = calcControls[hfx+'bakType']
1979    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1980    yc = np.zeros_like(yb)
1981    cw = np.diff(x)
1982    cw = np.append(cw,cw[-1])
1983       
1984    if 'C' in calcControls[hfx+'histType']:   
1985        shl = max(parmDict[hfx+'SH/L'],0.002)
1986        Ka2 = False
1987        if hfx+'Lam1' in parmDict.keys():
1988            wave = parmDict[hfx+'Lam1']
1989            Ka2 = True
1990            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1991            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1992        else:
1993            wave = parmDict[hfx+'Lam']
1994    for phase in Histogram['Reflection Lists']:
1995        refDict = Histogram['Reflection Lists'][phase]
1996        if phase not in Phases:     #skips deleted or renamed phases silently!
1997            continue
1998        Phase = Phases[phase]
1999        pId = Phase['pId']
2000        pfx = '%d::'%(pId)
2001        phfx = '%d:%d:'%(pId,hId)
2002        hfx = ':%d:'%(hId)
2003        SGData = Phase['General']['SGData']
2004        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2005        im = 0
2006        if Phase['General']['Type'] in ['modulated','magnetic']:
2007            SSGData = Phase['General']['SSGData']
2008            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2009            im = 1  #offset in SS reflection list
2010            #??
2011        Dij = GetDij(phfx,SGData,parmDict)
2012        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2013        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2014        if np.any(np.diag(G)<0.):
2015            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2016        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2017        Vst = np.sqrt(nl.det(G))    #V*
2018        if not Phase['General'].get('doPawley'):
2019            time0 = time.time()
2020            if im:
2021                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2022            else:
2023                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2024#            print 'sf calc time: %.3fs'%(time.time()-time0)
2025        time0 = time.time()
2026        badPeak = False
2027        for iref,refl in enumerate(refDict['RefList']):
2028            if 'C' in calcControls[hfx+'histType']:
2029                if im:
2030                    h,k,l,m = refl[:4]
2031                else:
2032                    h,k,l = refl[:3]
2033                Uniq = np.inner(refl[:3],SGMT)
2034                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2035                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2036                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2037                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2038                refl[11+im] *= Vst*Lorenz
2039                 
2040                if Phase['General'].get('doPawley'):
2041                    try:
2042                        if im:
2043                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2044                        else:
2045                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2046                        refl[9+im] = parmDict[pInd]
2047                    except KeyError:
2048#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2049                        continue
2050                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2051                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2052                iFin = np.searchsorted(x,refl[5+im]+fmax)
2053                if not iBeg+iFin:       #peak below low limit - skip peak
2054                    continue
2055                elif not iBeg-iFin:     #peak above high limit - done
2056                    break
2057                elif iBeg > iFin:   #bad peak coeff - skip
2058                    badPeak = True
2059                    continue
2060                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
2061                if Ka2:
2062                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2063                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2064                    iBeg = np.searchsorted(x,pos2-fmin)
2065                    iFin = np.searchsorted(x,pos2+fmax)
2066                    if not iBeg+iFin:       #peak below low limit - skip peak
2067                        continue
2068                    elif not iBeg-iFin:     #peak above high limit - done
2069                        return yc,yb
2070                    elif iBeg > iFin:   #bad peak coeff - skip
2071                        continue
2072                    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
2073            elif 'T' in calcControls[hfx+'histType']:
2074                h,k,l = refl[:3]
2075                Uniq = np.inner(refl[:3],SGMT)
2076                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2077                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2078#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2079                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2080                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2081                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2082                refl[11+im] *= Vst*Lorenz
2083                if Phase['General'].get('doPawley'):
2084                    try:
2085                        if im:
2086                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2087                        else:
2088                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2089                        refl[9+im] = parmDict[pInd]
2090                    except KeyError:
2091#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2092                        continue
2093                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2094                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2095                iFin = np.searchsorted(x,refl[5+im]+fmax)
2096                if not iBeg+iFin:       #peak below low limit - skip peak
2097                    continue
2098                elif not iBeg-iFin:     #peak above high limit - done
2099                    break
2100                elif iBeg > iFin:   #bad peak coeff - skip
2101                    badPeak = True
2102                    continue
2103                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]
2104#        print 'profile calc time: %.3fs'%(time.time()-time0)
2105    if badPeak:
2106        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2107    return yc,yb
2108   
2109def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2110    'Needs a doc string'
2111   
2112    def cellVaryDerv(pfx,SGData,dpdA): 
2113        if SGData['SGLaue'] in ['-1',]:
2114            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2115                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2116        elif SGData['SGLaue'] in ['2/m',]:
2117            if SGData['SGUniq'] == 'a':
2118                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2119            elif SGData['SGUniq'] == 'b':
2120                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2121            else:
2122                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2123        elif SGData['SGLaue'] in ['mmm',]:
2124            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2125        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2126            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2127        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2128            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2129        elif SGData['SGLaue'] in ['3R', '3mR']:
2130            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2131        elif SGData['SGLaue'] in ['m3m','m3']:
2132            return [[pfx+'A0',dpdA[0]]]
2133           
2134    # create a list of dependent variables and set up a dictionary to hold their derivatives
2135    dependentVars = G2mv.GetDependentVars()
2136    depDerivDict = {}
2137    for j in dependentVars:
2138        depDerivDict[j] = np.zeros(shape=(len(x)))
2139    #print 'dependent vars',dependentVars
2140    lenX = len(x)               
2141    hId = Histogram['hId']
2142    hfx = ':%d:'%(hId)
2143    bakType = calcControls[hfx+'bakType']
2144    dMdv = np.zeros(shape=(len(varylist),len(x)))
2145    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2146    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2147        bBpos =varylist.index(hfx+'Back;0')
2148        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2149    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2150    for name in varylist:
2151        if 'Debye' in name:
2152            id = int(name.split(';')[-1])
2153            parm = name[:int(name.rindex(';'))]
2154            ip = names.index(parm)
2155            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2156    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2157    for name in varylist:
2158        if 'BkPk' in name:
2159            parm,id = name.split(';')
2160            id = int(id)
2161            if parm in names:
2162                ip = names.index(parm)
2163                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2164    cw = np.diff(x)
2165    cw = np.append(cw,cw[-1])
2166    Ka2 = False #also for TOF!
2167    if 'C' in calcControls[hfx+'histType']:   
2168        shl = max(parmDict[hfx+'SH/L'],0.002)
2169        if hfx+'Lam1' in parmDict.keys():
2170            wave = parmDict[hfx+'Lam1']
2171            Ka2 = True
2172            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2173            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2174        else:
2175            wave = parmDict[hfx+'Lam']
2176    for phase in Histogram['Reflection Lists']:
2177        refDict = Histogram['Reflection Lists'][phase]
2178        if phase not in Phases:     #skips deleted or renamed phases silently!
2179            continue
2180        Phase = Phases[phase]
2181        SGData = Phase['General']['SGData']
2182        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2183        im = 0
2184        if Phase['General']['Type'] in ['modulated','magnetic']:
2185            SSGData = Phase['General']['SSGData']
2186            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2187            im = 1  #offset in SS reflection list
2188            #??
2189        pId = Phase['pId']
2190        pfx = '%d::'%(pId)
2191        phfx = '%d:%d:'%(pId,hId)
2192        Dij = GetDij(phfx,SGData,parmDict)
2193        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2194        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2195        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2196        if not Phase['General'].get('doPawley'):
2197            time0 = time.time()
2198            if im:
2199                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2200            else:
2201                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2202#            print 'sf-derv time %.3fs'%(time.time()-time0)
2203            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2204        time0 = time.time()
2205        for iref,refl in enumerate(refDict['RefList']):
2206            if im:
2207                h,k,l,m = refl[:4]
2208            else:
2209                h,k,l = refl[:3]
2210            Uniq = np.inner(refl[:3],SGMT)
2211            if 'T' in calcControls[hfx+'histType']:
2212                wave = refl[14+im]
2213            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2214            if 'C' in calcControls[hfx+'histType']:        #CW powder
2215                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2216            else: #'T'OF
2217                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2218            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2219            iFin = np.searchsorted(x,refl[5+im]+fmax)
2220            if not iBeg+iFin:       #peak below low limit - skip peak
2221                continue
2222            elif not iBeg-iFin:     #peak above high limit - done
2223                break
2224            pos = refl[5+im]
2225            if 'C' in calcControls[hfx+'histType']:
2226                tanth = tand(pos/2.0)
2227                costh = cosd(pos/2.0)
2228                lenBF = iFin-iBeg
2229                dMdpk = np.zeros(shape=(6,lenBF))
2230                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2231                for i in range(5):
2232                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2233                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2234                if Ka2:
2235                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2236                    iBeg2 = np.searchsorted(x,pos2-fmin)
2237                    iFin2 = np.searchsorted(x,pos2+fmax)
2238                    if iBeg2-iFin2:
2239                        lenBF2 = iFin2-iBeg2
2240                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2241                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2242                        for i in range(5):
2243                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2244                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2245                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2246            else:   #'T'OF
2247                lenBF = iFin-iBeg
2248                if lenBF < 0:   #bad peak coeff
2249                    break
2250                dMdpk = np.zeros(shape=(6,lenBF))
2251                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2252                for i in range(6):
2253                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2254                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2255            if Phase['General'].get('doPawley'):
2256                dMdpw = np.zeros(len(x))
2257                try:
2258                    if im:
2259                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2260                    else:
2261                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2262                    idx = varylist.index(pIdx)
2263                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2264                    if Ka2: #not for TOF either
2265                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2266                    dMdv[idx] = dMdpw
2267                except: # ValueError:
2268                    pass
2269            if 'C' in calcControls[hfx+'histType']:
2270                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2271                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2272                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2273                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2274                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2275                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2276                    hfx+'DisplaceY':[dpdY,'pos'],}
2277                if 'Bragg' in calcControls[hfx+'instType']:
2278                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2279                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2280                else:
2281                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2282            else:   #'T'OF
2283                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2284                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2285                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2286                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2287                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2288                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2289                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2290                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2291            for name in names:
2292                item = names[name]
2293                if name in varylist:
2294                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2295                    if Ka2 and iFin2-iBeg2:
2296                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2297                elif name in dependentVars:
2298                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2299                    if Ka2 and iFin2-iBeg2:
2300                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2301            for iPO in dIdPO:
2302                if iPO in varylist:
2303                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2304                    if Ka2 and iFin2-iBeg2:
2305                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2306                elif iPO in dependentVars:
2307                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2308                    if Ka2 and iFin2-iBeg2:
2309                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2310            for i,name in enumerate(['omega','chi','phi']):
2311                aname = pfx+'SH '+name
2312                if aname in varylist:
2313                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2314                    if Ka2 and iFin2-iBeg2:
2315                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2316                elif aname in dependentVars:
2317                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2318                    if Ka2 and iFin2-iBeg2:
2319                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2320            for iSH in dFdODF:
2321                if iSH in varylist:
2322                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2323                    if Ka2 and iFin2-iBeg2:
2324                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2325                elif iSH in dependentVars:
2326                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2327                    if Ka2 and iFin2-iBeg2:
2328                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2329            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2330            for name,dpdA in cellDervNames:
2331                if name in varylist:
2332                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2333                    if Ka2 and iFin2-iBeg2:
2334                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2335                elif name in dependentVars:
2336                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2337                    if Ka2 and iFin2-iBeg2:
2338                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2339            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2340            for name in dDijDict:
2341                if name in varylist:
2342                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2343                    if Ka2 and iFin2-iBeg2:
2344                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2345                elif name in dependentVars:
2346                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2347                    if Ka2 and iFin2-iBeg2:
2348                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2349            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2350                if name in varylist:
2351                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2352                    if Ka2 and iFin2-iBeg2:
2353                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2354                elif name in dependentVars:
2355                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2356                    if Ka2 and iFin2-iBeg2:
2357                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2358            if 'C' in calcControls[hfx+'histType']:
2359                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2360            else:   #'T'OF
2361                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2362            for name in gamDict:
2363                if name in varylist:
2364                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2365                    if Ka2 and iFin2-iBeg2:
2366                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2367                elif name in dependentVars:
2368                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2369                    if Ka2 and iFin2-iBeg2:
2370                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2371            for name in sigDict:
2372                if name in varylist:
2373                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2374                    if Ka2 and iFin2-iBeg2:
2375                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2376                elif name in dependentVars:
2377                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2378                    if Ka2 and iFin2-iBeg2:
2379                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2380            for name in ['BabA','BabU']:
2381                if refl[9+im]:
2382                    if phfx+name in varylist:
2383                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2384                        if Ka2 and iFin2-iBeg2:
2385                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2386                    elif phfx+name in dependentVars:                   
2387                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2388                        if Ka2 and iFin2-iBeg2:
2389                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2390            if not Phase['General'].get('doPawley'):
2391                #do atom derivatives -  for RB,F,X & U so far             
2392                corr = dervDict['int']/refl[9+im]
2393                if Ka2 and iFin2-iBeg2:
2394                    corr2 = dervDict2['int']/refl[9+im]
2395                for name in varylist+dependentVars:
2396                    if '::RBV;' in name:
2397                        pass
2398                    else:
2399                        try:
2400                            aname = name.split(pfx)[1][:2]
2401                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2402                        except IndexError:
2403                            continue
2404                    if name in varylist:
2405                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2406                        if Ka2 and iFin2-iBeg2:
2407                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2408                    elif name in dependentVars:
2409                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2410                        if Ka2 and iFin2-iBeg2:
2411                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2412    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2413    # now process derivatives in constraints
2414    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2415    return dMdv
2416   
2417def UserRejectHKL(ref,im,userReject):
2418    if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2419        return False
2420    elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2421        return False
2422    elif ref[11+im] < userReject['MinExt']:
2423        return False
2424    elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2425        return False
2426    return True
2427   
2428def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2429    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2430    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2431    are returned to either dervRefine or HessRefine.
2432
2433    :returns:
2434    '''
2435    nobs = Histogram['Residuals']['Nobs']
2436    hId = Histogram['hId']
2437    hfx = ':%d:'%(hId)
2438    pfx = '%d::'%(Phase['pId'])
2439    phfx = '%d:%d:'%(Phase['pId'],hId)
2440    SGData = Phase['General']['SGData']
2441    im = 0
2442    if Phase['General']['Type'] in ['modulated','magnetic']:
2443        SSGData = Phase['General']['SSGData']
2444        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2445        im = 1  #offset in SS reflection list
2446    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2447    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2448    refDict = Histogram['Data']
2449    if im:
2450        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2451    else:
2452        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2453    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2454    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2455    dependentVars = G2mv.GetDependentVars()
2456    depDerivDict = {}
2457    for j in dependentVars:
2458        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2459    wdf = np.zeros(len(refDict['RefList']))
2460    if calcControls['F**2']:
2461        for iref,ref in enumerate(refDict['RefList']):
2462            if ref[6+im] > 0:
2463                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2464                w = 1.0/ref[6+im]
2465                if ref[3+im] > 0:
2466                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2467                    for j,var in enumerate(varylist):
2468                        if var in dFdvDict:
2469                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2470                    for var in dependentVars:
2471                        if var in dFdvDict:
2472                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2473                    if phfx+'Scale' in varylist:
2474                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2475                    elif phfx+'Scale' in dependentVars:
2476                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK
2477                    for item in ['Ep','Es','Eg']:
2478                        if phfx+item in varylist and phfx+item in dervDict:
2479                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2480                        elif phfx+item in dependentVars and phfx+item in dervDict:
2481                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2482                    for item in ['BabA','BabU']:
2483                        if phfx+item in varylist:
2484                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2485                        elif phfx+item in dependentVars:
2486                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2487    else:   #F refinement
2488        for iref,ref in enumerate(refDict['RefList']):
2489            if ref[5+im] > 0.:
2490                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2491                Fo = np.sqrt(ref[5+im])
2492                Fc = np.sqrt(ref[7+im])
2493                w = 1.0/ref[6+im]
2494                if ref[3+im] > 0:
2495                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
2496                    for j,var in enumerate(varylist):
2497                        if var in dFdvDict:
2498                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2499                    for var in dependentVars:
2500                        if var in dFdvDict:
2501                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2502                    if phfx+'Scale' in varylist:
2503                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2504                    elif phfx+'Scale' in dependentVars:
2505                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK                   
2506                    for item in ['Ep','Es','Eg']:   #OK!
2507                        if phfx+item in varylist and phfx+item in dervDict:
2508                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] 
2509                        elif phfx+item in dependentVars and phfx+item in dervDict:
2510                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2511                    for item in ['BabA','BabU']:
2512                        if phfx+item in varylist:
2513                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2514                        elif phfx+item in dependentVars:
2515                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2516    return dMdvh,depDerivDict,wdf
2517   
2518
2519def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2520    '''Loop over histograms and compute derivatives of the fitting
2521    model (M) with respect to all parameters.  Results are returned in
2522    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2523    n is the number of parameters and m is the number of data
2524    points. This can exceed memory when m gets large. This routine is
2525    used when refinement derivatives are selected as "analtytic
2526    Jacobian" in Controls.
2527
2528    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2529    '''
2530    parmDict.update(zip(varylist,values))
2531    G2mv.Dict2Map(parmDict,varylist)
2532    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2533    nvar = len(varylist)
2534    dMdv = np.empty(0)
2535    histoList = Histograms.keys()
2536    histoList.sort()
2537    for histogram in histoList:
2538        if 'PWDR' in histogram[:4]:
2539            Histogram = Histograms[histogram]
2540            hId = Histogram['hId']
2541            hfx = ':%d:'%(hId)
2542            wtFactor = calcControls[hfx+'wtFactor']
2543            Limits = calcControls[hfx+'Limits']
2544            x,y,w,yc,yb,yd = Histogram['Data']
2545            xB = np.searchsorted(x,Limits[0])
2546            xF = np.searchsorted(x,Limits[1])
2547            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2548                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2549        elif 'HKLF' in histogram[:4]:
2550            Histogram = Histograms[histogram]
2551            phase = Histogram['Reflection Lists']
2552            Phase = Phases[phase]
2553            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2554            hfx = ':%d:'%(Histogram['hId'])
2555            wtFactor = calcControls[hfx+'wtFactor']
2556            # now process derivatives in constraints
2557            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2558        else:
2559            continue        #skip non-histogram entries
2560        if len(dMdv):
2561            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2562        else:
2563            dMdv = np.sqrt(wtFactor)*dMdvh
2564           
2565    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2566    if np.any(pVals):
2567        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2568        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2569       
2570    return dMdv
2571
2572def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2573    '''Loop over histograms and compute derivatives of the fitting
2574    model (M) with respect to all parameters.  For each histogram, the
2575    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2576    number of parameters and m is the number of data points *in the
2577    histogram*. The (n by n) Hessian is computed from each Jacobian
2578    and it is returned.  This routine is used when refinement
2579    derivatives are selected as "analtytic Hessian" in Controls.
2580
2581    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2582    '''
2583    parmDict.update(zip(varylist,values))
2584    G2mv.Dict2Map(parmDict,varylist)
2585    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2586    #fixup H atom positions here?
2587    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2588    nvar = len(varylist)
2589    Hess = np.empty(0)
2590    histoList = Histograms.keys()
2591    histoList.sort()
2592    for histogram in histoList:
2593        if 'PWDR' in histogram[:4]:
2594            Histogram = Histograms[histogram]
2595            hId = Histogram['hId']
2596            hfx = ':%d:'%(hId)
2597            wtFactor = calcControls[hfx+'wtFactor']
2598            Limits = calcControls[hfx+'Limits']
2599            x,y,w,yc,yb,yd = Histogram['Data']
2600            W = wtFactor*w
2601            dy = y-yc
2602            xB = np.searchsorted(x,Limits[0])
2603            xF = np.searchsorted(x,Limits[1])
2604            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2605                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2606            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2607            Dy = dy[xB:xF][np.newaxis,:]
2608            dMdvh *= Wt
2609            if dlg:
2610                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2611            if len(Hess):
2612                Hess += np.inner(dMdvh,dMdvh)
2613                dMdvh *= Wt*Dy
2614                Vec += np.sum(dMdvh,axis=1)
2615            else:
2616                Hess = np.inner(dMdvh,dMdvh)
2617                dMdvh *= Wt*Dy
2618                Vec = np.sum(dMdvh,axis=1)
2619        elif 'HKLF' in histogram[:4]:
2620            Histogram = Histograms[histogram]
2621            phase = Histogram['Reflection Lists']
2622            Phase = Phases[phase]
2623            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2624            hId = Histogram['hId']
2625            hfx = ':%d:'%(Histogram['hId'])
2626            wtFactor = calcControls[hfx+'wtFactor']
2627            # now process derivatives in constraints
2628            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2629#            print 'matrix build time: %.3f'%(time.time()-time0)
2630
2631            if dlg:
2632                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2633            if len(Hess):
2634                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2635                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2636            else:
2637                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2638                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2639        else:
2640            continue        #skip non-histogram entries
2641    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2642    if np.any(pVals):
2643        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2644        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2645        Hess += np.inner(dpdv*pWt,dpdv)
2646    return Vec,Hess
2647
2648def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2649    '''Computes the point-by-point discrepancies between every data point in every histogram
2650    and the observed value
2651   
2652    :returns: an np array of differences between observed and computed diffraction values.
2653    '''
2654    Values2Dict(parmDict, varylist, values)
2655    G2mv.Dict2Map(parmDict,varylist)
2656    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2657    M = np.empty(0)
2658    SumwYo = 0
2659    Nobs = 0
2660    Nrej = 0
2661    Next = 0
2662    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2663    #fixup Hatom positions here....
2664    histoList = Histograms.keys()
2665    histoList.sort()
2666    for histogram in histoList:
2667        if 'PWDR' in histogram[:4]:
2668            Histogram = Histograms[histogram]
2669            hId = Histogram['hId']
2670            hfx = ':%d:'%(hId)
2671            wtFactor = calcControls[hfx+'wtFactor']
2672            Limits = calcControls[hfx+'Limits']
2673            x,y,w,yc,yb,yd = Histogram['Data']
2674            yc *= 0.0                           #zero full calcd profiles
2675            yb *= 0.0
2676            yd *= 0.0
2677            xB = np.searchsorted(x,Limits[0])
2678            xF = np.searchsorted(x,Limits[1])
2679            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2680                varylist,Histogram,Phases,calcControls,pawleyLookup)
2681            yc[xB:xF] += yb[xB:xF]
2682            if not np.any(y):                   #fill dummy data
2683                rv = st.poisson(yc[xB:xF])
2684                y[xB:xF] = rv.rvs()
2685                Z = np.ones_like(yc[xB:xF])
2686                Z[1::2] *= -1
2687                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2688                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2689            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2690            W = wtFactor*w
2691            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2692            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2693            Nobs += Histogram['Residuals']['Nobs']
2694            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2695            SumwYo += Histogram['Residuals']['sumwYo']
2696            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2697            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2698            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2699            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2700            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.))
2701            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.))
2702            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2703            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2704            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2705            if dlg:
2706                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2707            M = np.concatenate((M,wdy))
2708#end of PWDR processing
2709        elif 'HKLF' in histogram[:4]:
2710            Histogram = Histograms[histogram]
2711            Histogram['Residuals'] = {}
2712            phase = Histogram['Reflection Lists']
2713            Phase = Phases[phase]
2714            hId = Histogram['hId']
2715            hfx = ':%d:'%(hId)
2716            wtFactor = calcControls[hfx+'wtFactor']
2717            pfx = '%d::'%(Phase['pId'])
2718            phfx = '%d:%d:'%(Phase['pId'],hId)
2719            SGData = Phase['General']['SGData']
2720            im = 0
2721            if Phase['General']['Type'] in ['modulated','magnetic']:
2722                SSGData = Phase['General']['SSGData']
2723                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2724                im = 1  #offset in SS reflection list
2725            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2726            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2727            refDict = Histogram['Data']
2728            time0 = time.time()
2729            if im:
2730                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2731            else:
2732                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2733#            print 'sf-calc time: %.3f'%(time.time()-time0)
2734            df = np.zeros(len(refDict['RefList']))
2735            sumwYo = 0
2736            sumFo = 0
2737            sumFo2 = 0
2738            sumdF = 0
2739            sumdF2 = 0
2740            nobs = 0
2741            nrej = 0
2742            next = 0
2743            if calcControls['F**2']:
2744                for i,ref in enumerate(refDict['RefList']):
2745                    if ref[6+im] > 0:
2746                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2747                        w = 1.0/ref[6+im]   # 1/sig(F^2)
2748                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]  #correct Fc^2 for extinction
2749                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2750                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2751                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2752                            Fo = np.sqrt(ref[5+im])
2753                            sumFo += Fo
2754                            sumFo2 += ref[5+im]
2755                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2756                            sumdF2 += abs(ref[5+im]-ref[7+im])
2757                            nobs += 1
2758                            df[i] = -w*(ref[5+im]-ref[7+im])
2759                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
2760                        else:
2761                            if ref[3+im]:
2762                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2763                                nrej += 1
2764                            else:   #sp.gp.extinct
2765                                next += 1
2766            else:
2767                for i,ref in enumerate(refDict['RefList']):
2768                    if ref[5+im] > 0.:
2769                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2770                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]    #correct Fc^2 for extinction
2771                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2772                        Fo = np.sqrt(ref[5+im])
2773                        Fc = np.sqrt(ref[7+im])
2774                        w = 2.0*Fo/ref[6+im]    # 1/sig(F)?
2775                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2776                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2777                            sumFo += Fo
2778                            sumFo2 += ref[5+im]
2779                            sumdF += abs(Fo-Fc)
2780                            sumdF2 += abs(ref[5+im]-ref[7+im])
2781                            nobs += 1
2782                            df[i] = -w*(Fo-Fc)
2783                            sumwYo += (w*Fo)**2
2784                        else:
2785                            if ref[3+im]:
2786                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2787                                nrej += 1
2788                            else:   #sp.gp.extinct
2789                                next += 1
2790            Histogram['Residuals']['Nobs'] = nobs
2791            Histogram['Residuals']['sumwYo'] = sumwYo
2792            SumwYo += sumwYo
2793            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/sumwYo)*100.)
2794            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2795            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2796            Histogram['Residuals'][phfx+'Nref'] = nobs
2797            Histogram['Residuals'][phfx+'Nrej'] = nrej
2798            Histogram['Residuals'][phfx+'Next'] = next
2799            Nobs += nobs
2800            Nrej += nrej
2801            Next += next
2802            if dlg:
2803                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2804            M = np.concatenate((M,wtFactor*df))
2805# end of HKLF processing
2806    Histograms['sumwYo'] = SumwYo
2807    Histograms['Nobs'] = Nobs
2808    Histograms['Nrej'] = Nrej
2809    Histograms['Next'] = Next
2810    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2811    if dlg:
2812        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2813        if not GoOn:
2814            parmDict['saved values'] = values
2815            dlg.Destroy()
2816            raise G2obj.G2Exception('User abort')         #Abort!!
2817    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2818    if len(pVals):
2819        pSum = np.sum(pWt*pVals**2)
2820        for name in pWsum:
2821            if pWsum[name]:
2822                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2823        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2824        Nobs += len(pVals)
2825        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2826    return M
2827                       
Note: See TracBrowser for help on using the repository browser.