source: trunk/GSASIIstrMath.py @ 2070

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

remove extraneous argument from G2mth.posFourier
work on ZigZag/Block? SS functions; function OK, derivatives not.

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