source: trunk/GSASIIstrMath.py @ 1781

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

fix bug from old Pref.Ori. not having penalty items
Implement colors for bad items in Reflection Lists; Prfo < 0 for PWDR
red Fc for abs(Fo-Fc)/sig > 10 & yellow if > 3.
Implement user selected rejection scheme for HKLF; set mul < 0. Get calculated but not used in least squares

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