source: trunk/GSASIIstrMath.py @ 1625

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

add a parameter to result from G2stIO.GetPhaseData?
add modulation functions to G2Math
add modulation names to G2obj
implement various wave types for modulations
plot position modulation wave function on contour plots
implement shift of modulation plot by +/-/0 keys
temporarily default G2spc.GetSSfxuinel to 1,-1 site symm.
move maxSSwave dict out of parmDict - now in controlDict
implement import of Sawtooth parms from J2K files.

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