source: trunk/GSASIIstrMath.py @ 1613

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

missing menu items in PDF controls
trap PDF plotting if no PDF to plot
some more on SS refinement constraints
Begin implementation of SStructureFactor & SStructureFactorDerv

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