source: trunk/GSASIIstrMath.py @ 1766

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

correct old LANSCE/HIPD import to give correct Omega, chi & phi from goniometer
correct 2-theta from TOF & use abs(Tth) as needed in Lorentz calc for TOF
some corrections to texture calcs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 134.6 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-03-27 19:46:30 +0000 (Fri, 27 Mar 2015) $
8# $Author: vondreele $
9# $Revision: 1766 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1766 2015-03-27 19:46:30Z vondreele $
12########### SVN repository information ###################
13import time
14import math
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.optimize as so
20import scipy.stats as st
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 1766 $")
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],abs(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],abs(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(abs(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    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1327    pi2 = np.sqrt(2./np.pi)
1328    if 'T' in calcControls[hfx+'histType']:
1329        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1330        wave = refl[14+im]
1331    else:   #'C'W
1332        sth2 = sind(refl[5+im]/2.)**2
1333        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1334    c2th = 1.-2.0*sth2
1335    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1336    if 'X' in calcControls[hfx+'histType']:
1337        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1338    xfac = flv2*parmDict[phfx+'Extinction']
1339    dbde = -500.*flv2
1340    if xfac > -1.:
1341        dbde = -flv2/(1.+xfac)**3
1342    dlde = 0.
1343    if 0 < xfac <= 1.:
1344        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1345        dlde = np.sum(xn*coef)
1346    elif xfac > 1.:
1347        xfac2 = 1./np.sqrt(xfac)
1348        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1349       
1350    return dbde*sth2+dlde*(1.-sth2)
1351
1352
1353#    delt = 0.01
1354#    parmDict[phfx+'Extinction'] += delt
1355#    plus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1356#    parmDict[phfx+'Extinction'] -= 2.*delt
1357#    minus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1358#    parmDict[phfx+'Extinction'] += delt
1359#    return (plus-minus)/(2.*delt)   
1360   
1361def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1362    'Needs a doc string'    #need powder extinction!
1363    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1364    if 'X' in parmDict[hfx+'Type']:
1365        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1366    POcorr = 1.0
1367    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture
1368        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1369    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1370        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1371    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1372        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1373    Icorr *= POcorr
1374    AbsCorr = 1.0
1375    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1376    Icorr *= AbsCorr
1377    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1378    Icorr *= ExtCorr
1379    return Icorr,POcorr,AbsCorr,ExtCorr
1380   
1381def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1382    'Needs a doc string'    #need powder extinction derivs!
1383    dIdsh = 1./parmDict[hfx+'Scale']
1384    dIdsp = 1./parmDict[phfx+'Scale']
1385    if 'X' in parmDict[hfx+'Type']:
1386        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1387        dIdPola /= pola
1388    else:       #'N'
1389        dIdPola = 0.0
1390    dFdODF = {}
1391    dFdSA = [0,0,0]
1392    dIdPO = {}
1393    if pfx+'SHorder' in parmDict:
1394        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1395        for iSH in dFdODF:
1396            dFdODF[iSH] /= odfCor
1397        for i in range(3):
1398            dFdSA[i] /= odfCor
1399    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1400        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1401        for iPO in dIdPO:
1402            dIdPO[iPO] /= POcorr
1403    if 'T' in parmDict[hfx+'Type']:
1404        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1405        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1406    else:
1407        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1408        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1409    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1410       
1411def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1412    'Needs a doc string'
1413    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1414        costh = cosd(refl[5+im]/2.)
1415        #crystallite size
1416        if calcControls[phfx+'SizeType'] == 'isotropic':
1417            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1418        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1419            H = np.array(refl[:3])
1420            P = np.array(calcControls[phfx+'SizeAxis'])
1421            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1422            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1423            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1424        else:           #ellipsoidal crystallites
1425            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1426            H = np.array(refl[:3])
1427            lenR = G2pwd.ellipseSize(H,Sij,GB)
1428            Sgam = 1.8*wave/(np.pi*costh*lenR)
1429        #microstrain               
1430        if calcControls[phfx+'MustrainType'] == 'isotropic':
1431            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1432        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1433            H = np.array(refl[:3])
1434            P = np.array(calcControls[phfx+'MustrainAxis'])
1435            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1436            Si = parmDict[phfx+'Mustrain;i']
1437            Sa = parmDict[phfx+'Mustrain;a']
1438            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1439        else:       #generalized - P.W. Stephens model
1440            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1441            Sum = 0
1442            for i,strm in enumerate(Strms):
1443                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1444            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1445    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1446        #crystallite size
1447        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1448            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1449        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1450            H = np.array(refl[:3])
1451            P = np.array(calcControls[phfx+'SizeAxis'])
1452            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1453            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1454            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1455        else:           #ellipsoidal crystallites   #OK
1456            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1457            H = np.array(refl[:3])
1458            lenR = G2pwd.ellipseSize(H,Sij,GB)
1459            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1460        #microstrain               
1461        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1462            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1463        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1464            H = np.array(refl[:3])
1465            P = np.array(calcControls[phfx+'MustrainAxis'])
1466            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1467            Si = parmDict[phfx+'Mustrain;i']
1468            Sa = parmDict[phfx+'Mustrain;a']
1469            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1470        else:       #generalized - P.W. Stephens model  OK
1471            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1472            Sum = 0
1473            for i,strm in enumerate(Strms):
1474                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1475            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1476           
1477    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1478    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1479    sig /= ateln2
1480    return sig,gam
1481       
1482def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1483    'Needs a doc string'
1484    gamDict = {}
1485    sigDict = {}
1486    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1487        costh = cosd(refl[5+im]/2.)
1488        tanth = tand(refl[5+im]/2.)
1489        #crystallite size derivatives
1490        if calcControls[phfx+'SizeType'] == 'isotropic':
1491            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1492            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1493            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1494        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1495            H = np.array(refl[:3])
1496            P = np.array(calcControls[phfx+'SizeAxis'])
1497            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1498            Si = parmDict[phfx+'Size;i']
1499            Sa = parmDict[phfx+'Size;a']
1500            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1501            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1502            Sgam = gami*sqtrm
1503            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1504            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1505            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1506            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1507            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1508            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1509        else:           #ellipsoidal crystallites
1510            const = 1.8*wave/(np.pi*costh)
1511            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1512            H = np.array(refl[:3])
1513            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1514            Sgam = const/lenR
1515            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1516                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1517                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1518        gamDict[phfx+'Size;mx'] = Sgam
1519        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1520               
1521        #microstrain derivatives               
1522        if calcControls[phfx+'MustrainType'] == 'isotropic':
1523            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1524            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1525            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1526        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1527            H = np.array(refl[:3])
1528            P = np.array(calcControls[phfx+'MustrainAxis'])
1529            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1530            Si = parmDict[phfx+'Mustrain;i']
1531            Sa = parmDict[phfx+'Mustrain;a']
1532            gami = 0.018*Si*Sa*tanth/np.pi
1533            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1534            Mgam = gami/sqtrm
1535            dsi = -gami*Si*cosP**2/sqtrm**3
1536            dsa = -gami*Sa*sinP**2/sqtrm**3
1537            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1538            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1539            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1540            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1541        else:       #generalized - P.W. Stephens model
1542            const = 0.018*refl[4+im]**2*tanth/np.pi
1543            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1544            Sum = 0
1545            for i,strm in enumerate(Strms):
1546                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1547                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1548                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1549            Mgam = const*np.sqrt(Sum)
1550            for i in range(len(Strms)):
1551                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1552                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1553        gamDict[phfx+'Mustrain;mx'] = Mgam
1554        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1555    else:   #'T'OF - All checked & OK
1556        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1557            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1558            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1559            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1560        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1561            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1562            H = np.array(refl[:3])
1563            P = np.array(calcControls[phfx+'SizeAxis'])
1564            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1565            Si = parmDict[phfx+'Size;i']
1566            Sa = parmDict[phfx+'Size;a']
1567            gami = const/(Si*Sa)
1568            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1569            Sgam = gami*sqtrm
1570            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1571            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1572            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1573            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1574            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1575            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1576        else:           #OK  ellipsoidal crystallites
1577            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1578            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1579            H = np.array(refl[:3])
1580            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1581            Sgam = const/lenR
1582            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1583                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1584                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1585        gamDict[phfx+'Size;mx'] = Sgam  #OK
1586        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1587               
1588        #microstrain derivatives               
1589        if calcControls[phfx+'MustrainType'] == 'isotropic':
1590            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1591            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1592            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1593        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1594            H = np.array(refl[:3])
1595            P = np.array(calcControls[phfx+'MustrainAxis'])
1596            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1597            Si = parmDict[phfx+'Mustrain;i']
1598            Sa = parmDict[phfx+'Mustrain;a']
1599            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1600            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1601            Mgam = gami/sqtrm
1602            dsi = -gami*Si*cosP**2/sqtrm**3
1603            dsa = -gami*Sa*sinP**2/sqtrm**3
1604            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1605            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1606            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1607            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1608        else:       #generalized - P.W. Stephens model OK
1609            pwrs = calcControls[phfx+'MuPwrs']
1610            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1611            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1612            Sum = 0
1613            for i,strm in enumerate(Strms):
1614                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1615                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1616                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1617            Mgam = const*np.sqrt(Sum)
1618            for i in range(len(Strms)):
1619                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1620                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1621        gamDict[phfx+'Mustrain;mx'] = Mgam
1622        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1623       
1624    return sigDict,gamDict
1625       
1626def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1627    'Needs a doc string'
1628    if im:
1629        h,k,l,m = refl[:4]
1630        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1631        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1632    else:
1633        h,k,l = refl[:3]
1634        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1635    refl[4+im] = d
1636    if 'C' in calcControls[hfx+'histType']:
1637        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1638        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1639        if 'Bragg' in calcControls[hfx+'instType']:
1640            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1641                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1642        else:               #Debye-Scherrer - simple but maybe not right
1643            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1644    elif 'T' in calcControls[hfx+'histType']:
1645        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1646        #do I need sample position effects - maybe?
1647    return pos
1648
1649def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1650    'Needs a doc string'
1651    dpr = 180./np.pi
1652    if im:
1653        h,k,l,m = refl[:4]
1654        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1655        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1656    else:
1657        m = 0
1658        h,k,l = refl[:3]       
1659        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1660    dst = np.sqrt(dstsq)
1661    dsp = 1./dst
1662    if 'C' in calcControls[hfx+'histType']:
1663        pos = refl[5+im]-parmDict[hfx+'Zero']
1664        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1665        dpdw = const*dst
1666        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1667        dpdZ = 1.0
1668        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1669            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1670        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1671        if 'Bragg' in calcControls[hfx+'instType']:
1672            dpdSh = -4.*shft*cosd(pos/2.0)
1673            dpdTr = -shft*sind(pos)*100.0
1674            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1675        else:               #Debye-Scherrer - simple but maybe not right
1676            dpdXd = -shft*cosd(pos)
1677            dpdYd = -shft*sind(pos)
1678            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1679    elif 'T' in calcControls[hfx+'histType']:
1680        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1681        dpdZ = 1.0
1682        dpdDC = dsp
1683        dpdDA = dsp**2
1684        dpdDB = 1./dsp
1685        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1686            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1687        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1688           
1689def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1690    'Needs a doc string'
1691    laue = SGData['SGLaue']
1692    uniq = SGData['SGUniq']
1693    h,k,l = refl[:3]
1694    if laue in ['m3','m3m']:
1695        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1696            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1697    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1698        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1699    elif laue in ['3R','3mR']:
1700        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1701    elif laue in ['4/m','4/mmm']:
1702        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1703    elif laue in ['mmm']:
1704        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1705    elif laue in ['2/m']:
1706        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1707        if uniq == 'a':
1708            Dij += parmDict[phfx+'D23']*k*l
1709        elif uniq == 'b':
1710            Dij += parmDict[phfx+'D13']*h*l
1711        elif uniq == 'c':
1712            Dij += parmDict[phfx+'D12']*h*k
1713    else:
1714        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1715            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1716    if 'C' in calcControls[hfx+'histType']:
1717        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1718    else:
1719        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1720           
1721def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1722    'Needs a doc string'
1723    laue = SGData['SGLaue']
1724    uniq = SGData['SGUniq']
1725    h,k,l = refl[:3]
1726    if laue in ['m3','m3m']:
1727        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1728            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1729    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1730        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1731    elif laue in ['3R','3mR']:
1732        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1733    elif laue in ['4/m','4/mmm']:
1734        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1735    elif laue in ['mmm']:
1736        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1737    elif laue in ['2/m']:
1738        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1739        if uniq == 'a':
1740            dDijDict[phfx+'D23'] = k*l
1741        elif uniq == 'b':
1742            dDijDict[phfx+'D13'] = h*l
1743        elif uniq == 'c':
1744            dDijDict[phfx+'D12'] = h*k
1745    else:
1746        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1747            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1748    if 'C' in calcControls[hfx+'histType']:
1749        for item in dDijDict:
1750            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1751    else:
1752        for item in dDijDict:
1753            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1754    return dDijDict
1755   
1756def GetDij(phfx,SGData,parmDict):
1757    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1758    return G2spc.HStrainVals(HSvals,SGData)
1759               
1760def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1761    'Needs a doc string'
1762    histoList = Histograms.keys()
1763    histoList.sort()
1764    for histogram in histoList:
1765        if 'PWDR' in histogram[:4]:
1766            Histogram = Histograms[histogram]
1767            hId = Histogram['hId']
1768            hfx = ':%d:'%(hId)
1769            Limits = calcControls[hfx+'Limits']
1770            if 'C' in calcControls[hfx+'histType']:
1771                shl = max(parmDict[hfx+'SH/L'],0.0005)
1772                Ka2 = False
1773                kRatio = 0.0
1774                if hfx+'Lam1' in parmDict.keys():
1775                    Ka2 = True
1776                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1777                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1778            x,y,w,yc,yb,yd = Histogram['Data']
1779            xB = np.searchsorted(x,Limits[0])
1780            xF = np.searchsorted(x,Limits[1])
1781            ymb = np.array(y-yb)
1782            ymb = np.where(ymb,ymb,1.0)
1783            ycmb = np.array(yc-yb)
1784            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1785            refLists = Histogram['Reflection Lists']
1786            for phase in refLists:
1787                if phase not in Phases:     #skips deleted or renamed phases silently!
1788                    continue
1789                Phase = Phases[phase]
1790                im = 0
1791                if Phase['General']['Type'] in ['modulated','magnetic']:
1792                    im = 1
1793                pId = Phase['pId']
1794                phfx = '%d:%d:'%(pId,hId)
1795                refDict = refLists[phase]
1796                sumFo = 0.0
1797                sumdF = 0.0
1798                sumFosq = 0.0
1799                sumdFsq = 0.0
1800                sumInt = 0.0
1801                for refl in refDict['RefList']:
1802                    if 'C' in calcControls[hfx+'histType']:
1803                        yp = np.zeros_like(yb)
1804                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1805                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1806                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1807                        iFin2 = iFin
1808                        if not iBeg+iFin:       #peak below low limit - skip peak
1809                            continue
1810                        elif not iBeg-iFin:     #peak above high limit - done
1811                            break
1812                        elif iBeg < iFin:
1813                            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
1814                            sumInt += refl[11+im]*refl[9+im]
1815                            if Ka2:
1816                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1817                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1818                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1819                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1820                                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
1821                                sumInt += refl[11+im]*refl[9+im]*kRatio
1822                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1823                               
1824                    elif 'T' in calcControls[hfx+'histType']:
1825                        yp = np.zeros_like(yb)
1826                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1827                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1828                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1829                        if iBeg < iFin:
1830                            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
1831                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1832                            sumInt += refl[11+im]*refl[9+im]
1833                    Fo = np.sqrt(np.abs(refl[8+im]))
1834                    Fc = np.sqrt(np.abs(refl[9]+im))
1835                    sumFo += Fo
1836                    sumFosq += refl[8+im]**2
1837                    sumdF += np.abs(Fo-Fc)
1838                    sumdFsq += (refl[8+im]-refl[9+im])**2
1839                if sumFo:
1840                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1841                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1842                else:
1843                    Histogram['Residuals'][phfx+'Rf'] = 100.
1844                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1845                Histogram['Residuals'][phfx+'sumInt'] = sumInt
1846                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1847                Histogram['Residuals']['hId'] = hId
1848        elif 'HKLF' in histogram[:4]:
1849            Histogram = Histograms[histogram]
1850            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1851               
1852def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1853    'Needs a doc string'
1854   
1855    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1856        U = parmDict[hfx+'U']
1857        V = parmDict[hfx+'V']
1858        W = parmDict[hfx+'W']
1859        X = parmDict[hfx+'X']
1860        Y = parmDict[hfx+'Y']
1861        tanPos = tand(refl[5+im]/2.0)
1862        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1863        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1864        sig = max(0.001,sig)
1865        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1866        gam = max(0.001,gam)
1867        return sig,gam
1868               
1869    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1870        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1871            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1872        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1873        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1874        sig += Ssig
1875        gam += Sgam
1876        return sig,gam
1877       
1878    def GetReflAlpBet(refl,im,hfx,parmDict):
1879        alp = parmDict[hfx+'alpha']/refl[4+im]
1880        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1881        return alp,bet
1882       
1883    hId = Histogram['hId']
1884    hfx = ':%d:'%(hId)
1885    bakType = calcControls[hfx+'bakType']
1886    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1887    yc = np.zeros_like(yb)
1888    cw = np.diff(x)
1889    cw = np.append(cw,cw[-1])
1890       
1891    if 'C' in calcControls[hfx+'histType']:   
1892        shl = max(parmDict[hfx+'SH/L'],0.002)
1893        Ka2 = False
1894        if hfx+'Lam1' in parmDict.keys():
1895            wave = parmDict[hfx+'Lam1']
1896            Ka2 = True
1897            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1898            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1899        else:
1900            wave = parmDict[hfx+'Lam']
1901    for phase in Histogram['Reflection Lists']:
1902        refDict = Histogram['Reflection Lists'][phase]
1903        if phase not in Phases:     #skips deleted or renamed phases silently!
1904            continue
1905        Phase = Phases[phase]
1906        pId = Phase['pId']
1907        pfx = '%d::'%(pId)
1908        phfx = '%d:%d:'%(pId,hId)
1909        hfx = ':%d:'%(hId)
1910        SGData = Phase['General']['SGData']
1911        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1912        im = 0
1913        if Phase['General']['Type'] in ['modulated','magnetic']:
1914            SSGData = Phase['General']['SSGData']
1915            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1916            im = 1  #offset in SS reflection list
1917            #??
1918        Dij = GetDij(phfx,SGData,parmDict)
1919        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
1920        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1921        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1922        Vst = np.sqrt(nl.det(G))    #V*
1923        if not Phase['General'].get('doPawley'):
1924            time0 = time.time()
1925            if im:
1926                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
1927            else:
1928                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1929#            print 'sf calc time: %.3fs'%(time.time()-time0)
1930        time0 = time.time()
1931        badPeak = False
1932        for iref,refl in enumerate(refDict['RefList']):
1933            if 'C' in calcControls[hfx+'histType']:
1934                if im:
1935                    h,k,l,m = refl[:4]
1936                else:
1937                    h,k,l = refl[:3]
1938                Uniq = np.inner(refl[:3],SGMT)
1939                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
1940                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
1941#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1942                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1943                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1944                refl[11+im] *= Vst*Lorenz
1945                 
1946                if Phase['General'].get('doPawley'):
1947                    try:
1948                        if im:
1949                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
1950                        else:
1951                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1952                        refl[9+im] = parmDict[pInd]
1953                    except KeyError:
1954#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1955                        continue
1956                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1957                iBeg = np.searchsorted(x,refl[5+im]-fmin)
1958                iFin = np.searchsorted(x,refl[5+im]+fmax)
1959                if not iBeg+iFin:       #peak below low limit - skip peak
1960                    continue
1961                elif not iBeg-iFin:     #peak above high limit - done
1962                    break
1963                elif iBeg > iFin:   #bad peak coeff - skip
1964                    badPeak = True
1965                    continue
1966                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
1967                if Ka2:
1968                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1969                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1970                    iBeg = np.searchsorted(x,pos2-fmin)
1971                    iFin = np.searchsorted(x,pos2+fmax)
1972                    if not iBeg+iFin:       #peak below low limit - skip peak
1973                        continue
1974                    elif not iBeg-iFin:     #peak above high limit - done
1975                        return yc,yb
1976                    elif iBeg > iFin:   #bad peak coeff - skip
1977                        continue
1978                    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
1979            elif 'T' in calcControls[hfx+'histType']:
1980                h,k,l = refl[:3]
1981                Uniq = np.inner(refl[:3],SGMT)
1982                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
1983                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
1984#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1985                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1986                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
1987                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1988                refl[11+im] *= Vst*Lorenz
1989                if Phase['General'].get('doPawley'):
1990                    try:
1991                        if im:
1992                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
1993                        else:
1994                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1995                        refl[9+im] = parmDict[pInd]
1996                    except KeyError:
1997#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1998                        continue
1999                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2000                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2001                iFin = np.searchsorted(x,refl[5+im]+fmax)
2002                if not iBeg+iFin:       #peak below low limit - skip peak
2003                    continue
2004                elif not iBeg-iFin:     #peak above high limit - done
2005                    break
2006                elif iBeg > iFin:   #bad peak coeff - skip
2007                    badPeak = True
2008                    continue
2009                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]
2010#        print 'profile calc time: %.3fs'%(time.time()-time0)
2011    if badPeak:
2012        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2013    return yc,yb
2014   
2015def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2016    'Needs a doc string'
2017   
2018    def cellVaryDerv(pfx,SGData,dpdA): 
2019        if SGData['SGLaue'] in ['-1',]:
2020            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2021                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2022        elif SGData['SGLaue'] in ['2/m',]:
2023            if SGData['SGUniq'] == 'a':
2024                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2025            elif SGData['SGUniq'] == 'b':
2026                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2027            else:
2028                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2029        elif SGData['SGLaue'] in ['mmm',]:
2030            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2031        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2032            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2033        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2034            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2035        elif SGData['SGLaue'] in ['3R', '3mR']:
2036            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2037        elif SGData['SGLaue'] in ['m3m','m3']:
2038            return [[pfx+'A0',dpdA[0]]]
2039           
2040    # create a list of dependent variables and set up a dictionary to hold their derivatives
2041    dependentVars = G2mv.GetDependentVars()
2042    depDerivDict = {}
2043    for j in dependentVars:
2044        depDerivDict[j] = np.zeros(shape=(len(x)))
2045    #print 'dependent vars',dependentVars
2046    lenX = len(x)               
2047    hId = Histogram['hId']
2048    hfx = ':%d:'%(hId)
2049    bakType = calcControls[hfx+'bakType']
2050    dMdv = np.zeros(shape=(len(varylist),len(x)))
2051    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2052    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2053        bBpos =varylist.index(hfx+'Back;0')
2054        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2055    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2056    for name in varylist:
2057        if 'Debye' in name:
2058            id = int(name.split(';')[-1])
2059            parm = name[:int(name.rindex(';'))]
2060            ip = names.index(parm)
2061            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2062    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2063    for name in varylist:
2064        if 'BkPk' in name:
2065            parm,id = name.split(';')
2066            id = int(id)
2067            if parm in names:
2068                ip = names.index(parm)
2069                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2070    cw = np.diff(x)
2071    cw = np.append(cw,cw[-1])
2072    Ka2 = False #also for TOF!
2073    if 'C' in calcControls[hfx+'histType']:   
2074        shl = max(parmDict[hfx+'SH/L'],0.002)
2075        if hfx+'Lam1' in parmDict.keys():
2076            wave = parmDict[hfx+'Lam1']
2077            Ka2 = True
2078            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2079            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2080        else:
2081            wave = parmDict[hfx+'Lam']
2082    for phase in Histogram['Reflection Lists']:
2083        refDict = Histogram['Reflection Lists'][phase]
2084        if phase not in Phases:     #skips deleted or renamed phases silently!
2085            continue
2086        Phase = Phases[phase]
2087        SGData = Phase['General']['SGData']
2088        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2089        im = 0
2090        if Phase['General']['Type'] in ['modulated','magnetic']:
2091            SSGData = Phase['General']['SSGData']
2092            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2093            im = 1  #offset in SS reflection list
2094            #??
2095        pId = Phase['pId']
2096        pfx = '%d::'%(pId)
2097        phfx = '%d:%d:'%(pId,hId)
2098        Dij = GetDij(phfx,SGData,parmDict)
2099        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2100        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2101        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2102        if not Phase['General'].get('doPawley'):
2103            time0 = time.time()
2104            if im:
2105                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2106            else:
2107                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2108#            print 'sf-derv time %.3fs'%(time.time()-time0)
2109            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2110        time0 = time.time()
2111        for iref,refl in enumerate(refDict['RefList']):
2112            if im:
2113                h,k,l,m = refl[:4]
2114            else:
2115                h,k,l = refl[:3]
2116            Uniq = np.inner(refl[:3],SGMT)
2117            if 'T' in calcControls[hfx+'histType']:
2118                wave = refl[14+im]
2119            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2120            if 'C' in calcControls[hfx+'histType']:        #CW powder
2121                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2122            else: #'T'OF
2123                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2124            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2125            iFin = np.searchsorted(x,refl[5+im]+fmax)
2126            if not iBeg+iFin:       #peak below low limit - skip peak
2127                continue
2128            elif not iBeg-iFin:     #peak above high limit - done
2129                break
2130            pos = refl[5+im]
2131            if 'C' in calcControls[hfx+'histType']:
2132                tanth = tand(pos/2.0)
2133                costh = cosd(pos/2.0)
2134                lenBF = iFin-iBeg
2135                dMdpk = np.zeros(shape=(6,lenBF))
2136                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2137                for i in range(5):
2138                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2139                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2140                if Ka2:
2141                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2142                    iBeg2 = np.searchsorted(x,pos2-fmin)
2143                    iFin2 = np.searchsorted(x,pos2+fmax)
2144                    if iBeg2-iFin2:
2145                        lenBF2 = iFin2-iBeg2
2146                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2147                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2148                        for i in range(5):
2149                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2150                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2151                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2152            else:   #'T'OF
2153                lenBF = iFin-iBeg
2154                if lenBF < 0:   #bad peak coeff
2155                    break
2156                dMdpk = np.zeros(shape=(6,lenBF))
2157                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2158                for i in range(6):
2159                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2160                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2161            if Phase['General'].get('doPawley'):
2162                dMdpw = np.zeros(len(x))
2163                try:
2164                    if im:
2165                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2166                    else:
2167                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2168                    idx = varylist.index(pIdx)
2169                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2170                    if Ka2: #not for TOF either
2171                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2172                    dMdv[idx] = dMdpw
2173                except: # ValueError:
2174                    pass
2175            if 'C' in calcControls[hfx+'histType']:
2176                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2177                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2178                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2179                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2180                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2181                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2182                    hfx+'DisplaceY':[dpdY,'pos'],}
2183                if 'Bragg' in calcControls[hfx+'instType']:
2184                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2185                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2186                else:
2187                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2188            else:   #'T'OF
2189                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2190                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2191                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2192                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2193                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2194                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2195                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2196                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2197            for name in names:
2198                item = names[name]
2199                if name in varylist:
2200                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2201                    if Ka2:
2202                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2203                elif name in dependentVars:
2204                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2205                    if Ka2:
2206                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2207            for iPO in dIdPO:
2208                if iPO in varylist:
2209                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2210                    if Ka2:
2211                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2212                elif iPO in dependentVars:
2213                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2214                    if Ka2:
2215                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2216            for i,name in enumerate(['omega','chi','phi']):
2217                aname = pfx+'SH '+name
2218                if aname in varylist:
2219                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2220                    if Ka2:
2221                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2222                elif aname in dependentVars:
2223                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2224                    if Ka2:
2225                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2226            for iSH in dFdODF:
2227                if iSH in varylist:
2228                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2229                    if Ka2:
2230                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2231                elif iSH in dependentVars:
2232                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2233                    if Ka2:
2234                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2235            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2236            for name,dpdA in cellDervNames:
2237                if name in varylist:
2238                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2239                    if Ka2:
2240                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2241                elif name in dependentVars:
2242                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2243                    if Ka2:
2244                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2245            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2246            for name in dDijDict:
2247                if name in varylist:
2248                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2249                    if Ka2:
2250                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2251                elif name in dependentVars:
2252                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2253                    if Ka2:
2254                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2255            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2256                if name in varylist:
2257                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2258                    if Ka2:
2259                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2260                elif name in dependentVars:
2261                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2262                    if Ka2:
2263                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2264            if 'C' in calcControls[hfx+'histType']:
2265                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2266            else:   #'T'OF
2267                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2268            for name in gamDict:
2269                if name in varylist:
2270                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2271                    if Ka2:
2272                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2273                elif name in dependentVars:
2274                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2275                    if Ka2:
2276                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2277            for name in sigDict:
2278                if name in varylist:
2279                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2280                    if Ka2:
2281                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2282                elif name in dependentVars:
2283                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2284                    if Ka2:
2285                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2286            for name in ['BabA','BabU']:
2287                if refl[9+im]:
2288                    if phfx+name in varylist:
2289                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2290                        if Ka2:
2291                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2292                    elif phfx+name in dependentVars:                   
2293                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2294                        if Ka2:
2295                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2296            if not Phase['General'].get('doPawley'):
2297                #do atom derivatives -  for RB,F,X & U so far             
2298                corr = dervDict['int']/refl[9+im]
2299                if Ka2:
2300                    corr2 = dervDict2['int']/refl[9+im]
2301                for name in varylist+dependentVars:
2302                    if '::RBV;' in name:
2303                        pass
2304                    else:
2305                        try:
2306                            aname = name.split(pfx)[1][:2]
2307                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2308                        except IndexError:
2309                            continue
2310                    if name in varylist:
2311                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2312                        if Ka2:
2313                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2314                    elif name in dependentVars:
2315                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2316                        if Ka2:
2317                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2318    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2319    # now process derivatives in constraints
2320    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2321    return dMdv
2322   
2323def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2324    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2325    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2326    are returned to either dervRefine or HessRefine.
2327
2328    :returns:
2329    '''
2330    nobs = Histogram['Residuals']['Nobs']
2331    hId = Histogram['hId']
2332    hfx = ':%d:'%(hId)
2333    pfx = '%d::'%(Phase['pId'])
2334    phfx = '%d:%d:'%(Phase['pId'],hId)
2335    SGData = Phase['General']['SGData']
2336    im = 0
2337    if Phase['General']['Type'] in ['modulated','magnetic']:
2338        SSGData = Phase['General']['SSGData']
2339        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2340        im = 1  #offset in SS reflection list
2341        #??
2342    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2343    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2344    refDict = Histogram['Data']
2345    if im:
2346        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2347    else:
2348        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2349    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2350    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2351    dependentVars = G2mv.GetDependentVars()
2352    depDerivDict = {}
2353    for j in dependentVars:
2354        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2355    wdf = np.zeros(len(refDict['RefList']))
2356    if calcControls['F**2']:
2357        for iref,ref in enumerate(refDict['RefList']):
2358            if ref[6+im] > 0:
2359                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1] 
2360                w = 1.0/ref[6+im]
2361                if w*ref[5+im] >= calcControls['minF/sig']:
2362                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2363                    for j,var in enumerate(varylist):
2364                        if var in dFdvDict:
2365                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2366                    for var in dependentVars:
2367                        if var in dFdvDict:
2368                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2369                    if phfx+'Scale' in varylist:
2370                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2371                    elif phfx+'Scale' in dependentVars:
2372                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]
2373                    for item in ['Ep','Es','Eg']:
2374                        if phfx+item in varylist and phfx+item in dervDict:
2375                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2376                        elif phfx+item in dependentVars and phfx+item in dervDict:
2377                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2378                    for item in ['BabA','BabU']:
2379                        if phfx+item in varylist:
2380                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2381                        elif phfx+item in dependentVars:
2382                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2383    else:
2384        for iref,ref in enumerate(refDict['RefList']):
2385            if ref[5+im] > 0.:
2386                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2387                Fo = np.sqrt(ref[5+im])
2388                Fc = np.sqrt(ref[7+im])
2389                w = 1.0/ref[6+im]
2390                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2391                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2392                    for j,var in enumerate(varylist):
2393                        if var in dFdvDict:
2394                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2395                    for var in dependentVars:
2396                        if var in dFdvDict:
2397                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2398                    if phfx+'Scale' in varylist:
2399                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
2400                    elif phfx+'Scale' in dependentVars:
2401                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]                           
2402                    for item in ['Ep','Es','Eg']:
2403                        if phfx+item in varylist and phfx+item in dervDict:
2404                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #correct
2405                        elif phfx+item in dependentVars and phfx+item in dervDict:
2406                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2407                    for item in ['BabA','BabU']:
2408                        if phfx+item in varylist:
2409                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2410                        elif phfx+item in dependentVars:
2411                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2412    return dMdvh,depDerivDict,wdf
2413   
2414
2415def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2416    '''Loop over histograms and compute derivatives of the fitting
2417    model (M) with respect to all parameters.  Results are returned in
2418    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2419    n is the number of parameters and m is the number of data
2420    points. This can exceed memory when m gets large. This routine is
2421    used when refinement derivatives are selected as "analtytic
2422    Jacobian" in Controls.
2423
2424    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2425    '''
2426    parmDict.update(zip(varylist,values))
2427    G2mv.Dict2Map(parmDict,varylist)
2428    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2429    nvar = len(varylist)
2430    dMdv = np.empty(0)
2431    histoList = Histograms.keys()
2432    histoList.sort()
2433    for histogram in histoList:
2434        if 'PWDR' in histogram[:4]:
2435            Histogram = Histograms[histogram]
2436            hId = Histogram['hId']
2437            hfx = ':%d:'%(hId)
2438            wtFactor = calcControls[hfx+'wtFactor']
2439            Limits = calcControls[hfx+'Limits']
2440            x,y,w,yc,yb,yd = Histogram['Data']
2441            xB = np.searchsorted(x,Limits[0])
2442            xF = np.searchsorted(x,Limits[1])
2443            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2444                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2445        elif 'HKLF' in histogram[:4]:
2446            Histogram = Histograms[histogram]
2447            phase = Histogram['Reflection Lists']
2448            Phase = Phases[phase]
2449            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2450            hfx = ':%d:'%(Histogram['hId'])
2451            wtFactor = calcControls[hfx+'wtFactor']
2452            # now process derivatives in constraints
2453            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2454        else:
2455            continue        #skip non-histogram entries
2456        if len(dMdv):
2457            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2458        else:
2459            dMdv = np.sqrt(wtFactor)*dMdvh
2460           
2461    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2462    if np.any(pVals):
2463        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2464        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2465       
2466    return dMdv
2467
2468def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2469    '''Loop over histograms and compute derivatives of the fitting
2470    model (M) with respect to all parameters.  For each histogram, the
2471    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2472    number of parameters and m is the number of data points *in the
2473    histogram*. The (n by n) Hessian is computed from each Jacobian
2474    and it is returned.  This routine is used when refinement
2475    derivatives are selected as "analtytic Hessian" in Controls.
2476
2477    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2478    '''
2479    parmDict.update(zip(varylist,values))
2480    G2mv.Dict2Map(parmDict,varylist)
2481    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2482    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2483    nvar = len(varylist)
2484    Hess = np.empty(0)
2485    histoList = Histograms.keys()
2486    histoList.sort()
2487    for histogram in histoList:
2488        if 'PWDR' in histogram[:4]:
2489            Histogram = Histograms[histogram]
2490            hId = Histogram['hId']
2491            hfx = ':%d:'%(hId)
2492            wtFactor = calcControls[hfx+'wtFactor']
2493            Limits = calcControls[hfx+'Limits']
2494            x,y,w,yc,yb,yd = Histogram['Data']
2495            W = wtFactor*w
2496            dy = y-yc
2497            xB = np.searchsorted(x,Limits[0])
2498            xF = np.searchsorted(x,Limits[1])
2499            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2500                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2501            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2502            Dy = dy[xB:xF][np.newaxis,:]
2503            dMdvh *= Wt
2504            if dlg:
2505                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2506            if len(Hess):
2507                Hess += np.inner(dMdvh,dMdvh)
2508                dMdvh *= Wt*Dy
2509                Vec += np.sum(dMdvh,axis=1)
2510            else:
2511                Hess = np.inner(dMdvh,dMdvh)
2512                dMdvh *= Wt*Dy
2513                Vec = np.sum(dMdvh,axis=1)
2514        elif 'HKLF' in histogram[:4]:
2515            Histogram = Histograms[histogram]
2516            phase = Histogram['Reflection Lists']
2517            Phase = Phases[phase]
2518            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2519            hId = Histogram['hId']
2520            hfx = ':%d:'%(Histogram['hId'])
2521            wtFactor = calcControls[hfx+'wtFactor']
2522            # now process derivatives in constraints
2523            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2524#            print 'matrix build time: %.3f'%(time.time()-time0)
2525
2526            if dlg:
2527                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2528            if len(Hess):
2529                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2530                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2531            else:
2532                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2533                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2534        else:
2535            continue        #skip non-histogram entries
2536    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2537    if np.any(pVals):
2538        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2539        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2540        Hess += np.inner(dpdv*pWt,dpdv)
2541    return Vec,Hess
2542
2543def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2544    '''Computes the point-by-point discrepancies between every data point in every histogram
2545    and the observed value
2546   
2547    :returns: an np array of differences between observed and computed diffraction values.
2548    '''
2549    Values2Dict(parmDict, varylist, values)
2550    G2mv.Dict2Map(parmDict,varylist)
2551    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2552    M = np.empty(0)
2553    SumwYo = 0
2554    Nobs = 0
2555    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2556    histoList = Histograms.keys()
2557    histoList.sort()
2558    for histogram in histoList:
2559        if 'PWDR' in histogram[:4]:
2560            Histogram = Histograms[histogram]
2561            hId = Histogram['hId']
2562            hfx = ':%d:'%(hId)
2563            wtFactor = calcControls[hfx+'wtFactor']
2564            Limits = calcControls[hfx+'Limits']
2565            x,y,w,yc,yb,yd = Histogram['Data']
2566            yc *= 0.0                           #zero full calcd profiles
2567            yb *= 0.0
2568            yd *= 0.0
2569            xB = np.searchsorted(x,Limits[0])
2570            xF = np.searchsorted(x,Limits[1])
2571            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2572                varylist,Histogram,Phases,calcControls,pawleyLookup)
2573            yc[xB:xF] += yb[xB:xF]
2574            if not np.any(y):                   #fill dummy data
2575                rv = st.poisson(yc[xB:xF])
2576                y[xB:xF] = rv.rvs()
2577                Z = np.ones_like(yc[xB:xF])
2578                Z[1::2] *= -1
2579                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2580                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2581            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2582            W = wtFactor*w
2583            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2584            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2585            Nobs += Histogram['Residuals']['Nobs']
2586            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2587            SumwYo += Histogram['Residuals']['sumwYo']
2588            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2589            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2590            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2591            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2592            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.))
2593            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.))
2594            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2595            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2596            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2597            if dlg:
2598                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2599            M = np.concatenate((M,wdy))
2600#end of PWDR processing
2601        elif 'HKLF' in histogram[:4]:
2602            Histogram = Histograms[histogram]
2603            Histogram['Residuals'] = {}
2604            phase = Histogram['Reflection Lists']
2605            Phase = Phases[phase]
2606            hId = Histogram['hId']
2607            hfx = ':%d:'%(hId)
2608            wtFactor = calcControls[hfx+'wtFactor']
2609            pfx = '%d::'%(Phase['pId'])
2610            phfx = '%d:%d:'%(Phase['pId'],hId)
2611            SGData = Phase['General']['SGData']
2612            im = 0
2613            if Phase['General']['Type'] in ['modulated','magnetic']:
2614                SSGData = Phase['General']['SSGData']
2615                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2616                im = 1  #offset in SS reflection list
2617                #??
2618            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2619            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2620            refDict = Histogram['Data']
2621            time0 = time.time()
2622            if im:
2623                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2624            else:
2625                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2626#            print 'sf-calc time: %.3f'%(time.time()-time0)
2627            df = np.zeros(len(refDict['RefList']))
2628            sumwYo = 0
2629            sumFo = 0
2630            sumFo2 = 0
2631            sumdF = 0
2632            sumdF2 = 0
2633            nobs = 0
2634            if calcControls['F**2']:
2635                for i,ref in enumerate(refDict['RefList']):
2636                    if ref[6+im] > 0:
2637                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2638                        w = 1.0/ref[6+im]
2639                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
2640                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2641                        if w*ref[5+im] >= calcControls['minF/sig']:
2642                            Fo = np.sqrt(ref[5+im])
2643                            sumFo += Fo
2644                            sumFo2 += ref[5+im]
2645                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2646                            sumdF2 += abs(ref[5+im]-ref[7+im])
2647                            nobs += 1
2648                            df[i] = -w*(ref[5+im]-ref[7+im])
2649                            sumwYo += (w*ref[5+im])**2
2650            else:
2651                for i,ref in enumerate(refDict['RefList']):
2652                    if ref[5+im] > 0.:
2653                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2654                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
2655                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2656                        Fo = np.sqrt(ref[5+im])
2657                        Fc = np.sqrt(ref[7+im])
2658                        w = 2.0*Fo/ref[6+im]
2659                        if w*Fo >= calcControls['minF/sig']:
2660                            sumFo += Fo
2661                            sumFo2 += ref[5+im]
2662                            sumdF += abs(Fo-Fc)
2663                            sumdF2 += abs(ref[5+im]-ref[7+im])
2664                            nobs += 1
2665                            df[i] = -w*(Fo-Fc)
2666                            sumwYo += (w*Fo)**2
2667            Histogram['Residuals']['Nobs'] = nobs
2668            Histogram['Residuals']['sumwYo'] = sumwYo
2669            SumwYo += sumwYo
2670            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2671            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2672            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2673            Histogram['Residuals'][phfx+'Nref'] = nobs
2674            Nobs += nobs
2675            if dlg:
2676                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2677            M = np.concatenate((M,wtFactor*df))
2678# end of HKLF processing
2679    Histograms['sumwYo'] = SumwYo
2680    Histograms['Nobs'] = Nobs
2681    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2682    if dlg:
2683        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2684        if not GoOn:
2685            parmDict['saved values'] = values
2686            dlg.Destroy()
2687            raise Exception         #Abort!!
2688    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2689    if len(pVals):
2690        pSum = np.sum(pWt*pVals**2)
2691        for name in pWsum:
2692            if pWsum:
2693                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2694        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2695        Nobs += len(pVals)
2696        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2697    return M
2698                       
Note: See TracBrowser for help on using the repository browser.