source: trunk/GSASIIstrMath.py @ 1860

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