source: trunk/GSASIIstrMath.py @ 1459

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

add instrument parameters (flight path & detector 2-theta) needed for TOF
rework reflection list for TOF
change default diff curve & reflection marker offsets
change weighting to instrument constants calibration to be 1/esd2 from peak fit positions - works a lot better
1st shot at TOF Rietveld refinement with derivatives - need to be checked for correctness (some are wrong)

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