source: trunk/GSASIIstrMath.py @ 1786

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

fix anomalous dispersion calcs.
fix SC weights on derivatives

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 138.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-04-10 18:14:03 +0000 (Fri, 10 Apr 2015) $
8# $Author: vondreele $
9# $Revision: 1786 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1786 2015-04-10 18:14:03Z vondreele $
12########### SVN repository information ###################
13import time
14import math
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.optimize as so
20import scipy.stats as st
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 1786 $")
23import GSASIIElem as G2el
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import GSASIIpwd as G2pwd
27import GSASIImapvars as G2mv
28import GSASIImath as G2mth
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*math.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        name = 'SH-Pref.Ori.'
401        pWsum[name] = 0.0
402        for hist in Phases[phase]['Histograms']:
403            if hist in Histograms:
404                hId = Histograms[hist]['hId']
405                phfx = '%d:%d:'%(pId,hId)
406                if calcControls[phfx+'poType'] == 'SH':
407                    toler = calcControls[phfx+'SHtoler']
408                    wt = 1./toler**2
409                    HKLs = calcControls[phfx+'SHhkl']
410                    SHnames = calcControls[phfx+'SHnames']
411                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
412                    for i,PH in enumerate(HKLs):
413                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
414                        SH3Coef = {}
415                        for item in SHcof:
416                            L,N = eval(item.strip('C'))
417                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
418                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
419                        X = np.linspace(0,90.0,26)
420                        Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
421                        IndY = ma.nonzero(Y)
422                        for ind in IndY[0]:
423                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
424                            pVals.append(Y[ind])
425                            pWt.append(wt)
426                            pWsum[name] += wt*(Y[ind])**2
427    pWsum['PWLref'] = 0.
428    for item in varyList:
429        if 'PWLref' in item and parmDict[item] < 0.:
430            pId = int(item.split(':')[0])
431            if negWt[pId]:
432                pNames.append(item)
433                pVals.append(-parmDict[item])
434                pWt.append(negWt[pId])
435                pWsum['PWLref'] += negWt[pId]*(-parmDict[item])**2
436    pVals = np.array(pVals)
437    pWt = np.array(pWt)         #should this be np.sqrt?
438    return pNames,pVals,pWt,pWsum
439   
440def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
441    'Needs a doc string'
442    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
443    pDerv = np.zeros((len(varyList),len(pVal)))
444    for phase in Phases:
445#        if phase not in restraintDict:
446#            continue
447        pId = Phases[phase]['pId']
448        General = Phases[phase]['General']
449        cx,ct,cs,cia = General['AtomPtrs']
450        SGData = General['SGData']
451        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
452        cell = General['Cell'][1:7]
453        Amat,Bmat = G2lat.cell2AB(cell)
454        textureData = General['SH Texture']
455
456        SHkeys = textureData['SH Coeff'][1].keys()
457        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
458        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
459        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
460        sam = SamSym[textureData['Model']]
461        phaseRest = restraintDict.get(phase,{})
462        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
463            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
464            'ChemComp':'Sites','Texture':'HKLs'}
465        lasthkl = np.array([0,0,0])
466        for ip,pName in enumerate(pNames):
467            pnames = pName.split(':')
468            if pId == int(pnames[0]):
469                name = pnames[1]
470                if 'PWL' in pName:
471                    pDerv[varyList.index(pName)][ip] += 1.
472                    continue
473                elif 'SH-' in pName:
474                    continue
475                id = int(pnames[2]) 
476                itemRest = phaseRest[name]
477                if name in ['Bond','Angle','Plane','Chiral']:
478                    indx,ops,obs,esd = itemRest[names[name]][id]
479                    dNames = []
480                    for ind in indx:
481                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
482                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
483                    if name == 'Bond':
484                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
485                    elif name == 'Angle':
486                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
487                    elif name == 'Plane':
488                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
489                    elif name == 'Chiral':
490                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
491                elif name in ['Torsion','Rama']:
492                    coffDict = itemRest['Coeff']
493                    indx,ops,cofName,esd = itemRest[names[name]][id]
494                    dNames = []
495                    for ind in indx:
496                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
497                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
498                    if name == 'Torsion':
499                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
500                    else:
501                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
502                elif name == 'ChemComp':
503                    indx,factors,obs,esd = itemRest[names[name]][id]
504                    dNames = []
505                    for ind in indx:
506                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
507                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
508                        deriv = mul*factors
509                elif 'Texture' in name:
510                    deriv = []
511                    dNames = []
512                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
513                    hkl = np.array(hkl)
514                    if np.any(lasthkl-hkl):
515                        PH = np.array(hkl)
516                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
517                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
518                        lasthkl = copy.copy(hkl)                       
519                    if 'unit' in name:
520                        pass
521                    else:
522                        gam = float(pnames[3])
523                        psi = float(pnames[4])
524                        for SHname in ODFln:
525                            l,m,n = eval(SHname[1:])
526                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
527                            dNames += [str(pId)+'::'+SHname]
528                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
529                for dName,drv in zip(dNames,deriv):
530                    try:
531                        ind = varyList.index(dName)
532                        pDerv[ind][ip] += drv
533                    except ValueError:
534                        pass
535       
536        lasthkl = np.array([0,0,0])
537        for ip,pName in enumerate(pNames):
538            deriv = []
539            dNames = []
540            pnames = pName.split(':')
541            if 'SH-' in pName and pId == int(pnames[0]):
542                hId = int(pnames[1])
543                phfx = '%d:%d:'%(pId,hId)
544                psi = float(pnames[4])
545                HKLs = calcControls[phfx+'SHhkl']
546                SHnames = calcControls[phfx+'SHnames']
547                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
548                hkl = np.array(HKLs[int(pnames[3])])     
549                if np.any(lasthkl-hkl):
550                    PH = np.array(hkl)
551                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
552                    SH3Coef = {}
553                    for item in SHcof:
554                        L,N = eval(item.strip('C'))
555                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
556                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
557                    lasthkl = copy.copy(hkl)                       
558                for SHname in SHnames:
559                    l,n = eval(SHname[1:])
560                    SH3name = 'C%d,0,%d'%(l,n)
561                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
562                    dNames += [phfx+SHname]
563                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
564            for dName,drv in zip(dNames,deriv):
565                try:
566                    ind = varyList.index(dName)
567                    pDerv[ind][ip] += drv
568                except ValueError:
569                    pass
570                   
571#    raise Exception
572    return pDerv
573
574################################################################################
575##### Function & derivative calculations
576################################################################################       
577                   
578def GetAtomFXU(pfx,calcControls,parmDict):
579    'Needs a doc string'
580    Natoms = calcControls['Natoms'][pfx]
581    Tdata = Natoms*[' ',]
582    Mdata = np.zeros(Natoms)
583    IAdata = Natoms*[' ',]
584    Fdata = np.zeros(Natoms)
585    FFdata = []
586    BLdata = []
587    Xdata = np.zeros((3,Natoms))
588    dXdata = np.zeros((3,Natoms))
589    Uisodata = np.zeros(Natoms)
590    Uijdata = np.zeros((6,Natoms))
591    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
592        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
593        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
594        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
595        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
596    for iatm in range(Natoms):
597        for key in keys:
598            parm = pfx+key+str(iatm)
599            if parm in parmDict:
600                keys[key][iatm] = parmDict[parm]
601    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.?
602    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
603   
604def GetAtomSSFXU(pfx,calcControls,parmDict):
605    'Needs a doc string'
606    Natoms = calcControls['Natoms'][pfx]
607    maxSSwave = calcControls['maxSSwave'][pfx]
608    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
609        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
610    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
611    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
612    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
613    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
614    waveTypes = []
615    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
616        'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
617        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
618        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
619        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
620        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
621    for iatm in range(Natoms):
622        waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
623        for key in keys:
624            for m in range(Nwave[key[0]]):
625                parm = pfx+key+str(iatm)+':%d'%(m)
626                if parm in parmDict:
627                    keys[key][m][iatm] = parmDict[parm]
628    return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()   
629   
630def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
631    ''' Not Used: here only for comparison the StructureFactor2 - faster version
632    Compute structure factors for all h,k,l for phase
633    puts the result, F^2, in each ref[8] in refList
634    input:
635   
636    :param dict refDict: where
637        'RefList' list where each ref = h,k,l,m,d,...
638        'FF' dict of form factors - filed in below
639    :param np.array G:      reciprocal metric tensor
640    :param str pfx:    phase id string
641    :param dict SGData: space group info. dictionary output from SpcGroup
642    :param dict calcControls:
643    :param dict ParmDict:
644
645    '''       
646    phfx = pfx.split(':')[0]+hfx
647    ast = np.sqrt(np.diag(G))
648    Mast = twopisq*np.multiply.outer(ast,ast)
649    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
650    SGT = np.array([ops[1] for ops in SGData['SGOps']])
651    FFtables = calcControls['FFtables']
652    BLtables = calcControls['BLtables']
653    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
654    FF = np.zeros(len(Tdata))
655    if 'NC' in calcControls[hfx+'histType']:
656        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
657    else:
658        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
659        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
660    Uij = np.array(G2lat.U6toUij(Uijdata))
661    bij = Mast*Uij.T
662    if not len(refDict['FF']):
663        if 'N' in calcControls[hfx+'histType']:
664            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
665        else:
666            dat = G2el.getFFvalues(FFtables,0.)       
667        refDict['FF']['El'] = dat.keys()
668        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
669    for iref,refl in enumerate(refDict['RefList']):
670        if 'NT' in calcControls[hfx+'histType']:
671            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14])
672        fbs = np.array([0,0])
673        H = refl[:3]
674        SQ = 1./(2.*refl[4])**2
675        SQfactor = 4.0*SQ*twopisq
676        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
677        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
678            if 'N' in calcControls[hfx+'histType']:
679                dat = G2el.getBLvalues(BLtables)
680                refDict['FF']['FF'][iref] = dat.values()
681            else:       #'X'
682                dat = G2el.getFFvalues(FFtables,SQ)
683                refDict['FF']['FF'][iref] = dat.values()
684        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
685        FF = refDict['FF']['FF'][iref][Tindx]
686        Uniq = np.inner(H,SGMT)
687        Phi = np.inner(H,SGT)
688        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
689        sinp = np.sin(phase)
690        cosp = np.cos(phase)
691        biso = -SQfactor*Uisodata
692        Tiso = np.where(biso<1.,np.exp(biso),1.0)
693        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
694        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
695        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
696        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
697        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
698        if not SGData['SGInv']:
699            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
700            fbs = np.sum(np.sum(fb,axis=1),axis=1)
701        fasq = fas**2
702        fbsq = fbs**2        #imaginary
703        refl[9] = np.sum(fasq)+np.sum(fbsq)
704        refl[10] = atan2d(fbs[0],fas[0])
705   
706def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
707    '''
708    Compute super structure factors for all h,k,l,m for phase
709    puts the result, F^2, in each ref[8+im] in refList
710    input:
711   
712    :param dict refDict: where
713        'RefList' list where each ref = h,k,l,m,d,...
714        'FF' dict of form factors - filed in below
715    :param np.array G:      reciprocal metric tensor
716    :param str pfx:    phase id string
717    :param dict SGData: space group info. dictionary output from SpcGroup
718    :param dict calcControls:
719    :param dict ParmDict:
720
721    '''       
722    phfx = pfx.split(':')[0]+hfx
723    ast = np.sqrt(np.diag(G))
724    Mast = twopisq*np.multiply.outer(ast,ast)
725    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
726    SGT = np.array([ops[1] for ops in SGData['SGOps']])
727    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
728    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
729    FFtables = calcControls['FFtables']
730    BLtables = calcControls['BLtables']
731    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
732    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
733    FF = np.zeros(len(Tdata))
734    if 'NC' in calcControls[hfx+'histType']:
735        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
736    else:
737        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
738        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
739    Uij = np.array(G2lat.U6toUij(Uijdata))
740    bij = Mast*Uij.T
741    if not len(refDict['FF']):
742        if 'N' in calcControls[hfx+'histType']:
743            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
744        else:
745            dat = G2el.getFFvalues(FFtables,0.)       
746        refDict['FF']['El'] = dat.keys()
747        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
748    for iref,refl in enumerate(refDict['RefList']):
749        if 'NT' in calcControls[hfx+'histType']:
750            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
751        fbs = np.array([0,0])
752        H = refl[:4]
753        SQ = 1./(2.*refl[4+im])**2
754        SQfactor = 4.0*SQ*twopisq
755        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
756        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
757            if 'N' in calcControls[hfx+'histType']:
758                dat = G2el.getBLvalues(BLtables)
759                refDict['FF']['FF'][iref] = dat.values()
760            else:       #'X'
761                dat = G2el.getFFvalues(FFtables,SQ)
762                refDict['FF']['FF'][iref] = dat.values()
763        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
764        FF = refDict['FF']['FF'][iref][Tindx]
765        Uniq = np.inner(H[:3],SGMT)
766        SSUniq = np.inner(H,SSGMT)
767        Phi = np.inner(H[:3],SGT)
768        SSPhi = np.inner(H,SSGT)
769        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)       
770        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
771        sinp = np.sin(phase)
772        cosp = np.cos(phase)
773        biso = -SQfactor*Uisodata
774        Tiso = np.where(biso<1.,np.exp(biso),1.0)
775        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
776        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
777        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
778        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
779        fb = np.zeros_like(fa)
780        if not SGData['SGInv']:
781            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
782        fa = fa*GfpuA-fb*GfpuB
783        fb = fb*GfpuA+fa*GfpuB
784        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
785        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
786           
787        fasq = fas**2
788        fbsq = fbs**2        #imaginary
789        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
790        refl[10+im] = atan2d(fbs[0],fas[0])
791   
792def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
793    ''' Compute structure factors for all h,k,l for phase
794    puts the result, F^2, in each ref[8] in refList
795    input:
796   
797    :param dict refDict: where
798        'RefList' list where each ref = h,k,l,m,d,...
799        'FF' dict of form factors - filed in below
800    :param np.array G:      reciprocal metric tensor
801    :param str pfx:    phase id string
802    :param dict SGData: space group info. dictionary output from SpcGroup
803    :param dict calcControls:
804    :param dict ParmDict:
805
806    '''       
807    phfx = pfx.split(':')[0]+hfx
808    ast = np.sqrt(np.diag(G))
809    Mast = twopisq*np.multiply.outer(ast,ast)
810    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
811    SGT = np.array([ops[1] for ops in SGData['SGOps']])
812    FFtables = calcControls['FFtables']
813    BLtables = calcControls['BLtables']
814    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
815    FF = np.zeros(len(Tdata))
816    if 'NC' in calcControls[hfx+'histType']:
817        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
818    elif 'X' in calcControls[hfx+'histType']:
819        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
820        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
821    Uij = np.array(G2lat.U6toUij(Uijdata))
822    bij = Mast*Uij.T
823    blkSize = 100       #no. of reflections in a block
824    nRef = refDict['RefList'].shape[0]
825    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
826        if 'N' in calcControls[hfx+'histType']:
827            dat = G2el.getBLvalues(BLtables)
828            refDict['FF']['El'] = dat.keys()
829            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
830        else:       #'X'
831            dat = G2el.getFFvalues(FFtables,0.)
832            refDict['FF']['El'] = dat.keys()
833            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
834            for iref,ref in enumerate(refDict['RefList']):
835                SQ = 1./(2.*ref[4])**2
836                dat = G2el.getFFvalues(FFtables,SQ)
837                refDict['FF']['FF'][iref] *= dat.values()
838#reflection processing begins here - big arrays!
839    iBeg = 0           
840    while iBeg < nRef:
841        iFin = min(iBeg+blkSize,nRef)
842        refl = refDict['RefList'][iBeg:iFin]
843        H = refl.T[:3]
844        SQ = 1./(2.*refl.T[4])**2
845        SQfactor = 4.0*SQ*twopisq
846        if 'T' in calcControls[hfx+'histType']:
847            if 'P' in calcControls[hfx+'histType']:
848                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
849            else:
850                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
851            FP = np.repeat(FP.T,len(SGT),axis=0)
852            FPP = np.repeat(FPP.T,len(SGT),axis=0)
853        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
854        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
855        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
856        Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
857        Phi = np.inner(H.T,SGT).flatten()
858        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
859        sinp = np.sin(phase)
860        cosp = np.cos(phase)
861        biso = -SQfactor*Uisodata[:,np.newaxis]
862        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
863        HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
864        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
865        Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
866        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
867        fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))   #real A,-b
868        fas = np.sum(np.sum(fa,axis=2),axis=2)        #real sum over atoms & unique hkl
869        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
870        fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))   #imag -B,+a
871        fbs = np.sum(np.sum(fb,axis=2),axis=2)  #imag sum over atoms & uniq hkl
872        if SGData['SGInv']: #centrosymmetric; B=0
873            fbs[0] *= 0.
874        if 'P' in calcControls[hfx+'histType']:
875            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
876        else:
877            refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
878        refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"?
879#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
880        iBeg += blkSize
881   
882def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
883    'Needs a doc string'
884    phfx = pfx.split(':')[0]+hfx
885    ast = np.sqrt(np.diag(G))
886    Mast = twopisq*np.multiply.outer(ast,ast)
887    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
888    SGT = np.array([ops[1] for ops in SGData['SGOps']])
889    FFtables = calcControls['FFtables']
890    BLtables = calcControls['BLtables']
891    nRef = len(refDict['RefList'])
892    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
893    mSize = len(Mdata)
894    FF = np.zeros(len(Tdata))
895    if 'NC' in calcControls[hfx+'histType']:
896        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
897    elif 'X' in calcControls[hfx+'histType']:
898        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
899        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
900    Uij = np.array(G2lat.U6toUij(Uijdata))
901    bij = Mast*Uij.T
902    dFdvDict = {}
903    dFdfr = np.zeros((nRef,mSize))
904    dFdx = np.zeros((nRef,mSize,3))
905    dFdui = np.zeros((nRef,mSize))
906    dFdua = np.zeros((nRef,mSize,6))
907    dFdbab = np.zeros((nRef,2))
908    for iref,refl in enumerate(refDict['RefList']):
909        if 'T' in calcControls[hfx+'histType']:
910            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
911        H = np.array(refl[:3])
912        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
913        SQfactor = 8.0*SQ*np.pi**2
914        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
915        Bab = parmDict[phfx+'BabA']*dBabdA
916        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
917        FF = refDict['FF']['FF'][iref].T[Tindx]
918        Uniq = np.inner(H,SGMT)
919        Phi = np.inner(H,SGT)
920        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
921        sinp = np.sin(phase)
922        cosp = np.cos(phase)
923        occ = Mdata*Fdata/len(Uniq)
924        biso = -SQfactor*Uisodata
925        Tiso = np.where(biso<1.,np.exp(biso),1.0)
926        HbH = -np.inner(H,np.inner(bij,H))
927        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
928        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
929        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
930        Tcorr = Tiso*Tuij
931        fot = (FF+FP-Bab)*occ*Tcorr
932        fotp = FPP*occ*Tcorr
933        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
934        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
935       
936        fas = np.sum(np.sum(fa,axis=1),axis=1)
937        fbs = np.sum(np.sum(fb,axis=1),axis=1)
938        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
939        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
940        #sum below is over Uniq
941        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
942        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
943        dfadui = np.sum(-SQfactor*fa,axis=2)
944        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
945        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
946        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
947        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
948        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
949        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
950        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
951        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
952        if not SGData['SGInv']:
953            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
954            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
955            dfbdui = np.sum(-SQfactor*fb,axis=2)
956            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
957            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
958            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
959            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
960            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
961            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
962            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
963        #loop over atoms - each dict entry is list of derivatives for all the reflections
964    for i in range(len(Mdata)):     
965        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
966        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
967        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
968        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
969        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
970        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
971        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
972        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
973        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
974        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
975        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
976    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
977    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
978    return dFdvDict
979   
980def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
981    'Needs a doc string'
982    phfx = pfx.split(':')[0]+hfx
983    ast = np.sqrt(np.diag(G))
984    Mast = twopisq*np.multiply.outer(ast,ast)
985    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
986    SGT = np.array([ops[1] for ops in SGData['SGOps']])
987    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
988    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
989    FFtables = calcControls['FFtables']
990    BLtables = calcControls['BLtables']
991    nRef = len(refDict['RefList'])
992    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
993    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
994    mSize = len(Mdata)
995    FF = np.zeros(len(Tdata))
996    if 'NC' in calcControls[hfx+'histType']:
997        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
998    elif 'X' in calcControls[hfx+'histType']:
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    dFdvDict = {}
1004    dFdfr = np.zeros((nRef,mSize))
1005    dFdx = np.zeros((nRef,mSize,3))
1006    dFdui = np.zeros((nRef,mSize))
1007    dFdua = np.zeros((nRef,mSize,6))
1008    dFdbab = np.zeros((nRef,2))
1009    for iref,refl in enumerate(refDict['RefList']):
1010        if 'T' in calcControls[hfx+'histType']:
1011            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1012        H = np.array(refl[:4])
1013        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1014        SQfactor = 8.0*SQ*np.pi**2
1015        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1016        Bab = parmDict[phfx+'BabA']*dBabdA
1017        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1018        FF = refDict['FF']['FF'][iref].T[Tindx]
1019        Uniq = np.inner(H[:3],SGMT)
1020        SSUniq = np.inner(H,SSGMT)
1021        Phi = np.inner(H[:3],SGT)
1022        SSPhi = np.inner(H,SSGT)
1023        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1024        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1025        #need ModulationDerv here dGAdXsin, etc 
1026        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
1027        sinp = np.sin(phase)
1028        cosp = np.cos(phase)
1029        occ = Mdata*Fdata/len(Uniq)
1030        biso = -SQfactor*Uisodata
1031        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1032        HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
1033        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1034        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1035        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1036        Tcorr = Tiso*Tuij
1037        fot = (FF+FP-Bab)*occ*Tcorr
1038        fotp = FPP*occ*Tcorr
1039        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1040        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1041        GfpuA = np.swapaxes(GfpuA,1,2)
1042        GfpuB = np.swapaxes(GfpuB,1,2)
1043        fa = fa*GfpuA-fb*GfpuB
1044        fb = fb*GfpuA+fa*GfpuB
1045       
1046        fas = np.sum(np.sum(fa,axis=1),axis=1)
1047        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1048        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1049        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1050        fax = fax*GfpuA-fbx*GfpuB
1051        fbx = fbx*GfpuA+fax*GfpuB
1052        #sum below is over Uniq
1053        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1054        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1055        dfadui = np.sum(-SQfactor*fa,axis=2)
1056        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1057        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
1058        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
1059        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1060        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1061        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1062        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1063        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1064        #need dFdXsin, etc for modulations...
1065        if not SGData['SGInv']:
1066            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1067            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
1068            dfbdui = np.sum(-SQfactor*fb,axis=2)
1069            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1070            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
1071            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1072            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1073            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1074            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1075            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1076        #need dFdXsin, etc for modulations...
1077        #loop over atoms - each dict entry is list of derivatives for all the reflections
1078    for i in range(len(Mdata)):     
1079        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1080        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1081        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1082        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1083        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1084        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1085        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1086        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1087        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1088        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1089        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1090        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1091    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
1092    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
1093    return dFdvDict
1094   
1095def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1096    ''' Single crystal extinction function; returns extinction & derivative
1097    '''
1098    extCor = 1.0
1099    dervDict = {}
1100    if calcControls[phfx+'EType'] != 'None':
1101        SQ = 1/(4.*ref[4+im]**2)
1102        if 'C' in parmDict[hfx+'Type']:           
1103            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1104        else:   #'T'
1105            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1106        if 'SXC' in parmDict[hfx+'Type']:
1107            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1108            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1109            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1110            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1111        elif 'SNT' in parmDict[hfx+'Type']:
1112            AV = 1.e7/parmDict[pfx+'Vol']**2
1113            PL = SQ
1114            PLZ = AV*ref[9+im]*ref[12+im]**2
1115        elif 'SNC' in parmDict[hfx+'Type']:
1116            AV = 1.e7/parmDict[pfx+'Vol']**2
1117            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1118            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1119           
1120        if 'Primary' in calcControls[phfx+'EType']:
1121            PLZ *= 1.5
1122        else:
1123            if 'C' in parmDict[hfx+'Type']:
1124                PLZ *= calcControls[phfx+'Tbar']
1125            else: #'T'
1126                PLZ *= ref[13+im]      #t-bar
1127        if 'Primary' in calcControls[phfx+'EType']:
1128            PLZ *= 1.5
1129            PSIG = parmDict[phfx+'Ep']
1130        elif 'I & II' in calcControls[phfx+'EType']:
1131            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1132        elif 'Type II' in calcControls[phfx+'EType']:
1133            PSIG = parmDict[phfx+'Es']
1134        else:       # 'Secondary Type I'
1135            PSIG = parmDict[phfx+'Eg']/PL
1136           
1137        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1138        AL = 0.025+0.285*cos2T
1139        BG = 0.02-0.025*cos2T
1140        BL = 0.15-0.2*(0.75-cos2T)**2
1141        if cos2T < 0.:
1142            BL = -0.45*cos2T
1143        CG = 2.
1144        CL = 2.
1145        PF = PLZ*PSIG
1146       
1147        if 'Gaussian' in calcControls[phfx+'EApprox']:
1148            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1149            extCor = np.sqrt(PF4)
1150            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1151        else:
1152            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1153            extCor = np.sqrt(PF4)
1154            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1155
1156        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1157            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1158        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1159            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1160        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1161            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1162               
1163    return 1./extCor,dervDict
1164   
1165def Dict2Values(parmdict, varylist):
1166    '''Use before call to leastsq to setup list of values for the parameters
1167    in parmdict, as selected by key in varylist'''
1168    return [parmdict[key] for key in varylist] 
1169   
1170def Values2Dict(parmdict, varylist, values):
1171    ''' Use after call to leastsq to update the parameter dictionary with
1172    values corresponding to keys in varylist'''
1173    parmdict.update(zip(varylist,values))
1174   
1175def GetNewCellParms(parmDict,varyList):
1176    'Needs a doc string'
1177    newCellDict = {}
1178    Anames = ['A'+str(i) for i in range(6)]
1179    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1180    for item in varyList:
1181        keys = item.split(':')
1182        if keys[2] in Ddict:
1183            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1184            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1185            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1186    return newCellDict          # is e.g. {'0::D11':A0-D11}
1187   
1188def ApplyXYZshifts(parmDict,varyList):
1189    '''
1190    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1191   
1192    :param dict parmDict: parameter dictionary
1193    :param list varyList: list of variables (not used!)
1194    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1195
1196    '''
1197    newAtomDict = {}
1198    for item in parmDict:
1199        if 'dA' in item:
1200            parm = ''.join(item.split('d'))
1201            parmDict[parm] += parmDict[item]
1202            newAtomDict[item] = [parm,parmDict[parm]]
1203    return newAtomDict
1204   
1205def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1206    'Spherical harmonics texture'
1207    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1208    if 'T' in calcControls[hfx+'histType']:
1209        tth = parmDict[hfx+'2-theta']
1210    else:
1211        tth = refl[5+im]
1212    odfCor = 1.0
1213    H = refl[:3]
1214    cell = G2lat.Gmat2cell(g)
1215    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1216    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1217    phi,beta = G2lat.CrsAng(H,cell,SGData)
1218    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1219    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1220    for item in SHnames:
1221        L,M,N = eval(item.strip('C'))
1222        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1223        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1224        Lnorm = G2lat.Lnorm(L)
1225        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1226    return odfCor
1227   
1228def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1229    'Spherical harmonics texture derivatives'
1230    if 'T' in calcControls[hfx+'histType']:
1231        tth = parmDict[hfx+'2-theta']
1232    else:
1233        tth = refl[5+im]
1234    FORPI = 4.0*np.pi
1235    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1236    odfCor = 1.0
1237    dFdODF = {}
1238    dFdSA = [0,0,0]
1239    H = refl[:3]
1240    cell = G2lat.Gmat2cell(g)
1241    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1242    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1243    phi,beta = G2lat.CrsAng(H,cell,SGData)
1244    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1245    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1246    for item in SHnames:
1247        L,M,N = eval(item.strip('C'))
1248        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1249        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1250        Lnorm = G2lat.Lnorm(L)
1251        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1252        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1253        for i in range(3):
1254            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1255    return odfCor,dFdODF,dFdSA
1256   
1257def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1258    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1259    if 'T' in calcControls[hfx+'histType']:
1260        tth = parmDict[hfx+'2-theta']
1261    else:
1262        tth = refl[5+im]
1263    odfCor = 1.0
1264    H = refl[:3]
1265    cell = G2lat.Gmat2cell(g)
1266    Sangl = [0.,0.,0.]
1267    if 'Bragg' in calcControls[hfx+'instType']:
1268        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1269        IFCoup = True
1270    else:
1271        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1272        IFCoup = False
1273    phi,beta = G2lat.CrsAng(H,cell,SGData)
1274    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1275    SHnames = calcControls[phfx+'SHnames']
1276    for item in SHnames:
1277        L,N = eval(item.strip('C'))
1278        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
1279        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1280    return np.squeeze(odfCor)
1281   
1282def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1283    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1284    if 'T' in calcControls[hfx+'histType']:
1285        tth = parmDict[hfx+'2-theta']
1286    else:
1287        tth = refl[5+im]
1288    FORPI = 12.5663706143592
1289    odfCor = 1.0
1290    dFdODF = {}
1291    H = refl[:3]
1292    cell = G2lat.Gmat2cell(g)
1293    Sangl = [0.,0.,0.]
1294    if 'Bragg' in calcControls[hfx+'instType']:
1295        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1296        IFCoup = True
1297    else:
1298        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1299        IFCoup = False
1300    phi,beta = G2lat.CrsAng(H,cell,SGData)
1301    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1302    SHnames = calcControls[phfx+'SHnames']
1303    for item in SHnames:
1304        L,N = eval(item.strip('C'))
1305        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1306        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1307        dFdODF[phfx+item] = Kcsl*Lnorm
1308    return odfCor,dFdODF
1309   
1310def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1311    'March-Dollase preferred orientation correction'
1312    POcorr = 1.0
1313    MD = parmDict[phfx+'MD']
1314    if MD != 1.0:
1315        MDAxis = calcControls[phfx+'MDAxis']
1316        sumMD = 0
1317        for H in uniq:           
1318            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1319            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1320            sumMD += A**3
1321        POcorr = sumMD/len(uniq)
1322    return POcorr
1323   
1324def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1325    'Needs a doc string'
1326    POcorr = 1.0
1327    POderv = {}
1328    if calcControls[phfx+'poType'] == 'MD':
1329        MD = parmDict[phfx+'MD']
1330        MDAxis = calcControls[phfx+'MDAxis']
1331        sumMD = 0
1332        sumdMD = 0
1333        for H in uniq:           
1334            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1335            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1336            sumMD += A**3
1337            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1338        POcorr = sumMD/len(uniq)
1339        POderv[phfx+'MD'] = sumdMD/len(uniq)
1340    else:   #spherical harmonics
1341        if calcControls[phfx+'SHord']:
1342            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1343    return POcorr,POderv
1344   
1345def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1346    'Needs a doc string'
1347    if 'Debye' in calcControls[hfx+'instType']:
1348        if 'T' in calcControls[hfx+'histType']:
1349            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1350        else:
1351            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1352    else:
1353        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1354   
1355def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1356    'Needs a doc string'
1357    if 'Debye' in calcControls[hfx+'instType']:
1358        if 'T' in calcControls[hfx+'histType']:
1359            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1360        else:
1361            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1362    else:
1363        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1364       
1365def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1366    'Needs a doc string'
1367    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1368    pi2 = np.sqrt(2./np.pi)
1369    if 'T' in calcControls[hfx+'histType']:
1370        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1371        wave = refl[14+im]
1372    else:   #'C'W
1373        sth2 = sind(refl[5+im]/2.)**2
1374        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1375    c2th = 1.-2.0*sth2
1376    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1377    if 'X' in calcControls[hfx+'histType']:
1378        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1379    xfac = flv2*parmDict[phfx+'Extinction']
1380    exb = 1.0
1381    if xfac > -1.:
1382        exb = 1./(1.+xfac)
1383    exl = 1.0
1384    if 0 < xfac <= 1.:
1385        xn = np.array([xfac**(i+1) for i in range(6)])
1386        exl += np.sum(xn*coef)
1387    elif xfac > 1.:
1388        xfac2 = 1./np.sqrt(xfac)
1389        exl = pi2*(1.-0.125/xfac)*xfac2
1390    return exb*sth2+exl*(1.-sth2)
1391   
1392def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1393    'Needs a doc string'
1394    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1395    pi2 = np.sqrt(2./np.pi)
1396    if 'T' in calcControls[hfx+'histType']:
1397        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1398        wave = refl[14+im]
1399    else:   #'C'W
1400        sth2 = sind(refl[5+im]/2.)**2
1401        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1402    c2th = 1.-2.0*sth2
1403    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1404    if 'X' in calcControls[hfx+'histType']:
1405        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1406    xfac = flv2*parmDict[phfx+'Extinction']
1407    dbde = -500.*flv2
1408    if xfac > -1.:
1409        dbde = -flv2/(1.+xfac)**3
1410    dlde = 0.
1411    if 0 < xfac <= 1.:
1412        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1413        dlde = np.sum(xn*coef)
1414    elif xfac > 1.:
1415        xfac2 = 1./np.sqrt(xfac)
1416        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1417       
1418    return dbde*sth2+dlde*(1.-sth2)
1419
1420
1421#    delt = 0.01
1422#    parmDict[phfx+'Extinction'] += delt
1423#    plus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1424#    parmDict[phfx+'Extinction'] -= 2.*delt
1425#    minus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1426#    parmDict[phfx+'Extinction'] += delt
1427#    return (plus-minus)/(2.*delt)   
1428   
1429def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1430    'Needs a doc string'    #need powder extinction!
1431    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1432    if 'X' in parmDict[hfx+'Type']:
1433        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1434    POcorr = 1.0
1435    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1436        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1437    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1438        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1439    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1440        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1441    Icorr *= POcorr
1442    AbsCorr = 1.0
1443    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1444    Icorr *= AbsCorr
1445    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1446    Icorr *= ExtCorr
1447    return Icorr,POcorr,AbsCorr,ExtCorr
1448   
1449def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1450    'Needs a doc string'    #need powder extinction derivs!
1451    dIdsh = 1./parmDict[hfx+'Scale']
1452    dIdsp = 1./parmDict[phfx+'Scale']
1453    if 'X' in parmDict[hfx+'Type']:
1454        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1455        dIdPola /= pola
1456    else:       #'N'
1457        dIdPola = 0.0
1458    dFdODF = {}
1459    dFdSA = [0,0,0]
1460    dIdPO = {}
1461    if pfx+'SHorder' in parmDict:
1462        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1463        for iSH in dFdODF:
1464            dFdODF[iSH] /= odfCor
1465        for i in range(3):
1466            dFdSA[i] /= odfCor
1467    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1468        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1469        for iPO in dIdPO:
1470            dIdPO[iPO] /= POcorr
1471    if 'T' in parmDict[hfx+'Type']:
1472        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1473        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1474    else:
1475        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1476        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1477    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1478       
1479def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1480    'Needs a doc string'
1481    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1482        costh = cosd(refl[5+im]/2.)
1483        #crystallite size
1484        if calcControls[phfx+'SizeType'] == 'isotropic':
1485            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1486        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1487            H = np.array(refl[:3])
1488            P = np.array(calcControls[phfx+'SizeAxis'])
1489            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1490            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1491            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1492        else:           #ellipsoidal crystallites
1493            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1494            H = np.array(refl[:3])
1495            lenR = G2pwd.ellipseSize(H,Sij,GB)
1496            Sgam = 1.8*wave/(np.pi*costh*lenR)
1497        #microstrain               
1498        if calcControls[phfx+'MustrainType'] == 'isotropic':
1499            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1500        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1501            H = np.array(refl[:3])
1502            P = np.array(calcControls[phfx+'MustrainAxis'])
1503            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1504            Si = parmDict[phfx+'Mustrain;i']
1505            Sa = parmDict[phfx+'Mustrain;a']
1506            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1507        else:       #generalized - P.W. Stephens model
1508            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1509            Sum = 0
1510            for i,strm in enumerate(Strms):
1511                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1512            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1513    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1514        #crystallite size
1515        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1516            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1517        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1518            H = np.array(refl[:3])
1519            P = np.array(calcControls[phfx+'SizeAxis'])
1520            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1521            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1522            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1523        else:           #ellipsoidal crystallites   #OK
1524            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1525            H = np.array(refl[:3])
1526            lenR = G2pwd.ellipseSize(H,Sij,GB)
1527            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1528        #microstrain               
1529        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1530            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1531        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1532            H = np.array(refl[:3])
1533            P = np.array(calcControls[phfx+'MustrainAxis'])
1534            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1535            Si = parmDict[phfx+'Mustrain;i']
1536            Sa = parmDict[phfx+'Mustrain;a']
1537            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1538        else:       #generalized - P.W. Stephens model  OK
1539            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1540            Sum = 0
1541            for i,strm in enumerate(Strms):
1542                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1543            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1544           
1545    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1546    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1547    sig /= ateln2
1548    return sig,gam
1549       
1550def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1551    'Needs a doc string'
1552    gamDict = {}
1553    sigDict = {}
1554    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1555        costh = cosd(refl[5+im]/2.)
1556        tanth = tand(refl[5+im]/2.)
1557        #crystallite size derivatives
1558        if calcControls[phfx+'SizeType'] == 'isotropic':
1559            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1560            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1561            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1562        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1563            H = np.array(refl[:3])
1564            P = np.array(calcControls[phfx+'SizeAxis'])
1565            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1566            Si = parmDict[phfx+'Size;i']
1567            Sa = parmDict[phfx+'Size;a']
1568            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1569            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1570            Sgam = gami*sqtrm
1571            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1572            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1573            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1574            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1575            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1576            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1577        else:           #ellipsoidal crystallites
1578            const = 1.8*wave/(np.pi*costh)
1579            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1580            H = np.array(refl[:3])
1581            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1582            Sgam = const/lenR
1583            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1584                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1585                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1586        gamDict[phfx+'Size;mx'] = Sgam
1587        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1588               
1589        #microstrain derivatives               
1590        if calcControls[phfx+'MustrainType'] == 'isotropic':
1591            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1592            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1593            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1594        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1595            H = np.array(refl[:3])
1596            P = np.array(calcControls[phfx+'MustrainAxis'])
1597            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1598            Si = parmDict[phfx+'Mustrain;i']
1599            Sa = parmDict[phfx+'Mustrain;a']
1600            gami = 0.018*Si*Sa*tanth/np.pi
1601            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1602            Mgam = gami/sqtrm
1603            dsi = -gami*Si*cosP**2/sqtrm**3
1604            dsa = -gami*Sa*sinP**2/sqtrm**3
1605            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1606            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1607            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1608            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1609        else:       #generalized - P.W. Stephens model
1610            const = 0.018*refl[4+im]**2*tanth/np.pi
1611            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1612            Sum = 0
1613            for i,strm in enumerate(Strms):
1614                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1615                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1616                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1617            Mgam = const*np.sqrt(Sum)
1618            for i in range(len(Strms)):
1619                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1620                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1621        gamDict[phfx+'Mustrain;mx'] = Mgam
1622        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1623    else:   #'T'OF - All checked & OK
1624        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1625            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1626            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1627            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1628        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1629            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1630            H = np.array(refl[:3])
1631            P = np.array(calcControls[phfx+'SizeAxis'])
1632            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1633            Si = parmDict[phfx+'Size;i']
1634            Sa = parmDict[phfx+'Size;a']
1635            gami = const/(Si*Sa)
1636            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1637            Sgam = gami*sqtrm
1638            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1639            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1640            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1641            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1642            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1643            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1644        else:           #OK  ellipsoidal crystallites
1645            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1646            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1647            H = np.array(refl[:3])
1648            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1649            Sgam = const/lenR
1650            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1651                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1652                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1653        gamDict[phfx+'Size;mx'] = Sgam  #OK
1654        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1655               
1656        #microstrain derivatives               
1657        if calcControls[phfx+'MustrainType'] == 'isotropic':
1658            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1659            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1660            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1661        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1662            H = np.array(refl[:3])
1663            P = np.array(calcControls[phfx+'MustrainAxis'])
1664            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1665            Si = parmDict[phfx+'Mustrain;i']
1666            Sa = parmDict[phfx+'Mustrain;a']
1667            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1668            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1669            Mgam = gami/sqtrm
1670            dsi = -gami*Si*cosP**2/sqtrm**3
1671            dsa = -gami*Sa*sinP**2/sqtrm**3
1672            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1673            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1674            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1675            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1676        else:       #generalized - P.W. Stephens model OK
1677            pwrs = calcControls[phfx+'MuPwrs']
1678            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1679            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1680            Sum = 0
1681            for i,strm in enumerate(Strms):
1682                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1683                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1684                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1685            Mgam = const*np.sqrt(Sum)
1686            for i in range(len(Strms)):
1687                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1688                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1689        gamDict[phfx+'Mustrain;mx'] = Mgam
1690        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1691       
1692    return sigDict,gamDict
1693       
1694def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1695    'Needs a doc string'
1696    if im:
1697        h,k,l,m = refl[:4]
1698        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1699        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1700    else:
1701        h,k,l = refl[:3]
1702        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1703    refl[4+im] = d
1704    if 'C' in calcControls[hfx+'histType']:
1705        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1706        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1707        if 'Bragg' in calcControls[hfx+'instType']:
1708            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1709                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1710        else:               #Debye-Scherrer - simple but maybe not right
1711            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1712    elif 'T' in calcControls[hfx+'histType']:
1713        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1714        #do I need sample position effects - maybe?
1715    return pos
1716
1717def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1718    'Needs a doc string'
1719    dpr = 180./np.pi
1720    if im:
1721        h,k,l,m = refl[:4]
1722        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1723        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1724    else:
1725        m = 0
1726        h,k,l = refl[:3]       
1727        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1728    dst = np.sqrt(dstsq)
1729    dsp = 1./dst
1730    if 'C' in calcControls[hfx+'histType']:
1731        pos = refl[5+im]-parmDict[hfx+'Zero']
1732        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1733        dpdw = const*dst
1734        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1735        dpdZ = 1.0
1736        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1737            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1738        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1739        if 'Bragg' in calcControls[hfx+'instType']:
1740            dpdSh = -4.*shft*cosd(pos/2.0)
1741            dpdTr = -shft*sind(pos)*100.0
1742            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1743        else:               #Debye-Scherrer - simple but maybe not right
1744            dpdXd = -shft*cosd(pos)
1745            dpdYd = -shft*sind(pos)
1746            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1747    elif 'T' in calcControls[hfx+'histType']:
1748        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1749        dpdZ = 1.0
1750        dpdDC = dsp
1751        dpdDA = dsp**2
1752        dpdDB = 1./dsp
1753        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1754            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1755        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1756           
1757def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1758    'Needs a doc string'
1759    laue = SGData['SGLaue']
1760    uniq = SGData['SGUniq']
1761    h,k,l = refl[:3]
1762    if laue in ['m3','m3m']:
1763        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1764            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1765    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1766        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1767    elif laue in ['3R','3mR']:
1768        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1769    elif laue in ['4/m','4/mmm']:
1770        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1771    elif laue in ['mmm']:
1772        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1773    elif laue in ['2/m']:
1774        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1775        if uniq == 'a':
1776            Dij += parmDict[phfx+'D23']*k*l
1777        elif uniq == 'b':
1778            Dij += parmDict[phfx+'D13']*h*l
1779        elif uniq == 'c':
1780            Dij += parmDict[phfx+'D12']*h*k
1781    else:
1782        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1783            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1784    if 'C' in calcControls[hfx+'histType']:
1785        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1786    else:
1787        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1788           
1789def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1790    'Needs a doc string'
1791    laue = SGData['SGLaue']
1792    uniq = SGData['SGUniq']
1793    h,k,l = refl[:3]
1794    if laue in ['m3','m3m']:
1795        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1796            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1797    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1798        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1799    elif laue in ['3R','3mR']:
1800        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1801    elif laue in ['4/m','4/mmm']:
1802        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1803    elif laue in ['mmm']:
1804        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1805    elif laue in ['2/m']:
1806        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1807        if uniq == 'a':
1808            dDijDict[phfx+'D23'] = k*l
1809        elif uniq == 'b':
1810            dDijDict[phfx+'D13'] = h*l
1811        elif uniq == 'c':
1812            dDijDict[phfx+'D12'] = h*k
1813    else:
1814        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1815            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1816    if 'C' in calcControls[hfx+'histType']:
1817        for item in dDijDict:
1818            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1819    else:
1820        for item in dDijDict:
1821            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1822    return dDijDict
1823   
1824def GetDij(phfx,SGData,parmDict):
1825    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1826    return G2spc.HStrainVals(HSvals,SGData)
1827               
1828def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1829    'Needs a doc string'
1830    histoList = Histograms.keys()
1831    histoList.sort()
1832    for histogram in histoList:
1833        if 'PWDR' in histogram[:4]:
1834            Histogram = Histograms[histogram]
1835            hId = Histogram['hId']
1836            hfx = ':%d:'%(hId)
1837            Limits = calcControls[hfx+'Limits']
1838            if 'C' in calcControls[hfx+'histType']:
1839                shl = max(parmDict[hfx+'SH/L'],0.0005)
1840                Ka2 = False
1841                kRatio = 0.0
1842                if hfx+'Lam1' in parmDict.keys():
1843                    Ka2 = True
1844                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1845                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1846            x,y,w,yc,yb,yd = Histogram['Data']
1847            xB = np.searchsorted(x,Limits[0])
1848            xF = np.searchsorted(x,Limits[1])
1849            ymb = np.array(y-yb)
1850            ymb = np.where(ymb,ymb,1.0)
1851            ycmb = np.array(yc-yb)
1852            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1853            refLists = Histogram['Reflection Lists']
1854            for phase in refLists:
1855                if phase not in Phases:     #skips deleted or renamed phases silently!
1856                    continue
1857                Phase = Phases[phase]
1858                im = 0
1859                if Phase['General']['Type'] in ['modulated','magnetic']:
1860                    im = 1
1861                pId = Phase['pId']
1862                phfx = '%d:%d:'%(pId,hId)
1863                refDict = refLists[phase]
1864                sumFo = 0.0
1865                sumdF = 0.0
1866                sumFosq = 0.0
1867                sumdFsq = 0.0
1868                sumInt = 0.0
1869                for refl in refDict['RefList']:
1870                    if 'C' in calcControls[hfx+'histType']:
1871                        yp = np.zeros_like(yb)
1872                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1873                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1874                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1875                        iFin2 = iFin
1876                        if not iBeg+iFin:       #peak below low limit - skip peak
1877                            continue
1878                        elif not iBeg-iFin:     #peak above high limit - done
1879                            break
1880                        elif iBeg < iFin:
1881                            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
1882                            sumInt += refl[11+im]*refl[9+im]
1883                            if Ka2:
1884                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1885                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1886                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1887                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1888                                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
1889                                sumInt += refl[11+im]*refl[9+im]*kRatio
1890                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1891                               
1892                    elif 'T' in calcControls[hfx+'histType']:
1893                        yp = np.zeros_like(yb)
1894                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1895                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1896                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1897                        if iBeg < iFin:
1898                            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
1899                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1900                            sumInt += refl[11+im]*refl[9+im]
1901                    Fo = np.sqrt(np.abs(refl[8+im]))
1902                    Fc = np.sqrt(np.abs(refl[9]+im))
1903                    sumFo += Fo
1904                    sumFosq += refl[8+im]**2
1905                    sumdF += np.abs(Fo-Fc)
1906                    sumdFsq += (refl[8+im]-refl[9+im])**2
1907                if sumFo:
1908                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1909                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1910                else:
1911                    Histogram['Residuals'][phfx+'Rf'] = 100.
1912                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1913                Histogram['Residuals'][phfx+'sumInt'] = sumInt
1914                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1915                Histogram['Residuals']['hId'] = hId
1916        elif 'HKLF' in histogram[:4]:
1917            Histogram = Histograms[histogram]
1918            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1919               
1920def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1921    'Needs a doc string'
1922   
1923    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1924        U = parmDict[hfx+'U']
1925        V = parmDict[hfx+'V']
1926        W = parmDict[hfx+'W']
1927        X = parmDict[hfx+'X']
1928        Y = parmDict[hfx+'Y']
1929        tanPos = tand(refl[5+im]/2.0)
1930        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1931        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1932        sig = max(0.001,sig)
1933        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1934        gam = max(0.001,gam)
1935        return sig,gam
1936               
1937    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1938        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1939            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1940        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1941        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1942        sig += Ssig
1943        gam += Sgam
1944        return sig,gam
1945       
1946    def GetReflAlpBet(refl,im,hfx,parmDict):
1947        alp = parmDict[hfx+'alpha']/refl[4+im]
1948        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1949        return alp,bet
1950       
1951    hId = Histogram['hId']
1952    hfx = ':%d:'%(hId)
1953    bakType = calcControls[hfx+'bakType']
1954    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1955    yc = np.zeros_like(yb)
1956    cw = np.diff(x)
1957    cw = np.append(cw,cw[-1])
1958       
1959    if 'C' in calcControls[hfx+'histType']:   
1960        shl = max(parmDict[hfx+'SH/L'],0.002)
1961        Ka2 = False
1962        if hfx+'Lam1' in parmDict.keys():
1963            wave = parmDict[hfx+'Lam1']
1964            Ka2 = True
1965            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1966            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1967        else:
1968            wave = parmDict[hfx+'Lam']
1969    for phase in Histogram['Reflection Lists']:
1970        refDict = Histogram['Reflection Lists'][phase]
1971        if phase not in Phases:     #skips deleted or renamed phases silently!
1972            continue
1973        Phase = Phases[phase]
1974        pId = Phase['pId']
1975        pfx = '%d::'%(pId)
1976        phfx = '%d:%d:'%(pId,hId)
1977        hfx = ':%d:'%(hId)
1978        SGData = Phase['General']['SGData']
1979        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1980        im = 0
1981        if Phase['General']['Type'] in ['modulated','magnetic']:
1982            SSGData = Phase['General']['SSGData']
1983            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1984            im = 1  #offset in SS reflection list
1985            #??
1986        Dij = GetDij(phfx,SGData,parmDict)
1987        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
1988        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1989        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1990        Vst = np.sqrt(nl.det(G))    #V*
1991        if not Phase['General'].get('doPawley'):
1992            time0 = time.time()
1993            if im:
1994                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
1995            else:
1996                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1997#            print 'sf calc time: %.3fs'%(time.time()-time0)
1998        time0 = time.time()
1999        badPeak = False
2000        for iref,refl in enumerate(refDict['RefList']):
2001            if 'C' in calcControls[hfx+'histType']:
2002                if im:
2003                    h,k,l,m = refl[:4]
2004                else:
2005                    h,k,l = refl[:3]
2006                Uniq = np.inner(refl[:3],SGMT)
2007                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2008                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2009                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2010                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2011                refl[11+im] *= Vst*Lorenz
2012                 
2013                if Phase['General'].get('doPawley'):
2014                    try:
2015                        if im:
2016                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2017                        else:
2018                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2019                        refl[9+im] = parmDict[pInd]
2020                    except KeyError:
2021#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2022                        continue
2023                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2024                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2025                iFin = np.searchsorted(x,refl[5+im]+fmax)
2026                if not iBeg+iFin:       #peak below low limit - skip peak
2027                    continue
2028                elif not iBeg-iFin:     #peak above high limit - done
2029                    break
2030                elif iBeg > iFin:   #bad peak coeff - skip
2031                    badPeak = True
2032                    continue
2033                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
2034                if Ka2:
2035                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2036                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2037                    iBeg = np.searchsorted(x,pos2-fmin)
2038                    iFin = np.searchsorted(x,pos2+fmax)
2039                    if not iBeg+iFin:       #peak below low limit - skip peak
2040                        continue
2041                    elif not iBeg-iFin:     #peak above high limit - done
2042                        return yc,yb
2043                    elif iBeg > iFin:   #bad peak coeff - skip
2044                        continue
2045                    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
2046            elif 'T' in calcControls[hfx+'histType']:
2047                h,k,l = refl[:3]
2048                Uniq = np.inner(refl[:3],SGMT)
2049                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2050                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2051#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2052                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2053                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2054                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2055                refl[11+im] *= Vst*Lorenz
2056                if Phase['General'].get('doPawley'):
2057                    try:
2058                        if im:
2059                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2060                        else:
2061                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2062                        refl[9+im] = parmDict[pInd]
2063                    except KeyError:
2064#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2065                        continue
2066                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2067                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2068                iFin = np.searchsorted(x,refl[5+im]+fmax)
2069                if not iBeg+iFin:       #peak below low limit - skip peak
2070                    continue
2071                elif not iBeg-iFin:     #peak above high limit - done
2072                    break
2073                elif iBeg > iFin:   #bad peak coeff - skip
2074                    badPeak = True
2075                    continue
2076                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]
2077#        print 'profile calc time: %.3fs'%(time.time()-time0)
2078    if badPeak:
2079        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2080    return yc,yb
2081   
2082def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2083    'Needs a doc string'
2084   
2085    def cellVaryDerv(pfx,SGData,dpdA): 
2086        if SGData['SGLaue'] in ['-1',]:
2087            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2088                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2089        elif SGData['SGLaue'] in ['2/m',]:
2090            if SGData['SGUniq'] == 'a':
2091                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2092            elif SGData['SGUniq'] == 'b':
2093                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2094            else:
2095                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2096        elif SGData['SGLaue'] in ['mmm',]:
2097            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2098        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2099            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2100        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2101            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2102        elif SGData['SGLaue'] in ['3R', '3mR']:
2103            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2104        elif SGData['SGLaue'] in ['m3m','m3']:
2105            return [[pfx+'A0',dpdA[0]]]
2106           
2107    # create a list of dependent variables and set up a dictionary to hold their derivatives
2108    dependentVars = G2mv.GetDependentVars()
2109    depDerivDict = {}
2110    for j in dependentVars:
2111        depDerivDict[j] = np.zeros(shape=(len(x)))
2112    #print 'dependent vars',dependentVars
2113    lenX = len(x)               
2114    hId = Histogram['hId']
2115    hfx = ':%d:'%(hId)
2116    bakType = calcControls[hfx+'bakType']
2117    dMdv = np.zeros(shape=(len(varylist),len(x)))
2118    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2119    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2120        bBpos =varylist.index(hfx+'Back;0')
2121        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2122    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2123    for name in varylist:
2124        if 'Debye' in name:
2125            id = int(name.split(';')[-1])
2126            parm = name[:int(name.rindex(';'))]
2127            ip = names.index(parm)
2128            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2129    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2130    for name in varylist:
2131        if 'BkPk' in name:
2132            parm,id = name.split(';')
2133            id = int(id)
2134            if parm in names:
2135                ip = names.index(parm)
2136                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2137    cw = np.diff(x)
2138    cw = np.append(cw,cw[-1])
2139    Ka2 = False #also for TOF!
2140    if 'C' in calcControls[hfx+'histType']:   
2141        shl = max(parmDict[hfx+'SH/L'],0.002)
2142        if hfx+'Lam1' in parmDict.keys():
2143            wave = parmDict[hfx+'Lam1']
2144            Ka2 = True
2145            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2146            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2147        else:
2148            wave = parmDict[hfx+'Lam']
2149    for phase in Histogram['Reflection Lists']:
2150        refDict = Histogram['Reflection Lists'][phase]
2151        if phase not in Phases:     #skips deleted or renamed phases silently!
2152            continue
2153        Phase = Phases[phase]
2154        SGData = Phase['General']['SGData']
2155        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2156        im = 0
2157        if Phase['General']['Type'] in ['modulated','magnetic']:
2158            SSGData = Phase['General']['SSGData']
2159            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2160            im = 1  #offset in SS reflection list
2161            #??
2162        pId = Phase['pId']
2163        pfx = '%d::'%(pId)
2164        phfx = '%d:%d:'%(pId,hId)
2165        Dij = GetDij(phfx,SGData,parmDict)
2166        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2167        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2168        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2169        if not Phase['General'].get('doPawley'):
2170            time0 = time.time()
2171            if im:
2172                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2173            else:
2174                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2175#            print 'sf-derv time %.3fs'%(time.time()-time0)
2176            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2177        time0 = time.time()
2178        for iref,refl in enumerate(refDict['RefList']):
2179            if im:
2180                h,k,l,m = refl[:4]
2181            else:
2182                h,k,l = refl[:3]
2183            Uniq = np.inner(refl[:3],SGMT)
2184            if 'T' in calcControls[hfx+'histType']:
2185                wave = refl[14+im]
2186            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2187            if 'C' in calcControls[hfx+'histType']:        #CW powder
2188                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2189            else: #'T'OF
2190                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2191            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2192            iFin = np.searchsorted(x,refl[5+im]+fmax)
2193            if not iBeg+iFin:       #peak below low limit - skip peak
2194                continue
2195            elif not iBeg-iFin:     #peak above high limit - done
2196                break
2197            pos = refl[5+im]
2198            if 'C' in calcControls[hfx+'histType']:
2199                tanth = tand(pos/2.0)
2200                costh = cosd(pos/2.0)
2201                lenBF = iFin-iBeg
2202                dMdpk = np.zeros(shape=(6,lenBF))
2203                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2204                for i in range(5):
2205                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2206                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2207                if Ka2:
2208                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2209                    iBeg2 = np.searchsorted(x,pos2-fmin)
2210                    iFin2 = np.searchsorted(x,pos2+fmax)
2211                    if iBeg2-iFin2:
2212                        lenBF2 = iFin2-iBeg2
2213                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2214                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2215                        for i in range(5):
2216                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2217                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2218                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2219            else:   #'T'OF
2220                lenBF = iFin-iBeg
2221                if lenBF < 0:   #bad peak coeff
2222                    break
2223                dMdpk = np.zeros(shape=(6,lenBF))
2224                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2225                for i in range(6):
2226                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2227                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2228            if Phase['General'].get('doPawley'):
2229                dMdpw = np.zeros(len(x))
2230                try:
2231                    if im:
2232                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2233                    else:
2234                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2235                    idx = varylist.index(pIdx)
2236                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2237                    if Ka2: #not for TOF either
2238                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2239                    dMdv[idx] = dMdpw
2240                except: # ValueError:
2241                    pass
2242            if 'C' in calcControls[hfx+'histType']:
2243                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2244                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2245                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2246                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2247                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2248                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2249                    hfx+'DisplaceY':[dpdY,'pos'],}
2250                if 'Bragg' in calcControls[hfx+'instType']:
2251                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2252                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2253                else:
2254                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2255            else:   #'T'OF
2256                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2257                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2258                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2259                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2260                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2261                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2262                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2263                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2264            for name in names:
2265                item = names[name]
2266                if name in varylist:
2267                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2268                    if Ka2:
2269                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2270                elif name in dependentVars:
2271                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2272                    if Ka2:
2273                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2274            for iPO in dIdPO:
2275                if iPO in varylist:
2276                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2277                    if Ka2:
2278                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2279                elif iPO in dependentVars:
2280                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2281                    if Ka2:
2282                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2283            for i,name in enumerate(['omega','chi','phi']):
2284                aname = pfx+'SH '+name
2285                if aname in varylist:
2286                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2287                    if Ka2:
2288                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2289                elif aname in dependentVars:
2290                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2291                    if Ka2:
2292                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2293            for iSH in dFdODF:
2294                if iSH in varylist:
2295                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2296                    if Ka2:
2297                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2298                elif iSH in dependentVars:
2299                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2300                    if Ka2:
2301                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2302            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2303            for name,dpdA in cellDervNames:
2304                if name in varylist:
2305                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2306                    if Ka2:
2307                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2308                elif name in dependentVars:
2309                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2310                    if Ka2:
2311                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2312            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2313            for name in dDijDict:
2314                if name in varylist:
2315                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2316                    if Ka2:
2317                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2318                elif name in dependentVars:
2319                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2320                    if Ka2:
2321                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2322            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2323                if name in varylist:
2324                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2325                    if Ka2:
2326                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2327                elif name in dependentVars:
2328                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2329                    if Ka2:
2330                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2331            if 'C' in calcControls[hfx+'histType']:
2332                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2333            else:   #'T'OF
2334                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2335            for name in gamDict:
2336                if name in varylist:
2337                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2338                    if Ka2:
2339                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2340                elif name in dependentVars:
2341                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2342                    if Ka2:
2343                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2344            for name in sigDict:
2345                if name in varylist:
2346                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2347                    if Ka2:
2348                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2349                elif name in dependentVars:
2350                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2351                    if Ka2:
2352                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2353            for name in ['BabA','BabU']:
2354                if refl[9+im]:
2355                    if phfx+name in varylist:
2356                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2357                        if Ka2:
2358                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2359                    elif phfx+name in dependentVars:                   
2360                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2361                        if Ka2:
2362                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2363            if not Phase['General'].get('doPawley'):
2364                #do atom derivatives -  for RB,F,X & U so far             
2365                corr = dervDict['int']/refl[9+im]
2366                if Ka2:
2367                    corr2 = dervDict2['int']/refl[9+im]
2368                for name in varylist+dependentVars:
2369                    if '::RBV;' in name:
2370                        pass
2371                    else:
2372                        try:
2373                            aname = name.split(pfx)[1][:2]
2374                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2375                        except IndexError:
2376                            continue
2377                    if name in varylist:
2378                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2379                        if Ka2:
2380                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2381                    elif name in dependentVars:
2382                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2383                        if Ka2:
2384                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2385    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2386    # now process derivatives in constraints
2387    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2388    return dMdv
2389   
2390def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2391    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2392    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2393    are returned to either dervRefine or HessRefine.
2394
2395    :returns:
2396    '''
2397    nobs = Histogram['Residuals']['Nobs']
2398    hId = Histogram['hId']
2399    hfx = ':%d:'%(hId)
2400    pfx = '%d::'%(Phase['pId'])
2401    phfx = '%d:%d:'%(Phase['pId'],hId)
2402    SGData = Phase['General']['SGData']
2403    im = 0
2404    if Phase['General']['Type'] in ['modulated','magnetic']:
2405        SSGData = Phase['General']['SSGData']
2406        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2407        im = 1  #offset in SS reflection list
2408        #??
2409    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2410    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2411    refDict = Histogram['Data']
2412    if im:
2413        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2414    else:
2415        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2416    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2417    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2418    dependentVars = G2mv.GetDependentVars()
2419    depDerivDict = {}
2420    for j in dependentVars:
2421        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2422    wdf = np.zeros(len(refDict['RefList']))
2423    if calcControls['F**2']:
2424        for iref,ref in enumerate(refDict['RefList']):
2425            if ref[6+im] > 0:
2426                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1] 
2427                w = 1.0/ref[6+im]
2428                if w*ref[5+im] >= calcControls['minF/sig'] and ref[3+im] > 0:
2429                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2430                    for j,var in enumerate(varylist):
2431                        if var in dFdvDict:
2432                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2433                    for var in dependentVars:
2434                        if var in dFdvDict:
2435                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2436                    if phfx+'Scale' in varylist:
2437                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2438                    elif phfx+'Scale' in dependentVars:
2439                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]
2440                    for item in ['Ep','Es','Eg']:
2441                        if phfx+item in varylist and phfx+item in dervDict:
2442                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2443                        elif phfx+item in dependentVars and phfx+item in dervDict:
2444                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2445                    for item in ['BabA','BabU']:
2446                        if phfx+item in varylist:
2447                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2448                        elif phfx+item in dependentVars:
2449                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2450    else:
2451        for iref,ref in enumerate(refDict['RefList']):
2452            if ref[5+im] > 0.:
2453                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2454                Fo = np.sqrt(ref[5+im])
2455                Fc = np.sqrt(ref[7+im])
2456                w = 1.0/ref[6+im]
2457                if ref[5+im]/ref[6+im] >= calcControls['minF/sig'] and ref[3+im] > 0:
2458                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2459                    for j,var in enumerate(varylist):
2460                        if var in dFdvDict:
2461                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2462                    for var in dependentVars:
2463                        if var in dFdvDict:
2464                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2465                    if phfx+'Scale' in varylist:
2466                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2467                    elif phfx+'Scale' in dependentVars:
2468                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]                           
2469                    for item in ['Ep','Es','Eg']:
2470                        if phfx+item in varylist and phfx+item in dervDict:
2471                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #correct
2472                        elif phfx+item in dependentVars and phfx+item in dervDict:
2473                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2474                    for item in ['BabA','BabU']:
2475                        if phfx+item in varylist:
2476                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2477                        elif phfx+item in dependentVars:
2478                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2479    return dMdvh,depDerivDict,wdf
2480   
2481
2482def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2483    '''Loop over histograms and compute derivatives of the fitting
2484    model (M) with respect to all parameters.  Results are returned in
2485    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2486    n is the number of parameters and m is the number of data
2487    points. This can exceed memory when m gets large. This routine is
2488    used when refinement derivatives are selected as "analtytic
2489    Jacobian" in Controls.
2490
2491    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2492    '''
2493    parmDict.update(zip(varylist,values))
2494    G2mv.Dict2Map(parmDict,varylist)
2495    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2496    nvar = len(varylist)
2497    dMdv = np.empty(0)
2498    histoList = Histograms.keys()
2499    histoList.sort()
2500    for histogram in histoList:
2501        if 'PWDR' in histogram[:4]:
2502            Histogram = Histograms[histogram]
2503            hId = Histogram['hId']
2504            hfx = ':%d:'%(hId)
2505            wtFactor = calcControls[hfx+'wtFactor']
2506            Limits = calcControls[hfx+'Limits']
2507            x,y,w,yc,yb,yd = Histogram['Data']
2508            xB = np.searchsorted(x,Limits[0])
2509            xF = np.searchsorted(x,Limits[1])
2510            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2511                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2512        elif 'HKLF' in histogram[:4]:
2513            Histogram = Histograms[histogram]
2514            phase = Histogram['Reflection Lists']
2515            Phase = Phases[phase]
2516            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2517            hfx = ':%d:'%(Histogram['hId'])
2518            wtFactor = calcControls[hfx+'wtFactor']
2519            # now process derivatives in constraints
2520            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2521        else:
2522            continue        #skip non-histogram entries
2523        if len(dMdv):
2524            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2525        else:
2526            dMdv = np.sqrt(wtFactor)*dMdvh
2527           
2528    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2529    if np.any(pVals):
2530        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2531        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2532       
2533    return dMdv
2534
2535def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2536    '''Loop over histograms and compute derivatives of the fitting
2537    model (M) with respect to all parameters.  For each histogram, the
2538    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2539    number of parameters and m is the number of data points *in the
2540    histogram*. The (n by n) Hessian is computed from each Jacobian
2541    and it is returned.  This routine is used when refinement
2542    derivatives are selected as "analtytic Hessian" in Controls.
2543
2544    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2545    '''
2546    parmDict.update(zip(varylist,values))
2547    G2mv.Dict2Map(parmDict,varylist)
2548    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2549    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2550    nvar = len(varylist)
2551    Hess = np.empty(0)
2552    histoList = Histograms.keys()
2553    histoList.sort()
2554    for histogram in histoList:
2555        if 'PWDR' in histogram[:4]:
2556            Histogram = Histograms[histogram]
2557            hId = Histogram['hId']
2558            hfx = ':%d:'%(hId)
2559            wtFactor = calcControls[hfx+'wtFactor']
2560            Limits = calcControls[hfx+'Limits']
2561            x,y,w,yc,yb,yd = Histogram['Data']
2562            W = wtFactor*w
2563            dy = y-yc
2564            xB = np.searchsorted(x,Limits[0])
2565            xF = np.searchsorted(x,Limits[1])
2566            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2567                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2568            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2569            Dy = dy[xB:xF][np.newaxis,:]
2570            dMdvh *= Wt
2571            if dlg:
2572                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2573            if len(Hess):
2574                Hess += np.inner(dMdvh,dMdvh)
2575                dMdvh *= Wt*Dy
2576                Vec += np.sum(dMdvh,axis=1)
2577            else:
2578                Hess = np.inner(dMdvh,dMdvh)
2579                dMdvh *= Wt*Dy
2580                Vec = np.sum(dMdvh,axis=1)
2581        elif 'HKLF' in histogram[:4]:
2582            Histogram = Histograms[histogram]
2583            phase = Histogram['Reflection Lists']
2584            Phase = Phases[phase]
2585            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2586            hId = Histogram['hId']
2587            hfx = ':%d:'%(Histogram['hId'])
2588            wtFactor = calcControls[hfx+'wtFactor']
2589            # now process derivatives in constraints
2590            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2591#            print 'matrix build time: %.3f'%(time.time()-time0)
2592
2593            if dlg:
2594                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2595            if len(Hess):
2596                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2597                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2598            else:
2599                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2600                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2601        else:
2602            continue        #skip non-histogram entries
2603    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2604    if np.any(pVals):
2605        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2606        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2607        Hess += np.inner(dpdv*pWt,dpdv)
2608    return Vec,Hess
2609
2610def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2611    '''Computes the point-by-point discrepancies between every data point in every histogram
2612    and the observed value
2613   
2614    :returns: an np array of differences between observed and computed diffraction values.
2615    '''
2616    Values2Dict(parmDict, varylist, values)
2617    G2mv.Dict2Map(parmDict,varylist)
2618    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2619    M = np.empty(0)
2620    SumwYo = 0
2621    Nobs = 0
2622    Nrej = 0
2623    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2624    histoList = Histograms.keys()
2625    histoList.sort()
2626    for histogram in histoList:
2627        if 'PWDR' in histogram[:4]:
2628            Histogram = Histograms[histogram]
2629            hId = Histogram['hId']
2630            hfx = ':%d:'%(hId)
2631            wtFactor = calcControls[hfx+'wtFactor']
2632            Limits = calcControls[hfx+'Limits']
2633            x,y,w,yc,yb,yd = Histogram['Data']
2634            yc *= 0.0                           #zero full calcd profiles
2635            yb *= 0.0
2636            yd *= 0.0
2637            xB = np.searchsorted(x,Limits[0])
2638            xF = np.searchsorted(x,Limits[1])
2639            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2640                varylist,Histogram,Phases,calcControls,pawleyLookup)
2641            yc[xB:xF] += yb[xB:xF]
2642            if not np.any(y):                   #fill dummy data
2643                rv = st.poisson(yc[xB:xF])
2644                y[xB:xF] = rv.rvs()
2645                Z = np.ones_like(yc[xB:xF])
2646                Z[1::2] *= -1
2647                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2648                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2649            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2650            W = wtFactor*w
2651            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2652            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2653            Nobs += Histogram['Residuals']['Nobs']
2654            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2655            SumwYo += Histogram['Residuals']['sumwYo']
2656            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2657            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2658            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2659            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2660            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.))
2661            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.))
2662            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2663            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2664            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2665            if dlg:
2666                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2667            M = np.concatenate((M,wdy))
2668#end of PWDR processing
2669        elif 'HKLF' in histogram[:4]:
2670            Histogram = Histograms[histogram]
2671            Histogram['Residuals'] = {}
2672            phase = Histogram['Reflection Lists']
2673            Phase = Phases[phase]
2674            hId = Histogram['hId']
2675            hfx = ':%d:'%(hId)
2676            wtFactor = calcControls[hfx+'wtFactor']
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                #??
2686            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2687            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2688            refDict = Histogram['Data']
2689            time0 = time.time()
2690            if im:
2691                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2692            else:
2693                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2694#            print 'sf-calc time: %.3f'%(time.time()-time0)
2695            df = np.zeros(len(refDict['RefList']))
2696            sumwYo = 0
2697            sumFo = 0
2698            sumFo2 = 0
2699            sumdF = 0
2700            sumdF2 = 0
2701            nobs = 0
2702            nrej = 0
2703            if calcControls['F**2']:
2704                for i,ref in enumerate(refDict['RefList']):
2705                    if ref[6+im] > 0:
2706                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2707                        w = 1.0/ref[6+im]   # 1/sig(F^2)
2708                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
2709                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2710                        if w*ref[5+im] >= calcControls['minF/sig'] and ref[3+im] > 0:  #min cutoff & user rejection
2711                            Fo = np.sqrt(ref[5+im])
2712                            sumFo += Fo
2713                            sumFo2 += ref[5+im]
2714                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2715                            sumdF2 += abs(ref[5+im]-ref[7+im])
2716                            nobs += 1
2717                            df[i] = -w*(ref[5+im]-ref[7+im])
2718                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
2719                        else:
2720                            nrej += 1
2721            else:
2722                for i,ref in enumerate(refDict['RefList']):
2723                    if ref[5+im] > 0.:
2724                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2725                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
2726                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2727                        Fo = np.sqrt(ref[5+im])
2728                        Fc = np.sqrt(ref[7+im])
2729                        w = 2.0*(Fo/ref[6+im])**2    # 1/sig(F)?
2730                        if ref[5+im]/ref[6+im] >= calcControls['minF/sig'] and ref[3+im] > 0:  #min cutoff & user rejection
2731                            sumFo += Fo
2732                            sumFo2 += ref[5+im]
2733                            sumdF += abs(Fo-Fc)
2734                            sumdF2 += abs(ref[5+im]-ref[7+im])
2735                            nobs += 1
2736                            df[i] = -w*(Fo-Fc)
2737                            sumwYo += (w*Fo)**2
2738                        else:
2739                            nrej += 1
2740            Histogram['Residuals']['Nobs'] = nobs
2741            Histogram['Residuals']['sumwYo'] = sumwYo
2742            SumwYo += sumwYo
2743            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/sumwYo)*100.)
2744            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2745            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2746            Histogram['Residuals'][phfx+'Nref'] = nobs
2747            Histogram['Residuals'][phfx+'Nrej'] = nrej
2748            Nobs += nobs
2749            Nrej += nrej
2750            if dlg:
2751                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2752            M = np.concatenate((M,wtFactor*df))
2753# end of HKLF processing
2754    Histograms['sumwYo'] = SumwYo
2755    Histograms['Nobs'] = Nobs
2756    Histograms['Nrej'] = Nrej
2757    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2758    if dlg:
2759        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2760        if not GoOn:
2761            parmDict['saved values'] = values
2762            dlg.Destroy()
2763            raise Exception         #Abort!!
2764    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2765    if len(pVals):
2766        pSum = np.sum(pWt*pVals**2)
2767        for name in pWsum:
2768            if pWsum[name]:
2769                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2770        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2771        Nobs += len(pVals)
2772        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2773    return M
2774                       
Note: See TracBrowser for help on using the repository browser.