source: trunk/GSASIIstrMath.py @ 1910

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

fix reading run-on numbers in Shelx HKLF 4 files
implement all I/O for nonmerohedral twin data from HKLF 5 files
Twin GUI stuff set for mero & nonmero twins

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