source: trunk/GSASIIstrMath.py @ 1474

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

1st MC/SA tutorial
various MC/SA fixes
fix to background peak fitting for CW & TOF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 115.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2014-08-21 18:31:59 +0000 (Thu, 21 Aug 2014) $
8# $Author: vondreele $
9# $Revision: 1474 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1474 2014-08-21 18:31:59Z 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: 1474 $")
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,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            pwrs = calcControls[phfx+'MuPwrs']
1189            sum = 0
1190            for i,pwr in enumerate(pwrs):
1191                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1192            Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*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']*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']*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1218        else:       #generalized - P.W. Stephens model
1219            pwrs = calcControls[phfx+'MuPwrs']
1220            sum = 0
1221            for i,pwr in enumerate(pwrs):
1222                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
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,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            pwrs = calcControls[phfx+'MuPwrs']
1292            const = 0.018*refl[4]**2*tanth
1293            sum = 0
1294            for i,pwr in enumerate(pwrs):
1295                term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1296                sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1297                gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1298                sigDict[phfx+'Mustrain:'+str(i)] = \
1299                    2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1300            Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1301            for i in range(len(pwrs)):
1302                sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1303        gamDict[phfx+'Mustrain;mx'] = Mgam
1304        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1305    else:   #'T'OF
1306        if calcControls[phfx+'SizeType'] == 'isotropic':
1307            Sgam = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;i']
1308            gamDict[phfx+'Size;i'] = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;mx']
1309            sigDict[phfx+'Size;i'] = 1.e-4*parmDict[hfx+'difC']*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1310        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1311            const = 1.e-4*parmDict[hfx+'difC']
1312            H = np.array(refl[:3])
1313            P = np.array(calcControls[phfx+'SizeAxis'])
1314            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1315            Si = parmDict[phfx+'Size;i']
1316            Sa = parmDict[phfx+'Size;a']
1317            gami = const*(Si*Sa)
1318            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1319            Sgam = gami*sqtrm
1320            dsi = gami*Si*cosP**2/sqtrm-gam/Si
1321            dsa = gami*Sa*sinP**2/sqtrm-gam/Sa
1322            gamDict[phfx+'Size;i'] = const*parmDict[phfx+'Size;mx']*Sa
1323            gamDict[phfx+'Size;a'] = const*parmDict[phfx+'Size;mx']*Si
1324            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1325            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1326        else:           #ellipsoidal crystallites
1327            const = 1.e-4*parmDict[hfx+'difC']
1328            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1329            H = np.array(refl[:3])
1330            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1331            Sgam = const/lenR
1332            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1333                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1334                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1335        gamDict[phfx+'Size;mx'] = Sgam
1336        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1337               
1338        #microstrain derivatives               
1339        if calcControls[phfx+'MustrainType'] == 'isotropic':
1340            Mgam = 1.e-6*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;i']
1341            gamDict[phfx+'Mustrain;i'] =  1.e-6*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']
1342            sigDict[phfx+'Mustrain;i'] =  2.e-6*parmDict[hfx+'difC']*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1343        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1344            H = np.array(refl[:3])
1345            P = np.array(calcControls[phfx+'MustrainAxis'])
1346            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1347            Si = parmDict[phfx+'Mustrain;i']
1348            Sa = parmDict[phfx+'Mustrain;a']
1349            gami = 1.e-6*parmDict[hfx+'difC']*Si*Sa
1350            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1351            Mgam = gami/sqtrm
1352            dsi = -gami*Si*cosP**2/sqtrm**3
1353            dsa = -gami*Sa*sinP**2/sqtrm**3
1354            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1355            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1356            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1357            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1358        else:       #generalized - P.W. Stephens model
1359            pwrs = calcControls[phfx+'MuPwrs']
1360            const = 1.e-6*parmDict[hfx+'difC']*refl[4]**2
1361            sum = 0
1362            for i,pwr in enumerate(pwrs):
1363                term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1364                sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1365                gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1366                sigDict[phfx+'Mustrain:'+str(i)] = \
1367                    2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1368            Mgam = const*sum
1369            for i in range(len(pwrs)):
1370                sigDict[phfx+'Mustrain:'+str(i)] *= Mgam       
1371        gamDict[phfx+'Mustrain;mx'] = Mgam
1372        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1373       
1374    return sigDict,gamDict
1375       
1376def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1377    'Needs a doc string'
1378    h,k,l = refl[:3]
1379    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1380    d = np.sqrt(dsq)
1381
1382    refl[4] = d
1383    if 'C' in calcControls[hfx+'histType']:
1384        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1385        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1386        if 'Bragg' in calcControls[hfx+'instType']:
1387            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1388                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1389        else:               #Debye-Scherrer - simple but maybe not right
1390            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1391    elif 'T' in calcControls[hfx+'histType']:
1392        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']*d**3+parmDict[hfx+'Zero']
1393        #do I need sample position effects - maybe?
1394    return pos
1395
1396def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1397    'Needs a doc string'
1398    dpr = 180./np.pi
1399    h,k,l = refl[:3]
1400    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1401    dst = np.sqrt(dstsq)
1402    dsp = 1./dst
1403    if 'C' in calcControls[hfx+'histType']:
1404        pos = refl[5]-parmDict[hfx+'Zero']
1405        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1406        dpdw = const*dst
1407        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1408        dpdA *= const*wave/(2.0*dst)
1409        dpdZ = 1.0
1410        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1411        if 'Bragg' in calcControls[hfx+'instType']:
1412            dpdSh = -4.*const*cosd(pos/2.0)
1413            dpdTr = -const*sind(pos)*100.0
1414            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1415        else:               #Debye-Scherrer - simple but maybe not right
1416            dpdXd = -const*cosd(pos)
1417            dpdYd = -const*sind(pos)
1418            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1419    elif 'T' in calcControls[hfx+'histType']:
1420        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1421        dpdZ = 1.0
1422        dpdDC = dsp
1423        dpdDA = dsp**2
1424        dpdDB = dsp**3
1425        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
1426           
1427def GetHStrainShift(refl,SGData,phfx,parmDict):
1428    'Needs a doc string'
1429    laue = SGData['SGLaue']
1430    uniq = SGData['SGUniq']
1431    h,k,l = refl[:3]
1432    if laue in ['m3','m3m']:
1433        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1434            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1435    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1436        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1437    elif laue in ['3R','3mR']:
1438        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1439    elif laue in ['4/m','4/mmm']:
1440        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1441    elif laue in ['mmm']:
1442        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1443    elif laue in ['2/m']:
1444        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1445        if uniq == 'a':
1446            Dij += parmDict[phfx+'D23']*k*l
1447        elif uniq == 'b':
1448            Dij += parmDict[phfx+'D13']*h*l
1449        elif uniq == 'c':
1450            Dij += parmDict[phfx+'D12']*h*k
1451    else:
1452        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1453            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1454    return -Dij*refl[4]**2*tand(refl[5]/2.0)
1455           
1456def GetHStrainShiftDerv(refl,SGData,phfx):
1457    'Needs a doc string'
1458    laue = SGData['SGLaue']
1459    uniq = SGData['SGUniq']
1460    h,k,l = refl[:3]
1461    if laue in ['m3','m3m']:
1462        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1463            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1464    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1465        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1466    elif laue in ['3R','3mR']:
1467        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1468    elif laue in ['4/m','4/mmm']:
1469        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1470    elif laue in ['mmm']:
1471        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1472    elif laue in ['2/m']:
1473        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1474        if uniq == 'a':
1475            dDijDict[phfx+'D23'] = k*l
1476        elif uniq == 'b':
1477            dDijDict[phfx+'D13'] = h*l
1478        elif uniq == 'c':
1479            dDijDict[phfx+'D12'] = h*k
1480            names.append()
1481    else:
1482        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1483            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1484    for item in dDijDict:
1485        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
1486    return dDijDict
1487   
1488def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1489    'Needs a doc string'
1490    histoList = Histograms.keys()
1491    histoList.sort()
1492    for histogram in histoList:
1493        if 'PWDR' in histogram[:4]:
1494            Histogram = Histograms[histogram]
1495            hId = Histogram['hId']
1496            hfx = ':%d:'%(hId)
1497            Limits = calcControls[hfx+'Limits']
1498            if 'C' in calcControls[hfx+'histType']:
1499                shl = max(parmDict[hfx+'SH/L'],0.0005)
1500                Ka2 = False
1501                kRatio = 0.0
1502                if hfx+'Lam1' in parmDict.keys():
1503                    Ka2 = True
1504                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1505                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1506            x,y,w,yc,yb,yd = Histogram['Data']
1507            xB = np.searchsorted(x,Limits[0])
1508            xF = np.searchsorted(x,Limits[1])
1509            ymb = np.array(y-yb)
1510            ymb = np.where(ymb,ymb,1.0)
1511            ycmb = np.array(yc-yb)
1512            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1513            refLists = Histogram['Reflection Lists']
1514            for phase in refLists:
1515                Phase = Phases[phase]
1516                pId = Phase['pId']
1517                phfx = '%d:%d:'%(pId,hId)
1518                refDict = refLists[phase]
1519                sumFo = 0.0
1520                sumdF = 0.0
1521                sumFosq = 0.0
1522                sumdFsq = 0.0
1523                for refl in refDict['RefList']:
1524                    if 'C' in calcControls[hfx+'histType']:
1525                        yp = np.zeros_like(yb)
1526                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1527                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1528                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1529                        iFin2 = iFin
1530                        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
1531                        if Ka2:
1532                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1533                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1534                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1535                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1536                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1537                                continue
1538                            elif not iBeg2-iFin2:     #peak above high limit - done
1539                                break
1540                            yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
1541                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0))
1542                    elif 'T' in calcControls[hfx+'histType']:
1543                        yp = np.zeros_like(yb)
1544                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1545                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1546                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1547                        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
1548                        refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
1549                    Fo = np.sqrt(np.abs(refl[8]))
1550                    Fc = np.sqrt(np.abs(refl[9]))
1551                    sumFo += Fo
1552                    sumFosq += refl[8]**2
1553                    sumdF += np.abs(Fo-Fc)
1554                    sumdFsq += (refl[8]-refl[9])**2
1555                Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1556                Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1557                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1558                Histogram['Residuals']['hId'] = hId
1559        elif 'HKLF' in histogram[:4]:
1560            Histogram = Histograms[histogram]
1561            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1562               
1563def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1564    'Needs a doc string'
1565   
1566    def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
1567        U = parmDict[hfx+'U']
1568        V = parmDict[hfx+'V']
1569        W = parmDict[hfx+'W']
1570        X = parmDict[hfx+'X']
1571        Y = parmDict[hfx+'Y']
1572        tanPos = tand(refl[5]/2.0)
1573        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)
1574        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1575        sig = max(0.001,sig)
1576        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1577        gam = max(0.001,gam)
1578        return sig,gam
1579               
1580    def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
1581        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4]**2+   \
1582            parmDict[hfx+'sig-2']*refl[4]**4+parmDict[hfx+'sig-q']/refl[4]**2
1583        gam = parmDict[hfx+'X']*refl[4]+parmDict[hfx+'Y']*refl[4]**2
1584        Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict)
1585        sig += Ssig     #save peak sigma
1586        sig = max(0.001,sig)
1587        gam += Sgam     #save peak gamma
1588        gam = max(0.001,gam)
1589        return sig,gam
1590       
1591    def GetReflAlpBet(refl,hfx,parmDict):
1592        alp = parmDict[hfx+'alpha']/refl[4]
1593        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]**4+parmDict[hfx+'beta-q']/refl[4]**2
1594        return alp,bet
1595               
1596    hId = Histogram['hId']
1597    hfx = ':%d:'%(hId)
1598    bakType = calcControls[hfx+'bakType']
1599    yb = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1600    yc = np.zeros_like(yb)
1601    cw = np.diff(x)
1602    cw = np.append(cw,cw[-1])
1603       
1604    if 'C' in calcControls[hfx+'histType']:   
1605        shl = max(parmDict[hfx+'SH/L'],0.002)
1606        Ka2 = False
1607        if hfx+'Lam1' in parmDict.keys():
1608            wave = parmDict[hfx+'Lam1']
1609            Ka2 = True
1610            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1611            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1612        else:
1613            wave = parmDict[hfx+'Lam']
1614    for phase in Histogram['Reflection Lists']:
1615        refDict = Histogram['Reflection Lists'][phase]
1616        Phase = Phases[phase]
1617        pId = Phase['pId']
1618        pfx = '%d::'%(pId)
1619        phfx = '%d:%d:'%(pId,hId)
1620        hfx = ':%d:'%(hId)
1621        SGData = Phase['General']['SGData']
1622        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1623        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #Do I want to modify by Dij?
1624        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1625        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1626        Vst = np.sqrt(nl.det(G))    #V*
1627        if not Phase['General'].get('doPawley'):
1628            time0 = time.time()
1629            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1630#            print 'sf calc time: %.3fs'%(time.time()-time0)
1631        time0 = time.time()
1632        for iref,refl in enumerate(refDict['RefList']):
1633            if 'C' in calcControls[hfx+'histType']:
1634                h,k,l = refl[:3]
1635                Uniq = np.inner(refl[:3],SGMT)
1636                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1637                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1638                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1639                refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1640                refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1641                refl[11] *= Vst*Lorenz
1642                 
1643                if Phase['General'].get('doPawley'):
1644                    try:
1645                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1646                        refl[9] = parmDict[pInd]
1647                    except KeyError:
1648#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1649                        continue
1650                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1651                iBeg = np.searchsorted(x,refl[5]-fmin)
1652                iFin = np.searchsorted(x,refl[5]+fmax)
1653                if not iBeg+iFin:       #peak below low limit - skip peak
1654                    continue
1655                elif not iBeg-iFin:     #peak above high limit - done
1656                    break
1657                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
1658                if Ka2:
1659                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1660                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1661                    iBeg = np.searchsorted(x,pos2-fmin)
1662                    iFin = np.searchsorted(x,pos2+fmax)
1663                    if not iBeg+iFin:       #peak below low limit - skip peak
1664                        continue
1665                    elif not iBeg-iFin:     #peak above high limit - done
1666                        return yc,yb
1667                    yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1668            elif 'T' in calcControls[hfx+'histType']:
1669                h,k,l = refl[:3]
1670                Uniq = np.inner(refl[:3],SGMT)
1671                refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict)         #corrected reflection position
1672                Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4]**4                                                #TOF Lorentz correction
1673                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1674                refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1675                refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
1676                refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1677                refl[11] *= Vst*Lorenz
1678                if Phase['General'].get('doPawley'):
1679                    try:
1680                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1681                        refl[9] = parmDict[pInd]
1682                    except KeyError:
1683#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1684                        continue
1685                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1686                iBeg = np.searchsorted(x,refl[5]-fmin)
1687                iFin = np.searchsorted(x,refl[5]+fmax)
1688                if not iBeg+iFin:       #peak below low limit - skip peak
1689                    continue
1690                elif not iBeg-iFin:     #peak above high limit - done
1691                    break
1692                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]
1693#        print 'profile calc time: %.3fs'%(time.time()-time0)
1694    return yc,yb
1695   
1696def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1697    'Needs a doc string'
1698   
1699    def cellVaryDerv(pfx,SGData,dpdA): 
1700        if SGData['SGLaue'] in ['-1',]:
1701            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1702                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1703        elif SGData['SGLaue'] in ['2/m',]:
1704            if SGData['SGUniq'] == 'a':
1705                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1706            elif SGData['SGUniq'] == 'b':
1707                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1708            else:
1709                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1710        elif SGData['SGLaue'] in ['mmm',]:
1711            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1712        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1713            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1714        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1715            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1716        elif SGData['SGLaue'] in ['3R', '3mR']:
1717            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1718        elif SGData['SGLaue'] in ['m3m','m3']:
1719            return [[pfx+'A0',dpdA[0]]]
1720           
1721    # create a list of dependent variables and set up a dictionary to hold their derivatives
1722    dependentVars = G2mv.GetDependentVars()
1723    depDerivDict = {}
1724    for j in dependentVars:
1725        depDerivDict[j] = np.zeros(shape=(len(x)))
1726    #print 'dependent vars',dependentVars
1727    lenX = len(x)               
1728    hId = Histogram['hId']
1729    hfx = ':%d:'%(hId)
1730    bakType = calcControls[hfx+'bakType']
1731    dMdv = np.zeros(shape=(len(varylist),len(x)))
1732    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1733    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1734        bBpos =varylist.index(hfx+'Back:0')
1735        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1736    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1737    for name in varylist:
1738        if 'Debye' in name:
1739            id = int(name.split(':')[-1])
1740            parm = name[:int(name.rindex(':'))]
1741            ip = names.index(parm)
1742            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1743    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1744    for name in varylist:
1745        if 'BkPk' in name:
1746            parm,id = name.split(';')
1747            id = int(id)
1748            if parm in names:
1749                ip = names.index(parm)
1750                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1751    cw = np.diff(x)
1752    cw = np.append(cw,cw[-1])
1753    Ka2 = False #also for TOF!
1754    if 'C' in calcControls[hfx+'histType']:   
1755        shl = max(parmDict[hfx+'SH/L'],0.002)
1756        if hfx+'Lam1' in parmDict.keys():
1757            wave = parmDict[hfx+'Lam1']
1758            Ka2 = True
1759            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1760            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1761        else:
1762            wave = parmDict[hfx+'Lam']
1763    for phase in Histogram['Reflection Lists']:
1764        refDict = Histogram['Reflection Lists'][phase]
1765        Phase = Phases[phase]
1766        SGData = Phase['General']['SGData']
1767        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1768        pId = Phase['pId']
1769        pfx = '%d::'%(pId)
1770        phfx = '%d:%d:'%(pId,hId)
1771        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #And modify here by Dij?
1772        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1773        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1774        if not Phase['General'].get('doPawley'):
1775            time0 = time.time()
1776            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1777#            print 'sf-derv time %.3fs'%(time.time()-time0)
1778            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1779        time0 = time.time()
1780        for iref,refl in enumerate(refDict['RefList']):
1781            h,k,l = refl[:3]
1782            Uniq = np.inner(refl[:3],SGMT)
1783            if 'T' in calcControls[hfx+'histType']:
1784                wave = refl[14]
1785            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1786            if 'C' in calcControls[hfx+'histType']:        #CW powder
1787                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1788            else: #'T'OF
1789                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1790            iBeg = np.searchsorted(x,refl[5]-fmin)
1791            iFin = np.searchsorted(x,refl[5]+fmax)
1792            if not iBeg+iFin:       #peak below low limit - skip peak
1793                continue
1794            elif not iBeg-iFin:     #peak above high limit - done
1795                break
1796            pos = refl[5]
1797            if 'C' in calcControls[hfx+'histType']:
1798                tanth = tand(pos/2.0)
1799                costh = cosd(pos/2.0)
1800                lenBF = iFin-iBeg
1801                dMdpk = np.zeros(shape=(6,lenBF))
1802                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1803                for i in range(5):
1804                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
1805                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1806                if Ka2:
1807                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1808                    iBeg2 = np.searchsorted(x,pos2-fmin)
1809                    iFin2 = np.searchsorted(x,pos2+fmax)
1810                    if iBeg2-iFin2:
1811                        lenBF2 = iFin2-iBeg2
1812                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1813                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1814                        for i in range(5):
1815                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
1816                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
1817                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1818            else:   #'T'OF
1819                lenBF = iFin-iBeg
1820                dMdpk = np.zeros(shape=(6,lenBF))
1821                dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1822                for i in range(6):
1823                    dMdpk[i] += refl[11]*refl[9]*dMdipk[i]      #cw[iBeg:iFin]*
1824                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
1825            if Phase['General'].get('doPawley'):
1826                dMdpw = np.zeros(len(x))
1827                try:
1828                    pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1829                    idx = varylist.index(pIdx)
1830                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1831                    if Ka2: #not for TOF either
1832                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1833                    dMdv[idx] = dMdpw
1834                except: # ValueError:
1835                    pass
1836            if 'C' in calcControls[hfx+'histType']:
1837                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1838                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1839                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1840                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1841                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1842                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1843                    hfx+'DisplaceY':[dpdY,'pos'],}
1844                if 'Bragg' in calcControls[hfx+'instType']:
1845                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1846                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
1847                else:
1848                    names.update({hfx+'Absorption':[dFdAb,'int'],})
1849            else:   #'T'OF
1850                dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1851                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1852                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1853                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1854                    hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1855                    hfx+'beta-q':[1./refl[4]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'],
1856                    hfx+'sig-2':[refl[4]**4,'sig'],hfx+'sig-q':[1./refl[4]**2,'sig'],
1857                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1858            for name in names:
1859                item = names[name]
1860                if name in varylist:
1861                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1862                    if Ka2:
1863                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1864                elif name in dependentVars:
1865                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1866                    if Ka2:
1867                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1868            for iPO in dIdPO:
1869                if iPO in varylist:
1870                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1871                    if Ka2:
1872                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1873                elif iPO in dependentVars:
1874                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1875                    if Ka2:
1876                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1877            for i,name in enumerate(['omega','chi','phi']):
1878                aname = pfx+'SH '+name
1879                if aname in varylist:
1880                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1881                    if Ka2:
1882                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1883                elif aname in dependentVars:
1884                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1885                    if Ka2:
1886                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1887            for iSH in dFdODF:
1888                if iSH in varylist:
1889                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1890                    if Ka2:
1891                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1892                elif iSH in dependentVars:
1893                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1894                    if Ka2:
1895                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1896            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1897            for name,dpdA in cellDervNames:
1898                if name in varylist:
1899                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1900                    if Ka2:
1901                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1902                elif name in dependentVars:
1903                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1904                    if Ka2:
1905                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1906            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1907            for name in dDijDict:
1908                if name in varylist:
1909                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1910                    if Ka2:
1911                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1912                elif name in dependentVars:
1913                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1914                    if Ka2:
1915                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1916            if 'C' in calcControls[hfx+'histType']:
1917                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)
1918            else:   #'T'OF
1919                sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict)
1920            for name in gamDict:
1921                if name in varylist:
1922                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1923                    if Ka2:
1924                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1925                elif name in dependentVars:
1926                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1927                    if Ka2:
1928                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1929            for name in sigDict:
1930                if name in varylist:
1931                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1932                    if Ka2:
1933                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1934                elif name in dependentVars:
1935                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1936                    if Ka2:
1937                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1938            for name in ['BabA','BabU']:
1939                if refl[9]:
1940                    if phfx+name in varylist:
1941                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1942                        if Ka2:
1943                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1944                    elif phfx+name in dependentVars:                   
1945                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1946                        if Ka2:
1947                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
1948            if not Phase['General'].get('doPawley'):
1949                #do atom derivatives -  for RB,F,X & U so far             
1950                corr = dervDict['int']/refl[9]
1951                if Ka2:
1952                    corr2 = dervDict2['int']/refl[9]
1953                for name in varylist+dependentVars:
1954                    if '::RBV;' in name:
1955                        pass
1956                    else:
1957                        try:
1958                            aname = name.split(pfx)[1][:2]
1959                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1960                        except IndexError:
1961                            continue
1962                    if name in varylist:
1963                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1964                        if Ka2:
1965                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1966                    elif name in dependentVars:
1967                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1968                        if Ka2:
1969                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1970    #        print 'profile derv time: %.3fs'%(time.time()-time0)
1971    # now process derivatives in constraints
1972    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1973    return dMdv
1974   
1975def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1976    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1977    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
1978    are returned to either dervRefine or HessRefine.
1979
1980    :returns:
1981    '''
1982    nobs = Histogram['Residuals']['Nobs']
1983    hId = Histogram['hId']
1984    hfx = ':%d:'%(hId)
1985    pfx = '%d::'%(Phase['pId'])
1986    phfx = '%d:%d:'%(Phase['pId'],hId)
1987    SGData = Phase['General']['SGData']
1988    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1989    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1990    refDict = Histogram['Data']
1991    dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1992    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1993    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1994    dependentVars = G2mv.GetDependentVars()
1995    depDerivDict = {}
1996    for j in dependentVars:
1997        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1998    wdf = np.zeros(len(refDict['RefList']))
1999    if calcControls['F**2']:
2000        for iref,ref in enumerate(refDict['RefList']):
2001            if ref[6] > 0:
2002                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:] 
2003                w = 1.0/ref[6]
2004                if w*ref[5] >= calcControls['minF/sig']:
2005                    wdf[iref] = w*(ref[5]-ref[7])
2006                    for j,var in enumerate(varylist):
2007                        if var in dFdvDict:
2008                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]   #*dervCor
2009                    for var in dependentVars:
2010                        if var in dFdvDict:
2011                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]  #*dervCor
2012                    if phfx+'Scale' in varylist:
2013                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
2014                    elif phfx+'Scale' in dependentVars:
2015                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
2016                    for item in ['Ep','Es','Eg']:
2017                        if phfx+item in varylist and dervDict:
2018                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor
2019                        elif phfx+item in dependentVars and dervDict:
2020                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
2021                    for item in ['BabA','BabU']:
2022                        if phfx+item in varylist:
2023                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2024                        elif phfx+item in dependentVars:
2025                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2026    else:
2027        for iref,ref in enumerate(refDict['RefList']):
2028            if ref[5] > 0.:
2029                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:]
2030                Fo = np.sqrt(ref[5])
2031                Fc = np.sqrt(ref[7])
2032                w = 1.0/ref[6]
2033                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2034                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2035                    for j,var in enumerate(varylist):
2036                        if var in dFdvDict:
2037                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]       #*dervCor
2038                    for var in dependentVars:
2039                        if var in dFdvDict:
2040                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]      #*dervCor
2041                    if phfx+'Scale' in varylist:
2042                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]    #*dervCor
2043                    elif phfx+'Scale' in dependentVars:
2044                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11] #*dervCor                           
2045                    for item in ['Ep','Es','Eg']:
2046                        if phfx+item in varylist and dervDict:
2047                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #correct
2048                        elif phfx+item in dependentVars and dervDict:
2049                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
2050                    for item in ['BabA','BabU']:
2051                        if phfx+item in varylist:
2052                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2053                        elif phfx+item in dependentVars:
2054                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2055    return dMdvh,depDerivDict,wdf
2056   
2057
2058def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2059    '''Loop over histograms and compute derivatives of the fitting
2060    model (M) with respect to all parameters.  Results are returned in
2061    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2062    n is the number of parameters and m is the number of data
2063    points. This can exceed memory when m gets large. This routine is
2064    used when refinement derivatives are selected as "analtytic
2065    Jacobian" in Controls.
2066
2067    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2068    '''
2069    parmDict.update(zip(varylist,values))
2070    G2mv.Dict2Map(parmDict,varylist)
2071    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2072    nvar = len(varylist)
2073    dMdv = np.empty(0)
2074    histoList = Histograms.keys()
2075    histoList.sort()
2076    for histogram in histoList:
2077        if 'PWDR' in histogram[:4]:
2078            Histogram = Histograms[histogram]
2079            hId = Histogram['hId']
2080            hfx = ':%d:'%(hId)
2081            wtFactor = calcControls[hfx+'wtFactor']
2082            Limits = calcControls[hfx+'Limits']
2083            x,y,w,yc,yb,yd = Histogram['Data']
2084            W = wtFactor*w
2085            xB = np.searchsorted(x,Limits[0])
2086            xF = np.searchsorted(x,Limits[1])
2087            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2088                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2089        elif 'HKLF' in histogram[:4]:
2090            Histogram = Histograms[histogram]
2091            phase = Histogram['Reflection Lists']
2092            Phase = Phases[phase]
2093            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2094            hfx = ':%d:'%(Histogram['hId'])
2095            wtFactor = calcControls[hfx+'wtFactor']
2096            # now process derivatives in constraints
2097            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2098        else:
2099            continue        #skip non-histogram entries
2100        if len(dMdv):
2101            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2102        else:
2103            dMdv = np.sqrt(wtFactor)*dMdvh
2104           
2105    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2106    if np.any(pVals):
2107        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2108        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2109       
2110    return dMdv
2111
2112def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2113    '''Loop over histograms and compute derivatives of the fitting
2114    model (M) with respect to all parameters.  For each histogram, the
2115    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2116    number of parameters and m is the number of data points *in the
2117    histogram*. The (n by n) Hessian is computed from each Jacobian
2118    and it is returned.  This routine is used when refinement
2119    derivatives are selected as "analtytic Hessian" in Controls.
2120
2121    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2122    '''
2123    parmDict.update(zip(varylist,values))
2124    G2mv.Dict2Map(parmDict,varylist)
2125    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2126    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2127    nvar = len(varylist)
2128    Hess = np.empty(0)
2129    histoList = Histograms.keys()
2130    histoList.sort()
2131    for histogram in histoList:
2132        if 'PWDR' in histogram[:4]:
2133            Histogram = Histograms[histogram]
2134            hId = Histogram['hId']
2135            hfx = ':%d:'%(hId)
2136            wtFactor = calcControls[hfx+'wtFactor']
2137            Limits = calcControls[hfx+'Limits']
2138            x,y,w,yc,yb,yd = Histogram['Data']
2139            W = wtFactor*w
2140            dy = y-yc
2141            xB = np.searchsorted(x,Limits[0])
2142            xF = np.searchsorted(x,Limits[1])
2143            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2144                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2145            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2146            Dy = dy[xB:xF][np.newaxis,:]
2147            dMdvh *= Wt
2148            if dlg:
2149                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2150            if len(Hess):
2151                Hess += np.inner(dMdvh,dMdvh)
2152                dMdvh *= Wt*Dy
2153                Vec += np.sum(dMdvh,axis=1)
2154            else:
2155                Hess = np.inner(dMdvh,dMdvh)
2156                dMdvh *= Wt*Dy
2157                Vec = np.sum(dMdvh,axis=1)
2158        elif 'HKLF' in histogram[:4]:
2159            Histogram = Histograms[histogram]
2160            phase = Histogram['Reflection Lists']
2161            Phase = Phases[phase]
2162            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2163            hId = Histogram['hId']
2164            hfx = ':%d:'%(Histogram['hId'])
2165            wtFactor = calcControls[hfx+'wtFactor']
2166            # now process derivatives in constraints
2167            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2168#            print 'matrix build time: %.3f'%(time.time()-time0)
2169
2170            if dlg:
2171                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2172            if len(Hess):
2173                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2174                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2175            else:
2176                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2177                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2178        else:
2179            continue        #skip non-histogram entries
2180    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2181    if np.any(pVals):
2182        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2183        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2184        Hess += np.inner(dpdv*pWt,dpdv)
2185    return Vec,Hess
2186
2187def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
2188    'Needs a doc string'
2189    Values2Dict(parmDict, varylist, values)
2190    G2mv.Dict2Map(parmDict,varylist)
2191    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2192    M = np.empty(0)
2193    SumwYo = 0
2194    Nobs = 0
2195    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2196    histoList = Histograms.keys()
2197    histoList.sort()
2198    for histogram in histoList:
2199        if 'PWDR' in histogram[:4]:
2200            Histogram = Histograms[histogram]
2201            hId = Histogram['hId']
2202            hfx = ':%d:'%(hId)
2203            wtFactor = calcControls[hfx+'wtFactor']
2204            Limits = calcControls[hfx+'Limits']
2205            x,y,w,yc,yb,yd = Histogram['Data']
2206            yc *= 0.0                           #zero full calcd profiles
2207            yb *= 0.0
2208            yd *= 0.0
2209            xB = np.searchsorted(x,Limits[0])
2210            xF = np.searchsorted(x,Limits[1])
2211            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2212                varylist,Histogram,Phases,calcControls,pawleyLookup)
2213            yc[xB:xF] += yb[xB:xF]
2214            if not np.any(y):                   #fill dummy data
2215                rv = st.poisson(yc[xB:xF])
2216                y[xB:xF] = rv.rvs()
2217                Z = np.ones_like(yc[xB:xF])
2218                Z[1::2] *= -1
2219                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2220                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2221            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2222            W = wtFactor*w
2223            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2224            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2225            Nobs += Histogram['Residuals']['Nobs']
2226            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2227            SumwYo += Histogram['Residuals']['sumwYo']
2228            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2229            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2230            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2231            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2232            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.))
2233            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.))
2234            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2235            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2236            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2237            if dlg:
2238                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2239            M = np.concatenate((M,wdy))
2240#end of PWDR processing
2241        elif 'HKLF' in histogram[:4]:
2242            Histogram = Histograms[histogram]
2243            Histogram['Residuals'] = {}
2244            phase = Histogram['Reflection Lists']
2245            Phase = Phases[phase]
2246            hId = Histogram['hId']
2247            hfx = ':%d:'%(hId)
2248            wtFactor = calcControls[hfx+'wtFactor']
2249            pfx = '%d::'%(Phase['pId'])
2250            phfx = '%d:%d:'%(Phase['pId'],hId)
2251            SGData = Phase['General']['SGData']
2252            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2253            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2254            refDict = Histogram['Data']
2255            time0 = time.time()
2256            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2257#            print 'sf-calc time: %.3f'%(time.time()-time0)
2258            df = np.zeros(len(refDict['RefList']))
2259            sumwYo = 0
2260            sumFo = 0
2261            sumFo2 = 0
2262            sumdF = 0
2263            sumdF2 = 0
2264            nobs = 0
2265            if calcControls['F**2']:
2266                for i,ref in enumerate(refDict['RefList']):
2267                    if ref[6] > 0:
2268                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2269                        w = 1.0/ref[6]
2270                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]  #correct Fc^2 for extinction
2271                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2272                        if w*ref[5] >= calcControls['minF/sig']:
2273                            Fo = np.sqrt(ref[5])
2274                            sumFo += Fo
2275                            sumFo2 += ref[5]
2276                            sumdF += abs(Fo-np.sqrt(ref[7]))
2277                            sumdF2 += abs(ref[5]-ref[7])
2278                            nobs += 1
2279                            df[i] = -w*(ref[5]-ref[7])
2280                            sumwYo += (w*ref[5])**2
2281            else:
2282                for i,ref in enumerate(refDict['RefList']):
2283                    if ref[5] > 0.:
2284                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2285                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]    #correct Fc^2 for extinction
2286                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2287                        Fo = np.sqrt(ref[5])
2288                        Fc = np.sqrt(ref[7])
2289                        w = 2.0*Fo/ref[6]
2290                        if w*Fo >= calcControls['minF/sig']:
2291                            sumFo += Fo
2292                            sumFo2 += ref[5]
2293                            sumdF += abs(Fo-Fc)
2294                            sumdF2 += abs(ref[5]-ref[7])
2295                            nobs += 1
2296                            df[i] = -w*(Fo-Fc)
2297                            sumwYo += (w*Fo)**2
2298            Histogram['Residuals']['Nobs'] = nobs
2299            Histogram['Residuals']['sumwYo'] = sumwYo
2300            SumwYo += sumwYo
2301            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2302            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2303            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2304            Histogram['Residuals'][phfx+'Nref'] = nobs
2305            Nobs += nobs
2306            if dlg:
2307                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2308            M = np.concatenate((M,wtFactor*df))
2309# end of HKLF processing
2310    Histograms['sumwYo'] = SumwYo
2311    Histograms['Nobs'] = Nobs
2312    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2313    if dlg:
2314        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2315        if not GoOn:
2316            parmDict['saved values'] = values
2317            dlg.Destroy()
2318            raise Exception         #Abort!!
2319    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2320    if len(pVals):
2321        pSum = np.sum(pWt*pVals**2)
2322        for name in pWsum:
2323            if pWsum:
2324                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2325        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2326        Nobs += len(pVals)
2327        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2328    return M
2329                       
Note: See TracBrowser for help on using the repository browser.