source: trunk/GSASIIstrMath.py @ 1478

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

Fix CW & TOF hydrostatic strain calculations & change disply formats
work on TOF size/mustrain calcs.

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