Context Navigation

source: trunk/GSASIIstrMath.py @ 1878

Last change on this file since 1878 was 1878, checked in by vondreele, 8 years ago
refactor DDataGUI - mostly move event routines to be inside respective sizer routines Enable Flack parameter - function OK; derivatives need work Allow inversion of noncentrosymmetric structures in SymOpDialog? - to test enantiomers some work to make pdf stuff neutron TOF friendly - not complete fix Debye background function - now works; refactor result printing for it
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 141.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: