source: trunk/GSASIIstrMath.py @ 1993

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

numerical derivatives for x modulations - work but slow
show density map slices as tau changed for SS structures

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