source: trunk/GSASIIstrMath.py @ 1498

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

try to implement drag/drop of tree items - commented out as it doesn't work
add 'ShowCell?' to Controls to avoid cell errors for sequential single peak fitting when normal sequential refinement has been done
fix a neg. peak width error

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