source: trunk/GSASIIstrMath.py @ 1562

Last change on this file since 1562 was 1562, checked in by vondreele, 8 years ago

wrong sign on Aij+Dij

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