source: trunk/GSASIIstrMath.py @ 1990

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

Structure factor derivs for SS continued...

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