source: trunk/GSASIIstrMath.py @ 1479

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

fix Stephens model for mustrain, fxn & CW derivs are correct

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