source: trunk/GSASIIstrMath.py @ 1496

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

Add sequential peak fitting.
Change Back:n to Back;n for background parameters
Also DebyeA, etc.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 116.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2014-09-15 00:48:39 +0000 (Mon, 15 Sep 2014) $
8# $Author: vondreele $
9# $Revision: 1496 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1496 2014-09-15 00:48:39Z 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: 1496 $")
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 for al2O3!   
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; returns extinction & derivative
802    '''
803    extCor = 1.0
804    dervDict = {}
805    if calcControls[phfx+'EType'] != 'None':
806        SQ = 1/(4.*ref[4]**2)
807        if 'C' in parmDict[hfx+'Type']:           
808            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
809        else:   #'T'
810            cos2T = 1.0-2.*SQ*ref[12]**2                       #cos(2theta)           
811        if 'SXC' in parmDict[hfx+'Type']:
812            AV = 7.9406e5/parmDict[pfx+'Vol']**2
813            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
814            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
815            PLZ = AV*P12*ref[7]*parmDict[hfx+'Lam']**2
816        elif 'SNT' in parmDict[hfx+'Type']:
817            AV = 1.e7/parmDict[pfx+'Vol']**2
818            PL = SQ
819            PLZ = AV*ref[7]*ref[12]**2
820        elif 'SNC' in parmDict[hfx+'Type']:
821            AV = 1.e7/parmDict[pfx+'Vol']**2
822            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
823            PLZ = AV*ref[9]*parmDict[hfx+'Lam']**2      #Fcsq as per GSAS, why not FcTsq (ref[9])?
824           
825        if 'Primary' in calcControls[phfx+'EType']:
826            PLZ *= 1.5
827        else:
828            if 'C' in parmDict[hfx+'Type']:
829                PLZ *= calcControls[phfx+'Tbar']
830            else: #'T'
831                PLZ *= ref[13]      #t-bar
832        if 'Primary' in calcControls[phfx+'EType']:
833            PLZ *= 1.5
834            PSIG = parmDict[phfx+'Ep']
835        elif 'I & II' in calcControls[phfx+'EType']:
836            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
837        elif 'Type II' in calcControls[phfx+'EType']:
838            PSIG = parmDict[phfx+'Es']
839        else:       # 'Secondary Type I'
840            PSIG = parmDict[phfx+'Eg']/PL
841           
842        AG = 0.58+0.48*cos2T+0.24*cos2T**2
843        AL = 0.025+0.285*cos2T
844        BG = 0.02-0.025*cos2T
845        BL = 0.15-0.2*(0.75-cos2T)**2
846        if cos2T < 0.:
847            BL = -0.45*cos2T
848        CG = 2.
849        CL = 2.
850        PF = PLZ*PSIG
851       
852        if 'Gaussian' in calcControls[phfx+'EApprox']:
853            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
854            extCor = np.sqrt(PF4)
855            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
856        else:
857            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
858            extCor = np.sqrt(PF4)
859            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
860
861        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
862            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
863        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
864            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
865        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
866            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
867               
868    return 1./extCor,dervDict
869   
870def Dict2Values(parmdict, varylist):
871    '''Use before call to leastsq to setup list of values for the parameters
872    in parmdict, as selected by key in varylist'''
873    return [parmdict[key] for key in varylist] 
874   
875def Values2Dict(parmdict, varylist, values):
876    ''' Use after call to leastsq to update the parameter dictionary with
877    values corresponding to keys in varylist'''
878    parmdict.update(zip(varylist,values))
879   
880def GetNewCellParms(parmDict,varyList):
881    'Needs a doc string'
882    newCellDict = {}
883    Anames = ['A'+str(i) for i in range(6)]
884    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
885    for item in varyList:
886        keys = item.split(':')
887        if keys[2] in Ddict:
888            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
889            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
890            newCellDict[parm] = [key,parmDict[key]-parmDict[item]]
891    return newCellDict          # is e.g. {'0::D11':A0-D11}
892   
893def ApplyXYZshifts(parmDict,varyList):
894    '''
895    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
896   
897    :param dict parmDict: parameter dictionary
898    :param list varyList: list of variables (not used!)
899    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
900
901    '''
902    newAtomDict = {}
903    for item in parmDict:
904        if 'dA' in item:
905            parm = ''.join(item.split('d'))
906            parmDict[parm] += parmDict[item]
907            newAtomDict[item] = [parm,parmDict[parm]]
908    return newAtomDict
909   
910def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
911    'Spherical harmonics texture'
912    IFCoup = 'Bragg' in calcControls[hfx+'instType']
913    if 'T' in calcControls[hfx+'histType']:
914        tth = parmDict[hfx+'2-theta']
915    else:
916        tth = refl[5]
917    odfCor = 1.0
918    H = refl[:3]
919    cell = G2lat.Gmat2cell(g)
920    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
921    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
922    phi,beta = G2lat.CrsAng(H,cell,SGData)
923    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
924    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
925    for item in SHnames:
926        L,M,N = eval(item.strip('C'))
927        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
928        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
929        Lnorm = G2lat.Lnorm(L)
930        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
931    return odfCor
932   
933def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
934    'Spherical harmonics texture derivatives'
935    if 'T' in calcControls[hfx+'histType']:
936        tth = parmDict[hfx+'2-theta']
937    else:
938        tth = refl[5]
939    FORPI = 4.0*np.pi
940    IFCoup = 'Bragg' in calcControls[hfx+'instType']
941    odfCor = 1.0
942    dFdODF = {}
943    dFdSA = [0,0,0]
944    H = refl[:3]
945    cell = G2lat.Gmat2cell(g)
946    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
947    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
948    phi,beta = G2lat.CrsAng(H,cell,SGData)
949    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
950    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
951    for item in SHnames:
952        L,M,N = eval(item.strip('C'))
953        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
954        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
955        Lnorm = G2lat.Lnorm(L)
956        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
957        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
958        for i in range(3):
959            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
960    return odfCor,dFdODF,dFdSA
961   
962def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
963    'spherical harmonics preferred orientation (cylindrical symmetry only)'
964    if 'T' in calcControls[hfx+'histType']:
965        tth = parmDict[hfx+'2-theta']
966    else:
967        tth = refl[5]
968    odfCor = 1.0
969    H = refl[:3]
970    cell = G2lat.Gmat2cell(g)
971    Sangl = [0.,0.,0.]
972    if 'Bragg' in calcControls[hfx+'instType']:
973        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
974        IFCoup = True
975    else:
976        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
977        IFCoup = False
978    phi,beta = G2lat.CrsAng(H,cell,SGData)
979    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
980    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
981    for item in SHnames:
982        L,N = eval(item.strip('C'))
983        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
984        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
985    return np.squeeze(odfCor)
986   
987def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
988    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
989    if 'T' in calcControls[hfx+'histType']:
990        tth = parmDict[hfx+'2-theta']
991    else:
992        tth = refl[5]
993    FORPI = 12.5663706143592
994    odfCor = 1.0
995    dFdODF = {}
996    H = refl[:3]
997    cell = G2lat.Gmat2cell(g)
998    Sangl = [0.,0.,0.]
999    if 'Bragg' in calcControls[hfx+'instType']:
1000        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1001        IFCoup = True
1002    else:
1003        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1004        IFCoup = False
1005    phi,beta = G2lat.CrsAng(H,cell,SGData)
1006    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1007    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1008    for item in SHnames:
1009        L,N = eval(item.strip('C'))
1010        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1011        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1012        dFdODF[phfx+item] = Kcsl*Lnorm
1013    return odfCor,dFdODF
1014   
1015def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1016    'March-Dollase preferred orientation correction'
1017    POcorr = 1.0
1018    MD = parmDict[phfx+'MD']
1019    if MD != 1.0:
1020        MDAxis = calcControls[phfx+'MDAxis']
1021        sumMD = 0
1022        for H in uniq:           
1023            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1024            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1025            sumMD += A**3
1026        POcorr = sumMD/len(uniq)
1027    return POcorr
1028   
1029def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1030    'Needs a doc string'
1031    POcorr = 1.0
1032    POderv = {}
1033    if calcControls[phfx+'poType'] == 'MD':
1034        MD = parmDict[phfx+'MD']
1035        MDAxis = calcControls[phfx+'MDAxis']
1036        sumMD = 0
1037        sumdMD = 0
1038        for H in uniq:           
1039            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1040            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1041            sumMD += A**3
1042            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1043        POcorr = sumMD/len(uniq)
1044        POderv[phfx+'MD'] = sumdMD/len(uniq)
1045    else:   #spherical harmonics
1046        if calcControls[phfx+'SHord']:
1047            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1048    return POcorr,POderv
1049   
1050def GetAbsorb(refl,hfx,calcControls,parmDict):
1051    'Needs a doc string'
1052    if 'Debye' in calcControls[hfx+'instType']:
1053        if 'T' in calcControls[hfx+'histType']:
1054            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1055        else:
1056            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1057    else:
1058        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1059   
1060def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1061    'Needs a doc string'
1062    if 'Debye' in calcControls[hfx+'instType']:
1063        if 'T' in calcControls[hfx+'histType']:
1064            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1065        else:
1066            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1067    else:
1068        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]))
1069       
1070def GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict):
1071    'Needs a doc string'
1072    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1073    pi2 = np.sqrt(2./np.pi)
1074    if 'T' in calcControls[hfx+'histType']:
1075        sth2 = sind(parmDict[hfx+'2-theta']/2.)**2
1076        wave = refl[14]
1077    else:   #'C'W
1078        sth2 = sind(refl[5]/2.)**2
1079        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1080    c2th = 1.-2.0*sth2
1081    flv2 = refl[9]*(wave/parmDict[pfx+'Vol'])**2
1082    if 'X' in calcControls[hfx+'histType']:
1083        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1084    xfac = flv2*parmDict[phfx+'Extinction']
1085    exb = 1.0
1086    if xfac > -1.:
1087        exb = 1./(1.+xfac)
1088    exl = 1.0
1089    if 0 < xfac <= 1.:
1090        xn = np.array([xfac**(i+1) for i in range(6)])
1091        exl = np.sum(xn*coef)
1092    elif xfac > 1.:
1093        xfac2 = 1./np.sqrt(xfac)
1094        exl = pi2*(1.-0.125/xfac)*xfac2
1095    return exb*sth2+exl*(1.-sth2)
1096   
1097def GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict):
1098    'Needs a doc string'
1099    delt = 0.001
1100    parmDict[phfx+'Extinction'] += delt
1101    plus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1102    parmDict[phfx+'Extinction'] -= 2.*delt
1103    minus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1104    parmDict[phfx+'Extinction'] += delt
1105    return (plus-minus)/(2.*delt)   
1106   
1107def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1108    'Needs a doc string'    #need powder extinction!
1109    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1110    if 'X' in parmDict[hfx+'Type']:
1111        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1112    POcorr = 1.0
1113    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture
1114        POcorr = SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1115    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1116        POcorr = GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1117    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1118        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1119    Icorr *= POcorr
1120    AbsCorr = 1.0
1121    AbsCorr = GetAbsorb(refl,hfx,calcControls,parmDict)
1122    Icorr *= AbsCorr
1123    ExtCorr = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1124    Icorr *= ExtCorr
1125    return Icorr,POcorr,AbsCorr,ExtCorr
1126   
1127def GetIntensityDerv(refl,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1128    'Needs a doc string'    #need powder extinction derivs!
1129    dIdsh = 1./parmDict[hfx+'Scale']
1130    dIdsp = 1./parmDict[phfx+'Scale']
1131    if 'X' in parmDict[hfx+'Type']:
1132        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1133        dIdPola /= pola
1134    else:       #'N'
1135        dIdPola = 0.0
1136    dFdODF = {}
1137    dFdSA = [0,0,0]
1138    dIdPO = {}
1139    if pfx+'SHorder' in parmDict:
1140        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1141        for iSH in dFdODF:
1142            dFdODF[iSH] /= odfCor
1143        for i in range(3):
1144            dFdSA[i] /= odfCor
1145    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1146        POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1147        for iPO in dIdPO:
1148            dIdPO[iPO] /= POcorr
1149    if 'T' in parmDict[hfx+'Type']:
1150        dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[16] #wave/abs corr
1151        dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[17]    #/ext corr
1152    else:
1153        dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[13] #wave/abs corr
1154        dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[14]    #/ext corr       
1155    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1156       
1157def GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1158    'Needs a doc string'
1159    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1160        costh = cosd(refl[5]/2.)
1161        #crystallite size
1162        if calcControls[phfx+'SizeType'] == 'isotropic':
1163            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1164        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1165            H = np.array(refl[:3])
1166            P = np.array(calcControls[phfx+'SizeAxis'])
1167            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1168            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1169            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1170        else:           #ellipsoidal crystallites
1171            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1172            H = np.array(refl[:3])
1173            lenR = G2pwd.ellipseSize(H,Sij,GB)
1174            Sgam = 1.8*wave/(np.pi*costh*lenR)
1175        #microstrain               
1176        if calcControls[phfx+'MustrainType'] == 'isotropic':
1177            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1178        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1179            H = np.array(refl[:3])
1180            P = np.array(calcControls[phfx+'MustrainAxis'])
1181            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1182            Si = parmDict[phfx+'Mustrain;i']
1183            Sa = parmDict[phfx+'Mustrain;a']
1184            Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1185        else:       #generalized - P.W. Stephens model
1186            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1187            Sum = 0
1188            for i,strm in enumerate(Strms):
1189                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1190            Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*np.sqrt(Sum)/np.pi
1191    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1192        #crystallite size
1193        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1194            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
1195        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1196            H = np.array(refl[:3])
1197            P = np.array(calcControls[phfx+'SizeAxis'])
1198            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1199            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1200            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1201        else:           #ellipsoidal crystallites   #OK
1202            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1203            H = np.array(refl[:3])
1204            lenR = G2pwd.ellipseSize(H,Sij,GB)
1205            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/lenR
1206        #microstrain               
1207        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1208            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
1209        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1210            H = np.array(refl[:3])
1211            P = np.array(calcControls[phfx+'MustrainAxis'])
1212            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1213            Si = parmDict[phfx+'Mustrain;i']
1214            Sa = parmDict[phfx+'Mustrain;a']
1215            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1216        else:       #generalized - P.W. Stephens model  OK
1217            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1218            Sum = 0
1219            for i,strm in enumerate(Strms):
1220                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1221            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4]**3
1222           
1223    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1224    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1225    sig /= ateln2
1226    return sig,gam
1227       
1228def GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1229    'Needs a doc string'
1230    gamDict = {}
1231    sigDict = {}
1232    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1233        costh = cosd(refl[5]/2.)
1234        tanth = tand(refl[5]/2.)
1235        #crystallite size derivatives
1236        if calcControls[phfx+'SizeType'] == 'isotropic':
1237            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1238            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1239            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1240        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1241            H = np.array(refl[:3])
1242            P = np.array(calcControls[phfx+'SizeAxis'])
1243            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1244            Si = parmDict[phfx+'Size;i']
1245            Sa = parmDict[phfx+'Size;a']
1246            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1247            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1248            Sgam = gami*sqtrm
1249            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1250            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1251            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1252            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1253            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1254            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1255        else:           #ellipsoidal crystallites
1256            const = 1.8*wave/(np.pi*costh)
1257            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1258            H = np.array(refl[:3])
1259            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1260            Sgam = const/lenR
1261            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1262                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1263                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1264        gamDict[phfx+'Size;mx'] = Sgam
1265        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1266               
1267        #microstrain derivatives               
1268        if calcControls[phfx+'MustrainType'] == 'isotropic':
1269            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1270            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1271            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1272        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1273            H = np.array(refl[:3])
1274            P = np.array(calcControls[phfx+'MustrainAxis'])
1275            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1276            Si = parmDict[phfx+'Mustrain;i']
1277            Sa = parmDict[phfx+'Mustrain;a']
1278            gami = 0.018*Si*Sa*tanth/np.pi
1279            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1280            Mgam = gami/sqtrm
1281            dsi = -gami*Si*cosP**2/sqtrm**3
1282            dsa = -gami*Sa*sinP**2/sqtrm**3
1283            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1284            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1285            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1286            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1287        else:       #generalized - P.W. Stephens model
1288            const = 0.018*refl[4]**2*tanth/np.pi
1289            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1290            Sum = 0
1291            for i,strm in enumerate(Strms):
1292                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1293                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1294                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1295            Mgam = const*np.sqrt(Sum)
1296            for i in range(len(Strms)):
1297                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1298                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1299        gamDict[phfx+'Mustrain;mx'] = Mgam
1300        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1301    else:   #'T'OF - All checked & OK
1302        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1303            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
1304            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1305            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1306        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1307            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
1308            H = np.array(refl[:3])
1309            P = np.array(calcControls[phfx+'SizeAxis'])
1310            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1311            Si = parmDict[phfx+'Size;i']
1312            Sa = parmDict[phfx+'Size;a']
1313            gami = const/(Si*Sa)
1314            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1315            Sgam = gami*sqtrm
1316            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1317            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1318            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1319            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1320            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1321            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1322        else:           #OK  ellipsoidal crystallites
1323            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
1324            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1325            H = np.array(refl[:3])
1326            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1327            Sgam = const/lenR
1328            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1329                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1330                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1331        gamDict[phfx+'Size;mx'] = Sgam  #OK
1332        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1333               
1334        #microstrain derivatives               
1335        if calcControls[phfx+'MustrainType'] == 'isotropic':
1336            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
1337            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1338            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1339        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1340            H = np.array(refl[:3])
1341            P = np.array(calcControls[phfx+'MustrainAxis'])
1342            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1343            Si = parmDict[phfx+'Mustrain;i']
1344            Sa = parmDict[phfx+'Mustrain;a']
1345            gami = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa
1346            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1347            Mgam = gami/sqtrm
1348            dsi = -gami*Si*cosP**2/sqtrm**3
1349            dsa = -gami*Sa*sinP**2/sqtrm**3
1350            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1351            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1352            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1353            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1354        else:       #generalized - P.W. Stephens model OK
1355            pwrs = calcControls[phfx+'MuPwrs']
1356            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1357            const = 1.e-6*parmDict[hfx+'difC']*refl[4]**3
1358            Sum = 0
1359            for i,strm in enumerate(Strms):
1360                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1361                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1362                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1363            Mgam = const*np.sqrt(Sum)
1364            for i in range(len(Strms)):
1365                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1366                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1367        gamDict[phfx+'Mustrain;mx'] = Mgam
1368        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1369       
1370    return sigDict,gamDict
1371       
1372def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1373    'Needs a doc string'
1374    h,k,l = refl[:3]
1375    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1376    d = np.sqrt(dsq)
1377
1378    refl[4] = d
1379    if 'C' in calcControls[hfx+'histType']:
1380        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1381        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1382        if 'Bragg' in calcControls[hfx+'instType']:
1383            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1384                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1385        else:               #Debye-Scherrer - simple but maybe not right
1386            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1387    elif 'T' in calcControls[hfx+'histType']:
1388        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1389        #do I need sample position effects - maybe?
1390    return pos
1391
1392def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1393    'Needs a doc string'
1394    dpr = 180./np.pi
1395    h,k,l = refl[:3]
1396    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1397    dst = np.sqrt(dstsq)
1398    dsp = 1./dst
1399    if 'C' in calcControls[hfx+'histType']:
1400        pos = refl[5]-parmDict[hfx+'Zero']
1401        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1402        dpdw = const*dst
1403        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1404        dpdA *= const*wave/(2.0*dst)
1405        dpdZ = 1.0
1406        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1407        if 'Bragg' in calcControls[hfx+'instType']:
1408            dpdSh = -4.*const*cosd(pos/2.0)
1409            dpdTr = -const*sind(pos)*100.0
1410            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1411        else:               #Debye-Scherrer - simple but maybe not right
1412            dpdXd = -const*cosd(pos)
1413            dpdYd = -const*sind(pos)
1414            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1415    elif 'T' in calcControls[hfx+'histType']:
1416        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1417        dpdZ = 1.0
1418        dpdDC = dsp
1419        dpdDA = dsp**2
1420        dpdDB = 1./dsp
1421        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
1422           
1423def GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict):
1424    'Needs a doc string'
1425    laue = SGData['SGLaue']
1426    uniq = SGData['SGUniq']
1427    h,k,l = refl[:3]
1428    if laue in ['m3','m3m']:
1429        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1430            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1431    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1432        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1433    elif laue in ['3R','3mR']:
1434        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1435    elif laue in ['4/m','4/mmm']:
1436        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1437    elif laue in ['mmm']:
1438        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1439    elif laue in ['2/m']:
1440        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1441        if uniq == 'a':
1442            Dij += parmDict[phfx+'D23']*k*l
1443        elif uniq == 'b':
1444            Dij += parmDict[phfx+'D13']*h*l
1445        elif uniq == 'c':
1446            Dij += parmDict[phfx+'D12']*h*k
1447    else:
1448        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1449            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1450    if 'C' in calcControls[hfx+'histType']:
1451        return -180.*Dij*refl[4]**2*tand(refl[5]/2.0)/np.pi
1452    else:
1453        return -Dij*parmDict[hfx+'difC']*refl[4]**2/2.
1454           
1455def GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict):
1456    'Needs a doc string'
1457    laue = SGData['SGLaue']
1458    uniq = SGData['SGUniq']
1459    h,k,l = refl[:3]
1460    if laue in ['m3','m3m']:
1461        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1462            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1463    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1464        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1465    elif laue in ['3R','3mR']:
1466        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1467    elif laue in ['4/m','4/mmm']:
1468        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1469    elif laue in ['mmm']:
1470        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1471    elif laue in ['2/m']:
1472        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1473        if uniq == 'a':
1474            dDijDict[phfx+'D23'] = k*l
1475        elif uniq == 'b':
1476            dDijDict[phfx+'D13'] = h*l
1477        elif uniq == 'c':
1478            dDijDict[phfx+'D12'] = h*k
1479            names.append()
1480    else:
1481        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1482            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1483    if 'C' in calcControls[hfx+'histType']:
1484        for item in dDijDict:
1485            dDijDict[item] *= -180.0*refl[4]**2*tand(refl[5]/2.0)/np.pi
1486    else:
1487        for item in dDijDict:
1488            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4]**2/2.
1489    return dDijDict
1490   
1491def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1492    'Needs a doc string'
1493    histoList = Histograms.keys()
1494    histoList.sort()
1495    for histogram in histoList:
1496        if 'PWDR' in histogram[:4]:
1497            Histogram = Histograms[histogram]
1498            hId = Histogram['hId']
1499            hfx = ':%d:'%(hId)
1500            Limits = calcControls[hfx+'Limits']
1501            if 'C' in calcControls[hfx+'histType']:
1502                shl = max(parmDict[hfx+'SH/L'],0.0005)
1503                Ka2 = False
1504                kRatio = 0.0
1505                if hfx+'Lam1' in parmDict.keys():
1506                    Ka2 = True
1507                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1508                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
1509            x,y,w,yc,yb,yd = Histogram['Data']
1510            xB = np.searchsorted(x,Limits[0])
1511            xF = np.searchsorted(x,Limits[1])
1512            ymb = np.array(y-yb)
1513            ymb = np.where(ymb,ymb,1.0)
1514            ycmb = np.array(yc-yb)
1515            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1516            refLists = Histogram['Reflection Lists']
1517            for phase in refLists:
1518                Phase = Phases[phase]
1519                pId = Phase['pId']
1520                phfx = '%d:%d:'%(pId,hId)
1521                refDict = refLists[phase]
1522                sumFo = 0.0
1523                sumdF = 0.0
1524                sumFosq = 0.0
1525                sumdFsq = 0.0
1526                for refl in refDict['RefList']:
1527                    if 'C' in calcControls[hfx+'histType']:
1528                        yp = np.zeros_like(yb)
1529                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1530                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1531                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1532                        iFin2 = iFin
1533                        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
1534                        if Ka2:
1535                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1536                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1537                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1538                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1539                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1540                                continue
1541                            elif not iBeg2-iFin2:     #peak above high limit - done
1542                                break
1543                            yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
1544                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0))
1545                    elif 'T' in calcControls[hfx+'histType']:
1546                        yp = np.zeros_like(yb)
1547                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1548                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1549                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1550                        if iBeg < iFin:
1551                            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
1552                            refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
1553                    Fo = np.sqrt(np.abs(refl[8]))
1554                    Fc = np.sqrt(np.abs(refl[9]))
1555                    sumFo += Fo
1556                    sumFosq += refl[8]**2
1557                    sumdF += np.abs(Fo-Fc)
1558                    sumdFsq += (refl[8]-refl[9])**2
1559                Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1560                Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1561                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1562                Histogram['Residuals']['hId'] = hId
1563        elif 'HKLF' in histogram[:4]:
1564            Histogram = Histograms[histogram]
1565            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1566               
1567def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1568    'Needs a doc string'
1569   
1570    def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
1571        U = parmDict[hfx+'U']
1572        V = parmDict[hfx+'V']
1573        W = parmDict[hfx+'W']
1574        X = parmDict[hfx+'X']
1575        Y = parmDict[hfx+'Y']
1576        tanPos = tand(refl[5]/2.0)
1577        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1578        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1579        sig = max(0.001,sig)
1580        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1581        gam = max(0.001,gam)
1582        return sig,gam
1583               
1584    def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
1585        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4]**2+   \
1586            parmDict[hfx+'sig-2']*refl[4]**4+parmDict[hfx+'sig-q']/refl[4]**2
1587        gam = parmDict[hfx+'X']*refl[4]+parmDict[hfx+'Y']*refl[4]**2
1588        Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1589        sig += Ssig
1590        gam += Sgam
1591        return sig,gam
1592       
1593    def GetReflAlpBet(refl,hfx,parmDict):
1594        alp = parmDict[hfx+'alpha']/refl[4]
1595        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]**4+parmDict[hfx+'beta-q']/refl[4]**2
1596        return alp,bet
1597               
1598    hId = Histogram['hId']
1599    hfx = ':%d:'%(hId)
1600    bakType = calcControls[hfx+'bakType']
1601    yb = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1602    yc = np.zeros_like(yb)
1603    cw = np.diff(x)
1604    cw = np.append(cw,cw[-1])
1605       
1606    if 'C' in calcControls[hfx+'histType']:   
1607        shl = max(parmDict[hfx+'SH/L'],0.002)
1608        Ka2 = False
1609        if hfx+'Lam1' in parmDict.keys():
1610            wave = parmDict[hfx+'Lam1']
1611            Ka2 = True
1612            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1613            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1614        else:
1615            wave = parmDict[hfx+'Lam']
1616    for phase in Histogram['Reflection Lists']:
1617        refDict = Histogram['Reflection Lists'][phase]
1618        Phase = Phases[phase]
1619        pId = Phase['pId']
1620        pfx = '%d::'%(pId)
1621        phfx = '%d:%d:'%(pId,hId)
1622        hfx = ':%d:'%(hId)
1623        SGData = Phase['General']['SGData']
1624        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1625        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #Do I want to modify by Dij?
1626        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1627        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1628        Vst = np.sqrt(nl.det(G))    #V*
1629        if not Phase['General'].get('doPawley'):
1630            time0 = time.time()
1631            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1632#            print 'sf calc time: %.3fs'%(time.time()-time0)
1633        time0 = time.time()
1634        badPeak = False
1635        for iref,refl in enumerate(refDict['RefList']):
1636            if 'C' in calcControls[hfx+'histType']:
1637                h,k,l = refl[:3]
1638                Uniq = np.inner(refl[:3],SGMT)
1639                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1640                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1641                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1642                refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1643                refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1644                refl[11] *= Vst*Lorenz
1645                 
1646                if Phase['General'].get('doPawley'):
1647                    try:
1648                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1649                        refl[9] = parmDict[pInd]
1650                    except KeyError:
1651#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1652                        continue
1653                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1654                iBeg = np.searchsorted(x,refl[5]-fmin)
1655                iFin = np.searchsorted(x,refl[5]+fmax)
1656                if not iBeg+iFin:       #peak below low limit - skip peak
1657                    continue
1658                elif not iBeg-iFin:     #peak above high limit - done
1659                    break
1660                elif iBeg > iFin:   #bad peak coeff - skip
1661                    badPeak = True
1662                    continue
1663                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
1664                if Ka2:
1665                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1666                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1667                    iBeg = np.searchsorted(x,pos2-fmin)
1668                    iFin = np.searchsorted(x,pos2+fmax)
1669                    if not iBeg+iFin:       #peak below low limit - skip peak
1670                        continue
1671                    elif not iBeg-iFin:     #peak above high limit - done
1672                        return yc,yb
1673                    elif iBeg > iFin:   #bad peak coeff - skip
1674                        continue
1675                    yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1676            elif 'T' in calcControls[hfx+'histType']:
1677                h,k,l = refl[:3]
1678                Uniq = np.inner(refl[:3],SGMT)
1679                refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict)         #corrected reflection position
1680                Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4]**4                                                #TOF Lorentz correction
1681                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
1682                refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
1683                refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
1684                refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1685                refl[11] *= Vst*Lorenz
1686                if Phase['General'].get('doPawley'):
1687                    try:
1688                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1689                        refl[9] = parmDict[pInd]
1690                    except KeyError:
1691#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1692                        continue
1693                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1694                iBeg = np.searchsorted(x,refl[5]-fmin)
1695                iFin = np.searchsorted(x,refl[5]+fmax)
1696                if not iBeg+iFin:       #peak below low limit - skip peak
1697                    continue
1698                elif not iBeg-iFin:     #peak above high limit - done
1699                    break
1700                elif iBeg > iFin:   #bad peak coeff - skip
1701                    badPeak = True
1702                    continue
1703                yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
1704#        print 'profile calc time: %.3fs'%(time.time()-time0)
1705    if badPeak:
1706        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
1707    return yc,yb
1708   
1709def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1710    'Needs a doc string'
1711   
1712    def cellVaryDerv(pfx,SGData,dpdA): 
1713        if SGData['SGLaue'] in ['-1',]:
1714            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1715                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1716        elif SGData['SGLaue'] in ['2/m',]:
1717            if SGData['SGUniq'] == 'a':
1718                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1719            elif SGData['SGUniq'] == 'b':
1720                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1721            else:
1722                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1723        elif SGData['SGLaue'] in ['mmm',]:
1724            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1725        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1726            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1727        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1728            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1729        elif SGData['SGLaue'] in ['3R', '3mR']:
1730            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1731        elif SGData['SGLaue'] in ['m3m','m3']:
1732            return [[pfx+'A0',dpdA[0]]]
1733           
1734    # create a list of dependent variables and set up a dictionary to hold their derivatives
1735    dependentVars = G2mv.GetDependentVars()
1736    depDerivDict = {}
1737    for j in dependentVars:
1738        depDerivDict[j] = np.zeros(shape=(len(x)))
1739    #print 'dependent vars',dependentVars
1740    lenX = len(x)               
1741    hId = Histogram['hId']
1742    hfx = ':%d:'%(hId)
1743    bakType = calcControls[hfx+'bakType']
1744    dMdv = np.zeros(shape=(len(varylist),len(x)))
1745    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1746    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
1747        bBpos =varylist.index(hfx+'Back;0')
1748        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1749    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1750    for name in varylist:
1751        if 'Debye' in name:
1752            id = int(name.split(';')[-1])
1753            parm = name[:int(name.rindex(';'))]
1754            ip = names.index(parm)
1755            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1756    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1757    for name in varylist:
1758        if 'BkPk' in name:
1759            parm,id = name.split(';')
1760            id = int(id)
1761            if parm in names:
1762                ip = names.index(parm)
1763                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1764    cw = np.diff(x)
1765    cw = np.append(cw,cw[-1])
1766    Ka2 = False #also for TOF!
1767    if 'C' in calcControls[hfx+'histType']:   
1768        shl = max(parmDict[hfx+'SH/L'],0.002)
1769        if hfx+'Lam1' in parmDict.keys():
1770            wave = parmDict[hfx+'Lam1']
1771            Ka2 = True
1772            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1773            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1774        else:
1775            wave = parmDict[hfx+'Lam']
1776    for phase in Histogram['Reflection Lists']:
1777        refDict = Histogram['Reflection Lists'][phase]
1778        Phase = Phases[phase]
1779        SGData = Phase['General']['SGData']
1780        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1781        pId = Phase['pId']
1782        pfx = '%d::'%(pId)
1783        phfx = '%d:%d:'%(pId,hId)
1784        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #And modify here by Dij? - no
1785        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1786        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1787        if not Phase['General'].get('doPawley'):
1788            time0 = time.time()
1789            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1790#            print 'sf-derv time %.3fs'%(time.time()-time0)
1791            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1792        time0 = time.time()
1793        for iref,refl in enumerate(refDict['RefList']):
1794            h,k,l = refl[:3]
1795            Uniq = np.inner(refl[:3],SGMT)
1796            if 'T' in calcControls[hfx+'histType']:
1797                wave = refl[14]
1798            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1799            if 'C' in calcControls[hfx+'histType']:        #CW powder
1800                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1801            else: #'T'OF
1802                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1803            iBeg = np.searchsorted(x,refl[5]-fmin)
1804            iFin = np.searchsorted(x,refl[5]+fmax)
1805            if not iBeg+iFin:       #peak below low limit - skip peak
1806                continue
1807            elif not iBeg-iFin:     #peak above high limit - done
1808                break
1809            pos = refl[5]
1810            if 'C' in calcControls[hfx+'histType']:
1811                tanth = tand(pos/2.0)
1812                costh = cosd(pos/2.0)
1813                lenBF = iFin-iBeg
1814                dMdpk = np.zeros(shape=(6,lenBF))
1815                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1816                for i in range(5):
1817                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
1818                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1819                if Ka2:
1820                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1821                    iBeg2 = np.searchsorted(x,pos2-fmin)
1822                    iFin2 = np.searchsorted(x,pos2+fmax)
1823                    if iBeg2-iFin2:
1824                        lenBF2 = iFin2-iBeg2
1825                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1826                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1827                        for i in range(5):
1828                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
1829                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
1830                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1831            else:   #'T'OF
1832                lenBF = iFin-iBeg
1833                if lenBF < 0:   #bad peak coeff
1834                    break
1835                dMdpk = np.zeros(shape=(6,lenBF))
1836                dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1837                for i in range(6):
1838                    dMdpk[i] += refl[11]*refl[9]*dMdipk[i]      #cw[iBeg:iFin]*
1839                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
1840            if Phase['General'].get('doPawley'):
1841                dMdpw = np.zeros(len(x))
1842                try:
1843                    pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1844                    idx = varylist.index(pIdx)
1845                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1846                    if Ka2: #not for TOF either
1847                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1848                    dMdv[idx] = dMdpw
1849                except: # ValueError:
1850                    pass
1851            if 'C' in calcControls[hfx+'histType']:
1852                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1853                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1854                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1855                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1856                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1857                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1858                    hfx+'DisplaceY':[dpdY,'pos'],}
1859                if 'Bragg' in calcControls[hfx+'instType']:
1860                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1861                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
1862                else:
1863                    names.update({hfx+'Absorption':[dFdAb,'int'],})
1864            else:   #'T'OF
1865                dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1866                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1867                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1868                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1869                    hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1870                    hfx+'beta-q':[1./refl[4]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'],
1871                    hfx+'sig-2':[refl[4]**4,'sig'],hfx+'sig-q':[1./refl[4]**2,'sig'],
1872                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1873            for name in names:
1874                item = names[name]
1875                if name in varylist:
1876                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1877                    if Ka2:
1878                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1879                elif name in dependentVars:
1880                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1881                    if Ka2:
1882                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1883            for iPO in dIdPO:
1884                if iPO in varylist:
1885                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1886                    if Ka2:
1887                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1888                elif iPO in dependentVars:
1889                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1890                    if Ka2:
1891                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1892            for i,name in enumerate(['omega','chi','phi']):
1893                aname = pfx+'SH '+name
1894                if aname in varylist:
1895                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1896                    if Ka2:
1897                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1898                elif aname in dependentVars:
1899                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1900                    if Ka2:
1901                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1902            for iSH in dFdODF:
1903                if iSH in varylist:
1904                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1905                    if Ka2:
1906                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1907                elif iSH in dependentVars:
1908                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1909                    if Ka2:
1910                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1911            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1912            for name,dpdA in cellDervNames:
1913                if name in varylist:
1914                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1915                    if Ka2:
1916                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1917                elif name in dependentVars:
1918                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1919                    if Ka2:
1920                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1921            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict)
1922            for name in dDijDict:
1923                if name in varylist:
1924                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1925                    if Ka2:
1926                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1927                elif name in dependentVars:
1928                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1929                    if Ka2:
1930                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1931            if 'C' in calcControls[hfx+'histType']:
1932                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1933            else:   #'T'OF
1934                sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1935            for name in gamDict:
1936                if name in varylist:
1937                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1938                    if Ka2:
1939                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1940                elif name in dependentVars:
1941                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1942                    if Ka2:
1943                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1944            for name in sigDict:
1945                if name in varylist:
1946                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1947                    if Ka2:
1948                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1949                elif name in dependentVars:
1950                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1951                    if Ka2:
1952                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1953            for name in ['BabA','BabU']:
1954                if refl[9]:
1955                    if phfx+name in varylist:
1956                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1957                        if Ka2:
1958                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1959                    elif phfx+name in dependentVars:                   
1960                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1961                        if Ka2:
1962                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
1963            if not Phase['General'].get('doPawley'):
1964                #do atom derivatives -  for RB,F,X & U so far             
1965                corr = dervDict['int']/refl[9]
1966                if Ka2:
1967                    corr2 = dervDict2['int']/refl[9]
1968                for name in varylist+dependentVars:
1969                    if '::RBV;' in name:
1970                        pass
1971                    else:
1972                        try:
1973                            aname = name.split(pfx)[1][:2]
1974                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1975                        except IndexError:
1976                            continue
1977                    if name in varylist:
1978                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1979                        if Ka2:
1980                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1981                    elif name in dependentVars:
1982                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1983                        if Ka2:
1984                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1985    #        print 'profile derv time: %.3fs'%(time.time()-time0)
1986    # now process derivatives in constraints
1987    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1988    return dMdv
1989   
1990def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1991    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1992    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
1993    are returned to either dervRefine or HessRefine.
1994
1995    :returns:
1996    '''
1997    nobs = Histogram['Residuals']['Nobs']
1998    hId = Histogram['hId']
1999    hfx = ':%d:'%(hId)
2000    pfx = '%d::'%(Phase['pId'])
2001    phfx = '%d:%d:'%(Phase['pId'],hId)
2002    SGData = Phase['General']['SGData']
2003    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2004    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2005    refDict = Histogram['Data']
2006    dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2007    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2008    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2009    dependentVars = G2mv.GetDependentVars()
2010    depDerivDict = {}
2011    for j in dependentVars:
2012        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2013    wdf = np.zeros(len(refDict['RefList']))
2014    if calcControls['F**2']:
2015        for iref,ref in enumerate(refDict['RefList']):
2016            if ref[6] > 0:
2017                dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1] 
2018                w = 1.0/ref[6]
2019                if w*ref[5] >= calcControls['minF/sig']:
2020                    wdf[iref] = w*(ref[5]-ref[7])
2021                    for j,var in enumerate(varylist):
2022                        if var in dFdvDict:
2023                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
2024                    for var in dependentVars:
2025                        if var in dFdvDict:
2026                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
2027                    if phfx+'Scale' in varylist:
2028                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]
2029                    elif phfx+'Scale' in dependentVars:
2030                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11]
2031                    for item in ['Ep','Es','Eg']:
2032                        if phfx+item in varylist and phfx+item in dervDict:
2033                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #OK
2034                        elif phfx+item in dependentVars and phfx+item in dervDict:
2035                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]  #OK
2036                    for item in ['BabA','BabU']:
2037                        if phfx+item in varylist:
2038                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
2039                        elif phfx+item in dependentVars:
2040                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
2041    else:
2042        for iref,ref in enumerate(refDict['RefList']):
2043            if ref[5] > 0.:
2044                dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2045                Fo = np.sqrt(ref[5])
2046                Fc = np.sqrt(ref[7])
2047                w = 1.0/ref[6]
2048                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
2049                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
2050                    for j,var in enumerate(varylist):
2051                        if var in dFdvDict:
2052                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
2053                    for var in dependentVars:
2054                        if var in dFdvDict:
2055                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
2056                    if phfx+'Scale' in varylist:
2057                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]
2058                    elif phfx+'Scale' in dependentVars:
2059                        depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11]                           
2060                    for item in ['Ep','Es','Eg']:
2061                        if phfx+item in varylist and phfx+item in dervDict:
2062                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #correct
2063                        elif phfx+item in dependentVars and phfx+item in dervDict:
2064                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
2065                    for item in ['BabA','BabU']:
2066                        if phfx+item in varylist:
2067                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
2068                        elif phfx+item in dependentVars:
2069                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
2070    return dMdvh,depDerivDict,wdf
2071   
2072
2073def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2074    '''Loop over histograms and compute derivatives of the fitting
2075    model (M) with respect to all parameters.  Results are returned in
2076    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2077    n is the number of parameters and m is the number of data
2078    points. This can exceed memory when m gets large. This routine is
2079    used when refinement derivatives are selected as "analtytic
2080    Jacobian" in Controls.
2081
2082    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2083    '''
2084    parmDict.update(zip(varylist,values))
2085    G2mv.Dict2Map(parmDict,varylist)
2086    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2087    nvar = len(varylist)
2088    dMdv = np.empty(0)
2089    histoList = Histograms.keys()
2090    histoList.sort()
2091    for histogram in histoList:
2092        if 'PWDR' in histogram[:4]:
2093            Histogram = Histograms[histogram]
2094            hId = Histogram['hId']
2095            hfx = ':%d:'%(hId)
2096            wtFactor = calcControls[hfx+'wtFactor']
2097            Limits = calcControls[hfx+'Limits']
2098            x,y,w,yc,yb,yd = Histogram['Data']
2099            xB = np.searchsorted(x,Limits[0])
2100            xF = np.searchsorted(x,Limits[1])
2101            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2102                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
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            hfx = ':%d:'%(Histogram['hId'])
2109            wtFactor = calcControls[hfx+'wtFactor']
2110            # now process derivatives in constraints
2111            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2112        else:
2113            continue        #skip non-histogram entries
2114        if len(dMdv):
2115            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2116        else:
2117            dMdv = np.sqrt(wtFactor)*dMdvh
2118           
2119    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2120    if np.any(pVals):
2121        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2122        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2123       
2124    return dMdv
2125
2126def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2127    '''Loop over histograms and compute derivatives of the fitting
2128    model (M) with respect to all parameters.  For each histogram, the
2129    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2130    number of parameters and m is the number of data points *in the
2131    histogram*. The (n by n) Hessian is computed from each Jacobian
2132    and it is returned.  This routine is used when refinement
2133    derivatives are selected as "analtytic Hessian" in Controls.
2134
2135    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2136    '''
2137    parmDict.update(zip(varylist,values))
2138    G2mv.Dict2Map(parmDict,varylist)
2139    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2140    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2141    nvar = len(varylist)
2142    Hess = np.empty(0)
2143    histoList = Histograms.keys()
2144    histoList.sort()
2145    for histogram in histoList:
2146        if 'PWDR' in histogram[:4]:
2147            Histogram = Histograms[histogram]
2148            hId = Histogram['hId']
2149            hfx = ':%d:'%(hId)
2150            wtFactor = calcControls[hfx+'wtFactor']
2151            Limits = calcControls[hfx+'Limits']
2152            x,y,w,yc,yb,yd = Histogram['Data']
2153            W = wtFactor*w
2154            dy = y-yc
2155            xB = np.searchsorted(x,Limits[0])
2156            xF = np.searchsorted(x,Limits[1])
2157            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2158                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2159            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2160            Dy = dy[xB:xF][np.newaxis,:]
2161            dMdvh *= Wt
2162            if dlg:
2163                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2164            if len(Hess):
2165                Hess += np.inner(dMdvh,dMdvh)
2166                dMdvh *= Wt*Dy
2167                Vec += np.sum(dMdvh,axis=1)
2168            else:
2169                Hess = np.inner(dMdvh,dMdvh)
2170                dMdvh *= Wt*Dy
2171                Vec = np.sum(dMdvh,axis=1)
2172        elif 'HKLF' in histogram[:4]:
2173            Histogram = Histograms[histogram]
2174            phase = Histogram['Reflection Lists']
2175            Phase = Phases[phase]
2176            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2177            hId = Histogram['hId']
2178            hfx = ':%d:'%(Histogram['hId'])
2179            wtFactor = calcControls[hfx+'wtFactor']
2180            # now process derivatives in constraints
2181            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2182#            print 'matrix build time: %.3f'%(time.time()-time0)
2183
2184            if dlg:
2185                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2186            if len(Hess):
2187                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2188                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2189            else:
2190                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2191                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2192        else:
2193            continue        #skip non-histogram entries
2194    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2195    if np.any(pVals):
2196        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2197        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2198        Hess += np.inner(dpdv*pWt,dpdv)
2199    return Vec,Hess
2200
2201def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
2202    'Needs a doc string'
2203    Values2Dict(parmDict, varylist, values)
2204    G2mv.Dict2Map(parmDict,varylist)
2205    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2206    M = np.empty(0)
2207    SumwYo = 0
2208    Nobs = 0
2209    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2210    histoList = Histograms.keys()
2211    histoList.sort()
2212    for histogram in histoList:
2213        if 'PWDR' in histogram[:4]:
2214            Histogram = Histograms[histogram]
2215            hId = Histogram['hId']
2216            hfx = ':%d:'%(hId)
2217            wtFactor = calcControls[hfx+'wtFactor']
2218            Limits = calcControls[hfx+'Limits']
2219            x,y,w,yc,yb,yd = Histogram['Data']
2220            yc *= 0.0                           #zero full calcd profiles
2221            yb *= 0.0
2222            yd *= 0.0
2223            xB = np.searchsorted(x,Limits[0])
2224            xF = np.searchsorted(x,Limits[1])
2225            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2226                varylist,Histogram,Phases,calcControls,pawleyLookup)
2227            yc[xB:xF] += yb[xB:xF]
2228            if not np.any(y):                   #fill dummy data
2229                rv = st.poisson(yc[xB:xF])
2230                y[xB:xF] = rv.rvs()
2231                Z = np.ones_like(yc[xB:xF])
2232                Z[1::2] *= -1
2233                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2234                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2235            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2236            W = wtFactor*w
2237            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2238            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2239            Nobs += Histogram['Residuals']['Nobs']
2240            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2241            SumwYo += Histogram['Residuals']['sumwYo']
2242            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2243            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2244            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2245            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2246            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.))
2247            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.))
2248            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2249            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2250            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
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,wdy))
2254#end of PWDR processing
2255        elif 'HKLF' in histogram[:4]:
2256            Histogram = Histograms[histogram]
2257            Histogram['Residuals'] = {}
2258            phase = Histogram['Reflection Lists']
2259            Phase = Phases[phase]
2260            hId = Histogram['hId']
2261            hfx = ':%d:'%(hId)
2262            wtFactor = calcControls[hfx+'wtFactor']
2263            pfx = '%d::'%(Phase['pId'])
2264            phfx = '%d:%d:'%(Phase['pId'],hId)
2265            SGData = Phase['General']['SGData']
2266            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2267            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2268            refDict = Histogram['Data']
2269            time0 = time.time()
2270            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2271#            print 'sf-calc time: %.3f'%(time.time()-time0)
2272            df = np.zeros(len(refDict['RefList']))
2273            sumwYo = 0
2274            sumFo = 0
2275            sumFo2 = 0
2276            sumdF = 0
2277            sumdF2 = 0
2278            nobs = 0
2279            if calcControls['F**2']:
2280                for i,ref in enumerate(refDict['RefList']):
2281                    if ref[6] > 0:
2282                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2283                        w = 1.0/ref[6]
2284                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]  #correct Fc^2 for extinction
2285                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2286                        if w*ref[5] >= calcControls['minF/sig']:
2287                            Fo = np.sqrt(ref[5])
2288                            sumFo += Fo
2289                            sumFo2 += ref[5]
2290                            sumdF += abs(Fo-np.sqrt(ref[7]))
2291                            sumdF2 += abs(ref[5]-ref[7])
2292                            nobs += 1
2293                            df[i] = -w*(ref[5]-ref[7])
2294                            sumwYo += (w*ref[5])**2
2295            else:
2296                for i,ref in enumerate(refDict['RefList']):
2297                    if ref[5] > 0.:
2298                        ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2299                        ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]    #correct Fc^2 for extinction
2300                        ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2301                        Fo = np.sqrt(ref[5])
2302                        Fc = np.sqrt(ref[7])
2303                        w = 2.0*Fo/ref[6]
2304                        if w*Fo >= calcControls['minF/sig']:
2305                            sumFo += Fo
2306                            sumFo2 += ref[5]
2307                            sumdF += abs(Fo-Fc)
2308                            sumdF2 += abs(ref[5]-ref[7])
2309                            nobs += 1
2310                            df[i] = -w*(Fo-Fc)
2311                            sumwYo += (w*Fo)**2
2312            Histogram['Residuals']['Nobs'] = nobs
2313            Histogram['Residuals']['sumwYo'] = sumwYo
2314            SumwYo += sumwYo
2315            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
2316            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2317            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2318            Histogram['Residuals'][phfx+'Nref'] = nobs
2319            Nobs += nobs
2320            if dlg:
2321                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2322            M = np.concatenate((M,wtFactor*df))
2323# end of HKLF processing
2324    Histograms['sumwYo'] = SumwYo
2325    Histograms['Nobs'] = Nobs
2326    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
2327    if dlg:
2328        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2329        if not GoOn:
2330            parmDict['saved values'] = values
2331            dlg.Destroy()
2332            raise Exception         #Abort!!
2333    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2334    if len(pVals):
2335        pSum = np.sum(pWt*pVals**2)
2336        for name in pWsum:
2337            if pWsum:
2338                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
2339        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2340        Nobs += len(pVals)
2341        M = np.concatenate((M,np.sqrt(pWt)*pVals))
2342    return M
2343                       
Note: See TracBrowser for help on using the repository browser.