source: trunk/GSASIIstrMath.py @ 1630

Last change on this file since 1630 was 1630, checked in by vondreele, 8 years ago

add export routine for tables of x,y0,y1,y2,... for multiple powder data sets - apparently useful for plotting in Excel, etc.
implement SS structure factor calcs for position modulation - Fourier only
working on printing wave results

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