source: trunk/GSASIIstrMath.py @ 1460

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

fix missing radius Factor for MC/SA plots
add Prfo, Abs & xt to reflection list items
fix most TOF derivatives (remove a multiply by 100!)

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