source: trunk/GSASIIstrMath.py @ 1890

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

fix Flack derivative & other adjustments to single crystal derivatives

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