source: trunk/GSASIIstrMath.py @ 1676

Last change on this file since 1676 was 1676, checked in by toby, 7 years ago

Ignore equivalences that are not in use; start on svn switch implementation for help

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 133.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-02-28 03:46:37 +0000 (Sat, 28 Feb 2015) $
8# $Author: toby $
9# $Revision: 1676 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1676 2015-02-28 03:46:37Z toby $
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: 1676 $")
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 = np.zeros_like(fa)
714        if not SGData['SGInv']:
715            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
716        fa = fa*GfpuA-fb*GfpuB
717        fb = fb*GfpuA+fa*GfpuB
718        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
719        fbs = np.real(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        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
946        SQfactor = 8.0*SQ*np.pi**2
947        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
948        Bab = parmDict[phfx+'BabA']*dBabdA
949        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
950        FF = refDict['FF']['FF'][iref].T[Tindx]
951        Uniq = np.inner(H[:3],SGMT)
952        SSUniq = np.inner(H,SSGMT)
953        Phi = np.inner(H[:3],SGT)
954        SSPhi = np.inner(H,SSGT)
955        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
956        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
957        #need ModulationDerv here dGAdXsin, etc 
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[:3],np.inner(bij,H[:3]))
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        GfpuA = np.swapaxes(GfpuA,1,2)
974        GfpuB = np.swapaxes(GfpuB,1,2)
975        fa = fa*GfpuA-fb*GfpuB
976        fb = fb*GfpuA+fa*GfpuB
977       
978        fas = np.sum(np.sum(fa,axis=1),axis=1)
979        fbs = np.sum(np.sum(fb,axis=1),axis=1)
980        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
981        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
982        fax = fax*GfpuA-fbx*GfpuB
983        fbx = fbx*GfpuA+fax*GfpuB
984        #sum below is over Uniq
985        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
986        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
987        dfadui = np.sum(-SQfactor*fa,axis=2)
988        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
989        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
990        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
991        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
992        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
993        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
994        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
995        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
996        #need dFdXsin, etc for modulations...
997        if not SGData['SGInv']:
998            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
999            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
1000            dfbdui = np.sum(-SQfactor*fb,axis=2)
1001            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1002            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
1003            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1004            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1005            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1006            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1007            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1008        #need dFdXsin, etc for modulations...
1009        #loop over atoms - each dict entry is list of derivatives for all the reflections
1010    for i in range(len(Mdata)):     
1011        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1012        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1013        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1014        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1015        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1016        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1017        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1018        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1019        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1020        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1021        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1022        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1023    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
1024    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
1025    return dFdvDict
1026   
1027def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1028    ''' Single crystal extinction function; returns extinction & derivative
1029    '''
1030    extCor = 1.0
1031    dervDict = {}
1032    if calcControls[phfx+'EType'] != 'None':
1033        SQ = 1/(4.*ref[4+im]**2)
1034        if 'C' in parmDict[hfx+'Type']:           
1035            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1036        else:   #'T'
1037            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1038        if 'SXC' in parmDict[hfx+'Type']:
1039            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1040            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1041            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1042            PLZ = AV*P12*ref[7+im]*parmDict[hfx+'Lam']**2
1043        elif 'SNT' in parmDict[hfx+'Type']:
1044            AV = 1.e7/parmDict[pfx+'Vol']**2
1045            PL = SQ
1046            PLZ = AV*ref[7+im]*ref[12+im]**2
1047        elif 'SNC' in parmDict[hfx+'Type']:
1048            AV = 1.e7/parmDict[pfx+'Vol']**2
1049            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1050            PLZ = AV*ref[9]*parmDict[hfx+'Lam']**2      #Fcsq as per GSAS, why not FcTsq (ref[9])?
1051           
1052        if 'Primary' in calcControls[phfx+'EType']:
1053            PLZ *= 1.5
1054        else:
1055            if 'C' in parmDict[hfx+'Type']:
1056                PLZ *= calcControls[phfx+'Tbar']
1057            else: #'T'
1058                PLZ *= ref[13+im]      #t-bar
1059        if 'Primary' in calcControls[phfx+'EType']:
1060            PLZ *= 1.5
1061            PSIG = parmDict[phfx+'Ep']
1062        elif 'I & II' in calcControls[phfx+'EType']:
1063            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1064        elif 'Type II' in calcControls[phfx+'EType']:
1065            PSIG = parmDict[phfx+'Es']
1066        else:       # 'Secondary Type I'
1067            PSIG = parmDict[phfx+'Eg']/PL
1068           
1069        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1070        AL = 0.025+0.285*cos2T
1071        BG = 0.02-0.025*cos2T
1072        BL = 0.15-0.2*(0.75-cos2T)**2
1073        if cos2T < 0.:
1074            BL = -0.45*cos2T
1075        CG = 2.
1076        CL = 2.
1077        PF = PLZ*PSIG
1078       
1079        if 'Gaussian' in calcControls[phfx+'EApprox']:
1080            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1081            extCor = np.sqrt(PF4)
1082            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1083        else:
1084            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1085            extCor = np.sqrt(PF4)
1086            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1087
1088        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1089            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1090        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1091            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1092        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1093            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1094               
1095    return 1./extCor,dervDict
1096   
1097def Dict2Values(parmdict, varylist):
1098    '''Use before call to leastsq to setup list of values for the parameters
1099    in parmdict, as selected by key in varylist'''
1100    return [parmdict[key] for key in varylist] 
1101   
1102def Values2Dict(parmdict, varylist, values):
1103    ''' Use after call to leastsq to update the parameter dictionary with
1104    values corresponding to keys in varylist'''
1105    parmdict.update(zip(varylist,values))
1106   
1107def GetNewCellParms(parmDict,varyList):
1108    'Needs a doc string'
1109    newCellDict = {}
1110    Anames = ['A'+str(i) for i in range(6)]
1111    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1112    for item in varyList:
1113        keys = item.split(':')
1114        if keys[2] in Ddict:
1115            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1116            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1117            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1118    return newCellDict          # is e.g. {'0::D11':A0-D11}
1119   
1120def ApplyXYZshifts(parmDict,varyList):
1121    '''
1122    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1123   
1124    :param dict parmDict: parameter dictionary
1125    :param list varyList: list of variables (not used!)
1126    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1127
1128    '''
1129    newAtomDict = {}
1130    for item in parmDict:
1131        if 'dA' in item:
1132            parm = ''.join(item.split('d'))
1133            parmDict[parm] += parmDict[item]
1134            newAtomDict[item] = [parm,parmDict[parm]]
1135    return newAtomDict
1136   
1137def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1138    'Spherical harmonics texture'
1139    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1140    if 'T' in calcControls[hfx+'histType']:
1141        tth = parmDict[hfx+'2-theta']
1142    else:
1143        tth = refl[5+im]
1144    odfCor = 1.0
1145    H = refl[:3]
1146    cell = G2lat.Gmat2cell(g)
1147    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1148    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1149    phi,beta = G2lat.CrsAng(H,cell,SGData)
1150    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1151    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1152    for item in SHnames:
1153        L,M,N = eval(item.strip('C'))
1154        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1155        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1156        Lnorm = G2lat.Lnorm(L)
1157        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1158    return odfCor
1159   
1160def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1161    'Spherical harmonics texture derivatives'
1162    if 'T' in calcControls[hfx+'histType']:
1163        tth = parmDict[hfx+'2-theta']
1164    else:
1165        tth = refl[5+im]
1166    FORPI = 4.0*np.pi
1167    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1168    odfCor = 1.0
1169    dFdODF = {}
1170    dFdSA = [0,0,0]
1171    H = refl[:3]
1172    cell = G2lat.Gmat2cell(g)
1173    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1174    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1175    phi,beta = G2lat.CrsAng(H,cell,SGData)
1176    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1177    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1178    for item in SHnames:
1179        L,M,N = eval(item.strip('C'))
1180        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1181        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1182        Lnorm = G2lat.Lnorm(L)
1183        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1184        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1185        for i in range(3):
1186            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1187    return odfCor,dFdODF,dFdSA
1188   
1189def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1190    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1191    if 'T' in calcControls[hfx+'histType']:
1192        tth = parmDict[hfx+'2-theta']
1193    else:
1194        tth = refl[5+im]
1195    odfCor = 1.0
1196    H = refl[:3]
1197    cell = G2lat.Gmat2cell(g)
1198    Sangl = [0.,0.,0.]
1199    if 'Bragg' in calcControls[hfx+'instType']:
1200        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1201        IFCoup = True
1202    else:
1203        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1204        IFCoup = False
1205    phi,beta = G2lat.CrsAng(H,cell,SGData)
1206    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1207    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1208    for item in SHnames:
1209        L,N = eval(item.strip('C'))
1210        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
1211        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1212    return np.squeeze(odfCor)
1213   
1214def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1215    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1216    if 'T' in calcControls[hfx+'histType']:
1217        tth = parmDict[hfx+'2-theta']
1218    else:
1219        tth = refl[5+im]
1220    FORPI = 12.5663706143592
1221    odfCor = 1.0
1222    dFdODF = {}
1223    H = refl[:3]
1224    cell = G2lat.Gmat2cell(g)
1225    Sangl = [0.,0.,0.]
1226    if 'Bragg' in calcControls[hfx+'instType']:
1227        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1228        IFCoup = True
1229    else:
1230        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1231        IFCoup = False
1232    phi,beta = G2lat.CrsAng(H,cell,SGData)
1233    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1234    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1235    for item in SHnames:
1236        L,N = eval(item.strip('C'))
1237        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1238        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1239        dFdODF[phfx+item] = Kcsl*Lnorm
1240    return odfCor,dFdODF
1241   
1242def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1243    'March-Dollase preferred orientation correction'
1244    POcorr = 1.0
1245    MD = parmDict[phfx+'MD']
1246    if MD != 1.0:
1247        MDAxis = calcControls[phfx+'MDAxis']
1248        sumMD = 0
1249        for H in uniq:           
1250            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1251            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1252            sumMD += A**3
1253        POcorr = sumMD/len(uniq)
1254    return POcorr
1255   
1256def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1257    'Needs a doc string'
1258    POcorr = 1.0
1259    POderv = {}
1260    if calcControls[phfx+'poType'] == 'MD':
1261        MD = parmDict[phfx+'MD']
1262        MDAxis = calcControls[phfx+'MDAxis']
1263        sumMD = 0
1264        sumdMD = 0
1265        for H in uniq:           
1266            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1267            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1268            sumMD += A**3
1269            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1270        POcorr = sumMD/len(uniq)
1271        POderv[phfx+'MD'] = sumdMD/len(uniq)
1272    else:   #spherical harmonics
1273        if calcControls[phfx+'SHord']:
1274            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1275    return POcorr,POderv
1276   
1277def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1278    'Needs a doc string'
1279    if 'Debye' in calcControls[hfx+'instType']:
1280        if 'T' in calcControls[hfx+'histType']:
1281            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],parmDict[hfx+'2-theta'],0,0)
1282        else:
1283            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1284    else:
1285        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1286   
1287def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1288    'Needs a doc string'
1289    if 'Debye' in calcControls[hfx+'instType']:
1290        if 'T' in calcControls[hfx+'histType']:
1291            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],parmDict[hfx+'2-theta'],0,0)
1292        else:
1293            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1294    else:
1295        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1296       
1297def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1298    'Needs a doc string'
1299    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1300    pi2 = np.sqrt(2./np.pi)
1301    if 'T' in calcControls[hfx+'histType']:
1302        sth2 = sind(parmDict[hfx+'2-theta']/2.)**2
1303        wave = refl[14+im]
1304    else:   #'C'W
1305        sth2 = sind(refl[5+im]/2.)**2
1306        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1307    c2th = 1.-2.0*sth2
1308    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1309    if 'X' in calcControls[hfx+'histType']:
1310        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1311    xfac = flv2*parmDict[phfx+'Extinction']
1312    exb = 1.0
1313    if xfac > -1.:
1314        exb = 1./(1.+xfac)
1315    exl = 1.0
1316    if 0 < xfac <= 1.:
1317        xn = np.array([xfac**(i+1) for i in range(6)])
1318        exl = np.sum(xn*coef)
1319    elif xfac > 1.:
1320        xfac2 = 1./np.sqrt(xfac)
1321        exl = pi2*(1.-0.125/xfac)*xfac2
1322    return exb*sth2+exl*(1.-sth2)
1323   
1324def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1325    'Needs a doc string'
1326    delt = 0.001
1327    parmDict[phfx+'Extinction'] += delt
1328    plus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1329    parmDict[phfx+'Extinction'] -= 2.*delt
1330    minus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1331    parmDict[phfx+'Extinction'] += delt
1332    return (plus-minus)/(2.*delt)   
1333   
1334def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1335    'Needs a doc string'    #need powder extinction!
1336    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1337    if 'X' in parmDict[hfx+'Type']:
1338        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1339    POcorr = 1.0
1340    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture
1341        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1342    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1343        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1344    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1345        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1346    Icorr *= POcorr
1347    AbsCorr = 1.0
1348    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1349    Icorr *= AbsCorr
1350    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1351    Icorr *= ExtCorr
1352    return Icorr,POcorr,AbsCorr,ExtCorr
1353   
1354def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1355    'Needs a doc string'    #need powder extinction derivs!
1356    dIdsh = 1./parmDict[hfx+'Scale']
1357    dIdsp = 1./parmDict[phfx+'Scale']
1358    if 'X' in parmDict[hfx+'Type']:
1359        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1360        dIdPola /= pola
1361    else:       #'N'
1362        dIdPola = 0.0
1363    dFdODF = {}
1364    dFdSA = [0,0,0]
1365    dIdPO = {}
1366    if pfx+'SHorder' in parmDict:
1367        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1368        for iSH in dFdODF:
1369            dFdODF[iSH] /= odfCor
1370        for i in range(3):
1371            dFdSA[i] /= odfCor
1372    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1373        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1374        for iPO in dIdPO:
1375            dIdPO[iPO] /= POcorr
1376    if 'T' in parmDict[hfx+'Type']:
1377        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1378        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1379    else:
1380        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1381        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1382    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1383       
1384def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1385    'Needs a doc string'
1386    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1387        costh = cosd(refl[5+im]/2.)
1388        #crystallite size
1389        if calcControls[phfx+'SizeType'] == 'isotropic':
1390            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1391        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1392            H = np.array(refl[:3])
1393            P = np.array(calcControls[phfx+'SizeAxis'])
1394            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1395            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1396            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1397        else:           #ellipsoidal crystallites
1398            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1399            H = np.array(refl[:3])
1400            lenR = G2pwd.ellipseSize(H,Sij,GB)
1401            Sgam = 1.8*wave/(np.pi*costh*lenR)
1402        #microstrain               
1403        if calcControls[phfx+'MustrainType'] == 'isotropic':
1404            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1405        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1406            H = np.array(refl[:3])
1407            P = np.array(calcControls[phfx+'MustrainAxis'])
1408            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1409            Si = parmDict[phfx+'Mustrain;i']
1410            Sa = parmDict[phfx+'Mustrain;a']
1411            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1412        else:       #generalized - P.W. Stephens model
1413            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1414            Sum = 0
1415            for i,strm in enumerate(Strms):
1416                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1417            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1418    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1419        #crystallite size
1420        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1421            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1422        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1423            H = np.array(refl[:3])
1424            P = np.array(calcControls[phfx+'SizeAxis'])
1425            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1426            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1427            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1428        else:           #ellipsoidal crystallites   #OK
1429            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1430            H = np.array(refl[:3])
1431            lenR = G2pwd.ellipseSize(H,Sij,GB)
1432            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1433        #microstrain               
1434        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1435            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1436        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1437            H = np.array(refl[:3])
1438            P = np.array(calcControls[phfx+'MustrainAxis'])
1439            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1440            Si = parmDict[phfx+'Mustrain;i']
1441            Sa = parmDict[phfx+'Mustrain;a']
1442            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1443        else:       #generalized - P.W. Stephens model  OK
1444            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1445            Sum = 0
1446            for i,strm in enumerate(Strms):
1447                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1448            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1449           
1450    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1451    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1452    sig /= ateln2
1453    return sig,gam
1454       
1455def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1456    'Needs a doc string'
1457    gamDict = {}
1458    sigDict = {}
1459    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1460        costh = cosd(refl[5+im]/2.)
1461        tanth = tand(refl[5+im]/2.)
1462        #crystallite size derivatives
1463        if calcControls[phfx+'SizeType'] == 'isotropic':
1464            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1465            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1466            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1467        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1468            H = np.array(refl[:3])
1469            P = np.array(calcControls[phfx+'SizeAxis'])
1470            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1471            Si = parmDict[phfx+'Size;i']
1472            Sa = parmDict[phfx+'Size;a']
1473            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1474            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1475            Sgam = gami*sqtrm
1476            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1477            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1478            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1479            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1480            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1481            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1482        else:           #ellipsoidal crystallites
1483            const = 1.8*wave/(np.pi*costh)
1484            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1485            H = np.array(refl[:3])
1486            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1487            Sgam = const/lenR
1488            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1489                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1490                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1491        gamDict[phfx+'Size;mx'] = Sgam
1492        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1493               
1494        #microstrain derivatives               
1495        if calcControls[phfx+'MustrainType'] == 'isotropic':
1496            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1497            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1498            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1499        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1500            H = np.array(refl[:3])
1501            P = np.array(calcControls[phfx+'MustrainAxis'])
1502            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1503            Si = parmDict[phfx+'Mustrain;i']
1504            Sa = parmDict[phfx+'Mustrain;a']
1505            gami = 0.018*Si*Sa*tanth/np.pi
1506            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1507            Mgam = gami/sqtrm
1508            dsi = -gami*Si*cosP**2/sqtrm**3
1509            dsa = -gami*Sa*sinP**2/sqtrm**3
1510            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1511            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1512            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1513            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1514        else:       #generalized - P.W. Stephens model
1515            const = 0.018*refl[4+im]**2*tanth/np.pi
1516            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1517            Sum = 0
1518            for i,strm in enumerate(Strms):
1519                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1520                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1521                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1522            Mgam = const*np.sqrt(Sum)
1523            for i in range(len(Strms)):
1524                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1525                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1526        gamDict[phfx+'Mustrain;mx'] = Mgam
1527        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1528    else:   #'T'OF - All checked & OK
1529        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1530            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1531            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1532            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1533        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1534            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1535            H = np.array(refl[:3])
1536            P = np.array(calcControls[phfx+'SizeAxis'])
1537            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1538            Si = parmDict[phfx+'Size;i']
1539            Sa = parmDict[phfx+'Size;a']
1540            gami = const/(Si*Sa)
1541            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1542            Sgam = gami*sqtrm
1543            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1544            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1545            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1546            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1547            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1548            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1549        else:           #OK  ellipsoidal crystallites
1550            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1551            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1552            H = np.array(refl[:3])
1553            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1554            Sgam = const/lenR
1555            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1556                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1557                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1558        gamDict[phfx+'Size;mx'] = Sgam  #OK
1559        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1560               
1561        #microstrain derivatives               
1562        if calcControls[phfx+'MustrainType'] == 'isotropic':
1563            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1564            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1565            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1566        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1567            H = np.array(refl[:3])
1568            P = np.array(calcControls[phfx+'MustrainAxis'])
1569            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1570            Si = parmDict[phfx+'Mustrain;i']
1571            Sa = parmDict[phfx+'Mustrain;a']
1572            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1573            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1574            Mgam = gami/sqtrm
1575            dsi = -gami*Si*cosP**2/sqtrm**3
1576            dsa = -gami*Sa*sinP**2/sqtrm**3
1577            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1578            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1579            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1580            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1581        else:       #generalized - P.W. Stephens model OK
1582            pwrs = calcControls[phfx+'MuPwrs']
1583            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1584            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1585            Sum = 0
1586            for i,strm in enumerate(Strms):
1587                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1588                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1589                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1590            Mgam = const*np.sqrt(Sum)
1591            for i in range(len(Strms)):
1592                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1593                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1594        gamDict[phfx+'Mustrain;mx'] = Mgam
1595        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1596       
1597    return sigDict,gamDict
1598       
1599def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1600    'Needs a doc string'
1601    if im:
1602        h,k,l,m = refl[:4]
1603        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1604        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1605    else:
1606        h,k,l = refl[:3]
1607        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1608    refl[4+im] = d
1609    if 'C' in calcControls[hfx+'histType']:
1610        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1611        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1612        if 'Bragg' in calcControls[hfx+'instType']:
1613            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1614                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1615        else:               #Debye-Scherrer - simple but maybe not right
1616            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1617    elif 'T' in calcControls[hfx+'histType']:
1618        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1619        #do I need sample position effects - maybe?
1620    return pos
1621
1622def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1623    'Needs a doc string'
1624    dpr = 180./np.pi
1625    if im:
1626        h,k,l,m = refl[:4]
1627        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1628        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1629    else:
1630        m = 0
1631        h,k,l = refl[:3]       
1632        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1633    dst = np.sqrt(dstsq)
1634    dsp = 1./dst
1635    if 'C' in calcControls[hfx+'histType']:
1636        pos = refl[5+im]-parmDict[hfx+'Zero']
1637        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1638        dpdw = const*dst
1639        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1640        dpdZ = 1.0
1641        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1642            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1643        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1644        if 'Bragg' in calcControls[hfx+'instType']:
1645            dpdSh = -4.*shft*cosd(pos/2.0)
1646            dpdTr = -shft*sind(pos)*100.0
1647            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1648        else:               #Debye-Scherrer - simple but maybe not right
1649            dpdXd = -shft*cosd(pos)
1650            dpdYd = -shft*sind(pos)
1651            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1652    elif 'T' in calcControls[hfx+'histType']:
1653        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1654        dpdZ = 1.0
1655        dpdDC = dsp
1656        dpdDA = dsp**2
1657        dpdDB = 1./dsp
1658        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1659            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1660        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1661           
1662def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1663    'Needs a doc string'
1664    laue = SGData['SGLaue']
1665    uniq = SGData['SGUniq']
1666    h,k,l = refl[:3]
1667    if laue in ['m3','m3m']:
1668        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1669            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1670    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1671        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1672    elif laue in ['3R','3mR']:
1673        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1674    elif laue in ['4/m','4/mmm']:
1675        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1676    elif laue in ['mmm']:
1677        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1678    elif laue in ['2/m']:
1679        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1680        if uniq == 'a':
1681            Dij += parmDict[phfx+'D23']*k*l
1682        elif uniq == 'b':
1683            Dij += parmDict[phfx+'D13']*h*l
1684        elif uniq == 'c':
1685            Dij += parmDict[phfx+'D12']*h*k
1686    else:
1687        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1688            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1689    if 'C' in calcControls[hfx+'histType']:
1690        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1691    else:
1692        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1693           
1694def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1695    'Needs a doc string'
1696    laue = SGData['SGLaue']
1697    uniq = SGData['SGUniq']
1698    h,k,l = refl[:3]
1699    if laue in ['m3','m3m']:
1700        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1701            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1702    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1703        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1704    elif laue in ['3R','3mR']:
1705        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1706    elif laue in ['4/m','4/mmm']:
1707        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1708    elif laue in ['mmm']:
1709        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1710    elif laue in ['2/m']:
1711        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1712        if uniq == 'a':
1713            dDijDict[phfx+'D23'] = k*l
1714        elif uniq == 'b':
1715            dDijDict[phfx+'D13'] = h*l
1716        elif uniq == 'c':
1717            dDijDict[phfx+'D12'] = h*k
1718            names.append()
1719    else:
1720        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1721            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1722    if 'C' in calcControls[hfx+'histType']:
1723        for item in dDijDict:
1724            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1725    else:
1726        for item in dDijDict:
1727            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1728    return dDijDict
1729   
1730def GetDij(phfx,SGData,parmDict):
1731    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1732    return G2spc.HStrainVals(HSvals,SGData)
1733               
1734def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1735    'Needs a doc string'
1736    histoList = Histograms.keys()
1737    histoList.sort()
1738    for histogram in histoList:
1739        if 'PWDR' in histogram[:4]:
1740            Histogram = Histograms[histogram]
1741            hId = Histogram['hId']
1742            hfx = ':%d:'%(hId)
1743            Limits = calcControls[hfx+'Limits']
1744            if 'C' in calcControls[hfx+'histType']:
1745                shl = max(parmDict[hfx+'SH/L'],0.0005)
1746                Ka2 = False
1747                kRatio = 0.0
1748                if hfx+'Lam1' in parmDict.keys():
1749                    Ka2 = True
1750                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1751                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1752            x,y,w,yc,yb,yd = Histogram['Data']
1753            xB = np.searchsorted(x,Limits[0])
1754            xF = np.searchsorted(x,Limits[1])
1755            ymb = np.array(y-yb)
1756            ymb = np.where(ymb,ymb,1.0)
1757            ycmb = np.array(yc-yb)
1758            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1759            refLists = Histogram['Reflection Lists']
1760            for phase in refLists:
1761                if phase not in Phases:     #skips deleted or renamed phases silently!
1762                    continue
1763                Phase = Phases[phase]
1764                im = 0
1765                if Phase['General']['Type'] in ['modulated','magnetic']:
1766                    im = 1
1767                pId = Phase['pId']
1768                phfx = '%d:%d:'%(pId,hId)
1769                refDict = refLists[phase]
1770                sumFo = 0.0
1771                sumdF = 0.0
1772                sumFosq = 0.0
1773                sumdFsq = 0.0
1774                for refl in refDict['RefList']:
1775                    if 'C' in calcControls[hfx+'histType']:
1776                        yp = np.zeros_like(yb)
1777                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1778                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1779                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1780                        iFin2 = iFin
1781                        if not iBeg+iFin:       #peak below low limit - skip peak
1782                            continue
1783                        elif not iBeg-iFin:     #peak above high limit - done
1784                            break
1785                        elif iBeg < iFin:
1786                            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
1787                            if Ka2:
1788                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1789                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1790                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1791                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1792                                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
1793                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1794                    elif 'T' in calcControls[hfx+'histType']:
1795                        yp = np.zeros_like(yb)
1796                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1797                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1798                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1799                        if iBeg < iFin:
1800                            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
1801                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1802                    Fo = np.sqrt(np.abs(refl[8+im]))
1803                    Fc = np.sqrt(np.abs(refl[9]+im))
1804                    sumFo += Fo
1805                    sumFosq += refl[8+im]**2
1806                    sumdF += np.abs(Fo-Fc)
1807                    sumdFsq += (refl[8+im]-refl[9+im])**2
1808                if sumFo:
1809                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1810                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1811                else:
1812                    Histogram['Residuals'][phfx+'Rf'] = 100.
1813                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1814                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1815                Histogram['Residuals']['hId'] = hId
1816        elif 'HKLF' in histogram[:4]:
1817            Histogram = Histograms[histogram]
1818            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1819               
1820def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1821    'Needs a doc string'
1822   
1823    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1824        U = parmDict[hfx+'U']
1825        V = parmDict[hfx+'V']
1826        W = parmDict[hfx+'W']
1827        X = parmDict[hfx+'X']
1828        Y = parmDict[hfx+'Y']
1829        tanPos = tand(refl[5+im]/2.0)
1830        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1831        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1832        sig = max(0.001,sig)
1833        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1834        gam = max(0.001,gam)
1835        return sig,gam
1836               
1837    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1838        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1839            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1840        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1841        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1842        sig += Ssig
1843        gam += Sgam
1844        return sig,gam
1845       
1846    def GetReflAlpBet(refl,im,hfx,parmDict):
1847        alp = parmDict[hfx+'alpha']/refl[4+im]
1848        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1849        return alp,bet
1850       
1851    hId = Histogram['hId']
1852    hfx = ':%d:'%(hId)
1853    bakType = calcControls[hfx+'bakType']
1854    yb = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1855    yc = np.zeros_like(yb)
1856    cw = np.diff(x)
1857    cw = np.append(cw,cw[-1])
1858       
1859    if 'C' in calcControls[hfx+'histType']:   
1860        shl = max(parmDict[hfx+'SH/L'],0.002)
1861        Ka2 = False
1862        if hfx+'Lam1' in parmDict.keys():
1863            wave = parmDict[hfx+'Lam1']
1864            Ka2 = True
1865            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1866            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1867        else:
1868            wave = parmDict[hfx+'Lam']
1869    for phase in Histogram['Reflection Lists']:
1870        refDict = Histogram['Reflection Lists'][phase]
1871        if phase not in Phases:     #skips deleted or renamed phases silently!
1872            continue
1873        Phase = Phases[phase]
1874        pId = Phase['pId']
1875        pfx = '%d::'%(pId)
1876        phfx = '%d:%d:'%(pId,hId)
1877        hfx = ':%d:'%(hId)
1878        SGData = Phase['General']['SGData']
1879        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1880        im = 0
1881        if Phase['General']['Type'] in ['modulated','magnetic']:
1882            SSGData = Phase['General']['SSGData']
1883            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1884            im = 1  #offset in SS reflection list
1885            #??
1886        Dij = GetDij(phfx,SGData,parmDict)
1887        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
1888        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1889        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1890        Vst = np.sqrt(nl.det(G))    #V*
1891        if not Phase['General'].get('doPawley'):
1892            time0 = time.time()
1893            if im:
1894                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
1895            else:
1896                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1897#            print 'sf calc time: %.3fs'%(time.time()-time0)
1898        time0 = time.time()
1899        badPeak = False
1900        for iref,refl in enumerate(refDict['RefList']):
1901            if 'C' in calcControls[hfx+'histType']:
1902                if im:
1903                    h,k,l,m = refl[:4]
1904                else:
1905                    h,k,l = refl[:3]
1906                Uniq = np.inner(refl[:3],SGMT)
1907                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
1908                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
1909#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1910                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1911                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1912                refl[11+im] *= Vst*Lorenz
1913                 
1914                if Phase['General'].get('doPawley'):
1915                    try:
1916                        if im:
1917                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
1918                        else:
1919                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1920                        refl[9+im] = parmDict[pInd]
1921                    except KeyError:
1922#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1923                        continue
1924                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1925                iBeg = np.searchsorted(x,refl[5+im]-fmin)
1926                iFin = np.searchsorted(x,refl[5+im]+fmax)
1927                if not iBeg+iFin:       #peak below low limit - skip peak
1928                    continue
1929                elif not iBeg-iFin:     #peak above high limit - done
1930                    break
1931                elif iBeg > iFin:   #bad peak coeff - skip
1932                    badPeak = True
1933                    continue
1934                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
1935                if Ka2:
1936                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1937                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1938                    iBeg = np.searchsorted(x,pos2-fmin)
1939                    iFin = np.searchsorted(x,pos2+fmax)
1940                    if not iBeg+iFin:       #peak below low limit - skip peak
1941                        continue
1942                    elif not iBeg-iFin:     #peak above high limit - done
1943                        return yc,yb
1944                    elif iBeg > iFin:   #bad peak coeff - skip
1945                        continue
1946                    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
1947            elif 'T' in calcControls[hfx+'histType']:
1948                h,k,l = refl[:3]
1949                Uniq = np.inner(refl[:3],SGMT)
1950                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
1951                Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4+im]**4                                                #TOF Lorentz correction
1952#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1953                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1954                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
1955                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1956                refl[11+im] *= Vst*Lorenz
1957                if Phase['General'].get('doPawley'):
1958                    try:
1959                        if im:
1960                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
1961                        else:
1962                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1963                        refl[9+im] = parmDict[pInd]
1964                    except KeyError:
1965#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1966                        continue
1967                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1968                iBeg = np.searchsorted(x,refl[5+im]-fmin)
1969                iFin = np.searchsorted(x,refl[5+im]+fmax)
1970                if not iBeg+iFin:       #peak below low limit - skip peak
1971                    continue
1972                elif not iBeg-iFin:     #peak above high limit - done
1973                    break
1974                elif iBeg > iFin:   #bad peak coeff - skip
1975                    badPeak = True
1976                    continue
1977                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]
1978#        print 'profile calc time: %.3fs'%(time.time()-time0)
1979    if badPeak:
1980        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
1981    return yc,yb
1982   
1983def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1984    'Needs a doc string'
1985   
1986    def cellVaryDerv(pfx,SGData,dpdA): 
1987        if SGData['SGLaue'] in ['-1',]:
1988            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1989                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1990        elif SGData['SGLaue'] in ['2/m',]:
1991            if SGData['SGUniq'] == 'a':
1992                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1993            elif SGData['SGUniq'] == 'b':
1994                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1995            else:
1996                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1997        elif SGData['SGLaue'] in ['mmm',]:
1998            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1999        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2000            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2001        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2002            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2003        elif SGData['SGLaue'] in ['3R', '3mR']:
2004            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2005        elif SGData['SGLaue'] in ['m3m','m3']:
2006            return [[pfx+'A0',dpdA[0]]]
2007           
2008    # create a list of dependent variables and set up a dictionary to hold their derivatives
2009    dependentVars = G2mv.GetDependentVars()
2010    depDerivDict = {}
2011    for j in dependentVars:
2012        depDerivDict[j] = np.zeros(shape=(len(x)))
2013    #print 'dependent vars',dependentVars
2014    lenX = len(x)               
2015    hId = Histogram['hId']
2016    hfx = ':%d:'%(hId)
2017    bakType = calcControls[hfx+'bakType']
2018    dMdv = np.zeros(shape=(len(varylist),len(x)))
2019    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2020    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2021        bBpos =varylist.index(hfx+'Back;0')
2022        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2023    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2024    for name in varylist:
2025        if 'Debye' in name:
2026            id = int(name.split(';')[-1])
2027            parm = name[:int(name.rindex(';'))]
2028            ip = names.index(parm)
2029            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2030    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2031    for name in varylist:
2032        if 'BkPk' in name:
2033            parm,id = name.split(';')
2034            id = int(id)
2035            if parm in names:
2036                ip = names.index(parm)
2037                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2038    cw = np.diff(x)
2039    cw = np.append(cw,cw[-1])
2040    Ka2 = False #also for TOF!
2041    if 'C' in calcControls[hfx+'histType']:   
2042        shl = max(parmDict[hfx+'SH/L'],0.002)
2043        if hfx+'Lam1' in parmDict.keys():
2044            wave = parmDict[hfx+'Lam1']
2045            Ka2 = True
2046            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2047            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2048        else:
2049            wave = parmDict[hfx+'Lam']
2050    for phase in Histogram['Reflection Lists']:
2051        refDict = Histogram['Reflection Lists'][phase]
2052        if phase not in Phases:     #skips deleted or renamed phases silently!
2053            continue
2054        Phase = Phases[phase]
2055        SGData = Phase['General']['SGData']
2056        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2057        im = 0
2058        if Phase['General']['Type'] in ['modulated','magnetic']:
2059            SSGData = Phase['General']['SSGData']
2060            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2061            im = 1  #offset in SS reflection list
2062            #??
2063        pId = Phase['pId']
2064        pfx = '%d::'%(pId)
2065        phfx = '%d:%d:'%(pId,hId)
2066        Dij = GetDij(phfx,SGData,parmDict)
2067        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2068        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2069        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2070        if not Phase['General'].get('doPawley'):
2071            time0 = time.time()
2072            if im:
2073                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2074            else:
2075                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2076#            print 'sf-derv time %.3fs'%(time.time()-time0)
2077            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2078        time0 = time.time()
2079        for iref,refl in enumerate(refDict['RefList']):
2080            if im:
2081                h,k,l,m = refl[:4]
2082            else:
2083                h,k,l = refl[:3]
2084            Uniq = np.inner(refl[:3],SGMT)
2085            if 'T' in calcControls[hfx+'histType']:
2086                wave = refl[14+im]
2087            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2088            if 'C' in calcControls[hfx+'histType']:        #CW powder
2089                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2090            else: #'T'OF
2091                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2092            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2093            iFin = np.searchsorted(x,refl[5+im]+fmax)
2094            if not iBeg+iFin:       #peak below low limit - skip peak
2095                continue
2096            elif not iBeg-iFin:     #peak above high limit - done
2097                break
2098            pos = refl[5+im]
2099            if 'C' in calcControls[hfx+'histType']:
2100                tanth = tand(pos/2.0)
2101                costh = cosd(pos/2.0)
2102                lenBF = iFin-iBeg
2103                dMdpk = np.zeros(shape=(6,lenBF))
2104                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2105                for i in range(5):
2106                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2107                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2108                if Ka2:
2109                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2110                    iBeg2 = np.searchsorted(x,pos2-fmin)
2111                    iFin2 = np.searchsorted(x,pos2+fmax)
2112                    if iBeg2-iFin2:
2113                        lenBF2 = iFin2-iBeg2
2114                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2115                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2116                        for i in range(5):
2117                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2118                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2119                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2120            else:   #'T'OF
2121                lenBF = iFin-iBeg
2122                if lenBF < 0:   #bad peak coeff
2123                    break
2124                dMdpk = np.zeros(shape=(6,lenBF))
2125                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2126                for i in range(6):
2127                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2128                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2129            if Phase['General'].get('doPawley'):
2130                dMdpw = np.zeros(len(x))
2131                try:
2132                    if im:
2133                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2134                    else:
2135                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2136                    idx = varylist.index(pIdx)
2137                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2138                    if Ka2: #not for TOF either
2139                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2140                    dMdv[idx] = dMdpw
2141                except: # ValueError:
2142                    pass
2143            if 'C' in calcControls[hfx+'histType']:
2144                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2145                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2146                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2147                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2148                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2149                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2150                    hfx+'DisplaceY':[dpdY,'pos'],}
2151                if 'Bragg' in calcControls[hfx+'instType']:
2152                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2153                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2154                else:
2155                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2156            else:   #'T'OF
2157                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2158                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2159                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2160                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2161                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2162                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2163                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2164                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2165            for name in names:
2166                item = names[name]
2167                if name in varylist:
2168                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2169                    if Ka2:
2170                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2171                elif name in dependentVars:
2172                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2173                    if Ka2:
2174                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2175            for iPO in dIdPO:
2176                if iPO in varylist:
2177                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2178                    if Ka2:
2179                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2180                elif iPO in dependentVars:
2181                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2182                    if Ka2:
2183                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2184            for i,name in enumerate(['omega','chi','phi']):
2185                aname = pfx+'SH '+name
2186                if aname in varylist:
2187                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2188                    if Ka2:
2189                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2190                elif aname in dependentVars:
2191                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2192                    if Ka2:
2193                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2194            for iSH in dFdODF:
2195                if iSH in varylist:
2196                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2197                    if Ka2:
2198                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2199                elif iSH in dependentVars:
2200                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2201                    if Ka2:
2202                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2203            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2204            for name,dpdA in cellDervNames:
2205                if name in varylist:
2206                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2207                    if Ka2:
2208                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2209                elif name in dependentVars:
2210                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2211                    if Ka2:
2212                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2213            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2214            for name in dDijDict:
2215                if name in varylist:
2216                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2217                    if Ka2:
2218                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2219                elif name in dependentVars:
2220                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2221                    if Ka2:
2222                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2223            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2224                if name in varylist:
2225                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2226                    if Ka2:
2227                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2228                elif name in dependentVars:
2229                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2230                    if Ka2:
2231                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2232            if 'C' in calcControls[hfx+'histType']:
2233                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2234            else:   #'T'OF
2235                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2236            for name in gamDict:
2237                if name in varylist:
2238                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2239                    if Ka2:
2240                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2241                elif name in dependentVars:
2242                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2243                    if Ka2:
2244                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2245            for name in sigDict:
2246                if name in varylist:
2247                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2248                    if Ka2:
2249                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2250                elif name in dependentVars:
2251                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2252                    if Ka2:
2253                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2254            for name in ['BabA','BabU']:
2255                if refl[9+im]:
2256                    if phfx+name in varylist:
2257                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2258                        if Ka2:
2259                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2260                    elif phfx+name in dependentVars:                   
2261                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2262                        if Ka2:
2263                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2264            if not Phase['General'].get('doPawley'):
2265                #do atom derivatives -  for RB,F,X & U so far             
2266                corr = dervDict['int']/refl[9+im]
2267                if Ka2:
2268                    corr2 = dervDict2['int']/refl[9+im]
2269                for name in varylist+dependentVars:
2270                    if '::RBV;' in name:
2271                        pass
2272                    else:
2273                        try:
2274                            aname = name.split(pfx)[1][:2]
2275                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2276                        except IndexError:
2277                            continue
2278                    if name in varylist:
2279                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2280                        if Ka2:
2281                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2282                    elif name in dependentVars:
2283                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2284                        if Ka2:
2285                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2286    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2287    # now process derivatives in constraints
2288    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2289    return dMdv
2290   
2291def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2292    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2293    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2294    are returned to either dervRefine or HessRefine.
2295
2296    :returns:
2297    '''
2298    nobs = Histogram['Residuals']['Nobs']
2299    hId = Histogram['hId']
2300    hfx = ':%d:'%(hId)
2301    pfx = '%d::'%(Phase['pId'])
2302    phfx = '%d:%d:'%(Phase['pId'],hId)
2303    SGData = Phase['General']['SGData']
2304    im = 0
2305    if Phase['General']['Type'] in ['modulated','magnetic']:
2306        SSGData = Phase['General']['SSGData']
2307        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2308        im = 1  #offset in SS reflection list
2309        #??
2310    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2311    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2312    refDict = Histogram['Data']
2313    if im:
2314        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2315    else:
2316        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2317    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2318    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2319    dependentVars = G2mv.GetDependentVars()
2320    depDerivDict = {}
2321    for j in dependentVars:
2322        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2323    wdf = np.zeros(len(refDict['RefList']))
2324    if calcControls['F**2']:
2325        for iref,ref in enumerate(refDict['RefList']):
2326            if ref[6+im] > 0:
2327                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1] 
2328                w = 1.0/ref[6+im]
2329                if w*ref[5+im] >= calcControls['minF/sig']:
2330                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2331                    for j,var in enumerate(varylist):
2332                        if var in dFdvDict:
2333                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2334                    for var in dependentVars:
2335                        if var in dFdvDict:
2336                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2337                    if phfx+'Scale' in varylist:
2338                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2339                    elif phfx+'Scale' in dependentVars:
2340                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]
2341                    for item in ['Ep','Es','Eg']:
2342                        if phfx+item in varylist and phfx+item in dervDict:
2343                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2344                        elif phfx+item in dependentVars and phfx+item in dervDict:
2345                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2346                    for item in ['BabA','BabU']:
2347                        if phfx+item in varylist:
2348                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2349                        elif phfx+item in dependentVars:
2350                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2351    else:
2352        for iref,ref in enumerate(refDict['RefList']):
2353            if ref[5+im] > 0.:
2354                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2355                Fo = np.sqrt(ref[5+im])
2356                Fc = np.sqrt(ref[7+im])
2357                w = 1.0/ref[6+im]
2358                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2359                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2360                    for j,var in enumerate(varylist):
2361                        if var in dFdvDict:
2362                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2363                    for var in dependentVars:
2364                        if var in dFdvDict:
2365                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2366                    if phfx+'Scale' in varylist:
2367                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2368                    elif phfx+'Scale' in dependentVars:
2369                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]                           
2370                    for item in ['Ep','Es','Eg']:
2371                        if phfx+item in varylist and phfx+item in dervDict:
2372                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #correct
2373                        elif phfx+item in dependentVars and phfx+item in dervDict:
2374                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2375                    for item in ['BabA','BabU']:
2376                        if phfx+item in varylist:
2377                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2378                        elif phfx+item in dependentVars:
2379                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2380    return dMdvh,depDerivDict,wdf
2381   
2382
2383def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2384    '''Loop over histograms and compute derivatives of the fitting
2385    model (M) with respect to all parameters.  Results are returned in
2386    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2387    n is the number of parameters and m is the number of data
2388    points. This can exceed memory when m gets large. This routine is
2389    used when refinement derivatives are selected as "analtytic
2390    Jacobian" in Controls.
2391
2392    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2393    '''
2394    parmDict.update(zip(varylist,values))
2395    G2mv.Dict2Map(parmDict,varylist)
2396    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2397    nvar = len(varylist)
2398    dMdv = np.empty(0)
2399    histoList = Histograms.keys()
2400    histoList.sort()
2401    for histogram in histoList:
2402        if 'PWDR' in histogram[:4]:
2403            Histogram = Histograms[histogram]
2404            hId = Histogram['hId']
2405            hfx = ':%d:'%(hId)
2406            wtFactor = calcControls[hfx+'wtFactor']
2407            Limits = calcControls[hfx+'Limits']
2408            x,y,w,yc,yb,yd = Histogram['Data']
2409            xB = np.searchsorted(x,Limits[0])
2410            xF = np.searchsorted(x,Limits[1])
2411            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2412                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2413        elif 'HKLF' in histogram[:4]:
2414            Histogram = Histograms[histogram]
2415            phase = Histogram['Reflection Lists']
2416            Phase = Phases[phase]
2417            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2418            hfx = ':%d:'%(Histogram['hId'])
2419            wtFactor = calcControls[hfx+'wtFactor']
2420            # now process derivatives in constraints
2421            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2422        else:
2423            continue        #skip non-histogram entries
2424        if len(dMdv):
2425            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2426        else:
2427            dMdv = np.sqrt(wtFactor)*dMdvh
2428           
2429    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2430    if np.any(pVals):
2431        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2432        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2433       
2434    return dMdv
2435
2436def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2437    '''Loop over histograms and compute derivatives of the fitting
2438    model (M) with respect to all parameters.  For each histogram, the
2439    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2440    number of parameters and m is the number of data points *in the
2441    histogram*. The (n by n) Hessian is computed from each Jacobian
2442    and it is returned.  This routine is used when refinement
2443    derivatives are selected as "analtytic Hessian" in Controls.
2444
2445    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2446    '''
2447    parmDict.update(zip(varylist,values))
2448    G2mv.Dict2Map(parmDict,varylist)
2449    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2450    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2451    nvar = len(varylist)
2452    Hess = np.empty(0)
2453    histoList = Histograms.keys()
2454    histoList.sort()
2455    for histogram in histoList:
2456        if 'PWDR' in histogram[:4]:
2457            Histogram = Histograms[histogram]
2458            hId = Histogram['hId']
2459            hfx = ':%d:'%(hId)
2460            wtFactor = calcControls[hfx+'wtFactor']
2461            Limits = calcControls[hfx+'Limits']
2462            x,y,w,yc,yb,yd = Histogram['Data']
2463            W = wtFactor*w
2464            dy = y-yc
2465            xB = np.searchsorted(x,Limits[0])
2466            xF = np.searchsorted(x,Limits[1])
2467            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2468                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2469            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2470            Dy = dy[xB:xF][np.newaxis,:]
2471            dMdvh *= Wt
2472            if dlg:
2473                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2474            if len(Hess):
2475                Hess += np.inner(dMdvh,dMdvh)
2476                dMdvh *= Wt*Dy
2477                Vec += np.sum(dMdvh,axis=1)
2478            else:
2479                Hess = np.inner(dMdvh,dMdvh)
2480                dMdvh *= Wt*Dy
2481                Vec = np.sum(dMdvh,axis=1)
2482        elif 'HKLF' in histogram[:4]:
2483            Histogram = Histograms[histogram]
2484            phase = Histogram['Reflection Lists']
2485            Phase = Phases[phase]
2486            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2487            hId = Histogram['hId']
2488            hfx = ':%d:'%(Histogram['hId'])
2489            wtFactor = calcControls[hfx+'wtFactor']
2490            # now process derivatives in constraints
2491            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2492#            print 'matrix build time: %.3f'%(time.time()-time0)
2493
2494            if dlg:
2495                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2496            if len(Hess):
2497                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2498                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2499            else:
2500                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2501                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2502        else:
2503            continue        #skip non-histogram entries
2504    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2505    if np.any(pVals):
2506        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2507        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2508        Hess += np.inner(dpdv*pWt,dpdv)
2509    return Vec,Hess
2510
2511def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2512    '''Computes the point-by-point discrepancies between every data point in every histogram
2513    and the observed value
2514   
2515    :returns: an np array of differences between observed and computed diffraction values.
2516    '''
2517    Values2Dict(parmDict, varylist, values)
2518    G2mv.Dict2Map(parmDict,varylist)
2519    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2520    M = np.empty(0)
2521    SumwYo = 0
2522    Nobs = 0
2523    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2524    histoList = Histograms.keys()
2525    histoList.sort()
2526    for histogram in histoList:
2527        if 'PWDR' in histogram[:4]:
2528            Histogram = Histograms[histogram]
2529            hId = Histogram['hId']
2530            hfx = ':%d:'%(hId)
2531            wtFactor = calcControls[hfx+'wtFactor']
2532            Limits = calcControls[hfx+'Limits']
2533            x,y,w,yc,yb,yd = Histogram['Data']
2534            yc *= 0.0                           #zero full calcd profiles
2535            yb *= 0.0
2536            yd *= 0.0
2537            xB = np.searchsorted(x,Limits[0])
2538            xF = np.searchsorted(x,Limits[1])
2539            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2540                varylist,Histogram,Phases,calcControls,pawleyLookup)
2541            yc[xB:xF] += yb[xB:xF]
2542            if not np.any(y):                   #fill dummy data
2543                rv = st.poisson(yc[xB:xF])
2544                y[xB:xF] = rv.rvs()
2545                Z = np.ones_like(yc[xB:xF])
2546                Z[1::2] *= -1
2547                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2548                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2549            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2550            W = wtFactor*w
2551            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2552            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2553            Nobs += Histogram['Residuals']['Nobs']
2554            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2555            SumwYo += Histogram['Residuals']['sumwYo']
2556            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2557            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2558            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2559            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2560            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.))
2561            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.))
2562            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2563            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2564            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2565            if dlg:
2566                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2567            M = np.concatenate((M,wdy))
2568#end of PWDR processing
2569        elif 'HKLF' in histogram[:4]:
2570            Histogram = Histograms[histogram]
2571            Histogram['Residuals'] = {}
2572            phase = Histogram['Reflection Lists']
2573            Phase = Phases[phase]
2574            hId = Histogram['hId']
2575            hfx = ':%d:'%(hId)
2576            wtFactor = calcControls[hfx+'wtFactor']
2577            pfx = '%d::'%(Phase['pId'])
2578            phfx = '%d:%d:'%(Phase['pId'],hId)
2579            SGData = Phase['General']['SGData']
2580            im = 0
2581            if Phase['General']['Type'] in ['modulated','magnetic']:
2582                SSGData = Phase['General']['SSGData']
2583                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2584                im = 1  #offset in SS reflection list
2585                #??
2586            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2587            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2588            refDict = Histogram['Data']
2589            time0 = time.time()
2590            if im:
2591                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2592            else:
2593                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2594#            print 'sf-calc time: %.3f'%(time.time()-time0)
2595            df = np.zeros(len(refDict['RefList']))
2596            sumwYo = 0
2597            sumFo = 0
2598            sumFo2 = 0
2599            sumdF = 0
2600            sumdF2 = 0
2601            nobs = 0
2602            if calcControls['F**2']:
2603                for i,ref in enumerate(refDict['RefList']):
2604                    if ref[6+im] > 0:
2605                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2606                        w = 1.0/ref[6+im]
2607                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
2608                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2609                        if w*ref[5+im] >= calcControls['minF/sig']:
2610                            Fo = np.sqrt(ref[5+im])
2611                            sumFo += Fo
2612                            sumFo2 += ref[5+im]
2613                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2614                            sumdF2 += abs(ref[5+im]-ref[7+im])
2615                            nobs += 1
2616                            df[i] = -w*(ref[5+im]-ref[7+im])
2617                            sumwYo += (w*ref[5+im])**2
2618            else:
2619                for i,ref in enumerate(refDict['RefList']):
2620                    if ref[5+im] > 0.:
2621                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2622                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
2623                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2624                        Fo = np.sqrt(ref[5+im])
2625                        Fc = np.sqrt(ref[7+im])
2626                        w = 2.0*Fo/ref[6+im]
2627                        if w*Fo >= calcControls['minF/sig']:
2628                            sumFo += Fo
2629                            sumFo2 += ref[5+im]
2630                            sumdF += abs(Fo-Fc)
2631                            sumdF2 += abs(ref[5+im]-ref[7+im])
2632                            nobs += 1
2633                            df[i] = -w*(Fo-Fc)
2634                            sumwYo += (w*Fo)**2
2635            Histogram['Residuals']['Nobs'] = nobs
2636            Histogram['Residuals']['sumwYo'] = sumwYo
2637            SumwYo += sumwYo
2638            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2639            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2640            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2641            Histogram['Residuals'][phfx+'Nref'] = nobs
2642            Nobs += nobs
2643            if dlg:
2644                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2645            M = np.concatenate((M,wtFactor*df))
2646# end of HKLF processing
2647    Histograms['sumwYo'] = SumwYo
2648    Histograms['Nobs'] = Nobs
2649    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2650    if dlg:
2651        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2652        if not GoOn:
2653            parmDict['saved values'] = values
2654            dlg.Destroy()
2655            raise Exception         #Abort!!
2656    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2657    if len(pVals):
2658        pSum = np.sum(pWt*pVals**2)
2659        for name in pWsum:
2660            if pWsum:
2661                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2662        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2663        Nobs += len(pVals)
2664        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2665    return M
2666                       
Note: See TracBrowser for help on using the repository browser.