source: trunk/GSASIIstrMath.py @ 1886

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

refactor StructureFactor2 & StructureFactorDerv? for twinning - deriv not complete

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