source: trunk/GSASIIstrMath.py @ 2076

Last change on this file since 2076 was 2076, checked in by vondreele, 6 years ago

fix TOF bugs in G2index & SStructureFactor

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