source: trunk/GSASIIstrMath.py @ 2006

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

minor work on SS structure factors

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