source: trunk/GSASIIstrMath.py @ 1976

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

fix bug in plotting 3D HKLs
SS str. Fctr. mods - add site fraction & Uij modulation math

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