source: trunk/GSASIIstrMath.py @ 1810

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

texture PO penalty fxn. now active

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