source: trunk/GSASIIstrMath.py @ 1622

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

add 4D charge flipping - works to give 3D structure; only enabled for 4D data
fix problems with hkl limits in Fourier calcs. - see adjHKLmax
fix reflection table/2D & 3D display issues
disable MC/SA for 4D data

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