source: trunk/GSASIIstrMath.py @ 2004

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

SS derivative work

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