source: trunk/GSASIIstrMath.py @ 1875

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

add Flack parameter to GUI (commented out for now)

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