Context Navigation

source: trunk/GSASIIstrMath.py @ 1482

Last change on this file since 1482 was 1482, checked in by vondreele, 9 years ago
change comment about sign of Dij application fix TOF derivatives for Ai = met. tensor elements fix TOF mustrain derivatives - all now OK work on TOF size derivatives
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 115.8 KB

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

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

Download in other formats: