source: trunk/GSASIIstrMath.py @ 1880

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

fix wx bug in SetMapPeaksText? - 'GridWindow?' vs 'grid window'
add characteristic source choice to single crystal instrument parameters -
select mean Ka = (2*Ka1+Ka2)/3
change default single crystal wavelength to MoKa?

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