source: trunk/GSASIIstrMath.py @ 1466

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

make extinction math use FcTsq? not Fcsq
fix 3D HKL point lighting
fix reflection listing for CW data

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 115.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2014-08-18 19:53:26 +0000 (Mon, 18 Aug 2014) $
8# $Author: vondreele $
9# $Revision: 1466 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1466 2014-08-18 19:53:26Z vondreele $
12########### SVN repository information ###################
13import time
14import math
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.optimize as so
20import scipy.stats as st
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 1466 $")
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,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,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            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1784            if 'C' in calcControls[hfx+'histType']:        #CW powder
1785                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1786            else: #'T'OF
1787                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1788            iBeg = np.searchsorted(x,refl[5]-fmin)
1789            iFin = np.searchsorted(x,refl[5]+fmax)
1790            if not iBeg+iFin:       #peak below low limit - skip peak
1791                continue
1792            elif not iBeg-iFin:     #peak above high limit - done
1793                break
1794            pos = refl[5]
1795            if 'C' in calcControls[hfx+'histType']:
1796                tanth = tand(pos/2.0)
1797                costh = cosd(pos/2.0)
1798                lenBF = iFin-iBeg
1799                dMdpk = np.zeros(shape=(6,lenBF))
1800                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1801                for i in range(5):
1802                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
1803                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1804                if Ka2:
1805                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1806                    iBeg2 = np.searchsorted(x,pos2-fmin)
1807                    iFin2 = np.searchsorted(x,pos2+fmax)
1808                    if iBeg2-iFin2:
1809                        lenBF2 = iFin2-iBeg2
1810                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1811                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1812                        for i in range(5):
1813                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
1814                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
1815                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1816            else:   #'T'OF
1817                lenBF = iFin-iBeg
1818                dMdpk = np.zeros(shape=(6,lenBF))
1819                dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1820                for i in range(6):
1821                    dMdpk[i] += refl[11]*refl[9]*dMdipk[i]      #cw[iBeg:iFin]*
1822                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
1823            if Phase['General'].get('doPawley'):
1824                dMdpw = np.zeros(len(x))
1825                try:
1826                    pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1827                    idx = varylist.index(pIdx)
1828                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1829                    if Ka2: #not for TOF either
1830                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1831                    dMdv[idx] = dMdpw
1832                except: # ValueError:
1833                    pass
1834            if 'C' in calcControls[hfx+'histType']:
1835                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1836                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1837                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1838                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1839                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1840                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1841                    hfx+'DisplaceY':[dpdY,'pos'],}
1842                if 'Bragg' in calcControls[hfx+'instType']:
1843                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1844                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
1845                else:
1846                    names.update({hfx+'Absorption':[dFdAb,'int'],})
1847            else:   #'T'OF
1848                dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1849                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1850                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1851                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1852                    hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1853                    hfx+'beta-q':[1./refl[4]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'],
1854                    hfx+'sig-2':[refl[4]**4,'sig'],hfx+'sig-q':[1./refl[4]**2,'sig'],
1855                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1856            for name in names:
1857                item = names[name]
1858                if name in varylist:
1859                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1860                    if Ka2:
1861                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1862                elif name in dependentVars:
1863                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1864                    if Ka2:
1865                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1866            for iPO in dIdPO:
1867                if iPO in varylist:
1868                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1869                    if Ka2:
1870                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1871                elif iPO in dependentVars:
1872                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1873                    if Ka2:
1874                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1875            for i,name in enumerate(['omega','chi','phi']):
1876                aname = pfx+'SH '+name
1877                if aname in varylist:
1878                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1879                    if Ka2:
1880                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1881                elif aname in dependentVars:
1882                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1883                    if Ka2:
1884                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1885            for iSH in dFdODF:
1886                if iSH in varylist:
1887                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1888                    if Ka2:
1889                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1890                elif iSH in dependentVars:
1891                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1892                    if Ka2:
1893                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1894            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1895            for name,dpdA in cellDervNames:
1896                if name in varylist:
1897                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1898                    if Ka2:
1899                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1900                elif name in dependentVars:
1901                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1902                    if Ka2:
1903                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1904            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1905            for name in dDijDict:
1906                if name in varylist:
1907                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1908                    if Ka2:
1909                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1910                elif name in dependentVars:
1911                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1912                    if Ka2:
1913                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1914            if 'C' in calcControls[hfx+'histType']:
1915                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)
1916            else:   #'T'OF
1917                sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict)
1918            for name in gamDict:
1919                if name in varylist:
1920                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1921                    if Ka2:
1922                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1923                elif name in dependentVars:
1924                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1925                    if Ka2:
1926                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1927            for name in sigDict:
1928                if name in varylist:
1929                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1930                    if Ka2:
1931                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1932                elif name in dependentVars:
1933                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1934                    if Ka2:
1935                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1936            for name in ['BabA','BabU']:
1937                if refl[9]:
1938                    if phfx+name in varylist:
1939                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1940                        if Ka2:
1941                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1942                    elif phfx+name in dependentVars:                   
1943                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1944                        if Ka2:
1945                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
1946            if not Phase['General'].get('doPawley'):
1947                #do atom derivatives -  for RB,F,X & U so far             
1948                corr = dervDict['int']/refl[9]
1949                if Ka2:
1950                    corr2 = dervDict2['int']/refl[9]
1951                for name in varylist+dependentVars:
1952                    if '::RBV;' in name:
1953                        pass
1954                    else:
1955                        try:
1956                            aname = name.split(pfx)[1][:2]
1957                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1958                        except IndexError:
1959                            continue
1960                    if name in varylist:
1961                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1962                        if Ka2:
1963                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1964                    elif name in dependentVars:
1965                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1966                        if Ka2:
1967                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1968    #        print 'profile derv time: %.3fs'%(time.time()-time0)
1969    # now process derivatives in constraints
1970    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1971    return dMdv
1972   
1973def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1974    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1975    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
1976    are returned to either dervRefine or HessRefine.
1977
1978    :returns:
1979    '''
1980    nobs = Histogram['Residuals']['Nobs']
1981    hId = Histogram['hId']
1982    hfx = ':%d:'%(hId)
1983    pfx = '%d::'%(Phase['pId'])
1984    phfx = '%d:%d:'%(Phase['pId'],hId)
1985    SGData = Phase['General']['SGData']
1986    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1987    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1988    refDict = Histogram['Data']
1989    dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1990    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1991    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1992    dependentVars = G2mv.GetDependentVars()
1993    depDerivDict = {}
1994    for j in dependentVars:
1995        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1996    wdf = np.zeros(len(refDict['RefList']))
1997    if calcControls['F**2']:
1998        for iref,ref in enumerate(refDict['RefList']):
1999            if ref[6] > 0:
2000                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:] 
2001                w = 1.0/ref[6]
2002                if w*ref[5] >= calcControls['minF/sig']:
2003                    wdf[iref] = w*(ref[5]-ref[7])
2004                    for j,var in enumerate(varylist):
2005                        if var in dFdvDict:
2006                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]   #*dervCor
2007                    for var in dependentVars:
2008                        if var in dFdvDict:
2009                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]  #*dervCor
2010                    if phfx+'Scale' in varylist:
2011                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
2012                    elif phfx+'Scale' in dependentVars:
2013                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
2014                    for item in ['Ep','Es','Eg']:
2015                        if phfx+item in varylist and dervDict:
2016                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor
2017                        elif phfx+item in dependentVars and dervDict:
2018                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
2019                    for item in ['BabA','BabU']:
2020                        if phfx+item in varylist:
2021                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2022                        elif phfx+item in dependentVars:
2023                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2024    else:
2025        for iref,ref in enumerate(refDict['RefList']):
2026            if ref[5] > 0.:
2027                dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:]
2028                Fo = np.sqrt(ref[5])
2029                Fc = np.sqrt(ref[7])
2030                w = 1.0/ref[6]
2031                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2032                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2033                    for j,var in enumerate(varylist):
2034                        if var in dFdvDict:
2035                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]       #*dervCor
2036                    for var in dependentVars:
2037                        if var in dFdvDict:
2038                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]      #*dervCor
2039                    if phfx+'Scale' in varylist:
2040                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]    #*dervCor
2041                    elif phfx+'Scale' in dependentVars:
2042                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11] #*dervCor                           
2043                    for item in ['Ep','Es','Eg']:
2044                        if phfx+item in varylist and dervDict:
2045                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #correct
2046                        elif phfx+item in dependentVars and dervDict:
2047                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
2048                    for item in ['BabA','BabU']:
2049                        if phfx+item in varylist:
2050                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2051                        elif phfx+item in dependentVars:
2052                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
2053    return dMdvh,depDerivDict,wdf
2054   
2055
2056def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2057    '''Loop over histograms and compute derivatives of the fitting
2058    model (M) with respect to all parameters.  Results are returned in
2059    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2060    n is the number of parameters and m is the number of data
2061    points. This can exceed memory when m gets large. This routine is
2062    used when refinement derivatives are selected as "analtytic
2063    Jacobian" in Controls.
2064
2065    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2066    '''
2067    parmDict.update(zip(varylist,values))
2068    G2mv.Dict2Map(parmDict,varylist)
2069    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2070    nvar = len(varylist)
2071    dMdv = np.empty(0)
2072    histoList = Histograms.keys()
2073    histoList.sort()
2074    for histogram in histoList:
2075        if 'PWDR' in histogram[:4]:
2076            Histogram = Histograms[histogram]
2077            hId = Histogram['hId']
2078            hfx = ':%d:'%(hId)
2079            wtFactor = calcControls[hfx+'wtFactor']
2080            Limits = calcControls[hfx+'Limits']
2081            x,y,w,yc,yb,yd = Histogram['Data']
2082            W = wtFactor*w
2083            xB = np.searchsorted(x,Limits[0])
2084            xF = np.searchsorted(x,Limits[1])
2085            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2086                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2087        elif 'HKLF' in histogram[:4]:
2088            Histogram = Histograms[histogram]
2089            phase = Histogram['Reflection Lists']
2090            Phase = Phases[phase]
2091            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2092            hfx = ':%d:'%(Histogram['hId'])
2093            wtFactor = calcControls[hfx+'wtFactor']
2094            # now process derivatives in constraints
2095            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2096        else:
2097            continue        #skip non-histogram entries
2098        if len(dMdv):
2099            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2100        else:
2101            dMdv = np.sqrt(wtFactor)*dMdvh
2102           
2103    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2104    if np.any(pVals):
2105        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2106        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2107       
2108    return dMdv
2109
2110def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2111    '''Loop over histograms and compute derivatives of the fitting
2112    model (M) with respect to all parameters.  For each histogram, the
2113    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2114    number of parameters and m is the number of data points *in the
2115    histogram*. The (n by n) Hessian is computed from each Jacobian
2116    and it is returned.  This routine is used when refinement
2117    derivatives are selected as "analtytic Hessian" in Controls.
2118
2119    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2120    '''
2121    parmDict.update(zip(varylist,values))
2122    G2mv.Dict2Map(parmDict,varylist)
2123    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2124    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2125    nvar = len(varylist)
2126    Hess = np.empty(0)
2127    histoList = Histograms.keys()
2128    histoList.sort()
2129    for histogram in histoList:
2130        if 'PWDR' in histogram[:4]:
2131            Histogram = Histograms[histogram]
2132            hId = Histogram['hId']
2133            hfx = ':%d:'%(hId)
2134            wtFactor = calcControls[hfx+'wtFactor']
2135            Limits = calcControls[hfx+'Limits']
2136            x,y,w,yc,yb,yd = Histogram['Data']
2137            W = wtFactor*w
2138            dy = y-yc
2139            xB = np.searchsorted(x,Limits[0])
2140            xF = np.searchsorted(x,Limits[1])
2141            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2142                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2143            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2144            Dy = dy[xB:xF][np.newaxis,:]
2145            dMdvh *= Wt
2146            if dlg:
2147                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2148            if len(Hess):
2149                Hess += np.inner(dMdvh,dMdvh)
2150                dMdvh *= Wt*Dy
2151                Vec += np.sum(dMdvh,axis=1)
2152            else:
2153                Hess = np.inner(dMdvh,dMdvh)
2154                dMdvh *= Wt*Dy
2155                Vec = np.sum(dMdvh,axis=1)
2156        elif 'HKLF' in histogram[:4]:
2157            Histogram = Histograms[histogram]
2158            phase = Histogram['Reflection Lists']
2159            Phase = Phases[phase]
2160            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2161            hId = Histogram['hId']
2162            hfx = ':%d:'%(Histogram['hId'])
2163            wtFactor = calcControls[hfx+'wtFactor']
2164            # now process derivatives in constraints
2165            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2166#            print 'matrix build time: %.3f'%(time.time()-time0)
2167
2168            if dlg:
2169                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2170            if len(Hess):
2171                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2172                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2173            else:
2174                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2175                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2176        else:
2177            continue        #skip non-histogram entries
2178    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2179    if np.any(pVals):
2180        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2181        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2182        Hess += np.inner(dpdv*pWt,dpdv)
2183    return Vec,Hess
2184
2185def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
2186    'Needs a doc string'
2187    Values2Dict(parmDict, varylist, values)
2188    G2mv.Dict2Map(parmDict,varylist)
2189    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2190    M = np.empty(0)
2191    SumwYo = 0
2192    Nobs = 0
2193    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2194    histoList = Histograms.keys()
2195    histoList.sort()
2196    for histogram in histoList:
2197        if 'PWDR' in histogram[:4]:
2198            Histogram = Histograms[histogram]
2199            hId = Histogram['hId']
2200            hfx = ':%d:'%(hId)
2201            wtFactor = calcControls[hfx+'wtFactor']
2202            Limits = calcControls[hfx+'Limits']
2203            x,y,w,yc,yb,yd = Histogram['Data']
2204            yc *= 0.0                           #zero full calcd profiles
2205            yb *= 0.0
2206            yd *= 0.0
2207            xB = np.searchsorted(x,Limits[0])
2208            xF = np.searchsorted(x,Limits[1])
2209            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2210                varylist,Histogram,Phases,calcControls,pawleyLookup)
2211            yc[xB:xF] += yb[xB:xF]
2212            if not np.any(y):                   #fill dummy data
2213                rv = st.poisson(yc[xB:xF])
2214                y[xB:xF] = rv.rvs()
2215                Z = np.ones_like(yc[xB:xF])
2216                Z[1::2] *= -1
2217                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2218                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2219            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2220            W = wtFactor*w
2221            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2222            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2223            Nobs += Histogram['Residuals']['Nobs']
2224            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2225            SumwYo += Histogram['Residuals']['sumwYo']
2226            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2227            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2228            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2229            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2230            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.))
2231            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.))
2232            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2233            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2234            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2235            if dlg:
2236                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2237            M = np.concatenate((M,wdy))
2238#end of PWDR processing
2239        elif 'HKLF' in histogram[:4]:
2240            Histogram = Histograms[histogram]
2241            Histogram['Residuals'] = {}
2242            phase = Histogram['Reflection Lists']
2243            Phase = Phases[phase]
2244            hId = Histogram['hId']
2245            hfx = ':%d:'%(hId)
2246            wtFactor = calcControls[hfx+'wtFactor']
2247            pfx = '%d::'%(Phase['pId'])
2248            phfx = '%d:%d:'%(Phase['pId'],hId)
2249            SGData = Phase['General']['SGData']
2250            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2251            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2252            refDict = Histogram['Data']
2253            time0 = time.time()
2254            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2255#            print 'sf-calc time: %.3f'%(time.time()-time0)
2256            df = np.zeros(len(refDict['RefList']))
2257            sumwYo = 0
2258            sumFo = 0
2259            sumFo2 = 0
2260            sumdF = 0
2261            sumdF2 = 0
2262            nobs = 0
2263            if calcControls['F**2']:
2264                for i,ref in enumerate(refDict['RefList']):
2265                    if ref[6] > 0:
2266                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2267                        w = 1.0/ref[6]
2268                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]  #correct Fc^2 for extinction
2269                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2270                        if w*ref[5] >= calcControls['minF/sig']:
2271                            Fo = np.sqrt(ref[5])
2272                            sumFo += Fo
2273                            sumFo2 += ref[5]
2274                            sumdF += abs(Fo-np.sqrt(ref[7]))
2275                            sumdF2 += abs(ref[5]-ref[7])
2276                            nobs += 1
2277                            df[i] = -w*(ref[5]-ref[7])
2278                            sumwYo += (w*ref[5])**2
2279            else:
2280                for i,ref in enumerate(refDict['RefList']):
2281                    if ref[5] > 0.:
2282                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2283                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]    #correct Fc^2 for extinction
2284                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2285                        Fo = np.sqrt(ref[5])
2286                        Fc = np.sqrt(ref[7])
2287                        w = 2.0*Fo/ref[6]
2288                        if w*Fo >= calcControls['minF/sig']:
2289                            sumFo += Fo
2290                            sumFo2 += ref[5]
2291                            sumdF += abs(Fo-Fc)
2292                            sumdF2 += abs(ref[5]-ref[7])
2293                            nobs += 1
2294                            df[i] = -w*(Fo-Fc)
2295                            sumwYo += (w*Fo)**2
2296            Histogram['Residuals']['Nobs'] = nobs
2297            Histogram['Residuals']['sumwYo'] = sumwYo
2298            SumwYo += sumwYo
2299            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2300            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2301            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2302            Histogram['Residuals'][phfx+'Nref'] = nobs
2303            Nobs += nobs
2304            if dlg:
2305                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2306            M = np.concatenate((M,wtFactor*df))
2307# end of HKLF processing
2308    Histograms['sumwYo'] = SumwYo
2309    Histograms['Nobs'] = Nobs
2310    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2311    if dlg:
2312        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2313        if not GoOn:
2314            parmDict['saved values'] = values
2315            dlg.Destroy()
2316            raise Exception         #Abort!!
2317    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2318    if len(pVals):
2319        pSum = np.sum(pWt*pVals**2)
2320        for name in pWsum:
2321            if pWsum:
2322                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2323        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2324        Nobs += len(pVals)
2325        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2326    return M
2327                       
Note: See TracBrowser for help on using the repository browser.