source: trunk/GSASIIstrMath.py @ 1874

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

minor formatting change for TOF peak fitting
fix powder extinction & derivatives

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 140.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-06-01 16:31:55 +0000 (Mon, 01 Jun 2015) $
8# $Author: vondreele $
9# $Revision: 1874 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1874 2015-06-01 16:31:55Z 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: 1874 $")
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        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
859        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
860        Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
861        Phi = np.inner(H.T,SGT).flatten()
862        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
863        sinp = np.sin(phase)
864        cosp = np.cos(phase)
865        biso = -SQfactor*Uisodata[:,np.newaxis]
866        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
867        HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
868        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
869        Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
870        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
871        fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))   #real A,-b
872        fas = np.sum(np.sum(fa,axis=2),axis=2)        #real sum over atoms & unique hkl
873        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
874        fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))   #imag -B,+a
875        fbs = np.sum(np.sum(fb,axis=2),axis=2)  #imag sum over atoms & uniq hkl
876        if SGData['SGInv']: #centrosymmetric; B=0
877            fbs[0] *= 0.
878        if 'P' in calcControls[hfx+'histType']:
879            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
880        else:
881            refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
882        refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"?
883#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
884        iBeg += blkSize
885   
886def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
887    'Needs a doc string'
888    phfx = pfx.split(':')[0]+hfx
889    ast = np.sqrt(np.diag(G))
890    Mast = twopisq*np.multiply.outer(ast,ast)
891    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
892    SGT = np.array([ops[1] for ops in SGData['SGOps']])
893    FFtables = calcControls['FFtables']
894    BLtables = calcControls['BLtables']
895    nRef = len(refDict['RefList'])
896    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
897    mSize = len(Mdata)
898    FF = np.zeros(len(Tdata))
899    if 'NC' in calcControls[hfx+'histType']:
900        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
901    elif 'X' in calcControls[hfx+'histType']:
902        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
903        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
904    Uij = np.array(G2lat.U6toUij(Uijdata))
905    bij = Mast*Uij.T
906    dFdvDict = {}
907    dFdfr = np.zeros((nRef,mSize))
908    dFdx = np.zeros((nRef,mSize,3))
909    dFdui = np.zeros((nRef,mSize))
910    dFdua = np.zeros((nRef,mSize,6))
911    dFdbab = np.zeros((nRef,2))
912    for iref,refl in enumerate(refDict['RefList']):
913        if 'T' in calcControls[hfx+'histType']:
914            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
915        H = np.array(refl[:3])
916        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
917        SQfactor = 8.0*SQ*np.pi**2
918        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
919        Bab = parmDict[phfx+'BabA']*dBabdA
920        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
921        FF = refDict['FF']['FF'][iref].T[Tindx]
922        Uniq = np.inner(H,SGMT)
923        Phi = np.inner(H,SGT)
924        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
925        sinp = np.sin(phase)
926        cosp = np.cos(phase)
927        occ = Mdata*Fdata/len(Uniq)
928        biso = -SQfactor*Uisodata
929        Tiso = np.where(biso<1.,np.exp(biso),1.0)
930        HbH = -np.inner(H,np.inner(bij,H))
931        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
932        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
933        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
934        Tcorr = Tiso*Tuij
935        fot = (FF+FP-Bab)*occ*Tcorr
936        fotp = FPP*occ*Tcorr
937        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
938        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
939       
940        fas = np.sum(np.sum(fa,axis=1),axis=1)      #real sum over atoms & unique hkl
941        fbs = np.sum(np.sum(fb,axis=1),axis=1)      #imag sum over atoms & uniq hkl
942        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
943        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
944        #sum below is over Uniq
945        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
946        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
947        dfadui = np.sum(-SQfactor*fa,axis=2)
948        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
949        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
950        if not SGData['SGInv']:
951            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
952            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
953            dfbdui = np.sum(-SQfactor*fb,axis=2)
954            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
955            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
956        else:
957            dfbdfr = np.zeros_like(dfadfr)
958            dfbdx = np.zeros_like(dfadx)
959            dfbdui = np.zeros_like(dfadui)
960            dfbdua = np.zeros_like(dfadua)
961            dfbdba = np.zeros_like(dfadba)
962        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
963        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
964            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
965                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
966            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
967                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
968            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
969                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
970            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
971                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
972        else:
973            SA = fas[0]-fbs[1]
974            SB = fbs[0]+fas[1]
975            dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
976                2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
977            dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])
978            dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])
979            dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])
980        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
981            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
982        #loop over atoms - each dict entry is list of derivatives for all the reflections
983    for i in range(len(Mdata)):
984        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
985        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
986        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
987        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
988        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
989        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
990        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
991        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
992        dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
993        dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
994        dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
995    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
996    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
997    return dFdvDict
998   
999def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1000    'Needs a doc string'
1001    phfx = pfx.split(':')[0]+hfx
1002    ast = np.sqrt(np.diag(G))
1003    Mast = twopisq*np.multiply.outer(ast,ast)
1004    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1005    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1006    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1007    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1008    FFtables = calcControls['FFtables']
1009    BLtables = calcControls['BLtables']
1010    nRef = len(refDict['RefList'])
1011    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1012    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1013    mSize = len(Mdata)
1014    FF = np.zeros(len(Tdata))
1015    if 'NC' in calcControls[hfx+'histType']:
1016        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1017    elif 'X' in calcControls[hfx+'histType']:
1018        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1019        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1020    Uij = np.array(G2lat.U6toUij(Uijdata))
1021    bij = Mast*Uij.T
1022    dFdvDict = {}
1023    dFdfr = np.zeros((nRef,mSize))
1024    dFdx = np.zeros((nRef,mSize,3))
1025    dFdui = np.zeros((nRef,mSize))
1026    dFdua = np.zeros((nRef,mSize,6))
1027    dFdbab = np.zeros((nRef,2))
1028    for iref,refl in enumerate(refDict['RefList']):
1029        if 'T' in calcControls[hfx+'histType']:
1030            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1031        H = np.array(refl[:4])
1032        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1033        SQfactor = 8.0*SQ*np.pi**2
1034        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1035        Bab = parmDict[phfx+'BabA']*dBabdA
1036        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1037        FF = refDict['FF']['FF'][iref].T[Tindx]
1038        Uniq = np.inner(H[:3],SGMT)
1039        SSUniq = np.inner(H,SSGMT)
1040        Phi = np.inner(H[:3],SGT)
1041        SSPhi = np.inner(H,SSGT)
1042        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1043        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
1044        #need ModulationDerv here dGAdXsin, etc 
1045        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
1046        sinp = np.sin(phase)
1047        cosp = np.cos(phase)
1048        occ = Mdata*Fdata/len(Uniq)
1049        biso = -SQfactor*Uisodata
1050        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1051        HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
1052        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1053        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1054        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1055        Tcorr = Tiso*Tuij
1056        fot = (FF+FP-Bab)*occ*Tcorr
1057        fotp = FPP*occ*Tcorr
1058        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1059        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1060        GfpuA = np.swapaxes(GfpuA,1,2)
1061        GfpuB = np.swapaxes(GfpuB,1,2)
1062        fa = fa*GfpuA-fb*GfpuB
1063        fb = fb*GfpuA+fa*GfpuB
1064       
1065        fas = np.sum(np.sum(fa,axis=1),axis=1)
1066        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1067        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1068        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1069        fax = fax*GfpuA-fbx*GfpuB
1070        fbx = fbx*GfpuA+fax*GfpuB
1071        #sum below is over Uniq
1072        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1073        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1074        dfadui = np.sum(-SQfactor*fa,axis=2)
1075        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1076        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
1077        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
1078        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1079        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1080        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1081        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1082        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1083        #need dFdXsin, etc for modulations...
1084        if not SGData['SGInv']:
1085            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
1086            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
1087            dfbdui = np.sum(-SQfactor*fb,axis=2)
1088            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1089            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
1090            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1091            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1092            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1093            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1094            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1095        #need dFdXsin, etc for modulations...
1096        #loop over atoms - each dict entry is list of derivatives for all the reflections
1097    for i in range(len(Mdata)):     
1098        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1099        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1100        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1101        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1102        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1103        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1104        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1105        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1106        dFdvDict[pfx+'AU12:'+str(i)] = .5*dFdua.T[3][i]
1107        dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
1108        dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
1109        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1110    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
1111    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
1112    return dFdvDict
1113   
1114def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1115    ''' Single crystal extinction function; returns extinction & derivative
1116    '''
1117    extCor = 1.0
1118    dervDict = {}
1119    dervCor = 1.0
1120    if calcControls[phfx+'EType'] != 'None':
1121        SQ = 1/(4.*ref[4+im]**2)
1122        if 'C' in parmDict[hfx+'Type']:           
1123            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1124        else:   #'T'
1125            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1126        if 'SXC' in parmDict[hfx+'Type']:
1127            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1128            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1129            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1130            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1131        elif 'SNT' in parmDict[hfx+'Type']:
1132            AV = 1.e7/parmDict[pfx+'Vol']**2
1133            PL = SQ
1134            PLZ = AV*ref[9+im]*ref[12+im]**2
1135        elif 'SNC' in parmDict[hfx+'Type']:
1136            AV = 1.e7/parmDict[pfx+'Vol']**2
1137            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1138            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1139           
1140        if 'Primary' in calcControls[phfx+'EType']:
1141            PLZ *= 1.5
1142        else:
1143            if 'C' in parmDict[hfx+'Type']:
1144                PLZ *= calcControls[phfx+'Tbar']
1145            else: #'T'
1146                PLZ *= ref[13+im]      #t-bar
1147        if 'Primary' in calcControls[phfx+'EType']:
1148            PLZ *= 1.5
1149            PSIG = parmDict[phfx+'Ep']
1150        elif 'I & II' in calcControls[phfx+'EType']:
1151            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1152        elif 'Type II' in calcControls[phfx+'EType']:
1153            PSIG = parmDict[phfx+'Es']
1154        else:       # 'Secondary Type I'
1155            PSIG = parmDict[phfx+'Eg']/PL
1156           
1157        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1158        AL = 0.025+0.285*cos2T
1159        BG = 0.02-0.025*cos2T
1160        BL = 0.15-0.2*(0.75-cos2T)**2
1161        if cos2T < 0.:
1162            BL = -0.45*cos2T
1163        CG = 2.
1164        CL = 2.
1165        PF = PLZ*PSIG
1166       
1167        if 'Gaussian' in calcControls[phfx+'EApprox']:
1168            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1169            extCor = np.sqrt(PF4)
1170            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1171        else:
1172            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1173            extCor = np.sqrt(PF4)
1174            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1175
1176        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
1177        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1178            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1179        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1180            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1181        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1182            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1183               
1184    return 1./extCor,dervDict,dervCor
1185   
1186def Dict2Values(parmdict, varylist):
1187    '''Use before call to leastsq to setup list of values for the parameters
1188    in parmdict, as selected by key in varylist'''
1189    return [parmdict[key] for key in varylist] 
1190   
1191def Values2Dict(parmdict, varylist, values):
1192    ''' Use after call to leastsq to update the parameter dictionary with
1193    values corresponding to keys in varylist'''
1194    parmdict.update(zip(varylist,values))
1195   
1196def GetNewCellParms(parmDict,varyList):
1197    'Needs a doc string'
1198    newCellDict = {}
1199    Anames = ['A'+str(i) for i in range(6)]
1200    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1201    for item in varyList:
1202        keys = item.split(':')
1203        if keys[2] in Ddict:
1204            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1205            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1206            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1207    return newCellDict          # is e.g. {'0::D11':A0-D11}
1208   
1209def ApplyXYZshifts(parmDict,varyList):
1210    '''
1211    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1212   
1213    :param dict parmDict: parameter dictionary
1214    :param list varyList: list of variables (not used!)
1215    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1216
1217    '''
1218    newAtomDict = {}
1219    for item in parmDict:
1220        if 'dA' in item:
1221            parm = ''.join(item.split('d'))
1222            parmDict[parm] += parmDict[item]
1223            newAtomDict[item] = [parm,parmDict[parm]]
1224    return newAtomDict
1225   
1226def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1227    'Spherical harmonics texture'
1228    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1229    if 'T' in calcControls[hfx+'histType']:
1230        tth = parmDict[hfx+'2-theta']
1231    else:
1232        tth = refl[5+im]
1233    odfCor = 1.0
1234    H = refl[:3]
1235    cell = G2lat.Gmat2cell(g)
1236    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1237    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1238    phi,beta = G2lat.CrsAng(H,cell,SGData)
1239    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1240    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1241    for item in SHnames:
1242        L,M,N = eval(item.strip('C'))
1243        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1244        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1245        Lnorm = G2lat.Lnorm(L)
1246        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1247    return odfCor
1248   
1249def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1250    'Spherical harmonics texture derivatives'
1251    if 'T' in calcControls[hfx+'histType']:
1252        tth = parmDict[hfx+'2-theta']
1253    else:
1254        tth = refl[5+im]
1255    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1256    odfCor = 1.0
1257    dFdODF = {}
1258    dFdSA = [0,0,0]
1259    H = refl[:3]
1260    cell = G2lat.Gmat2cell(g)
1261    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1262    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1263    phi,beta = G2lat.CrsAng(H,cell,SGData)
1264    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1265    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1266    for item in SHnames:
1267        L,M,N = eval(item.strip('C'))
1268        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1269        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1270        Lnorm = G2lat.Lnorm(L)
1271        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1272        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1273        for i in range(3):
1274            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1275    return odfCor,dFdODF,dFdSA
1276   
1277def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1278    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1279    if 'T' in calcControls[hfx+'histType']:
1280        tth = parmDict[hfx+'2-theta']
1281    else:
1282        tth = refl[5+im]
1283    odfCor = 1.0
1284    H = refl[:3]
1285    cell = G2lat.Gmat2cell(g)
1286    Sangls = [0.,0.,0.]
1287    if 'Bragg' in calcControls[hfx+'instType']:
1288        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1289        IFCoup = True
1290    else:
1291        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1292        IFCoup = False
1293    phi,beta = G2lat.CrsAng(H,cell,SGData)
1294    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1295    SHnames = calcControls[phfx+'SHnames']
1296    for item in SHnames:
1297        L,N = eval(item.strip('C'))
1298        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1299        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1300        Lnorm = G2lat.Lnorm(L)
1301        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1302    return np.squeeze(odfCor)
1303   
1304def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1305    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1306    if 'T' in calcControls[hfx+'histType']:
1307        tth = parmDict[hfx+'2-theta']
1308    else:
1309        tth = refl[5+im]
1310    odfCor = 1.0
1311    dFdODF = {}
1312    H = refl[:3]
1313    cell = G2lat.Gmat2cell(g)
1314    Sangls = [0.,0.,0.]
1315    if 'Bragg' in calcControls[hfx+'instType']:
1316        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1317        IFCoup = True
1318    else:
1319        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1320        IFCoup = False
1321    phi,beta = G2lat.CrsAng(H,cell,SGData)
1322    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1323    SHnames = calcControls[phfx+'SHnames']
1324    for item in SHnames:
1325        L,N = eval(item.strip('C'))
1326        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1327        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1328        Lnorm = G2lat.Lnorm(L)
1329        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1330        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
1331    return odfCor,dFdODF
1332   
1333def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1334    'March-Dollase preferred orientation correction'
1335    POcorr = 1.0
1336    MD = parmDict[phfx+'MD']
1337    if MD != 1.0:
1338        MDAxis = calcControls[phfx+'MDAxis']
1339        sumMD = 0
1340        for H in uniq:           
1341            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1342            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1343            sumMD += A**3
1344        POcorr = sumMD/len(uniq)
1345    return POcorr
1346   
1347def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1348    'Needs a doc string'
1349    POcorr = 1.0
1350    POderv = {}
1351    if calcControls[phfx+'poType'] == 'MD':
1352        MD = parmDict[phfx+'MD']
1353        MDAxis = calcControls[phfx+'MDAxis']
1354        sumMD = 0
1355        sumdMD = 0
1356        for H in uniq:           
1357            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1358            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1359            sumMD += A**3
1360            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1361        POcorr = sumMD/len(uniq)
1362        POderv[phfx+'MD'] = sumdMD/len(uniq)
1363    else:   #spherical harmonics
1364        if calcControls[phfx+'SHord']:
1365            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1366    return POcorr,POderv
1367   
1368def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1369    'Needs a doc string'
1370    if 'Debye' in calcControls[hfx+'instType']:
1371        if 'T' in calcControls[hfx+'histType']:
1372            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1373        else:
1374            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1375    else:
1376        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1377   
1378def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1379    'Needs a doc string'
1380    if 'Debye' in calcControls[hfx+'instType']:
1381        if 'T' in calcControls[hfx+'histType']:
1382            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1383        else:
1384            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1385    else:
1386        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1387       
1388def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1389    'Needs a doc string'
1390    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1391    pi2 = np.sqrt(2./np.pi)
1392    if 'T' in calcControls[hfx+'histType']:
1393        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1394        wave = refl[14+im]
1395    else:   #'C'W
1396        sth2 = sind(refl[5+im]/2.)**2
1397        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1398    c2th = 1.-2.0*sth2
1399    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1400    if 'X' in calcControls[hfx+'histType']:
1401        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1402    xfac = flv2*parmDict[phfx+'Extinction']
1403    exb = 1.0
1404    if xfac > -1.:
1405        exb = 1./np.sqrt(1.+xfac)
1406    exl = 1.0
1407    if 0 < xfac <= 1.:
1408        xn = np.array([xfac**(i+1) for i in range(6)])
1409        exl += np.sum(xn*coef)
1410    elif xfac > 1.:
1411        xfac2 = 1./np.sqrt(xfac)
1412        exl = pi2*(1.-0.125/xfac)*xfac2
1413    return exb*sth2+exl*(1.-sth2)
1414   
1415def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1416    'Needs a doc string'
1417    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1418    pi2 = np.sqrt(2./np.pi)
1419    if 'T' in calcControls[hfx+'histType']:
1420        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1421        wave = refl[14+im]
1422    else:   #'C'W
1423        sth2 = sind(refl[5+im]/2.)**2
1424        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1425    c2th = 1.-2.0*sth2
1426    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1427    if 'X' in calcControls[hfx+'histType']:
1428        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1429    xfac = flv2*parmDict[phfx+'Extinction']
1430    dbde = -500.*flv2
1431    if xfac > -1.:
1432        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
1433    dlde = 0.
1434    if 0 < xfac <= 1.:
1435        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1436        dlde = np.sum(xn*coef)
1437    elif xfac > 1.:
1438        xfac2 = 1./np.sqrt(xfac)
1439        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1440       
1441    return dbde*sth2+dlde*(1.-sth2)
1442   
1443def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1444    'Needs a doc string'    #need powder extinction!
1445    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1446    if 'X' in parmDict[hfx+'Type']:
1447        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1448    POcorr = 1.0
1449    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1450        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1451    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1452        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1453    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1454        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1455    Icorr *= POcorr
1456    AbsCorr = 1.0
1457    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1458    Icorr *= AbsCorr
1459    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1460    Icorr *= ExtCorr
1461    return Icorr,POcorr,AbsCorr,ExtCorr
1462   
1463def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1464    'Needs a doc string'    #need powder extinction derivs!
1465    dIdsh = 1./parmDict[hfx+'Scale']
1466    dIdsp = 1./parmDict[phfx+'Scale']
1467    if 'X' in parmDict[hfx+'Type']:
1468        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1469        dIdPola /= pola
1470    else:       #'N'
1471        dIdPola = 0.0
1472    dFdODF = {}
1473    dFdSA = [0,0,0]
1474    dIdPO = {}
1475    if pfx+'SHorder' in parmDict:
1476        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1477        for iSH in dFdODF:
1478            dFdODF[iSH] /= odfCor
1479        for i in range(3):
1480            dFdSA[i] /= odfCor
1481    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1482        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1483        for iPO in dIdPO:
1484            dIdPO[iPO] /= POcorr
1485    if 'T' in parmDict[hfx+'Type']:
1486        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1487        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1488    else:
1489        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1490        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1491    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1492       
1493def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1494    'Needs a doc string'
1495    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1496        costh = cosd(refl[5+im]/2.)
1497        #crystallite size
1498        if calcControls[phfx+'SizeType'] == 'isotropic':
1499            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1500        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1501            H = np.array(refl[:3])
1502            P = np.array(calcControls[phfx+'SizeAxis'])
1503            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1504            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1505            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1506        else:           #ellipsoidal crystallites
1507            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1508            H = np.array(refl[:3])
1509            lenR = G2pwd.ellipseSize(H,Sij,GB)
1510            Sgam = 1.8*wave/(np.pi*costh*lenR)
1511        #microstrain               
1512        if calcControls[phfx+'MustrainType'] == 'isotropic':
1513            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1514        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1515            H = np.array(refl[:3])
1516            P = np.array(calcControls[phfx+'MustrainAxis'])
1517            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1518            Si = parmDict[phfx+'Mustrain;i']
1519            Sa = parmDict[phfx+'Mustrain;a']
1520            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1521        else:       #generalized - P.W. Stephens model
1522            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1523            Sum = 0
1524            for i,strm in enumerate(Strms):
1525                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1526            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1527    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1528        #crystallite size
1529        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1530            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1531        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1532            H = np.array(refl[:3])
1533            P = np.array(calcControls[phfx+'SizeAxis'])
1534            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1535            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1536            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1537        else:           #ellipsoidal crystallites   #OK
1538            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1539            H = np.array(refl[:3])
1540            lenR = G2pwd.ellipseSize(H,Sij,GB)
1541            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1542        #microstrain               
1543        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1544            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1545        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1546            H = np.array(refl[:3])
1547            P = np.array(calcControls[phfx+'MustrainAxis'])
1548            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1549            Si = parmDict[phfx+'Mustrain;i']
1550            Sa = parmDict[phfx+'Mustrain;a']
1551            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1552        else:       #generalized - P.W. Stephens model  OK
1553            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1554            Sum = 0
1555            for i,strm in enumerate(Strms):
1556                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1557            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1558           
1559    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1560    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1561    sig /= ateln2
1562    return sig,gam
1563       
1564def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1565    'Needs a doc string'
1566    gamDict = {}
1567    sigDict = {}
1568    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1569        costh = cosd(refl[5+im]/2.)
1570        tanth = tand(refl[5+im]/2.)
1571        #crystallite size derivatives
1572        if calcControls[phfx+'SizeType'] == 'isotropic':
1573            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1574            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1575            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1576        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1577            H = np.array(refl[:3])
1578            P = np.array(calcControls[phfx+'SizeAxis'])
1579            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1580            Si = parmDict[phfx+'Size;i']
1581            Sa = parmDict[phfx+'Size;a']
1582            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1583            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1584            Sgam = gami*sqtrm
1585            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1586            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1587            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1588            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1589            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1590            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1591        else:           #ellipsoidal crystallites
1592            const = 1.8*wave/(np.pi*costh)
1593            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1594            H = np.array(refl[:3])
1595            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1596            Sgam = const/lenR
1597            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1598                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1599                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1600        gamDict[phfx+'Size;mx'] = Sgam
1601        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1602               
1603        #microstrain derivatives               
1604        if calcControls[phfx+'MustrainType'] == 'isotropic':
1605            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1606            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1607            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1608        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1609            H = np.array(refl[:3])
1610            P = np.array(calcControls[phfx+'MustrainAxis'])
1611            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1612            Si = parmDict[phfx+'Mustrain;i']
1613            Sa = parmDict[phfx+'Mustrain;a']
1614            gami = 0.018*Si*Sa*tanth/np.pi
1615            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1616            Mgam = gami/sqtrm
1617            dsi = -gami*Si*cosP**2/sqtrm**3
1618            dsa = -gami*Sa*sinP**2/sqtrm**3
1619            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1620            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1621            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1622            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1623        else:       #generalized - P.W. Stephens model
1624            const = 0.018*refl[4+im]**2*tanth/np.pi
1625            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1626            Sum = 0
1627            for i,strm in enumerate(Strms):
1628                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1629                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1630                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1631            Mgam = const*np.sqrt(Sum)
1632            for i in range(len(Strms)):
1633                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1634                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1635        gamDict[phfx+'Mustrain;mx'] = Mgam
1636        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1637    else:   #'T'OF - All checked & OK
1638        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1639            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1640            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1641            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1642        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1643            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1644            H = np.array(refl[:3])
1645            P = np.array(calcControls[phfx+'SizeAxis'])
1646            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1647            Si = parmDict[phfx+'Size;i']
1648            Sa = parmDict[phfx+'Size;a']
1649            gami = const/(Si*Sa)
1650            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1651            Sgam = gami*sqtrm
1652            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1653            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1654            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1655            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1656            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1657            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1658        else:           #OK  ellipsoidal crystallites
1659            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1660            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1661            H = np.array(refl[:3])
1662            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1663            Sgam = const/lenR
1664            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1665                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1666                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1667        gamDict[phfx+'Size;mx'] = Sgam  #OK
1668        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1669               
1670        #microstrain derivatives               
1671        if calcControls[phfx+'MustrainType'] == 'isotropic':
1672            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1673            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1674            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1675        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1676            H = np.array(refl[:3])
1677            P = np.array(calcControls[phfx+'MustrainAxis'])
1678            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1679            Si = parmDict[phfx+'Mustrain;i']
1680            Sa = parmDict[phfx+'Mustrain;a']
1681            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1682            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1683            Mgam = gami/sqtrm
1684            dsi = -gami*Si*cosP**2/sqtrm**3
1685            dsa = -gami*Sa*sinP**2/sqtrm**3
1686            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1687            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1688            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1689            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1690        else:       #generalized - P.W. Stephens model OK
1691            pwrs = calcControls[phfx+'MuPwrs']
1692            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1693            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1694            Sum = 0
1695            for i,strm in enumerate(Strms):
1696                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1697                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1698                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1699            Mgam = const*np.sqrt(Sum)
1700            for i in range(len(Strms)):
1701                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1702                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1703        gamDict[phfx+'Mustrain;mx'] = Mgam
1704        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1705       
1706    return sigDict,gamDict
1707       
1708def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1709    'Needs a doc string'
1710    if im:
1711        h,k,l,m = refl[:4]
1712        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1713        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1714    else:
1715        h,k,l = refl[:3]
1716        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1717    refl[4+im] = d
1718    if 'C' in calcControls[hfx+'histType']:
1719        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1720        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1721        if 'Bragg' in calcControls[hfx+'instType']:
1722            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1723                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1724        else:               #Debye-Scherrer - simple but maybe not right
1725            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1726    elif 'T' in calcControls[hfx+'histType']:
1727        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1728        #do I need sample position effects - maybe?
1729    return pos
1730
1731def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1732    'Needs a doc string'
1733    dpr = 180./np.pi
1734    if im:
1735        h,k,l,m = refl[:4]
1736        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1737        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1738    else:
1739        m = 0
1740        h,k,l = refl[:3]       
1741        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1742    dst = np.sqrt(dstsq)
1743    dsp = 1./dst
1744    if 'C' in calcControls[hfx+'histType']:
1745        pos = refl[5+im]-parmDict[hfx+'Zero']
1746        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1747        dpdw = const*dst
1748        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1749        dpdZ = 1.0
1750        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1751            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1752        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1753        if 'Bragg' in calcControls[hfx+'instType']:
1754            dpdSh = -4.*shft*cosd(pos/2.0)
1755            dpdTr = -shft*sind(pos)*100.0
1756            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1757        else:               #Debye-Scherrer - simple but maybe not right
1758            dpdXd = -shft*cosd(pos)
1759            dpdYd = -shft*sind(pos)
1760            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1761    elif 'T' in calcControls[hfx+'histType']:
1762        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1763        dpdZ = 1.0
1764        dpdDC = dsp
1765        dpdDA = dsp**2
1766        dpdDB = 1./dsp
1767        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1768            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1769        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1770           
1771def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1772    'Needs a doc string'
1773    laue = SGData['SGLaue']
1774    uniq = SGData['SGUniq']
1775    h,k,l = refl[:3]
1776    if laue in ['m3','m3m']:
1777        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1778            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1779    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1780        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1781    elif laue in ['3R','3mR']:
1782        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1783    elif laue in ['4/m','4/mmm']:
1784        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1785    elif laue in ['mmm']:
1786        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1787    elif laue in ['2/m']:
1788        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1789        if uniq == 'a':
1790            Dij += parmDict[phfx+'D23']*k*l
1791        elif uniq == 'b':
1792            Dij += parmDict[phfx+'D13']*h*l
1793        elif uniq == 'c':
1794            Dij += parmDict[phfx+'D12']*h*k
1795    else:
1796        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1797            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1798    if 'C' in calcControls[hfx+'histType']:
1799        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1800    else:
1801        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1802           
1803def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1804    'Needs a doc string'
1805    laue = SGData['SGLaue']
1806    uniq = SGData['SGUniq']
1807    h,k,l = refl[:3]
1808    if laue in ['m3','m3m']:
1809        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1810            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1811    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1812        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1813    elif laue in ['3R','3mR']:
1814        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1815    elif laue in ['4/m','4/mmm']:
1816        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1817    elif laue in ['mmm']:
1818        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1819    elif laue in ['2/m']:
1820        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1821        if uniq == 'a':
1822            dDijDict[phfx+'D23'] = k*l
1823        elif uniq == 'b':
1824            dDijDict[phfx+'D13'] = h*l
1825        elif uniq == 'c':
1826            dDijDict[phfx+'D12'] = h*k
1827    else:
1828        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1829            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1830    if 'C' in calcControls[hfx+'histType']:
1831        for item in dDijDict:
1832            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1833    else:
1834        for item in dDijDict:
1835            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
1836    return dDijDict
1837   
1838def GetDij(phfx,SGData,parmDict):
1839    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
1840    return G2spc.HStrainVals(HSvals,SGData)
1841               
1842def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1843    'Needs a doc string'
1844    histoList = Histograms.keys()
1845    histoList.sort()
1846    for histogram in histoList:
1847        if 'PWDR' in histogram[:4]:
1848            Histogram = Histograms[histogram]
1849            hId = Histogram['hId']
1850            hfx = ':%d:'%(hId)
1851            Limits = calcControls[hfx+'Limits']
1852            if 'C' in calcControls[hfx+'histType']:
1853                shl = max(parmDict[hfx+'SH/L'],0.0005)
1854                Ka2 = False
1855                kRatio = 0.0
1856                if hfx+'Lam1' in parmDict.keys():
1857                    Ka2 = True
1858                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1859                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1860            x,y,w,yc,yb,yd = Histogram['Data']
1861            xB = np.searchsorted(x,Limits[0])
1862            xF = np.searchsorted(x,Limits[1])
1863            ymb = np.array(y-yb)
1864            ymb = np.where(ymb,ymb,1.0)
1865            ycmb = np.array(yc-yb)
1866            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1867            refLists = Histogram['Reflection Lists']
1868            for phase in refLists:
1869                if phase not in Phases:     #skips deleted or renamed phases silently!
1870                    continue
1871                Phase = Phases[phase]
1872                im = 0
1873                if Phase['General']['Type'] in ['modulated','magnetic']:
1874                    im = 1
1875                pId = Phase['pId']
1876                phfx = '%d:%d:'%(pId,hId)
1877                refDict = refLists[phase]
1878                sumFo = 0.0
1879                sumdF = 0.0
1880                sumFosq = 0.0
1881                sumdFsq = 0.0
1882                sumInt = 0.0
1883                for refl in refDict['RefList']:
1884                    if 'C' in calcControls[hfx+'histType']:
1885                        yp = np.zeros_like(yb)
1886                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
1887                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1888                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1889                        iFin2 = iFin
1890                        if not iBeg+iFin:       #peak below low limit - skip peak
1891                            continue
1892                        elif not iBeg-iFin:     #peak above high limit - done
1893                            break
1894                        elif iBeg < iFin:
1895                            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
1896                            sumInt += refl[11+im]*refl[9+im]
1897                            if Ka2:
1898                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1899                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
1900                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1901                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1902                                if iFin2 > iBeg2: 
1903                                    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
1904                                    sumInt += refl[11+im]*refl[9+im]*kRatio
1905                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
1906                               
1907                    elif 'T' in calcControls[hfx+'histType']:
1908                        yp = np.zeros_like(yb)
1909                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
1910                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
1911                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
1912                        if iBeg < iFin:
1913                            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
1914                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
1915                            sumInt += refl[11+im]*refl[9+im]
1916                    Fo = np.sqrt(np.abs(refl[8+im]))
1917                    Fc = np.sqrt(np.abs(refl[9]+im))
1918                    sumFo += Fo
1919                    sumFosq += refl[8+im]**2
1920                    sumdF += np.abs(Fo-Fc)
1921                    sumdFsq += (refl[8+im]-refl[9+im])**2
1922                if sumFo:
1923                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1924                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1925                else:
1926                    Histogram['Residuals'][phfx+'Rf'] = 100.
1927                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
1928                Histogram['Residuals'][phfx+'sumInt'] = sumInt
1929                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1930                Histogram['Residuals']['hId'] = hId
1931        elif 'HKLF' in histogram[:4]:
1932            Histogram = Histograms[histogram]
1933            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1934               
1935def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1936    'Needs a doc string'
1937   
1938    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
1939        U = parmDict[hfx+'U']
1940        V = parmDict[hfx+'V']
1941        W = parmDict[hfx+'W']
1942        X = parmDict[hfx+'X']
1943        Y = parmDict[hfx+'Y']
1944        tanPos = tand(refl[5+im]/2.0)
1945        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1946        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1947        sig = max(0.001,sig)
1948        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
1949        gam = max(0.001,gam)
1950        return sig,gam
1951               
1952    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
1953        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
1954            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
1955        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
1956        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1957        sig += Ssig
1958        gam += Sgam
1959        return sig,gam
1960       
1961    def GetReflAlpBet(refl,im,hfx,parmDict):
1962        alp = parmDict[hfx+'alpha']/refl[4+im]
1963        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
1964        return alp,bet
1965       
1966    hId = Histogram['hId']
1967    hfx = ':%d:'%(hId)
1968    bakType = calcControls[hfx+'bakType']
1969    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1970    yc = np.zeros_like(yb)
1971    cw = np.diff(x)
1972    cw = np.append(cw,cw[-1])
1973       
1974    if 'C' in calcControls[hfx+'histType']:   
1975        shl = max(parmDict[hfx+'SH/L'],0.002)
1976        Ka2 = False
1977        if hfx+'Lam1' in parmDict.keys():
1978            wave = parmDict[hfx+'Lam1']
1979            Ka2 = True
1980            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1981            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1982        else:
1983            wave = parmDict[hfx+'Lam']
1984    for phase in Histogram['Reflection Lists']:
1985        refDict = Histogram['Reflection Lists'][phase]
1986        if phase not in Phases:     #skips deleted or renamed phases silently!
1987            continue
1988        Phase = Phases[phase]
1989        pId = Phase['pId']
1990        pfx = '%d::'%(pId)
1991        phfx = '%d:%d:'%(pId,hId)
1992        hfx = ':%d:'%(hId)
1993        SGData = Phase['General']['SGData']
1994        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1995        im = 0
1996        if Phase['General']['Type'] in ['modulated','magnetic']:
1997            SSGData = Phase['General']['SSGData']
1998            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1999            im = 1  #offset in SS reflection list
2000            #??
2001        Dij = GetDij(phfx,SGData,parmDict)
2002        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2003        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2004        if np.any(np.diag(G)<0.):
2005            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2006        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2007        Vst = np.sqrt(nl.det(G))    #V*
2008        if not Phase['General'].get('doPawley'):
2009            time0 = time.time()
2010            if im:
2011                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2012            else:
2013                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2014#            print 'sf calc time: %.3fs'%(time.time()-time0)
2015        time0 = time.time()
2016        badPeak = False
2017        for iref,refl in enumerate(refDict['RefList']):
2018            if 'C' in calcControls[hfx+'histType']:
2019                if im:
2020                    h,k,l,m = refl[:4]
2021                else:
2022                    h,k,l = refl[:3]
2023                Uniq = np.inner(refl[:3],SGMT)
2024                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2025                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2026                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2027                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2028                refl[11+im] *= Vst*Lorenz
2029                 
2030                if Phase['General'].get('doPawley'):
2031                    try:
2032                        if im:
2033                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2034                        else:
2035                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2036                        refl[9+im] = parmDict[pInd]
2037                    except KeyError:
2038#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2039                        continue
2040                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2041                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2042                iFin = np.searchsorted(x,refl[5+im]+fmax)
2043                if not iBeg+iFin:       #peak below low limit - skip peak
2044                    continue
2045                elif not iBeg-iFin:     #peak above high limit - done
2046                    break
2047                elif iBeg > iFin:   #bad peak coeff - skip
2048                    badPeak = True
2049                    continue
2050                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
2051                if Ka2:
2052                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2053                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2054                    iBeg = np.searchsorted(x,pos2-fmin)
2055                    iFin = np.searchsorted(x,pos2+fmax)
2056                    if not iBeg+iFin:       #peak below low limit - skip peak
2057                        continue
2058                    elif not iBeg-iFin:     #peak above high limit - done
2059                        return yc,yb
2060                    elif iBeg > iFin:   #bad peak coeff - skip
2061                        continue
2062                    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
2063            elif 'T' in calcControls[hfx+'histType']:
2064                h,k,l = refl[:3]
2065                Uniq = np.inner(refl[:3],SGMT)
2066                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2067                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2068#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2069                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2070                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2071                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2072                refl[11+im] *= Vst*Lorenz
2073                if Phase['General'].get('doPawley'):
2074                    try:
2075                        if im:
2076                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2077                        else:
2078                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2079                        refl[9+im] = parmDict[pInd]
2080                    except KeyError:
2081#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2082                        continue
2083                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2084                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2085                iFin = np.searchsorted(x,refl[5+im]+fmax)
2086                if not iBeg+iFin:       #peak below low limit - skip peak
2087                    continue
2088                elif not iBeg-iFin:     #peak above high limit - done
2089                    break
2090                elif iBeg > iFin:   #bad peak coeff - skip
2091                    badPeak = True
2092                    continue
2093                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]
2094#        print 'profile calc time: %.3fs'%(time.time()-time0)
2095    if badPeak:
2096        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2097    return yc,yb
2098   
2099def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2100    'Needs a doc string'
2101   
2102    def cellVaryDerv(pfx,SGData,dpdA): 
2103        if SGData['SGLaue'] in ['-1',]:
2104            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2105                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2106        elif SGData['SGLaue'] in ['2/m',]:
2107            if SGData['SGUniq'] == 'a':
2108                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2109            elif SGData['SGUniq'] == 'b':
2110                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2111            else:
2112                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2113        elif SGData['SGLaue'] in ['mmm',]:
2114            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2115        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2116            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2117        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2118            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2119        elif SGData['SGLaue'] in ['3R', '3mR']:
2120            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2121        elif SGData['SGLaue'] in ['m3m','m3']:
2122            return [[pfx+'A0',dpdA[0]]]
2123           
2124    # create a list of dependent variables and set up a dictionary to hold their derivatives
2125    dependentVars = G2mv.GetDependentVars()
2126    depDerivDict = {}
2127    for j in dependentVars:
2128        depDerivDict[j] = np.zeros(shape=(len(x)))
2129    #print 'dependent vars',dependentVars
2130    lenX = len(x)               
2131    hId = Histogram['hId']
2132    hfx = ':%d:'%(hId)
2133    bakType = calcControls[hfx+'bakType']
2134    dMdv = np.zeros(shape=(len(varylist),len(x)))
2135    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2136    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2137        bBpos =varylist.index(hfx+'Back;0')
2138        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2139    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2140    for name in varylist:
2141        if 'Debye' in name:
2142            id = int(name.split(';')[-1])
2143            parm = name[:int(name.rindex(';'))]
2144            ip = names.index(parm)
2145            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2146    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2147    for name in varylist:
2148        if 'BkPk' in name:
2149            parm,id = name.split(';')
2150            id = int(id)
2151            if parm in names:
2152                ip = names.index(parm)
2153                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2154    cw = np.diff(x)
2155    cw = np.append(cw,cw[-1])
2156    Ka2 = False #also for TOF!
2157    if 'C' in calcControls[hfx+'histType']:   
2158        shl = max(parmDict[hfx+'SH/L'],0.002)
2159        if hfx+'Lam1' in parmDict.keys():
2160            wave = parmDict[hfx+'Lam1']
2161            Ka2 = True
2162            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2163            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2164        else:
2165            wave = parmDict[hfx+'Lam']
2166    for phase in Histogram['Reflection Lists']:
2167        refDict = Histogram['Reflection Lists'][phase]
2168        if phase not in Phases:     #skips deleted or renamed phases silently!
2169            continue
2170        Phase = Phases[phase]
2171        SGData = Phase['General']['SGData']
2172        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2173        im = 0
2174        if Phase['General']['Type'] in ['modulated','magnetic']:
2175            SSGData = Phase['General']['SSGData']
2176            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2177            im = 1  #offset in SS reflection list
2178            #??
2179        pId = Phase['pId']
2180        pfx = '%d::'%(pId)
2181        phfx = '%d:%d:'%(pId,hId)
2182        Dij = GetDij(phfx,SGData,parmDict)
2183        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2184        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2185        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2186        if not Phase['General'].get('doPawley'):
2187            time0 = time.time()
2188            if im:
2189                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2190            else:
2191                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2192#            print 'sf-derv time %.3fs'%(time.time()-time0)
2193            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2194        time0 = time.time()
2195        for iref,refl in enumerate(refDict['RefList']):
2196            if im:
2197                h,k,l,m = refl[:4]
2198            else:
2199                h,k,l = refl[:3]
2200            Uniq = np.inner(refl[:3],SGMT)
2201            if 'T' in calcControls[hfx+'histType']:
2202                wave = refl[14+im]
2203            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2204            if 'C' in calcControls[hfx+'histType']:        #CW powder
2205                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2206            else: #'T'OF
2207                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2208            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2209            iFin = np.searchsorted(x,refl[5+im]+fmax)
2210            if not iBeg+iFin:       #peak below low limit - skip peak
2211                continue
2212            elif not iBeg-iFin:     #peak above high limit - done
2213                break
2214            pos = refl[5+im]
2215            if 'C' in calcControls[hfx+'histType']:
2216                tanth = tand(pos/2.0)
2217                costh = cosd(pos/2.0)
2218                lenBF = iFin-iBeg
2219                dMdpk = np.zeros(shape=(6,lenBF))
2220                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2221                for i in range(5):
2222                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2223                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2224                if Ka2:
2225                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2226                    iBeg2 = np.searchsorted(x,pos2-fmin)
2227                    iFin2 = np.searchsorted(x,pos2+fmax)
2228                    if iBeg2-iFin2:
2229                        lenBF2 = iFin2-iBeg2
2230                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2231                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2232                        for i in range(5):
2233                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2234                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2235                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2236            else:   #'T'OF
2237                lenBF = iFin-iBeg
2238                if lenBF < 0:   #bad peak coeff
2239                    break
2240                dMdpk = np.zeros(shape=(6,lenBF))
2241                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2242                for i in range(6):
2243                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2244                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2245            if Phase['General'].get('doPawley'):
2246                dMdpw = np.zeros(len(x))
2247                try:
2248                    if im:
2249                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2250                    else:
2251                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2252                    idx = varylist.index(pIdx)
2253                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2254                    if Ka2: #not for TOF either
2255                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2256                    dMdv[idx] = dMdpw
2257                except: # ValueError:
2258                    pass
2259            if 'C' in calcControls[hfx+'histType']:
2260                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2261                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2262                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2263                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2264                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2265                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2266                    hfx+'DisplaceY':[dpdY,'pos'],}
2267                if 'Bragg' in calcControls[hfx+'instType']:
2268                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2269                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2270                else:
2271                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2272            else:   #'T'OF
2273                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2274                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2275                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2276                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2277                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2278                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2279                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2280                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2281            for name in names:
2282                item = names[name]
2283                if name in varylist:
2284                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2285                    if Ka2 and iFin2-iBeg2:
2286                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2287                elif name in dependentVars:
2288                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2289                    if Ka2 and iFin2-iBeg2:
2290                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2291            for iPO in dIdPO:
2292                if iPO in varylist:
2293                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2294                    if Ka2 and iFin2-iBeg2:
2295                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2296                elif iPO in dependentVars:
2297                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2298                    if Ka2 and iFin2-iBeg2:
2299                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2300            for i,name in enumerate(['omega','chi','phi']):
2301                aname = pfx+'SH '+name
2302                if aname in varylist:
2303                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2304                    if Ka2 and iFin2-iBeg2:
2305                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2306                elif aname in dependentVars:
2307                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2308                    if Ka2 and iFin2-iBeg2:
2309                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2310            for iSH in dFdODF:
2311                if iSH in varylist:
2312                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2313                    if Ka2 and iFin2-iBeg2:
2314                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2315                elif iSH in dependentVars:
2316                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2317                    if Ka2 and iFin2-iBeg2:
2318                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2319            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2320            for name,dpdA in cellDervNames:
2321                if name in varylist:
2322                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2323                    if Ka2 and iFin2-iBeg2:
2324                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2325                elif name in dependentVars:
2326                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2327                    if Ka2 and iFin2-iBeg2:
2328                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2329            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2330            for name in dDijDict:
2331                if name in varylist:
2332                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2333                    if Ka2 and iFin2-iBeg2:
2334                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2335                elif name in dependentVars:
2336                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2337                    if Ka2 and iFin2-iBeg2:
2338                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2339            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2340                if name in varylist:
2341                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2342                    if Ka2 and iFin2-iBeg2:
2343                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2344                elif name in dependentVars:
2345                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2346                    if Ka2 and iFin2-iBeg2:
2347                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2348            if 'C' in calcControls[hfx+'histType']:
2349                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2350            else:   #'T'OF
2351                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2352            for name in gamDict:
2353                if name in varylist:
2354                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2355                    if Ka2 and iFin2-iBeg2:
2356                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2357                elif name in dependentVars:
2358                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2359                    if Ka2 and iFin2-iBeg2:
2360                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2361            for name in sigDict:
2362                if name in varylist:
2363                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2364                    if Ka2 and iFin2-iBeg2:
2365                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2366                elif name in dependentVars:
2367                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2368                    if Ka2 and iFin2-iBeg2:
2369                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2370            for name in ['BabA','BabU']:
2371                if refl[9+im]:
2372                    if phfx+name in varylist:
2373                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2374                        if Ka2 and iFin2-iBeg2:
2375                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
2376                    elif phfx+name in dependentVars:                   
2377                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
2378                        if Ka2 and iFin2-iBeg2:
2379                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
2380            if not Phase['General'].get('doPawley'):
2381                #do atom derivatives -  for RB,F,X & U so far             
2382                corr = dervDict['int']/refl[9+im]
2383                if Ka2 and iFin2-iBeg2:
2384                    corr2 = dervDict2['int']/refl[9+im]
2385                for name in varylist+dependentVars:
2386                    if '::RBV;' in name:
2387                        pass
2388                    else:
2389                        try:
2390                            aname = name.split(pfx)[1][:2]
2391                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2392                        except IndexError:
2393                            continue
2394                    if name in varylist:
2395                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2396                        if Ka2 and iFin2-iBeg2:
2397                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2398                    elif name in dependentVars:
2399                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2400                        if Ka2 and iFin2-iBeg2:
2401                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2402    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2403    # now process derivatives in constraints
2404    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2405    return dMdv
2406   
2407def UserRejectHKL(ref,im,userReject):
2408    if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2409        return False
2410    elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2411        return False
2412    elif ref[11+im] < userReject['MinExt']:
2413        return False
2414    elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2415        return False
2416    return True
2417   
2418def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2419    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2420    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2421    are returned to either dervRefine or HessRefine.
2422
2423    :returns:
2424    '''
2425    nobs = Histogram['Residuals']['Nobs']
2426    hId = Histogram['hId']
2427    hfx = ':%d:'%(hId)
2428    pfx = '%d::'%(Phase['pId'])
2429    phfx = '%d:%d:'%(Phase['pId'],hId)
2430    SGData = Phase['General']['SGData']
2431    im = 0
2432    if Phase['General']['Type'] in ['modulated','magnetic']:
2433        SSGData = Phase['General']['SSGData']
2434        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2435        im = 1  #offset in SS reflection list
2436    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2437    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2438    refDict = Histogram['Data']
2439    if im:
2440        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2441    else:
2442        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2443    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2444    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2445    dependentVars = G2mv.GetDependentVars()
2446    depDerivDict = {}
2447    for j in dependentVars:
2448        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2449    wdf = np.zeros(len(refDict['RefList']))
2450    if calcControls['F**2']:
2451        for iref,ref in enumerate(refDict['RefList']):
2452            if ref[6+im] > 0:
2453                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2454                w = 1.0/ref[6+im]
2455                if ref[3+im] > 0:
2456                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2457                    for j,var in enumerate(varylist):
2458                        if var in dFdvDict:
2459                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2460                    for var in dependentVars:
2461                        if var in dFdvDict:
2462                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2463                    if phfx+'Scale' in varylist:
2464                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2465                    elif phfx+'Scale' in dependentVars:
2466                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK
2467                    for item in ['Ep','Es','Eg']:
2468                        if phfx+item in varylist and phfx+item in dervDict:
2469                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2470                        elif phfx+item in dependentVars and phfx+item in dervDict:
2471                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2472                    for item in ['BabA','BabU']:
2473                        if phfx+item in varylist:
2474                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2475                        elif phfx+item in dependentVars:
2476                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2477    else:   #F refinement
2478        for iref,ref in enumerate(refDict['RefList']):
2479            if ref[5+im] > 0.:
2480                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2481                Fo = np.sqrt(ref[5+im])
2482                Fc = np.sqrt(ref[7+im])
2483                w = 1.0/ref[6+im]
2484                if ref[3+im] > 0:
2485                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
2486                    for j,var in enumerate(varylist):
2487                        if var in dFdvDict:
2488                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2489                    for var in dependentVars:
2490                        if var in dFdvDict:
2491                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2492                    if phfx+'Scale' in varylist:
2493                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]  #OK
2494                    elif phfx+'Scale' in dependentVars:
2495                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]   #OK                   
2496                    for item in ['Ep','Es','Eg']:   #OK!
2497                        if phfx+item in varylist and phfx+item in dervDict:
2498                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] 
2499                        elif phfx+item in dependentVars and phfx+item in dervDict:
2500                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2501                    for item in ['BabA','BabU']:
2502                        if phfx+item in varylist:
2503                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2504                        elif phfx+item in dependentVars:
2505                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2506    return dMdvh,depDerivDict,wdf
2507   
2508
2509def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2510    '''Loop over histograms and compute derivatives of the fitting
2511    model (M) with respect to all parameters.  Results are returned in
2512    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2513    n is the number of parameters and m is the number of data
2514    points. This can exceed memory when m gets large. This routine is
2515    used when refinement derivatives are selected as "analtytic
2516    Jacobian" in Controls.
2517
2518    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2519    '''
2520    parmDict.update(zip(varylist,values))
2521    G2mv.Dict2Map(parmDict,varylist)
2522    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2523    nvar = len(varylist)
2524    dMdv = np.empty(0)
2525    histoList = Histograms.keys()
2526    histoList.sort()
2527    for histogram in histoList:
2528        if 'PWDR' in histogram[:4]:
2529            Histogram = Histograms[histogram]
2530            hId = Histogram['hId']
2531            hfx = ':%d:'%(hId)
2532            wtFactor = calcControls[hfx+'wtFactor']
2533            Limits = calcControls[hfx+'Limits']
2534            x,y,w,yc,yb,yd = Histogram['Data']
2535            xB = np.searchsorted(x,Limits[0])
2536            xF = np.searchsorted(x,Limits[1])
2537            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2538                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2539        elif 'HKLF' in histogram[:4]:
2540            Histogram = Histograms[histogram]
2541            phase = Histogram['Reflection Lists']
2542            Phase = Phases[phase]
2543            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2544            hfx = ':%d:'%(Histogram['hId'])
2545            wtFactor = calcControls[hfx+'wtFactor']
2546            # now process derivatives in constraints
2547            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2548        else:
2549            continue        #skip non-histogram entries
2550        if len(dMdv):
2551            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2552        else:
2553            dMdv = np.sqrt(wtFactor)*dMdvh
2554           
2555    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2556    if np.any(pVals):
2557        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2558        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2559       
2560    return dMdv
2561
2562def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2563    '''Loop over histograms and compute derivatives of the fitting
2564    model (M) with respect to all parameters.  For each histogram, the
2565    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2566    number of parameters and m is the number of data points *in the
2567    histogram*. The (n by n) Hessian is computed from each Jacobian
2568    and it is returned.  This routine is used when refinement
2569    derivatives are selected as "analtytic Hessian" in Controls.
2570
2571    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2572    '''
2573    parmDict.update(zip(varylist,values))
2574    G2mv.Dict2Map(parmDict,varylist)
2575    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2576    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2577    nvar = len(varylist)
2578    Hess = np.empty(0)
2579    histoList = Histograms.keys()
2580    histoList.sort()
2581    for histogram in histoList:
2582        if 'PWDR' in histogram[:4]:
2583            Histogram = Histograms[histogram]
2584            hId = Histogram['hId']
2585            hfx = ':%d:'%(hId)
2586            wtFactor = calcControls[hfx+'wtFactor']
2587            Limits = calcControls[hfx+'Limits']
2588            x,y,w,yc,yb,yd = Histogram['Data']
2589            W = wtFactor*w
2590            dy = y-yc
2591            xB = np.searchsorted(x,Limits[0])
2592            xF = np.searchsorted(x,Limits[1])
2593            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2594                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2595            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2596            Dy = dy[xB:xF][np.newaxis,:]
2597            dMdvh *= Wt
2598            if dlg:
2599                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2600            if len(Hess):
2601                Hess += np.inner(dMdvh,dMdvh)
2602                dMdvh *= Wt*Dy
2603                Vec += np.sum(dMdvh,axis=1)
2604            else:
2605                Hess = np.inner(dMdvh,dMdvh)
2606                dMdvh *= Wt*Dy
2607                Vec = np.sum(dMdvh,axis=1)
2608        elif 'HKLF' in histogram[:4]:
2609            Histogram = Histograms[histogram]
2610            phase = Histogram['Reflection Lists']
2611            Phase = Phases[phase]
2612            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2613            hId = Histogram['hId']
2614            hfx = ':%d:'%(Histogram['hId'])
2615            wtFactor = calcControls[hfx+'wtFactor']
2616            # now process derivatives in constraints
2617            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2618#            print 'matrix build time: %.3f'%(time.time()-time0)
2619
2620            if dlg:
2621                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2622            if len(Hess):
2623                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2624                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2625            else:
2626                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2627                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2628        else:
2629            continue        #skip non-histogram entries
2630    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2631    if np.any(pVals):
2632        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2633        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2634        Hess += np.inner(dpdv*pWt,dpdv)
2635    return Vec,Hess
2636
2637def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2638    '''Computes the point-by-point discrepancies between every data point in every histogram
2639    and the observed value
2640   
2641    :returns: an np array of differences between observed and computed diffraction values.
2642    '''
2643    Values2Dict(parmDict, varylist, values)
2644    G2mv.Dict2Map(parmDict,varylist)
2645    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2646    M = np.empty(0)
2647    SumwYo = 0
2648    Nobs = 0
2649    Nrej = 0
2650    Next = 0
2651    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2652    histoList = Histograms.keys()
2653    histoList.sort()
2654    for histogram in histoList:
2655        if 'PWDR' in histogram[:4]:
2656            Histogram = Histograms[histogram]
2657            hId = Histogram['hId']
2658            hfx = ':%d:'%(hId)
2659            wtFactor = calcControls[hfx+'wtFactor']
2660            Limits = calcControls[hfx+'Limits']
2661            x,y,w,yc,yb,yd = Histogram['Data']
2662            yc *= 0.0                           #zero full calcd profiles
2663            yb *= 0.0
2664            yd *= 0.0
2665            xB = np.searchsorted(x,Limits[0])
2666            xF = np.searchsorted(x,Limits[1])
2667            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2668                varylist,Histogram,Phases,calcControls,pawleyLookup)
2669            yc[xB:xF] += yb[xB:xF]
2670            if not np.any(y):                   #fill dummy data
2671                rv = st.poisson(yc[xB:xF])
2672                y[xB:xF] = rv.rvs()
2673                Z = np.ones_like(yc[xB:xF])
2674                Z[1::2] *= -1
2675                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2676                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2677            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2678            W = wtFactor*w
2679            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2680            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2681            Nobs += Histogram['Residuals']['Nobs']
2682            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2683            SumwYo += Histogram['Residuals']['sumwYo']
2684            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2685            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2686            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2687            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2688            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.))
2689            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.))
2690            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2691            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2692            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2693            if dlg:
2694                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2695            M = np.concatenate((M,wdy))
2696#end of PWDR processing
2697        elif 'HKLF' in histogram[:4]:
2698            Histogram = Histograms[histogram]
2699            Histogram['Residuals'] = {}
2700            phase = Histogram['Reflection Lists']
2701            Phase = Phases[phase]
2702            hId = Histogram['hId']
2703            hfx = ':%d:'%(hId)
2704            wtFactor = calcControls[hfx+'wtFactor']
2705            pfx = '%d::'%(Phase['pId'])
2706            phfx = '%d:%d:'%(Phase['pId'],hId)
2707            SGData = Phase['General']['SGData']
2708            im = 0
2709            if Phase['General']['Type'] in ['modulated','magnetic']:
2710                SSGData = Phase['General']['SSGData']
2711                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2712                im = 1  #offset in SS reflection list
2713            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2714            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2715            refDict = Histogram['Data']
2716            time0 = time.time()
2717            if im:
2718                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2719            else:
2720                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2721#            print 'sf-calc time: %.3f'%(time.time()-time0)
2722            df = np.zeros(len(refDict['RefList']))
2723            sumwYo = 0
2724            sumFo = 0
2725            sumFo2 = 0
2726            sumdF = 0
2727            sumdF2 = 0
2728            nobs = 0
2729            nrej = 0
2730            next = 0
2731            if calcControls['F**2']:
2732                for i,ref in enumerate(refDict['RefList']):
2733                    if ref[6+im] > 0:
2734                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2735                        w = 1.0/ref[6+im]   # 1/sig(F^2)
2736                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
2737                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2738                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2739                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2740                            Fo = np.sqrt(ref[5+im])
2741                            sumFo += Fo
2742                            sumFo2 += ref[5+im]
2743                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2744                            sumdF2 += abs(ref[5+im]-ref[7+im])
2745                            nobs += 1
2746                            df[i] = -w*(ref[5+im]-ref[7+im])
2747                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
2748                        else:
2749                            if ref[3+im]:
2750                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2751                                nrej += 1
2752                            else:   #sp.gp.extinct
2753                                next += 1
2754            else:
2755                for i,ref in enumerate(refDict['RefList']):
2756                    if ref[5+im] > 0.:
2757                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2758                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
2759                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2760                        Fo = np.sqrt(ref[5+im])
2761                        Fc = np.sqrt(ref[7+im])
2762                        w = 2.0*Fo/ref[6+im]    # 1/sig(F)?
2763                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2764                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2765                            sumFo += Fo
2766                            sumFo2 += ref[5+im]
2767                            sumdF += abs(Fo-Fc)
2768                            sumdF2 += abs(ref[5+im]-ref[7+im])
2769                            nobs += 1
2770                            df[i] = -w*(Fo-Fc)
2771                            sumwYo += (w*Fo)**2
2772                        else:
2773                            if ref[3+im]:
2774                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2775                                nrej += 1
2776                            else:   #sp.gp.extinct
2777                                next += 1
2778            Histogram['Residuals']['Nobs'] = nobs
2779            Histogram['Residuals']['sumwYo'] = sumwYo
2780            SumwYo += sumwYo
2781            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/sumwYo)*100.)
2782            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2783            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2784            Histogram['Residuals'][phfx+'Nref'] = nobs
2785            Histogram['Residuals'][phfx+'Nrej'] = nrej
2786            Histogram['Residuals'][phfx+'Next'] = next
2787            Nobs += nobs
2788            Nrej += nrej
2789            Next += next
2790            if dlg:
2791                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2792            M = np.concatenate((M,wtFactor*df))
2793# end of HKLF processing
2794    Histograms['sumwYo'] = SumwYo
2795    Histograms['Nobs'] = Nobs
2796    Histograms['Nrej'] = Nrej
2797    Histograms['Next'] = Next
2798    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2799    if dlg:
2800        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2801        if not GoOn:
2802            parmDict['saved values'] = values
2803            dlg.Destroy()
2804            raise G2obj.G2Exception('User abort')         #Abort!!
2805    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2806    if len(pVals):
2807        pSum = np.sum(pWt*pVals**2)
2808        for name in pWsum:
2809            if pWsum[name]:
2810                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2811        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2812        Nobs += len(pVals)
2813        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2814    return M
2815                       
Note: See TracBrowser for help on using the repository browser.