source: trunk/GSASIIstrMath.py @ 2041

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

enter space group (R -3 c) for Alumina image calibrant - gives better set of lines
some revision to SS derivatives
add ImageJ format (both big & little endian) to tiff importer

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