source: trunk/GSASIIstrMath.py @ 1618

Last change on this file since 1618 was 1618, checked in by vondreele, 8 years ago

fix space group print in shelx export
fix errors when data is deleted
add HKLshow to change in space group in index page
more on supersymmetry refinement

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