source: trunk/GSASIIstrMath.py @ 1812

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

make new ranId for phases imported from gpx files
eliminate all raise Exceptions from Refine & seqRefine - now gives a ErrorMessage? popup upon e.g. user abort.

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