source: trunk/GSASIIstrMath.py @ 1490

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

fix TOF t-bar reference
"fix" TOF single xtal extinction?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 115.9 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2014-09-09 21:33:26 +0000 (Tue, 09 Sep 2014) $
8# $Author: vondreele $
9# $Revision: 1490 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1490 2014-09-09 21:33:26Z 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: 1490 $")
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 for al2O3!   
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; returns extinction & derivative
802    '''
803    extCor = 1.0
804    dervCor = 1.0
805    dervDict = {}
806    if calcControls[phfx+'EType'] != 'None':
807        SQ = 1/(4.*ref[4]**2)
808        if 'C' in parmDict[hfx+'Type']:           
809            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
810        else:   #'T'
811            cos2T = 1.0-2.*SQ*ref[12]**2                       #cos(2theta)           
812        if 'SXC' in parmDict[hfx+'Type']:
813            AV = 7.9406e5/parmDict[pfx+'Vol']**2
814            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
815            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
816            PLZ = AV*P12*ref[7]*parmDict[hfx+'Lam']**2
817        elif 'SNT' in parmDict[hfx+'Type']:
818            AV = 1.e7/parmDict[pfx+'Vol']**2
819            PL = SQ
820            PLZ = AV*ref[7]*ref[12]**2
821        elif 'SNC' in parmDict[hfx+'Type']:
822            AV = 1.e7/parmDict[pfx+'Vol']**2
823            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
824            PLZ = AV*ref[9]*parmDict[hfx+'Lam']**2      #Fcsq as per GSAS, why not FcTsq (ref[9])?
825           
826        if 'Primary' in calcControls[phfx+'EType']:
827            PLZ *= 1.5
828        else:
829            if 'C' in parmDict[hfx+'Type']:
830                PLZ *= calcControls[phfx+'Tbar']
831            else: #'T'
832                PLZ *= ref[13]      #t-bar
833        if 'Primary' in calcControls[phfx+'EType']:
834            PLZ *= 1.5
835            PSIG = parmDict[phfx+'Ep']
836        elif 'I & II' in calcControls[phfx+'EType']:
837            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
838        elif 'Type II' in calcControls[phfx+'EType']:
839            PSIG = parmDict[phfx+'Es']
840        else:       # 'Secondary Type I'
841            PSIG = parmDict[phfx+'Eg']/PL
842           
843        AG = 0.58+0.48*cos2T+0.24*cos2T**2
844        AL = 0.025+0.285*cos2T
845        BG = 0.02-0.025*cos2T
846        BL = 0.15-0.2*(0.75-cos2T)**2
847        if cos2T < 0.:
848            BL = -0.45*cos2T
849        CG = 2.
850        CL = 2.
851        PF = PLZ*PSIG
852       
853        if 'Gaussian' in calcControls[phfx+'EApprox']:
854            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
855            extCor = np.sqrt(PF4)
856            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
857        else:
858            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
859            extCor = np.sqrt(PF4)
860            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
861
862        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
863        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
864            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
865        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
866            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
867        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
868            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
869               
870    return 1./extCor,dervCor,dervDict
871   
872def Dict2Values(parmdict, varylist):
873    '''Use before call to leastsq to setup list of values for the parameters
874    in parmdict, as selected by key in varylist'''
875    return [parmdict[key] for key in varylist] 
876   
877def Values2Dict(parmdict, varylist, values):
878    ''' Use after call to leastsq to update the parameter dictionary with
879    values corresponding to keys in varylist'''
880    parmdict.update(zip(varylist,values))
881   
882def GetNewCellParms(parmDict,varyList):
883    'Needs a doc string'
884    newCellDict = {}
885    Anames = ['A'+str(i) for i in range(6)]
886    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
887    for item in varyList:
888        keys = item.split(':')
889        if keys[2] in Ddict:
890            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
891            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
892            newCellDict[parm] = [key,parmDict[key]-parmDict[item]]
893    return newCellDict          # is e.g. {'0::D11':A0-D11}
894   
895def ApplyXYZshifts(parmDict,varyList):
896    '''
897    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
898   
899    :param dict parmDict: parameter dictionary
900    :param list varyList: list of variables (not used!)
901    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
902
903    '''
904    newAtomDict = {}
905    for item in parmDict:
906        if 'dA' in item:
907            parm = ''.join(item.split('d'))
908            parmDict[parm] += parmDict[item]
909            newAtomDict[item] = [parm,parmDict[parm]]
910    return newAtomDict
911   
912def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
913    'Spherical harmonics texture'
914    IFCoup = 'Bragg' in calcControls[hfx+'instType']
915    if 'T' in calcControls[hfx+'histType']:
916        tth = parmDict[hfx+'2-theta']
917    else:
918        tth = refl[5]
919    odfCor = 1.0
920    H = refl[:3]
921    cell = G2lat.Gmat2cell(g)
922    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
923    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
924    phi,beta = G2lat.CrsAng(H,cell,SGData)
925    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
926    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
927    for item in SHnames:
928        L,M,N = eval(item.strip('C'))
929        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
930        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
931        Lnorm = G2lat.Lnorm(L)
932        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
933    return odfCor
934   
935def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
936    'Spherical harmonics texture derivatives'
937    if 'T' in calcControls[hfx+'histType']:
938        tth = parmDict[hfx+'2-theta']
939    else:
940        tth = refl[5]
941    FORPI = 4.0*np.pi
942    IFCoup = 'Bragg' in calcControls[hfx+'instType']
943    odfCor = 1.0
944    dFdODF = {}
945    dFdSA = [0,0,0]
946    H = refl[:3]
947    cell = G2lat.Gmat2cell(g)
948    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
949    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
950    phi,beta = G2lat.CrsAng(H,cell,SGData)
951    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
952    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
953    for item in SHnames:
954        L,M,N = eval(item.strip('C'))
955        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
956        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
957        Lnorm = G2lat.Lnorm(L)
958        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
959        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
960        for i in range(3):
961            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
962    return odfCor,dFdODF,dFdSA
963   
964def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
965    'spherical harmonics preferred orientation (cylindrical symmetry only)'
966    if 'T' in calcControls[hfx+'histType']:
967        tth = parmDict[hfx+'2-theta']
968    else:
969        tth = refl[5]
970    odfCor = 1.0
971    H = refl[:3]
972    cell = G2lat.Gmat2cell(g)
973    Sangl = [0.,0.,0.]
974    if 'Bragg' in calcControls[hfx+'instType']:
975        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
976        IFCoup = True
977    else:
978        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
979        IFCoup = False
980    phi,beta = G2lat.CrsAng(H,cell,SGData)
981    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
982    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
983    for item in SHnames:
984        L,N = eval(item.strip('C'))
985        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
986        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
987    return np.squeeze(odfCor)
988   
989def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
990    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
991    if 'T' in calcControls[hfx+'histType']:
992        tth = parmDict[hfx+'2-theta']
993    else:
994        tth = refl[5]
995    FORPI = 12.5663706143592
996    odfCor = 1.0
997    dFdODF = {}
998    H = refl[:3]
999    cell = G2lat.Gmat2cell(g)
1000    Sangl = [0.,0.,0.]
1001    if 'Bragg' in calcControls[hfx+'instType']:
1002        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1003        IFCoup = True
1004    else:
1005        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1006        IFCoup = False
1007    phi,beta = G2lat.CrsAng(H,cell,SGData)
1008    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1009    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1010    for item in SHnames:
1011        L,N = eval(item.strip('C'))
1012        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1013        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1014        dFdODF[phfx+item] = Kcsl*Lnorm
1015    return odfCor,dFdODF
1016   
1017def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1018    'March-Dollase preferred orientation correction'
1019    POcorr = 1.0
1020    MD = parmDict[phfx+'MD']
1021    if MD != 1.0:
1022        MDAxis = calcControls[phfx+'MDAxis']
1023        sumMD = 0
1024        for H in uniq:           
1025            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1026            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1027            sumMD += A**3
1028        POcorr = sumMD/len(uniq)
1029    return POcorr
1030   
1031def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1032    'Needs a doc string'
1033    POcorr = 1.0
1034    POderv = {}
1035    if calcControls[phfx+'poType'] == 'MD':
1036        MD = parmDict[phfx+'MD']
1037        MDAxis = calcControls[phfx+'MDAxis']
1038        sumMD = 0
1039        sumdMD = 0
1040        for H in uniq:           
1041            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1042            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1043            sumMD += A**3
1044            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1045        POcorr = sumMD/len(uniq)
1046        POderv[phfx+'MD'] = sumdMD/len(uniq)
1047    else:   #spherical harmonics
1048        if calcControls[phfx+'SHord']:
1049            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1050    return POcorr,POderv
1051   
1052def GetAbsorb(refl,hfx,calcControls,parmDict):
1053    'Needs a doc string'
1054    if 'Debye' in calcControls[hfx+'instType']:
1055        if 'T' in calcControls[hfx+'histType']:
1056            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1057        else:
1058            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1059    else:
1060        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1061   
1062def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1063    'Needs a doc string'
1064    if 'Debye' in calcControls[hfx+'instType']:
1065        if 'T' in calcControls[hfx+'histType']:
1066            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1067        else:
1068            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1069    else:
1070        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]))
1071       
1072def GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict):
1073    'Needs a doc string'
1074    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1075    pi2 = np.sqrt(2./np.pi)
1076    if 'T' in calcControls[hfx+'histType']:
1077        sth2 = sind(parmDict[hfx+'2-theta']/2.)**2
1078        wave = refl[14]
1079    else:   #'C'W
1080        sth2 = sind(refl[5]/2.)**2
1081        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1082    c2th = 1.-2.0*sth2
1083    flv2 = refl[9]*(wave/parmDict[pfx+'Vol'])**2
1084    if 'X' in calcControls[hfx+'histType']:
1085        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1086    xfac = flv2*parmDict[phfx+'Extinction']
1087    exb = 1.0
1088    if xfac > -1.:
1089        exb = 1./(1.+xfac)
1090    exl = 1.0
1091    if 0 < xfac <= 1.:
1092        xn = np.array([xfac**(i+1) for i in range(6)])
1093        exl = np.sum(xn*coef)
1094    elif xfac > 1.:
1095        xfac2 = 1./np.sqrt(xfac)
1096        exl = pi2*(1.-0.125/xfac)*xfac2
1097    return exb*sth2+exl*(1.-sth2)
1098   
1099def GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict):
1100    'Needs a doc string'
1101    delt = 0.001
1102    parmDict[phfx+'Extinction'] += delt
1103    plus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1104    parmDict[phfx+'Extinction'] -= 2.*delt
1105    minus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1106    parmDict[phfx+'Extinction'] += delt
1107    return (plus-minus)/(2.*delt)   
1108   
1109def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1110    'Needs a doc string'    #need powder extinction!
1111    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1112    if 'X' in parmDict[hfx+'Type']:
1113        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1114    POcorr = 1.0
1115    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture
1116        POcorr = SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1117    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1118        POcorr = GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1119    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1120        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1121    Icorr *= POcorr
1122    AbsCorr = 1.0
1123    AbsCorr = GetAbsorb(refl,hfx,calcControls,parmDict)
1124    Icorr *= AbsCorr
1125    ExtCorr = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1126    Icorr *= ExtCorr
1127    return Icorr,POcorr,AbsCorr,ExtCorr
1128   
1129def GetIntensityDerv(refl,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1130    'Needs a doc string'    #need powder extinction derivs!
1131    dIdsh = 1./parmDict[hfx+'Scale']
1132    dIdsp = 1./parmDict[phfx+'Scale']
1133    if 'X' in parmDict[hfx+'Type']:
1134        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1135        dIdPola /= pola
1136    else:       #'N'
1137        dIdPola = 0.0
1138    dFdODF = {}
1139    dFdSA = [0,0,0]
1140    dIdPO = {}
1141    if pfx+'SHorder' in parmDict:
1142        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1143        for iSH in dFdODF:
1144            dFdODF[iSH] /= odfCor
1145        for i in range(3):
1146            dFdSA[i] /= odfCor
1147    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1148        POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1149        for iPO in dIdPO:
1150            dIdPO[iPO] /= POcorr
1151    if 'T' in parmDict[hfx+'Type']:
1152        dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[16] #wave/abs corr
1153        dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[17]    #/ext corr
1154    else:
1155        dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[13] #wave/abs corr
1156        dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[14]    #/ext corr       
1157    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1158       
1159def GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1160    'Needs a doc string'
1161    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1162        costh = cosd(refl[5]/2.)
1163        #crystallite size
1164        if calcControls[phfx+'SizeType'] == 'isotropic':
1165            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1166        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1167            H = np.array(refl[:3])
1168            P = np.array(calcControls[phfx+'SizeAxis'])
1169            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1170            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1171            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1172        else:           #ellipsoidal crystallites
1173            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1174            H = np.array(refl[:3])
1175            lenR = G2pwd.ellipseSize(H,Sij,GB)
1176            Sgam = 1.8*wave/(np.pi*costh*lenR)
1177        #microstrain               
1178        if calcControls[phfx+'MustrainType'] == 'isotropic':
1179            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1180        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1181            H = np.array(refl[:3])
1182            P = np.array(calcControls[phfx+'MustrainAxis'])
1183            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1184            Si = parmDict[phfx+'Mustrain;i']
1185            Sa = parmDict[phfx+'Mustrain;a']
1186            Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1187        else:       #generalized - P.W. Stephens model
1188            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1189            Sum = 0
1190            for i,strm in enumerate(Strms):
1191                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1192            Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*np.sqrt(Sum)/np.pi
1193    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1194        #crystallite size
1195        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1196            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
1197        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1198            H = np.array(refl[:3])
1199            P = np.array(calcControls[phfx+'SizeAxis'])
1200            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1201            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1202            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1203        else:           #ellipsoidal crystallites   #OK
1204            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1205            H = np.array(refl[:3])
1206            lenR = G2pwd.ellipseSize(H,Sij,GB)
1207            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/lenR
1208        #microstrain               
1209        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1210            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
1211        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1212            H = np.array(refl[:3])
1213            P = np.array(calcControls[phfx+'MustrainAxis'])
1214            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1215            Si = parmDict[phfx+'Mustrain;i']
1216            Sa = parmDict[phfx+'Mustrain;a']
1217            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1218        else:       #generalized - P.W. Stephens model  OK
1219            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1220            Sum = 0
1221            for i,strm in enumerate(Strms):
1222                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1223            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4]**3
1224           
1225    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1226    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1227    sig /= ateln2
1228    return sig,gam
1229       
1230def GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1231    'Needs a doc string'
1232    gamDict = {}
1233    sigDict = {}
1234    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1235        costh = cosd(refl[5]/2.)
1236        tanth = tand(refl[5]/2.)
1237        #crystallite size derivatives
1238        if calcControls[phfx+'SizeType'] == 'isotropic':
1239            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1240            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1241            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1242        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1243            H = np.array(refl[:3])
1244            P = np.array(calcControls[phfx+'SizeAxis'])
1245            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1246            Si = parmDict[phfx+'Size;i']
1247            Sa = parmDict[phfx+'Size;a']
1248            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1249            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1250            Sgam = gami*sqtrm
1251            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1252            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1253            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1254            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1255            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1256            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1257        else:           #ellipsoidal crystallites
1258            const = 1.8*wave/(np.pi*costh)
1259            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1260            H = np.array(refl[:3])
1261            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1262            Sgam = const/lenR
1263            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1264                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1265                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1266        gamDict[phfx+'Size;mx'] = Sgam
1267        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1268               
1269        #microstrain derivatives               
1270        if calcControls[phfx+'MustrainType'] == 'isotropic':
1271            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1272            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1273            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1274        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1275            H = np.array(refl[:3])
1276            P = np.array(calcControls[phfx+'MustrainAxis'])
1277            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1278            Si = parmDict[phfx+'Mustrain;i']
1279            Sa = parmDict[phfx+'Mustrain;a']
1280            gami = 0.018*Si*Sa*tanth/np.pi
1281            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1282            Mgam = gami/sqtrm
1283            dsi = -gami*Si*cosP**2/sqtrm**3
1284            dsa = -gami*Sa*sinP**2/sqtrm**3
1285            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1286            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1287            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1288            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1289        else:       #generalized - P.W. Stephens model
1290            const = 0.018*refl[4]**2*tanth/np.pi
1291            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1292            Sum = 0
1293            for i,strm in enumerate(Strms):
1294                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1295                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1296                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1297            Mgam = const*np.sqrt(Sum)
1298            for i in range(len(Strms)):
1299                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1300                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1301        gamDict[phfx+'Mustrain;mx'] = Mgam
1302        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1303    else:   #'T'OF - All checked & OK
1304        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1305            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
1306            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1307            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1308        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1309            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
1310            H = np.array(refl[:3])
1311            P = np.array(calcControls[phfx+'SizeAxis'])
1312            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1313            Si = parmDict[phfx+'Size;i']
1314            Sa = parmDict[phfx+'Size;a']
1315            gami = const/(Si*Sa)
1316            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1317            Sgam = gami*sqtrm
1318            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1319            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1320            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1321            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1322            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1323            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1324        else:           #OK  ellipsoidal crystallites
1325            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
1326            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1327            H = np.array(refl[:3])
1328            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1329            Sgam = const/lenR
1330            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1331                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1332                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1333        gamDict[phfx+'Size;mx'] = Sgam  #OK
1334        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1335               
1336        #microstrain derivatives               
1337        if calcControls[phfx+'MustrainType'] == 'isotropic':
1338            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
1339            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1340            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1341        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1342            H = np.array(refl[:3])
1343            P = np.array(calcControls[phfx+'MustrainAxis'])
1344            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1345            Si = parmDict[phfx+'Mustrain;i']
1346            Sa = parmDict[phfx+'Mustrain;a']
1347            gami = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa
1348            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1349            Mgam = gami/sqtrm
1350            dsi = -gami*Si*cosP**2/sqtrm**3
1351            dsa = -gami*Sa*sinP**2/sqtrm**3
1352            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1353            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1354            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1355            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1356        else:       #generalized - P.W. Stephens model OK
1357            pwrs = calcControls[phfx+'MuPwrs']
1358            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1359            const = 1.e-6*parmDict[hfx+'difC']*refl[4]**3
1360            Sum = 0
1361            for i,strm in enumerate(Strms):
1362                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1363                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1364                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1365            Mgam = const*np.sqrt(Sum)
1366            for i in range(len(Strms)):
1367                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1368                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1369        gamDict[phfx+'Mustrain;mx'] = Mgam
1370        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1371       
1372    return sigDict,gamDict
1373       
1374def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1375    'Needs a doc string'
1376    h,k,l = refl[:3]
1377    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1378    d = np.sqrt(dsq)
1379
1380    refl[4] = d
1381    if 'C' in calcControls[hfx+'histType']:
1382        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1383        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1384        if 'Bragg' in calcControls[hfx+'instType']:
1385            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1386                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1387        else:               #Debye-Scherrer - simple but maybe not right
1388            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1389    elif 'T' in calcControls[hfx+'histType']:
1390        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1391        #do I need sample position effects - maybe?
1392    return pos
1393
1394def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1395    'Needs a doc string'
1396    dpr = 180./np.pi
1397    h,k,l = refl[:3]
1398    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1399    dst = np.sqrt(dstsq)
1400    dsp = 1./dst
1401    if 'C' in calcControls[hfx+'histType']:
1402        pos = refl[5]-parmDict[hfx+'Zero']
1403        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1404        dpdw = const*dst
1405        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1406        dpdA *= const*wave/(2.0*dst)
1407        dpdZ = 1.0
1408        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1409        if 'Bragg' in calcControls[hfx+'instType']:
1410            dpdSh = -4.*const*cosd(pos/2.0)
1411            dpdTr = -const*sind(pos)*100.0
1412            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1413        else:               #Debye-Scherrer - simple but maybe not right
1414            dpdXd = -const*cosd(pos)
1415            dpdYd = -const*sind(pos)
1416            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1417    elif 'T' in calcControls[hfx+'histType']:
1418        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1419        dpdZ = 1.0
1420        dpdDC = dsp
1421        dpdDA = dsp**2
1422        dpdDB = 1./dsp
1423        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
1424           
1425def GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict):
1426    'Needs a doc string'
1427    laue = SGData['SGLaue']
1428    uniq = SGData['SGUniq']
1429    h,k,l = refl[:3]
1430    if laue in ['m3','m3m']:
1431        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1432            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1433    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1434        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1435    elif laue in ['3R','3mR']:
1436        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1437    elif laue in ['4/m','4/mmm']:
1438        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1439    elif laue in ['mmm']:
1440        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1441    elif laue in ['2/m']:
1442        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1443        if uniq == 'a':
1444            Dij += parmDict[phfx+'D23']*k*l
1445        elif uniq == 'b':
1446            Dij += parmDict[phfx+'D13']*h*l
1447        elif uniq == 'c':
1448            Dij += parmDict[phfx+'D12']*h*k
1449    else:
1450        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1451            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1452    if 'C' in calcControls[hfx+'histType']:
1453        return -180.*Dij*refl[4]**2*tand(refl[5]/2.0)/np.pi
1454    else:
1455        return -Dij*parmDict[hfx+'difC']*refl[4]**2/2.
1456           
1457def GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict):
1458    'Needs a doc string'
1459    laue = SGData['SGLaue']
1460    uniq = SGData['SGUniq']
1461    h,k,l = refl[:3]
1462    if laue in ['m3','m3m']:
1463        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1464            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1465    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1466        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1467    elif laue in ['3R','3mR']:
1468        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1469    elif laue in ['4/m','4/mmm']:
1470        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1471    elif laue in ['mmm']:
1472        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1473    elif laue in ['2/m']:
1474        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1475        if uniq == 'a':
1476            dDijDict[phfx+'D23'] = k*l
1477        elif uniq == 'b':
1478            dDijDict[phfx+'D13'] = h*l
1479        elif uniq == 'c':
1480            dDijDict[phfx+'D12'] = h*k
1481            names.append()
1482    else:
1483        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1484            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1485    if 'C' in calcControls[hfx+'histType']:
1486        for item in dDijDict:
1487            dDijDict[item] *= -180.0*refl[4]**2*tand(refl[5]/2.0)/np.pi
1488    else:
1489        for item in dDijDict:
1490            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4]**2/2.
1491    return dDijDict
1492   
1493def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1494    'Needs a doc string'
1495    histoList = Histograms.keys()
1496    histoList.sort()
1497    for histogram in histoList:
1498        if 'PWDR' in histogram[:4]:
1499            Histogram = Histograms[histogram]
1500            hId = Histogram['hId']
1501            hfx = ':%d:'%(hId)
1502            Limits = calcControls[hfx+'Limits']
1503            if 'C' in calcControls[hfx+'histType']:
1504                shl = max(parmDict[hfx+'SH/L'],0.0005)
1505                Ka2 = False
1506                kRatio = 0.0
1507                if hfx+'Lam1' in parmDict.keys():
1508                    Ka2 = True
1509                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1510                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1511            x,y,w,yc,yb,yd = Histogram['Data']
1512            xB = np.searchsorted(x,Limits[0])
1513            xF = np.searchsorted(x,Limits[1])
1514            ymb = np.array(y-yb)
1515            ymb = np.where(ymb,ymb,1.0)
1516            ycmb = np.array(yc-yb)
1517            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1518            refLists = Histogram['Reflection Lists']
1519            for phase in refLists:
1520                Phase = Phases[phase]
1521                pId = Phase['pId']
1522                phfx = '%d:%d:'%(pId,hId)
1523                refDict = refLists[phase]
1524                sumFo = 0.0
1525                sumdF = 0.0
1526                sumFosq = 0.0
1527                sumdFsq = 0.0
1528                for refl in refDict['RefList']:
1529                    if 'C' in calcControls[hfx+'histType']:
1530                        yp = np.zeros_like(yb)
1531                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1532                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1533                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1534                        iFin2 = iFin
1535                        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
1536                        if Ka2:
1537                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1538                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1539                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1540                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1541                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1542                                continue
1543                            elif not iBeg2-iFin2:     #peak above high limit - done
1544                                break
1545                            yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
1546                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0))
1547                    elif 'T' in calcControls[hfx+'histType']:
1548                        yp = np.zeros_like(yb)
1549                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1550                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1551                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1552                        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
1553                        refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
1554                    Fo = np.sqrt(np.abs(refl[8]))
1555                    Fc = np.sqrt(np.abs(refl[9]))
1556                    sumFo += Fo
1557                    sumFosq += refl[8]**2
1558                    sumdF += np.abs(Fo-Fc)
1559                    sumdFsq += (refl[8]-refl[9])**2
1560                Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1561                Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1562                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1563                Histogram['Residuals']['hId'] = hId
1564        elif 'HKLF' in histogram[:4]:
1565            Histogram = Histograms[histogram]
1566            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1567               
1568def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1569    'Needs a doc string'
1570   
1571    def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
1572        U = parmDict[hfx+'U']
1573        V = parmDict[hfx+'V']
1574        W = parmDict[hfx+'W']
1575        X = parmDict[hfx+'X']
1576        Y = parmDict[hfx+'Y']
1577        tanPos = tand(refl[5]/2.0)
1578        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1579        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1580        sig = max(0.001,sig)
1581        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1582        gam = max(0.001,gam)
1583        return sig,gam
1584               
1585    def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
1586        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4]**2+   \
1587            parmDict[hfx+'sig-2']*refl[4]**4+parmDict[hfx+'sig-q']/refl[4]**2
1588        gam = parmDict[hfx+'X']*refl[4]+parmDict[hfx+'Y']*refl[4]**2
1589        Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1590        sig += Ssig
1591        gam += Sgam
1592        return sig,gam
1593       
1594    def GetReflAlpBet(refl,hfx,parmDict):
1595        alp = parmDict[hfx+'alpha']/refl[4]
1596        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]**4+parmDict[hfx+'beta-q']/refl[4]**2
1597        return alp,bet
1598               
1599    hId = Histogram['hId']
1600    hfx = ':%d:'%(hId)
1601    bakType = calcControls[hfx+'bakType']
1602    yb = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1603    yc = np.zeros_like(yb)
1604    cw = np.diff(x)
1605    cw = np.append(cw,cw[-1])
1606       
1607    if 'C' in calcControls[hfx+'histType']:   
1608        shl = max(parmDict[hfx+'SH/L'],0.002)
1609        Ka2 = False
1610        if hfx+'Lam1' in parmDict.keys():
1611            wave = parmDict[hfx+'Lam1']
1612            Ka2 = True
1613            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1614            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1615        else:
1616            wave = parmDict[hfx+'Lam']
1617    for phase in Histogram['Reflection Lists']:
1618        refDict = Histogram['Reflection Lists'][phase]
1619        Phase = Phases[phase]
1620        pId = Phase['pId']
1621        pfx = '%d::'%(pId)
1622        phfx = '%d:%d:'%(pId,hId)
1623        hfx = ':%d:'%(hId)
1624        SGData = Phase['General']['SGData']
1625        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1626        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #Do I want to modify by Dij?
1627        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1628        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1629        Vst = np.sqrt(nl.det(G))    #V*
1630        if not Phase['General'].get('doPawley'):
1631            time0 = time.time()
1632            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1633#            print 'sf calc time: %.3fs'%(time.time()-time0)
1634        time0 = time.time()
1635        for iref,refl in enumerate(refDict['RefList']):
1636            if 'C' in calcControls[hfx+'histType']:
1637                h,k,l = refl[:3]
1638                Uniq = np.inner(refl[:3],SGMT)
1639                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1640                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1641                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1642                refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1643                refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1644                refl[11] *= Vst*Lorenz
1645                 
1646                if Phase['General'].get('doPawley'):
1647                    try:
1648                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1649                        refl[9] = parmDict[pInd]
1650                    except KeyError:
1651#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1652                        continue
1653                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1654                iBeg = np.searchsorted(x,refl[5]-fmin)
1655                iFin = np.searchsorted(x,refl[5]+fmax)
1656                if not iBeg+iFin:       #peak below low limit - skip peak
1657                    continue
1658                elif not iBeg-iFin:     #peak above high limit - done
1659                    break
1660                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
1661                if Ka2:
1662                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1663                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1664                    iBeg = np.searchsorted(x,pos2-fmin)
1665                    iFin = np.searchsorted(x,pos2+fmax)
1666                    if not iBeg+iFin:       #peak below low limit - skip peak
1667                        continue
1668                    elif not iBeg-iFin:     #peak above high limit - done
1669                        return yc,yb
1670                    yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1671            elif 'T' in calcControls[hfx+'histType']:
1672                h,k,l = refl[:3]
1673                Uniq = np.inner(refl[:3],SGMT)
1674                refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict)         #corrected reflection position
1675                Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4]**4                                                #TOF Lorentz correction
1676                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1677                refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1678                refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
1679                refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1680                refl[11] *= Vst*Lorenz
1681                if Phase['General'].get('doPawley'):
1682                    try:
1683                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1684                        refl[9] = parmDict[pInd]
1685                    except KeyError:
1686#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1687                        continue
1688                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1689                iBeg = np.searchsorted(x,refl[5]-fmin)
1690                iFin = np.searchsorted(x,refl[5]+fmax)
1691                if not iBeg+iFin:       #peak below low limit - skip peak
1692                    continue
1693                elif not iBeg-iFin:     #peak above high limit - done
1694                    break
1695                yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
1696#        print 'profile calc time: %.3fs'%(time.time()-time0)
1697    return yc,yb
1698   
1699def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1700    'Needs a doc string'
1701   
1702    def cellVaryDerv(pfx,SGData,dpdA): 
1703        if SGData['SGLaue'] in ['-1',]:
1704            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1705                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1706        elif SGData['SGLaue'] in ['2/m',]:
1707            if SGData['SGUniq'] == 'a':
1708                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1709            elif SGData['SGUniq'] == 'b':
1710                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1711            else:
1712                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1713        elif SGData['SGLaue'] in ['mmm',]:
1714            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1715        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1716            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1717        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1718            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1719        elif SGData['SGLaue'] in ['3R', '3mR']:
1720            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1721        elif SGData['SGLaue'] in ['m3m','m3']:
1722            return [[pfx+'A0',dpdA[0]]]
1723           
1724    # create a list of dependent variables and set up a dictionary to hold their derivatives
1725    dependentVars = G2mv.GetDependentVars()
1726    depDerivDict = {}
1727    for j in dependentVars:
1728        depDerivDict[j] = np.zeros(shape=(len(x)))
1729    #print 'dependent vars',dependentVars
1730    lenX = len(x)               
1731    hId = Histogram['hId']
1732    hfx = ':%d:'%(hId)
1733    bakType = calcControls[hfx+'bakType']
1734    dMdv = np.zeros(shape=(len(varylist),len(x)))
1735    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1736    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1737        bBpos =varylist.index(hfx+'Back:0')
1738        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1739    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1740    for name in varylist:
1741        if 'Debye' in name:
1742            id = int(name.split(':')[-1])
1743            parm = name[:int(name.rindex(':'))]
1744            ip = names.index(parm)
1745            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1746    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1747    for name in varylist:
1748        if 'BkPk' in name:
1749            parm,id = name.split(';')
1750            id = int(id)
1751            if parm in names:
1752                ip = names.index(parm)
1753                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1754    cw = np.diff(x)
1755    cw = np.append(cw,cw[-1])
1756    Ka2 = False #also for TOF!
1757    if 'C' in calcControls[hfx+'histType']:   
1758        shl = max(parmDict[hfx+'SH/L'],0.002)
1759        if hfx+'Lam1' in parmDict.keys():
1760            wave = parmDict[hfx+'Lam1']
1761            Ka2 = True
1762            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1763            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1764        else:
1765            wave = parmDict[hfx+'Lam']
1766    for phase in Histogram['Reflection Lists']:
1767        refDict = Histogram['Reflection Lists'][phase]
1768        Phase = Phases[phase]
1769        SGData = Phase['General']['SGData']
1770        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1771        pId = Phase['pId']
1772        pfx = '%d::'%(pId)
1773        phfx = '%d:%d:'%(pId,hId)
1774        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #And modify here by Dij? - no
1775        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1776        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1777        if not Phase['General'].get('doPawley'):
1778            time0 = time.time()
1779            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1780#            print 'sf-derv time %.3fs'%(time.time()-time0)
1781            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1782        time0 = time.time()
1783        for iref,refl in enumerate(refDict['RefList']):
1784            h,k,l = refl[:3]
1785            Uniq = np.inner(refl[:3],SGMT)
1786            if 'T' in calcControls[hfx+'histType']:
1787                wave = refl[14]
1788            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1789            if 'C' in calcControls[hfx+'histType']:        #CW powder
1790                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1791            else: #'T'OF
1792                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1793            iBeg = np.searchsorted(x,refl[5]-fmin)
1794            iFin = np.searchsorted(x,refl[5]+fmax)
1795            if not iBeg+iFin:       #peak below low limit - skip peak
1796                continue
1797            elif not iBeg-iFin:     #peak above high limit - done
1798                break
1799            pos = refl[5]
1800            if 'C' in calcControls[hfx+'histType']:
1801                tanth = tand(pos/2.0)
1802                costh = cosd(pos/2.0)
1803                lenBF = iFin-iBeg
1804                dMdpk = np.zeros(shape=(6,lenBF))
1805                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1806                for i in range(5):
1807                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
1808                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1809                if Ka2:
1810                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1811                    iBeg2 = np.searchsorted(x,pos2-fmin)
1812                    iFin2 = np.searchsorted(x,pos2+fmax)
1813                    if iBeg2-iFin2:
1814                        lenBF2 = iFin2-iBeg2
1815                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1816                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1817                        for i in range(5):
1818                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
1819                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
1820                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1821            else:   #'T'OF
1822                lenBF = iFin-iBeg
1823                dMdpk = np.zeros(shape=(6,lenBF))
1824                dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1825                for i in range(6):
1826                    dMdpk[i] += refl[11]*refl[9]*dMdipk[i]      #cw[iBeg:iFin]*
1827                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
1828            if Phase['General'].get('doPawley'):
1829                dMdpw = np.zeros(len(x))
1830                try:
1831                    pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1832                    idx = varylist.index(pIdx)
1833                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1834                    if Ka2: #not for TOF either
1835                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1836                    dMdv[idx] = dMdpw
1837                except: # ValueError:
1838                    pass
1839            if 'C' in calcControls[hfx+'histType']:
1840                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1841                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1842                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1843                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1844                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1845                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1846                    hfx+'DisplaceY':[dpdY,'pos'],}
1847                if 'Bragg' in calcControls[hfx+'instType']:
1848                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1849                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
1850                else:
1851                    names.update({hfx+'Absorption':[dFdAb,'int'],})
1852            else:   #'T'OF
1853                dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1854                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1855                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1856                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1857                    hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1858                    hfx+'beta-q':[1./refl[4]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'],
1859                    hfx+'sig-2':[refl[4]**4,'sig'],hfx+'sig-q':[1./refl[4]**2,'sig'],
1860                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1861            for name in names:
1862                item = names[name]
1863                if name in varylist:
1864                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1865                    if Ka2:
1866                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1867                elif name in dependentVars:
1868                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1869                    if Ka2:
1870                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1871            for iPO in dIdPO:
1872                if iPO in varylist:
1873                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1874                    if Ka2:
1875                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1876                elif iPO in dependentVars:
1877                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1878                    if Ka2:
1879                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1880            for i,name in enumerate(['omega','chi','phi']):
1881                aname = pfx+'SH '+name
1882                if aname in varylist:
1883                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1884                    if Ka2:
1885                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1886                elif aname in dependentVars:
1887                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1888                    if Ka2:
1889                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1890            for iSH in dFdODF:
1891                if iSH in varylist:
1892                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1893                    if Ka2:
1894                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1895                elif iSH in dependentVars:
1896                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1897                    if Ka2:
1898                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1899            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1900            for name,dpdA in cellDervNames:
1901                if name in varylist:
1902                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1903                    if Ka2:
1904                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1905                elif name in dependentVars:
1906                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1907                    if Ka2:
1908                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1909            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict)
1910            for name in dDijDict:
1911                if name in varylist:
1912                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1913                    if Ka2:
1914                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1915                elif name in dependentVars:
1916                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1917                    if Ka2:
1918                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1919            if 'C' in calcControls[hfx+'histType']:
1920                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1921            else:   #'T'OF
1922                sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1923            for name in gamDict:
1924                if name in varylist:
1925                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1926                    if Ka2:
1927                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1928                elif name in dependentVars:
1929                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1930                    if Ka2:
1931                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1932            for name in sigDict:
1933                if name in varylist:
1934                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1935                    if Ka2:
1936                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1937                elif name in dependentVars:
1938                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1939                    if Ka2:
1940                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1941            for name in ['BabA','BabU']:
1942                if refl[9]:
1943                    if phfx+name in varylist:
1944                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1945                        if Ka2:
1946                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1947                    elif phfx+name in dependentVars:                   
1948                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1949                        if Ka2:
1950                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
1951            if not Phase['General'].get('doPawley'):
1952                #do atom derivatives -  for RB,F,X & U so far             
1953                corr = dervDict['int']/refl[9]
1954                if Ka2:
1955                    corr2 = dervDict2['int']/refl[9]
1956                for name in varylist+dependentVars:
1957                    if '::RBV;' in name:
1958                        pass
1959                    else:
1960                        try:
1961                            aname = name.split(pfx)[1][:2]
1962                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1963                        except IndexError:
1964                            continue
1965                    if name in varylist:
1966                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1967                        if Ka2:
1968                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1969                    elif name in dependentVars:
1970                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1971                        if Ka2:
1972                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1973    #        print 'profile derv time: %.3fs'%(time.time()-time0)
1974    # now process derivatives in constraints
1975    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1976    return dMdv
1977   
1978def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1979    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1980    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
1981    are returned to either dervRefine or HessRefine.
1982
1983    :returns:
1984    '''
1985    nobs = Histogram['Residuals']['Nobs']
1986    hId = Histogram['hId']
1987    hfx = ':%d:'%(hId)
1988    pfx = '%d::'%(Phase['pId'])
1989    phfx = '%d:%d:'%(Phase['pId'],hId)
1990    SGData = Phase['General']['SGData']
1991    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1992    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1993    refDict = Histogram['Data']
1994    dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1995    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1996    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1997    dependentVars = G2mv.GetDependentVars()
1998    depDerivDict = {}
1999    for j in dependentVars:
2000        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2001    wdf = np.zeros(len(refDict['RefList']))
2002    if calcControls['F**2']:
2003        for iref,ref in enumerate(refDict['RefList']):
2004            if ref[6] > 0:
2005                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:] 
2006                w = 1.0/ref[6]
2007                if w*ref[5] >= calcControls['minF/sig']:
2008                    wdf[iref] = w*(ref[5]-ref[7])
2009                    for j,var in enumerate(varylist):
2010                        if var in dFdvDict:
2011                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]   #*dervCor
2012                    for var in dependentVars:
2013                        if var in dFdvDict:
2014                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]  #*dervCor
2015                    if phfx+'Scale' in varylist:
2016                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]  #OK? *dervCor
2017                    elif phfx+'Scale' in dependentVars:
2018                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11]  #OK? *dervCor
2019                    for item in ['Ep','Es','Eg']:
2020                        if phfx+item in varylist and dervDict:
2021                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #OK? /dervCor
2022                        elif phfx+item in dependentVars and dervDict:
2023                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]  #OK? /dervCor
2024                    for item in ['BabA','BabU']:
2025                        if phfx+item in varylist:
2026                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2027                        elif phfx+item in dependentVars:
2028                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2029    else:
2030        for iref,ref in enumerate(refDict['RefList']):
2031            if ref[5] > 0.:
2032                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:]
2033                Fo = np.sqrt(ref[5])
2034                Fc = np.sqrt(ref[7])
2035                w = 1.0/ref[6]
2036                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2037                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2038                    for j,var in enumerate(varylist):
2039                        if var in dFdvDict:
2040                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor       #*ref[11]
2041                    for var in dependentVars:
2042                        if var in dFdvDict:
2043                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor      #*ref[11]
2044                    if phfx+'Scale' in varylist:
2045                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]    #*dervCor
2046                    elif phfx+'Scale' in dependentVars:
2047                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11] #*dervCor                           
2048                    for item in ['Ep','Es','Eg']:
2049                        if phfx+item in varylist and dervDict:
2050                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #correct
2051                        elif phfx+item in dependentVars and dervDict:
2052                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
2053                    for item in ['BabA','BabU']:
2054                        if phfx+item in varylist:
2055                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2056                        elif phfx+item in dependentVars:
2057                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2058    return dMdvh,depDerivDict,wdf
2059   
2060
2061def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2062    '''Loop over histograms and compute derivatives of the fitting
2063    model (M) with respect to all parameters.  Results are returned in
2064    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2065    n is the number of parameters and m is the number of data
2066    points. This can exceed memory when m gets large. This routine is
2067    used when refinement derivatives are selected as "analtytic
2068    Jacobian" in Controls.
2069
2070    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2071    '''
2072    parmDict.update(zip(varylist,values))
2073    G2mv.Dict2Map(parmDict,varylist)
2074    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2075    nvar = len(varylist)
2076    dMdv = np.empty(0)
2077    histoList = Histograms.keys()
2078    histoList.sort()
2079    for histogram in histoList:
2080        if 'PWDR' in histogram[:4]:
2081            Histogram = Histograms[histogram]
2082            hId = Histogram['hId']
2083            hfx = ':%d:'%(hId)
2084            wtFactor = calcControls[hfx+'wtFactor']
2085            Limits = calcControls[hfx+'Limits']
2086            x,y,w,yc,yb,yd = Histogram['Data']
2087            xB = np.searchsorted(x,Limits[0])
2088            xF = np.searchsorted(x,Limits[1])
2089            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2090                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2091        elif 'HKLF' in histogram[:4]:
2092            Histogram = Histograms[histogram]
2093            phase = Histogram['Reflection Lists']
2094            Phase = Phases[phase]
2095            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2096            hfx = ':%d:'%(Histogram['hId'])
2097            wtFactor = calcControls[hfx+'wtFactor']
2098            # now process derivatives in constraints
2099            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2100        else:
2101            continue        #skip non-histogram entries
2102        if len(dMdv):
2103            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2104        else:
2105            dMdv = np.sqrt(wtFactor)*dMdvh
2106           
2107    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2108    if np.any(pVals):
2109        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2110        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2111       
2112    return dMdv
2113
2114def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2115    '''Loop over histograms and compute derivatives of the fitting
2116    model (M) with respect to all parameters.  For each histogram, the
2117    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2118    number of parameters and m is the number of data points *in the
2119    histogram*. The (n by n) Hessian is computed from each Jacobian
2120    and it is returned.  This routine is used when refinement
2121    derivatives are selected as "analtytic Hessian" in Controls.
2122
2123    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2124    '''
2125    parmDict.update(zip(varylist,values))
2126    G2mv.Dict2Map(parmDict,varylist)
2127    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2128    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2129    nvar = len(varylist)
2130    Hess = np.empty(0)
2131    histoList = Histograms.keys()
2132    histoList.sort()
2133    for histogram in histoList:
2134        if 'PWDR' in histogram[:4]:
2135            Histogram = Histograms[histogram]
2136            hId = Histogram['hId']
2137            hfx = ':%d:'%(hId)
2138            wtFactor = calcControls[hfx+'wtFactor']
2139            Limits = calcControls[hfx+'Limits']
2140            x,y,w,yc,yb,yd = Histogram['Data']
2141            W = wtFactor*w
2142            dy = y-yc
2143            xB = np.searchsorted(x,Limits[0])
2144            xF = np.searchsorted(x,Limits[1])
2145            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2146                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2147            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2148            Dy = dy[xB:xF][np.newaxis,:]
2149            dMdvh *= Wt
2150            if dlg:
2151                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2152            if len(Hess):
2153                Hess += np.inner(dMdvh,dMdvh)
2154                dMdvh *= Wt*Dy
2155                Vec += np.sum(dMdvh,axis=1)
2156            else:
2157                Hess = np.inner(dMdvh,dMdvh)
2158                dMdvh *= Wt*Dy
2159                Vec = np.sum(dMdvh,axis=1)
2160        elif 'HKLF' in histogram[:4]:
2161            Histogram = Histograms[histogram]
2162            phase = Histogram['Reflection Lists']
2163            Phase = Phases[phase]
2164            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2165            hId = Histogram['hId']
2166            hfx = ':%d:'%(Histogram['hId'])
2167            wtFactor = calcControls[hfx+'wtFactor']
2168            # now process derivatives in constraints
2169            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2170#            print 'matrix build time: %.3f'%(time.time()-time0)
2171
2172            if dlg:
2173                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2174            if len(Hess):
2175                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2176                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2177            else:
2178                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2179                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2180        else:
2181            continue        #skip non-histogram entries
2182    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2183    if np.any(pVals):
2184        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2185        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2186        Hess += np.inner(dpdv*pWt,dpdv)
2187    return Vec,Hess
2188
2189def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
2190    'Needs a doc string'
2191    Values2Dict(parmDict, varylist, values)
2192    G2mv.Dict2Map(parmDict,varylist)
2193    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2194    M = np.empty(0)
2195    SumwYo = 0
2196    Nobs = 0
2197    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2198    histoList = Histograms.keys()
2199    histoList.sort()
2200    for histogram in histoList:
2201        if 'PWDR' in histogram[:4]:
2202            Histogram = Histograms[histogram]
2203            hId = Histogram['hId']
2204            hfx = ':%d:'%(hId)
2205            wtFactor = calcControls[hfx+'wtFactor']
2206            Limits = calcControls[hfx+'Limits']
2207            x,y,w,yc,yb,yd = Histogram['Data']
2208            yc *= 0.0                           #zero full calcd profiles
2209            yb *= 0.0
2210            yd *= 0.0
2211            xB = np.searchsorted(x,Limits[0])
2212            xF = np.searchsorted(x,Limits[1])
2213            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2214                varylist,Histogram,Phases,calcControls,pawleyLookup)
2215            yc[xB:xF] += yb[xB:xF]
2216            if not np.any(y):                   #fill dummy data
2217                rv = st.poisson(yc[xB:xF])
2218                y[xB:xF] = rv.rvs()
2219                Z = np.ones_like(yc[xB:xF])
2220                Z[1::2] *= -1
2221                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2222                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2223            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2224            W = wtFactor*w
2225            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2226            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2227            Nobs += Histogram['Residuals']['Nobs']
2228            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2229            SumwYo += Histogram['Residuals']['sumwYo']
2230            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2231            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2232            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2233            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2234            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.))
2235            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.))
2236            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2237            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2238            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2239            if dlg:
2240                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2241            M = np.concatenate((M,wdy))
2242#end of PWDR processing
2243        elif 'HKLF' in histogram[:4]:
2244            Histogram = Histograms[histogram]
2245            Histogram['Residuals'] = {}
2246            phase = Histogram['Reflection Lists']
2247            Phase = Phases[phase]
2248            hId = Histogram['hId']
2249            hfx = ':%d:'%(hId)
2250            wtFactor = calcControls[hfx+'wtFactor']
2251            pfx = '%d::'%(Phase['pId'])
2252            phfx = '%d:%d:'%(Phase['pId'],hId)
2253            SGData = Phase['General']['SGData']
2254            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2255            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2256            refDict = Histogram['Data']
2257            time0 = time.time()
2258            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2259#            print 'sf-calc time: %.3f'%(time.time()-time0)
2260            df = np.zeros(len(refDict['RefList']))
2261            sumwYo = 0
2262            sumFo = 0
2263            sumFo2 = 0
2264            sumdF = 0
2265            sumdF2 = 0
2266            nobs = 0
2267            if calcControls['F**2']:
2268                for i,ref in enumerate(refDict['RefList']):
2269                    if ref[6] > 0:
2270                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2271                        w = 1.0/ref[6]
2272                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]  #correct Fc^2 for extinction
2273                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2274                        if w*ref[5] >= calcControls['minF/sig']:
2275                            Fo = np.sqrt(ref[5])
2276                            sumFo += Fo
2277                            sumFo2 += ref[5]
2278                            sumdF += abs(Fo-np.sqrt(ref[7]))
2279                            sumdF2 += abs(ref[5]-ref[7])
2280                            nobs += 1
2281                            df[i] = -w*(ref[5]-ref[7])
2282                            sumwYo += (w*ref[5])**2
2283            else:
2284                for i,ref in enumerate(refDict['RefList']):
2285                    if ref[5] > 0.:
2286                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2287                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]    #correct Fc^2 for extinction
2288                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2289                        Fo = np.sqrt(ref[5])
2290                        Fc = np.sqrt(ref[7])
2291                        w = 2.0*Fo/ref[6]
2292                        if w*Fo >= calcControls['minF/sig']:
2293                            sumFo += Fo
2294                            sumFo2 += ref[5]
2295                            sumdF += abs(Fo-Fc)
2296                            sumdF2 += abs(ref[5]-ref[7])
2297                            nobs += 1
2298                            df[i] = -w*(Fo-Fc)
2299                            sumwYo += (w*Fo)**2
2300            Histogram['Residuals']['Nobs'] = nobs
2301            Histogram['Residuals']['sumwYo'] = sumwYo
2302            SumwYo += sumwYo
2303            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2304            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2305            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2306            Histogram['Residuals'][phfx+'Nref'] = nobs
2307            Nobs += nobs
2308            if dlg:
2309                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2310            M = np.concatenate((M,wtFactor*df))
2311# end of HKLF processing
2312    Histograms['sumwYo'] = SumwYo
2313    Histograms['Nobs'] = Nobs
2314    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2315    if dlg:
2316        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2317        if not GoOn:
2318            parmDict['saved values'] = values
2319            dlg.Destroy()
2320            raise Exception         #Abort!!
2321    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2322    if len(pVals):
2323        pSum = np.sum(pWt*pVals**2)
2324        for name in pWsum:
2325            if pWsum:
2326                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2327        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2328        Nobs += len(pVals)
2329        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2330    return M
2331                       
Note: See TracBrowser for help on using the repository browser.