source: trunk/GSASIIstrMath.py @ 1776

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

remove Copy, Copy flags buttons from Data display; add them as menu items under edit
Clean up DDataGUI
use the IRAD value on the old IPARM file to set teh x-ray source, e.g. CuKa?
fix SC extinction - now more closely matches old GSAS

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