source: trunk/GSASIIstrMath.py @ 1902

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

suppress editing of 1st twin element; it is always the identity matrix and frac is 1-Sum(other twin fr)
make FWHM consistent with measured FWHM for all lists & exports of peak lists.

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