source: trunk/GSASIIstrMath.py @ 1970

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

Add pygauleg to pypowder.for - recompile Win2.7 & Win2.7-64
work on SS structure factor

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