source: trunk/GSASIIstrMath.py @ 1913

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

fix sfact importers so twin id=1 by default
GUI details for twins fixed

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