Context Navigation

source: trunk/GSASIIstrMath.py @ 1498

Last change on this file since 1498 was 1498, checked in by vondreele, 9 years ago
try to implement drag/drop of tree items - commented out as it doesn't work add 'ShowCell?' to Controls to avoid cell errors for sequential single peak fitting when normal sequential refinement has been done fix a neg. peak width error
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 116.4 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2014-09-16 21:20:57 +0000 (Tue, 16 Sep 2014) $
8	# $Author: vondreele $
9	# $Revision: 1498 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 1498 2014-09-16 21:20:57Z 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: 1498 $")
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 for al2O3!
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	dervDict = {}
805	if calcControls[phfx+'EType'] != 'None':
806	SQ = 1/(4.ref[4]*2)
807	if 'C' in parmDict[hfx+'Type']:
808	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
809	else: #'T'
810	cos2T = 1.0-2.SQref[12]**2 #cos(2theta)
811	if 'SXC' in parmDict[hfx+'Type']:
812	AV = 7.9406e5/parmDict[pfx+'Vol']**2
813	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
814	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
815	PLZ = AVP12ref[7]parmDict[hfx+'Lam']*2
816	elif 'SNT' in parmDict[hfx+'Type']:
817	AV = 1.e7/parmDict[pfx+'Vol']**2
818	PL = SQ
819	PLZ = AVref[7]ref[12]**2
820	elif 'SNC' in parmDict[hfx+'Type']:
821	AV = 1.e7/parmDict[pfx+'Vol']**2
822	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
823	PLZ = AVref[9]parmDict[hfx+'Lam']**2 #Fcsq as per GSAS, why not FcTsq (ref[9])?
824
825	if 'Primary' in calcControls[phfx+'EType']:
826	PLZ *= 1.5
827	else:
828	if 'C' in parmDict[hfx+'Type']:
829	PLZ *= calcControls[phfx+'Tbar']
830	else: #'T'
831	PLZ *= ref[13] #t-bar
832	if 'Primary' in calcControls[phfx+'EType']:
833	PLZ *= 1.5
834	PSIG = parmDict[phfx+'Ep']
835	elif 'I & II' in calcControls[phfx+'EType']:
836	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
837	elif 'Type II' in calcControls[phfx+'EType']:
838	PSIG = parmDict[phfx+'Es']
839	else: # 'Secondary Type I'
840	PSIG = parmDict[phfx+'Eg']/PL
841
842	AG = 0.58+0.48cos2T+0.24cos2T**2
843	AL = 0.025+0.285*cos2T
844	BG = 0.02-0.025*cos2T
845	BL = 0.15-0.2(0.75-cos2T)*2
846	if cos2T < 0.:
847	BL = -0.45*cos2T
848	CG = 2.
849	CL = 2.
850	PF = PLZ*PSIG
851
852	if 'Gaussian' in calcControls[phfx+'EApprox']:
853	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
854	extCor = np.sqrt(PF4)
855	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
856	else:
857	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
858	extCor = np.sqrt(PF4)
859	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
860
861	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
862	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
863	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
864	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
865	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
866	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
867
868	return 1./extCor,dervDict
869
870	def Dict2Values(parmdict, varylist):
871	'''Use before call to leastsq to setup list of values for the parameters
872	in parmdict, as selected by key in varylist'''
873	return [parmdict[key] for key in varylist]
874
875	def Values2Dict(parmdict, varylist, values):
876	''' Use after call to leastsq to update the parameter dictionary with
877	values corresponding to keys in varylist'''
878	parmdict.update(zip(varylist,values))
879
880	def GetNewCellParms(parmDict,varyList):
881	'Needs a doc string'
882	newCellDict = {}
883	Anames = ['A'+str(i) for i in range(6)]
884	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
885	for item in varyList:
886	keys = item.split(':')
887	if keys[2] in Ddict:
888	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
889	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
890	newCellDict[parm] = [key,parmDict[key]-parmDict[item]]
891	return newCellDict # is e.g. {'0::D11':A0-D11}
892
893	def ApplyXYZshifts(parmDict,varyList):
894	'''
895	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
896
897	:param dict parmDict: parameter dictionary
898	:param list varyList: list of variables (not used!)
899	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
900
901	'''
902	newAtomDict = {}
903	for item in parmDict:
904	if 'dA' in item:
905	parm = ''.join(item.split('d'))
906	parmDict[parm] += parmDict[item]
907	newAtomDict[item] = [parm,parmDict[parm]]
908	return newAtomDict
909
910	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
911	'Spherical harmonics texture'
912	IFCoup = 'Bragg' in calcControls[hfx+'instType']
913	if 'T' in calcControls[hfx+'histType']:
914	tth = parmDict[hfx+'2-theta']
915	else:
916	tth = refl[5]
917	odfCor = 1.0
918	H = refl[:3]
919	cell = G2lat.Gmat2cell(g)
920	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
921	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
922	phi,beta = G2lat.CrsAng(H,cell,SGData)
923	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
924	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
925	for item in SHnames:
926	L,M,N = eval(item.strip('C'))
927	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
928	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
929	Lnorm = G2lat.Lnorm(L)
930	odfCor += parmDict[pfx+item]LnormKcl*Ksl
931	return odfCor
932
933	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
934	'Spherical harmonics texture derivatives'
935	if 'T' in calcControls[hfx+'histType']:
936	tth = parmDict[hfx+'2-theta']
937	else:
938	tth = refl[5]
939	FORPI = 4.0*np.pi
940	IFCoup = 'Bragg' in calcControls[hfx+'instType']
941	odfCor = 1.0
942	dFdODF = {}
943	dFdSA = [0,0,0]
944	H = refl[:3]
945	cell = G2lat.Gmat2cell(g)
946	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
947	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
948	phi,beta = G2lat.CrsAng(H,cell,SGData)
949	psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
950	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
951	for item in SHnames:
952	L,M,N = eval(item.strip('C'))
953	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
954	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
955	Lnorm = G2lat.Lnorm(L)
956	odfCor += parmDict[pfx+item]LnormKcl*Ksl
957	dFdODF[pfx+item] = LnormKclKsl
958	for i in range(3):
959	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
960	return odfCor,dFdODF,dFdSA
961
962	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
963	'spherical harmonics preferred orientation (cylindrical symmetry only)'
964	if 'T' in calcControls[hfx+'histType']:
965	tth = parmDict[hfx+'2-theta']
966	else:
967	tth = refl[5]
968	odfCor = 1.0
969	H = refl[:3]
970	cell = G2lat.Gmat2cell(g)
971	Sangl = [0.,0.,0.]
972	if 'Bragg' in calcControls[hfx+'instType']:
973	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
974	IFCoup = True
975	else:
976	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
977	IFCoup = False
978	phi,beta = G2lat.CrsAng(H,cell,SGData)
979	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
980	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
981	for item in SHnames:
982	L,N = eval(item.strip('C'))
983	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
984	odfCor += parmDict[phfx+item]LnormKcsl
985	return np.squeeze(odfCor)
986
987	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
988	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
989	if 'T' in calcControls[hfx+'histType']:
990	tth = parmDict[hfx+'2-theta']
991	else:
992	tth = refl[5]
993	FORPI = 12.5663706143592
994	odfCor = 1.0
995	dFdODF = {}
996	H = refl[:3]
997	cell = G2lat.Gmat2cell(g)
998	Sangl = [0.,0.,0.]
999	if 'Bragg' in calcControls[hfx+'instType']:
1000	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1001	IFCoup = True
1002	else:
1003	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1004	IFCoup = False
1005	phi,beta = G2lat.CrsAng(H,cell,SGData)
1006	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1007	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1008	for item in SHnames:
1009	L,N = eval(item.strip('C'))
1010	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
1011	odfCor += parmDict[phfx+item]LnormKcsl
1012	dFdODF[phfx+item] = Kcsl*Lnorm
1013	return odfCor,dFdODF
1014
1015	def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1016	'March-Dollase preferred orientation correction'
1017	POcorr = 1.0
1018	MD = parmDict[phfx+'MD']
1019	if MD != 1.0:
1020	MDAxis = calcControls[phfx+'MDAxis']
1021	sumMD = 0
1022	for H in uniq:
1023	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1024	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1025	sumMD += A**3
1026	POcorr = sumMD/len(uniq)
1027	return POcorr
1028
1029	def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1030	'Needs a doc string'
1031	POcorr = 1.0
1032	POderv = {}
1033	if calcControls[phfx+'poType'] == 'MD':
1034	MD = parmDict[phfx+'MD']
1035	MDAxis = calcControls[phfx+'MDAxis']
1036	sumMD = 0
1037	sumdMD = 0
1038	for H in uniq:
1039	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1040	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1041	sumMD += A**3
1042	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1043	POcorr = sumMD/len(uniq)
1044	POderv[phfx+'MD'] = sumdMD/len(uniq)
1045	else: #spherical harmonics
1046	if calcControls[phfx+'SHord']:
1047	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1048	return POcorr,POderv
1049
1050	def GetAbsorb(refl,hfx,calcControls,parmDict):
1051	'Needs a doc string'
1052	if 'Debye' in calcControls[hfx+'instType']:
1053	if 'T' in calcControls[hfx+'histType']:
1054	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1055	else:
1056	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1057	else:
1058	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1059
1060	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1061	'Needs a doc string'
1062	if 'Debye' in calcControls[hfx+'instType']:
1063	if 'T' in calcControls[hfx+'histType']:
1064	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
1065	else:
1066	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1067	else:
1068	return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]))
1069
1070	def GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict):
1071	'Needs a doc string'
1072	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1073	pi2 = np.sqrt(2./np.pi)
1074	if 'T' in calcControls[hfx+'histType']:
1075	sth2 = sind(parmDict[hfx+'2-theta']/2.)**2
1076	wave = refl[14]
1077	else: #'C'W
1078	sth2 = sind(refl[5]/2.)**2
1079	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1080	c2th = 1.-2.0*sth2
1081	flv2 = refl[9](wave/parmDict[pfx+'Vol'])*2
1082	if 'X' in calcControls[hfx+'histType']:
1083	flv2 = 0.079411(1.0+c2th**2)/2.0
1084	xfac = flv2*parmDict[phfx+'Extinction']
1085	exb = 1.0
1086	if xfac > -1.:
1087	exb = 1./(1.+xfac)
1088	exl = 1.0
1089	if 0 < xfac <= 1.:
1090	xn = np.array([xfac**(i+1) for i in range(6)])
1091	exl = np.sum(xn*coef)
1092	elif xfac > 1.:
1093	xfac2 = 1./np.sqrt(xfac)
1094	exl = pi2(1.-0.125/xfac)xfac2
1095	return exbsth2+exl(1.-sth2)
1096
1097	def GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict):
1098	'Needs a doc string'
1099	delt = 0.001
1100	parmDict[phfx+'Extinction'] += delt
1101	plus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1102	parmDict[phfx+'Extinction'] -= 2.*delt
1103	minus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1104	parmDict[phfx+'Extinction'] += delt
1105	return (plus-minus)/(2.*delt)
1106
1107	def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1108	'Needs a doc string' #need powder extinction!
1109	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
1110	if 'X' in parmDict[hfx+'Type']:
1111	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1112	POcorr = 1.0
1113	if pfx+'SHorder' in parmDict: #generalized spherical harmonics texture
1114	POcorr = SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1115	elif calcControls[phfx+'poType'] == 'MD': #March-Dollase
1116	POcorr = GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1117	elif calcControls[phfx+'SHord']: #cylindrical spherical harmonics
1118	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1119	Icorr *= POcorr
1120	AbsCorr = 1.0
1121	AbsCorr = GetAbsorb(refl,hfx,calcControls,parmDict)
1122	Icorr *= AbsCorr
1123	ExtCorr = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
1124	Icorr *= ExtCorr
1125	return Icorr,POcorr,AbsCorr,ExtCorr
1126
1127	def GetIntensityDerv(refl,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1128	'Needs a doc string' #need powder extinction derivs!
1129	dIdsh = 1./parmDict[hfx+'Scale']
1130	dIdsp = 1./parmDict[phfx+'Scale']
1131	if 'X' in parmDict[hfx+'Type']:
1132	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1133	dIdPola /= pola
1134	else: #'N'
1135	dIdPola = 0.0
1136	dFdODF = {}
1137	dFdSA = [0,0,0]
1138	dIdPO = {}
1139	if pfx+'SHorder' in parmDict:
1140	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1141	for iSH in dFdODF:
1142	dFdODF[iSH] /= odfCor
1143	for i in range(3):
1144	dFdSA[i] /= odfCor
1145	elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1146	POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1147	for iPO in dIdPO:
1148	dIdPO[iPO] /= POcorr
1149	if 'T' in parmDict[hfx+'Type']:
1150	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[16] #wave/abs corr
1151	dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[17] #/ext corr
1152	else:
1153	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[13] #wave/abs corr
1154	dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[14] #/ext corr
1155	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1156
1157	def GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1158	'Needs a doc string'
1159	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1160	costh = cosd(refl[5]/2.)
1161	#crystallite size
1162	if calcControls[phfx+'SizeType'] == 'isotropic':
1163	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1164	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1165	H = np.array(refl[:3])
1166	P = np.array(calcControls[phfx+'SizeAxis'])
1167	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1168	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1169	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1170	else: #ellipsoidal crystallites
1171	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1172	H = np.array(refl[:3])
1173	lenR = G2pwd.ellipseSize(H,Sij,GB)
1174	Sgam = 1.8wave/(np.picosth*lenR)
1175	#microstrain
1176	if calcControls[phfx+'MustrainType'] == 'isotropic':
1177	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1178	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1179	H = np.array(refl[:3])
1180	P = np.array(calcControls[phfx+'MustrainAxis'])
1181	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1182	Si = parmDict[phfx+'Mustrain;i']
1183	Sa = parmDict[phfx+'Mustrain;a']
1184	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1185	else: #generalized - P.W. Stephens model
1186	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1187	Sum = 0
1188	for i,strm in enumerate(Strms):
1189	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1190	Mgam = 0.018refl[4]2tand(refl[5]/2.)*np.sqrt(Sum)/np.pi
1191	elif 'T' in calcControls[hfx+'histType']: #All checked & OK
1192	#crystallite size
1193	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1194	Sgam = 1.e-4parmDict[hfx+'difC']refl[4]**2/parmDict[phfx+'Size;i']
1195	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1196	H = np.array(refl[:3])
1197	P = np.array(calcControls[phfx+'SizeAxis'])
1198	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1199	Sgam = 1.e-4parmDict[hfx+'difC']refl[4]*2/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a'])
1200	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1201	else: #ellipsoidal crystallites #OK
1202	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1203	H = np.array(refl[:3])
1204	lenR = G2pwd.ellipseSize(H,Sij,GB)
1205	Sgam = 1.e-4parmDict[hfx+'difC']refl[4]**2/lenR
1206	#microstrain
1207	if calcControls[phfx+'MustrainType'] == 'isotropic': #OK
1208	Mgam = 1.e-6parmDict[hfx+'difC']refl[4]*parmDict[phfx+'Mustrain;i']
1209	elif calcControls[phfx+'MustrainType'] == 'uniaxial': #OK
1210	H = np.array(refl[:3])
1211	P = np.array(calcControls[phfx+'MustrainAxis'])
1212	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1213	Si = parmDict[phfx+'Mustrain;i']
1214	Sa = parmDict[phfx+'Mustrain;a']
1215	Mgam = 1.e-6parmDict[hfx+'difC']refl[4]SiSa/np.sqrt((SicosP)2+(SasinP)**2)
1216	else: #generalized - P.W. Stephens model OK
1217	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1218	Sum = 0
1219	for i,strm in enumerate(Strms):
1220	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1221	Mgam = 1.e-6parmDict[hfx+'difC']np.sqrt(Sum)refl[4]*3
1222
1223	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1224	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1225	sig /= ateln2
1226	return sig,gam
1227
1228	def GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1229	'Needs a doc string'
1230	gamDict = {}
1231	sigDict = {}
1232	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1233	costh = cosd(refl[5]/2.)
1234	tanth = tand(refl[5]/2.)
1235	#crystallite size derivatives
1236	if calcControls[phfx+'SizeType'] == 'isotropic':
1237	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1238	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1239	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1240	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1241	H = np.array(refl[:3])
1242	P = np.array(calcControls[phfx+'SizeAxis'])
1243	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1244	Si = parmDict[phfx+'Size;i']
1245	Sa = parmDict[phfx+'Size;a']
1246	gami = 1.8wave/(costhnp.piSiSa)
1247	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1248	Sgam = gami*sqtrm
1249	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1250	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1251	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1252	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1253	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1254	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1255	else: #ellipsoidal crystallites
1256	const = 1.8wave/(np.picosth)
1257	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1258	H = np.array(refl[:3])
1259	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1260	Sgam = const/lenR
1261	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1262	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1263	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1264	gamDict[phfx+'Size;mx'] = Sgam
1265	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1266
1267	#microstrain derivatives
1268	if calcControls[phfx+'MustrainType'] == 'isotropic':
1269	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1270	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1271	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1272	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1273	H = np.array(refl[:3])
1274	P = np.array(calcControls[phfx+'MustrainAxis'])
1275	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1276	Si = parmDict[phfx+'Mustrain;i']
1277	Sa = parmDict[phfx+'Mustrain;a']
1278	gami = 0.018SiSa*tanth/np.pi
1279	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1280	Mgam = gami/sqtrm
1281	dsi = -gamiSicosP2/sqtrm3
1282	dsa = -gamiSasinP2/sqtrm3
1283	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1284	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1285	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1286	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1287	else: #generalized - P.W. Stephens model
1288	const = 0.018refl[4]2tanth/np.pi
1289	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1290	Sum = 0
1291	for i,strm in enumerate(Strms):
1292	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1293	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1294	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1295	Mgam = const*np.sqrt(Sum)
1296	for i in range(len(Strms)):
1297	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1298	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1299	gamDict[phfx+'Mustrain;mx'] = Mgam
1300	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1301	else: #'T'OF - All checked & OK
1302	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1303	Sgam = 1.e-4parmDict[hfx+'difC']refl[4]**2/parmDict[phfx+'Size;i']
1304	gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1305	sigDict[phfx+'Size;i'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])*2/(ateln2parmDict[phfx+'Size;i'])
1306	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1307	const = 1.e-4parmDict[hfx+'difC']refl[4]**2
1308	H = np.array(refl[:3])
1309	P = np.array(calcControls[phfx+'SizeAxis'])
1310	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1311	Si = parmDict[phfx+'Size;i']
1312	Sa = parmDict[phfx+'Size;a']
1313	gami = const/(Si*Sa)
1314	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1315	Sgam = gami*sqtrm
1316	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1317	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1318	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1319	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1320	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1321	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1322	else: #OK ellipsoidal crystallites
1323	const = 1.e-4parmDict[hfx+'difC']refl[4]**2
1324	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1325	H = np.array(refl[:3])
1326	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1327	Sgam = const/lenR
1328	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1329	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1330	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1331	gamDict[phfx+'Size;mx'] = Sgam #OK
1332	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2 #OK
1333
1334	#microstrain derivatives
1335	if calcControls[phfx+'MustrainType'] == 'isotropic':
1336	Mgam = 1.e-6parmDict[hfx+'difC']refl[4]*parmDict[phfx+'Mustrain;i']
1337	gamDict[phfx+'Mustrain;i'] = 1.e-6refl[4]parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] #OK
1338	sigDict[phfx+'Mustrain;i'] = 2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])*2/(ateln2parmDict[phfx+'Mustrain;i'])
1339	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1340	H = np.array(refl[:3])
1341	P = np.array(calcControls[phfx+'MustrainAxis'])
1342	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1343	Si = parmDict[phfx+'Mustrain;i']
1344	Sa = parmDict[phfx+'Mustrain;a']
1345	gami = 1.e-6parmDict[hfx+'difC']refl[4]SiSa
1346	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1347	Mgam = gami/sqtrm
1348	dsi = -gamiSicosP2/sqtrm3
1349	dsa = -gamiSasinP2/sqtrm3
1350	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1351	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1352	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1353	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1354	else: #generalized - P.W. Stephens model OK
1355	pwrs = calcControls[phfx+'MuPwrs']
1356	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1357	const = 1.e-6parmDict[hfx+'difC']refl[4]**3
1358	Sum = 0
1359	for i,strm in enumerate(Strms):
1360	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1361	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1362	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1363	Mgam = const*np.sqrt(Sum)
1364	for i in range(len(Strms)):
1365	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1366	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1367	gamDict[phfx+'Mustrain;mx'] = Mgam
1368	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1369
1370	return sigDict,gamDict
1371
1372	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1373	'Needs a doc string'
1374	h,k,l = refl[:3]
1375	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1376	d = np.sqrt(dsq)
1377
1378	refl[4] = d
1379	if 'C' in calcControls[hfx+'histType']:
1380	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1381	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1382	if 'Bragg' in calcControls[hfx+'instType']:
1383	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1384	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1385	else: #Debye-Scherrer - simple but maybe not right
1386	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1387	elif 'T' in calcControls[hfx+'histType']:
1388	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1389	#do I need sample position effects - maybe?
1390	return pos
1391
1392	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1393	'Needs a doc string'
1394	dpr = 180./np.pi
1395	h,k,l = refl[:3]
1396	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1397	dst = np.sqrt(dstsq)
1398	dsp = 1./dst
1399	if 'C' in calcControls[hfx+'histType']:
1400	pos = refl[5]-parmDict[hfx+'Zero']
1401	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1402	dpdw = const*dst
1403	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1404	dpdA = constwave/(2.0*dst)
1405	dpdZ = 1.0
1406	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1407	if 'Bragg' in calcControls[hfx+'instType']:
1408	dpdSh = -4.constcosd(pos/2.0)
1409	dpdTr = -constsind(pos)100.0
1410	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1411	else: #Debye-Scherrer - simple but maybe not right
1412	dpdXd = -const*cosd(pos)
1413	dpdYd = -const*sind(pos)
1414	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1415	elif 'T' in calcControls[hfx+'histType']:
1416	dpdA = -np.array([h2,k2,l*2,hk,hl,kl])parmDict[hfx+'difC']dsp**3/2.
1417	dpdZ = 1.0
1418	dpdDC = dsp
1419	dpdDA = dsp**2
1420	dpdDB = 1./dsp
1421	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
1422
1423	def GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict):
1424	'Needs a doc string'
1425	laue = SGData['SGLaue']
1426	uniq = SGData['SGUniq']
1427	h,k,l = refl[:3]
1428	if laue in ['m3','m3m']:
1429	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
1430	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
1431	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1432	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
1433	elif laue in ['3R','3mR']:
1434	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
1435	elif laue in ['4/m','4/mmm']:
1436	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
1437	elif laue in ['mmm']:
1438	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1439	elif laue in ['2/m']:
1440	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1441	if uniq == 'a':
1442	Dij += parmDict[phfx+'D23']kl
1443	elif uniq == 'b':
1444	Dij += parmDict[phfx+'D13']hl
1445	elif uniq == 'c':
1446	Dij += parmDict[phfx+'D12']hk
1447	else:
1448	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
1449	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
1450	if 'C' in calcControls[hfx+'histType']:
1451	return -180.Dijrefl[4]*2tand(refl[5]/2.0)/np.pi
1452	else:
1453	return -DijparmDict[hfx+'difC']refl[4]**2/2.
1454
1455	def GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict):
1456	'Needs a doc string'
1457	laue = SGData['SGLaue']
1458	uniq = SGData['SGUniq']
1459	h,k,l = refl[:3]
1460	if laue in ['m3','m3m']:
1461	dDijDict = {phfx+'D11':h2+k2+l**2,
1462	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
1463	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1464	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
1465	elif laue in ['3R','3mR']:
1466	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
1467	elif laue in ['4/m','4/mmm']:
1468	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
1469	elif laue in ['mmm']:
1470	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1471	elif laue in ['2/m']:
1472	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1473	if uniq == 'a':
1474	dDijDict[phfx+'D23'] = k*l
1475	elif uniq == 'b':
1476	dDijDict[phfx+'D13'] = h*l
1477	elif uniq == 'c':
1478	dDijDict[phfx+'D12'] = h*k
1479	names.append()
1480	else:
1481	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
1482	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
1483	if 'C' in calcControls[hfx+'histType']:
1484	for item in dDijDict:
1485	dDijDict[item] = -180.0refl[4]*2tand(refl[5]/2.0)/np.pi
1486	else:
1487	for item in dDijDict:
1488	dDijDict[item] = -parmDict[hfx+'difC']refl[4]**2/2.
1489	return dDijDict
1490
1491	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1492	'Needs a doc string'
1493	histoList = Histograms.keys()
1494	histoList.sort()
1495	for histogram in histoList:
1496	if 'PWDR' in histogram[:4]:
1497	Histogram = Histograms[histogram]
1498	hId = Histogram['hId']
1499	hfx = ':%d:'%(hId)
1500	Limits = calcControls[hfx+'Limits']
1501	if 'C' in calcControls[hfx+'histType']:
1502	shl = max(parmDict[hfx+'SH/L'],0.0005)
1503	Ka2 = False
1504	kRatio = 0.0
1505	if hfx+'Lam1' in parmDict.keys():
1506	Ka2 = True
1507	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1508	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1509	x,y,w,yc,yb,yd = Histogram['Data']
1510	xB = np.searchsorted(x,Limits[0])
1511	xF = np.searchsorted(x,Limits[1])
1512	ymb = np.array(y-yb)
1513	ymb = np.where(ymb,ymb,1.0)
1514	ycmb = np.array(yc-yb)
1515	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
1516	refLists = Histogram['Reflection Lists']
1517	for phase in refLists:
1518	Phase = Phases[phase]
1519	pId = Phase['pId']
1520	phfx = '%d:%d:'%(pId,hId)
1521	refDict = refLists[phase]
1522	sumFo = 0.0
1523	sumdF = 0.0
1524	sumFosq = 0.0
1525	sumdFsq = 0.0
1526	for refl in refDict['RefList']:
1527	if 'C' in calcControls[hfx+'histType']:
1528	yp = np.zeros_like(yb)
1529	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1530	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1531	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1532	iFin2 = iFin
1533	if not iBeg+iFin: #peak below low limit - skip peak
1534	continue
1535	elif not iBeg-iFin: #peak above high limit - done
1536	break
1537	elif iBeg < iFin:
1538	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
1539	if Ka2:
1540	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1541	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1542	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1543	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1544	yp[iBeg2:iFin2] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2])) #and here
1545	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11](1.+kRatio)),0.0))
1546	elif 'T' in calcControls[hfx+'histType']:
1547	yp = np.zeros_like(yb)
1548	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1549	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1550	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1551	if iBeg < iFin:
1552	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
1553	refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
1554	Fo = np.sqrt(np.abs(refl[8]))
1555	Fc = np.sqrt(np.abs(refl[9]))
1556	sumFo += Fo
1557	sumFosq += refl[8]**2
1558	sumdF += np.abs(Fo-Fc)
1559	sumdFsq += (refl[8]-refl[9])**2
1560	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1561	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1562	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1563	Histogram['Residuals']['hId'] = hId
1564	elif 'HKLF' in histogram[:4]:
1565	Histogram = Histograms[histogram]
1566	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1567
1568	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1569	'Needs a doc string'
1570
1571	def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
1572	U = parmDict[hfx+'U']
1573	V = parmDict[hfx+'V']
1574	W = parmDict[hfx+'W']
1575	X = parmDict[hfx+'X']
1576	Y = parmDict[hfx+'Y']
1577	tanPos = tand(refl[5]/2.0)
1578	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1579	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
1580	sig = max(0.001,sig)
1581	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
1582	gam = max(0.001,gam)
1583	return sig,gam
1584
1585	def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
1586	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4]*2+ \
1587	parmDict[hfx+'sig-2']refl[4]4+parmDict[hfx+'sig-q']/refl[4]*2
1588	gam = parmDict[hfx+'X']refl[4]+parmDict[hfx+'Y']refl[4]**2
1589	Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1590	sig += Ssig
1591	gam += Sgam
1592	return sig,gam
1593
1594	def GetReflAlpBet(refl,hfx,parmDict):
1595	alp = parmDict[hfx+'alpha']/refl[4]
1596	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]4+parmDict[hfx+'beta-q']/refl[4]2
1597	return alp,bet
1598
1599	hId = Histogram['hId']
1600	hfx = ':%d:'%(hId)
1601	bakType = calcControls[hfx+'bakType']
1602	yb = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1603	yc = np.zeros_like(yb)
1604	cw = np.diff(x)
1605	cw = np.append(cw,cw[-1])
1606
1607	if 'C' in calcControls[hfx+'histType']:
1608	shl = max(parmDict[hfx+'SH/L'],0.002)
1609	Ka2 = False
1610	if hfx+'Lam1' in parmDict.keys():
1611	wave = parmDict[hfx+'Lam1']
1612	Ka2 = True
1613	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1614	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1615	else:
1616	wave = parmDict[hfx+'Lam']
1617	for phase in Histogram['Reflection Lists']:
1618	refDict = Histogram['Reflection Lists'][phase]
1619	Phase = Phases[phase]
1620	pId = Phase['pId']
1621	pfx = '%d::'%(pId)
1622	phfx = '%d:%d:'%(pId,hId)
1623	hfx = ':%d:'%(hId)
1624	SGData = Phase['General']['SGData']
1625	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1626	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #Do I want to modify by Dij?
1627	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1628	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1629	Vst = np.sqrt(nl.det(G)) #V*
1630	if not Phase['General'].get('doPawley'):
1631	time0 = time.time()
1632	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1633	# print 'sf calc time: %.3fs'%(time.time()-time0)
1634	time0 = time.time()
1635	badPeak = False
1636	for iref,refl in enumerate(refDict['RefList']):
1637	if 'C' in calcControls[hfx+'histType']:
1638	h,k,l = refl[:3]
1639	Uniq = np.inner(refl[:3],SGMT)
1640	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
1641	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
1642	refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
1643	refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
1644	refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1645	refl[11] = VstLorenz
1646
1647	if Phase['General'].get('doPawley'):
1648	try:
1649	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1650	refl[9] = parmDict[pInd]
1651	except KeyError:
1652	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1653	continue
1654	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1655	iBeg = np.searchsorted(x,refl[5]-fmin)
1656	iFin = np.searchsorted(x,refl[5]+fmax)
1657	if not iBeg+iFin: #peak below low limit - skip peak
1658	continue
1659	elif not iBeg-iFin: #peak above high limit - done
1660	break
1661	elif iBeg > iFin: #bad peak coeff - skip
1662	badPeak = True
1663	continue
1664	yc[iBeg:iFin] += refl[11]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
1665	if Ka2:
1666	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1667	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1668	iBeg = np.searchsorted(x,pos2-fmin)
1669	iFin = np.searchsorted(x,pos2+fmax)
1670	if not iBeg+iFin: #peak below low limit - skip peak
1671	continue
1672	elif not iBeg-iFin: #peak above high limit - done
1673	return yc,yb
1674	elif iBeg > iFin: #bad peak coeff - skip
1675	continue
1676	yc[iBeg:iFin] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here
1677	elif 'T' in calcControls[hfx+'histType']:
1678	h,k,l = refl[:3]
1679	Uniq = np.inner(refl[:3],SGMT)
1680	refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict) #corrected reflection position
1681	Lorenz = sind(parmDict[hfx+'2-theta']/2)refl[4]*4 #TOF Lorentz correction
1682	refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
1683	refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict) #peak sig & gam
1684	refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
1685	refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1686	refl[11] = VstLorenz
1687	if Phase['General'].get('doPawley'):
1688	try:
1689	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1690	refl[9] = parmDict[pInd]
1691	except KeyError:
1692	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1693	continue
1694	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1695	iBeg = np.searchsorted(x,refl[5]-fmin)
1696	iFin = np.searchsorted(x,refl[5]+fmax)
1697	if not iBeg+iFin: #peak below low limit - skip peak
1698	continue
1699	elif not iBeg-iFin: #peak above high limit - done
1700	break
1701	elif iBeg > iFin: #bad peak coeff - skip
1702	badPeak = True
1703	continue
1704	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]
1705	# print 'profile calc time: %.3fs'%(time.time()-time0)
1706	if badPeak:
1707	print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
1708	return yc,yb
1709
1710	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1711	'Needs a doc string'
1712
1713	def cellVaryDerv(pfx,SGData,dpdA):
1714	if SGData['SGLaue'] in ['-1',]:
1715	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1716	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1717	elif SGData['SGLaue'] in ['2/m',]:
1718	if SGData['SGUniq'] == 'a':
1719	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1720	elif SGData['SGUniq'] == 'b':
1721	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1722	else:
1723	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1724	elif SGData['SGLaue'] in ['mmm',]:
1725	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1726	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1727	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1728	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1729	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1730	elif SGData['SGLaue'] in ['3R', '3mR']:
1731	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
1732	elif SGData['SGLaue'] in ['m3m','m3']:
1733	return [[pfx+'A0',dpdA[0]]]
1734
1735	# create a list of dependent variables and set up a dictionary to hold their derivatives
1736	dependentVars = G2mv.GetDependentVars()
1737	depDerivDict = {}
1738	for j in dependentVars:
1739	depDerivDict[j] = np.zeros(shape=(len(x)))
1740	#print 'dependent vars',dependentVars
1741	lenX = len(x)
1742	hId = Histogram['hId']
1743	hfx = ':%d:'%(hId)
1744	bakType = calcControls[hfx+'bakType']
1745	dMdv = np.zeros(shape=(len(varylist),len(x)))
1746	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
1747	if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
1748	bBpos =varylist.index(hfx+'Back;0')
1749	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1750	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1751	for name in varylist:
1752	if 'Debye' in name:
1753	id = int(name.split(';')[-1])
1754	parm = name[:int(name.rindex(';'))]
1755	ip = names.index(parm)
1756	dMdv[varylist.index(name)] = dMddb[3*id+ip]
1757	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1758	for name in varylist:
1759	if 'BkPk' in name:
1760	parm,id = name.split(';')
1761	id = int(id)
1762	if parm in names:
1763	ip = names.index(parm)
1764	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1765	cw = np.diff(x)
1766	cw = np.append(cw,cw[-1])
1767	Ka2 = False #also for TOF!
1768	if 'C' in calcControls[hfx+'histType']:
1769	shl = max(parmDict[hfx+'SH/L'],0.002)
1770	if hfx+'Lam1' in parmDict.keys():
1771	wave = parmDict[hfx+'Lam1']
1772	Ka2 = True
1773	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1774	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1775	else:
1776	wave = parmDict[hfx+'Lam']
1777	for phase in Histogram['Reflection Lists']:
1778	refDict = Histogram['Reflection Lists'][phase]
1779	Phase = Phases[phase]
1780	SGData = Phase['General']['SGData']
1781	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1782	pId = Phase['pId']
1783	pfx = '%d::'%(pId)
1784	phfx = '%d:%d:'%(pId,hId)
1785	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #And modify here by Dij? - no
1786	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1787	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1788	if not Phase['General'].get('doPawley'):
1789	time0 = time.time()
1790	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1791	# print 'sf-derv time %.3fs'%(time.time()-time0)
1792	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1793	time0 = time.time()
1794	for iref,refl in enumerate(refDict['RefList']):
1795	h,k,l = refl[:3]
1796	Uniq = np.inner(refl[:3],SGMT)
1797	if 'T' in calcControls[hfx+'histType']:
1798	wave = refl[14]
1799	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1800	if 'C' in calcControls[hfx+'histType']: #CW powder
1801	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1802	else: #'T'OF
1803	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1804	iBeg = np.searchsorted(x,refl[5]-fmin)
1805	iFin = np.searchsorted(x,refl[5]+fmax)
1806	if not iBeg+iFin: #peak below low limit - skip peak
1807	continue
1808	elif not iBeg-iFin: #peak above high limit - done
1809	break
1810	pos = refl[5]
1811	if 'C' in calcControls[hfx+'histType']:
1812	tanth = tand(pos/2.0)
1813	costh = cosd(pos/2.0)
1814	lenBF = iFin-iBeg
1815	dMdpk = np.zeros(shape=(6,lenBF))
1816	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1817	for i in range(5):
1818	dMdpk[i] += 100.cw[iBeg:iFin]refl[11]refl[9]dMdipk[i]
1819	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1820	if Ka2:
1821	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
1822	iBeg2 = np.searchsorted(x,pos2-fmin)
1823	iFin2 = np.searchsorted(x,pos2+fmax)
1824	if iBeg2-iFin2:
1825	lenBF2 = iFin2-iBeg2
1826	dMdpk2 = np.zeros(shape=(6,lenBF2))
1827	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1828	for i in range(5):
1829	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11]refl[9]kRatio*dMdipk2[i]
1830	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11]*dMdipk2[0]
1831	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1832	else: #'T'OF
1833	lenBF = iFin-iBeg
1834	if lenBF < 0: #bad peak coeff
1835	break
1836	dMdpk = np.zeros(shape=(6,lenBF))
1837	dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1838	for i in range(6):
1839	dMdpk[i] += refl[11]refl[9]dMdipk[i] #cw[iBeg:iFin]*
1840	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1841	if Phase['General'].get('doPawley'):
1842	dMdpw = np.zeros(len(x))
1843	try:
1844	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1845	idx = varylist.index(pIdx)
1846	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1847	if Ka2: #not for TOF either
1848	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1849	dMdv[idx] = dMdpw
1850	except: # ValueError:
1851	pass
1852	if 'C' in calcControls[hfx+'histType']:
1853	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1854	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1855	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1856	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1857	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1858	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1859	hfx+'DisplaceY':[dpdY,'pos'],}
1860	if 'Bragg' in calcControls[hfx+'instType']:
1861	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1862	hfx+'SurfRoughB':[dFdAb[1],'int'],})
1863	else:
1864	names.update({hfx+'Absorption':[dFdAb,'int'],})
1865	else: #'T'OF
1866	dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1867	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1868	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1869	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1870	hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1871	hfx+'beta-q':[1./refl[4]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]2,'sig'],
1872	hfx+'sig-2':[refl[4]4,'sig'],hfx+'sig-q':[1./refl[4]2,'sig'],
1873	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1874	for name in names:
1875	item = names[name]
1876	if name in varylist:
1877	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1878	if Ka2:
1879	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1880	elif name in dependentVars:
1881	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1882	if Ka2:
1883	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1884	for iPO in dIdPO:
1885	if iPO in varylist:
1886	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1887	if Ka2:
1888	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1889	elif iPO in dependentVars:
1890	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1891	if Ka2:
1892	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1893	for i,name in enumerate(['omega','chi','phi']):
1894	aname = pfx+'SH '+name
1895	if aname in varylist:
1896	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1897	if Ka2:
1898	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1899	elif aname in dependentVars:
1900	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1901	if Ka2:
1902	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1903	for iSH in dFdODF:
1904	if iSH in varylist:
1905	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1906	if Ka2:
1907	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1908	elif iSH in dependentVars:
1909	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1910	if Ka2:
1911	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1912	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1913	for name,dpdA in cellDervNames:
1914	if name in varylist:
1915	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1916	if Ka2:
1917	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1918	elif name in dependentVars:
1919	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1920	if Ka2:
1921	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1922	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict)
1923	for name in dDijDict:
1924	if name in varylist:
1925	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1926	if Ka2:
1927	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1928	elif name in dependentVars:
1929	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1930	if Ka2:
1931	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1932	if 'C' in calcControls[hfx+'histType']:
1933	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1934	else: #'T'OF
1935	sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
1936	for name in gamDict:
1937	if name in varylist:
1938	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1939	if Ka2:
1940	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1941	elif name in dependentVars:
1942	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1943	if Ka2:
1944	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1945	for name in sigDict:
1946	if name in varylist:
1947	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1948	if Ka2:
1949	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1950	elif name in dependentVars:
1951	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1952	if Ka2:
1953	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1954	for name in ['BabA','BabU']:
1955	if refl[9]:
1956	if phfx+name in varylist:
1957	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1958	if Ka2:
1959	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1960	elif phfx+name in dependentVars:
1961	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1962	if Ka2:
1963	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1964	if not Phase['General'].get('doPawley'):
1965	#do atom derivatives - for RB,F,X & U so far
1966	corr = dervDict['int']/refl[9]
1967	if Ka2:
1968	corr2 = dervDict2['int']/refl[9]
1969	for name in varylist+dependentVars:
1970	if '::RBV;' in name:
1971	pass
1972	else:
1973	try:
1974	aname = name.split(pfx)[1][:2]
1975	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1976	except IndexError:
1977	continue
1978	if name in varylist:
1979	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1980	if Ka2:
1981	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1982	elif name in dependentVars:
1983	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1984	if Ka2:
1985	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1986	# print 'profile derv time: %.3fs'%(time.time()-time0)
1987	# now process derivatives in constraints
1988	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1989	return dMdv
1990
1991	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1992	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1993	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
1994	are returned to either dervRefine or HessRefine.
1995
1996	:returns:
1997	'''
1998	nobs = Histogram['Residuals']['Nobs']
1999	hId = Histogram['hId']
2000	hfx = ':%d:'%(hId)
2001	pfx = '%d::'%(Phase['pId'])
2002	phfx = '%d:%d:'%(Phase['pId'],hId)
2003	SGData = Phase['General']['SGData']
2004	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2005	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2006	refDict = Histogram['Data']
2007	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2008	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2009	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2010	dependentVars = G2mv.GetDependentVars()
2011	depDerivDict = {}
2012	for j in dependentVars:
2013	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2014	wdf = np.zeros(len(refDict['RefList']))
2015	if calcControls['F**2']:
2016	for iref,ref in enumerate(refDict['RefList']):
2017	if ref[6] > 0:
2018	dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2019	w = 1.0/ref[6]
2020	if w*ref[5] >= calcControls['minF/sig']:
2021	wdf[iref] = w*(ref[5]-ref[7])
2022	for j,var in enumerate(varylist):
2023	if var in dFdvDict:
2024	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11]
2025	for var in dependentVars:
2026	if var in dFdvDict:
2027	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11]
2028	if phfx+'Scale' in varylist:
2029	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]ref[11]
2030	elif phfx+'Scale' in dependentVars:
2031	depDerivDict[phfx+'Scale'][iref] = wref[9]ref[11]
2032	for item in ['Ep','Es','Eg']:
2033	if phfx+item in varylist and phfx+item in dervDict:
2034	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11] #OK
2035	elif phfx+item in dependentVars and phfx+item in dervDict:
2036	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11] #OK
2037	for item in ['BabA','BabU']:
2038	if phfx+item in varylist:
2039	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11]
2040	elif phfx+item in dependentVars:
2041	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11]
2042	else:
2043	for iref,ref in enumerate(refDict['RefList']):
2044	if ref[5] > 0.:
2045	dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
2046	Fo = np.sqrt(ref[5])
2047	Fc = np.sqrt(ref[7])
2048	w = 1.0/ref[6]
2049	if 2.0Fow*Fo >= calcControls['minF/sig']:
2050	wdf[iref] = 2.0Fow*(Fo-Fc)
2051	for j,var in enumerate(varylist):
2052	if var in dFdvDict:
2053	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11]
2054	for var in dependentVars:
2055	if var in dFdvDict:
2056	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11]
2057	if phfx+'Scale' in varylist:
2058	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]ref[11]
2059	elif phfx+'Scale' in dependentVars:
2060	depDerivDict[phfx+'Scale'][iref] = wref[9]ref[11]
2061	for item in ['Ep','Es','Eg']:
2062	if phfx+item in varylist and phfx+item in dervDict:
2063	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11] #correct
2064	elif phfx+item in dependentVars and phfx+item in dervDict:
2065	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
2066	for item in ['BabA','BabU']:
2067	if phfx+item in varylist:
2068	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11]
2069	elif phfx+item in dependentVars:
2070	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11]
2071	return dMdvh,depDerivDict,wdf
2072
2073
2074	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2075	'''Loop over histograms and compute derivatives of the fitting
2076	model (M) with respect to all parameters. Results are returned in
2077	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2078	n is the number of parameters and m is the number of data
2079	points. This can exceed memory when m gets large. This routine is
2080	used when refinement derivatives are selected as "analtytic
2081	Jacobian" in Controls.
2082
2083	:returns: Jacobian numpy.array dMdv for all histograms concatinated
2084	'''
2085	parmDict.update(zip(varylist,values))
2086	G2mv.Dict2Map(parmDict,varylist)
2087	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2088	nvar = len(varylist)
2089	dMdv = np.empty(0)
2090	histoList = Histograms.keys()
2091	histoList.sort()
2092	for histogram in histoList:
2093	if 'PWDR' in histogram[:4]:
2094	Histogram = Histograms[histogram]
2095	hId = Histogram['hId']
2096	hfx = ':%d:'%(hId)
2097	wtFactor = calcControls[hfx+'wtFactor']
2098	Limits = calcControls[hfx+'Limits']
2099	x,y,w,yc,yb,yd = Histogram['Data']
2100	xB = np.searchsorted(x,Limits[0])
2101	xF = np.searchsorted(x,Limits[1])
2102	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2103	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2104	elif 'HKLF' in histogram[:4]:
2105	Histogram = Histograms[histogram]
2106	phase = Histogram['Reflection Lists']
2107	Phase = Phases[phase]
2108	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2109	hfx = ':%d:'%(Histogram['hId'])
2110	wtFactor = calcControls[hfx+'wtFactor']
2111	# now process derivatives in constraints
2112	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2113	else:
2114	continue #skip non-histogram entries
2115	if len(dMdv):
2116	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2117	else:
2118	dMdv = np.sqrt(wtFactor)*dMdvh
2119
2120	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2121	if np.any(pVals):
2122	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2123	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2124
2125	return dMdv
2126
2127	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2128	'''Loop over histograms and compute derivatives of the fitting
2129	model (M) with respect to all parameters. For each histogram, the
2130	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2131	number of parameters and m is the number of data points *in the
2132	histogram*. The (n by n) Hessian is computed from each Jacobian
2133	and it is returned. This routine is used when refinement
2134	derivatives are selected as "analtytic Hessian" in Controls.
2135
2136	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2137	'''
2138	parmDict.update(zip(varylist,values))
2139	G2mv.Dict2Map(parmDict,varylist)
2140	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2141	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
2142	nvar = len(varylist)
2143	Hess = np.empty(0)
2144	histoList = Histograms.keys()
2145	histoList.sort()
2146	for histogram in histoList:
2147	if 'PWDR' in histogram[:4]:
2148	Histogram = Histograms[histogram]
2149	hId = Histogram['hId']
2150	hfx = ':%d:'%(hId)
2151	wtFactor = calcControls[hfx+'wtFactor']
2152	Limits = calcControls[hfx+'Limits']
2153	x,y,w,yc,yb,yd = Histogram['Data']
2154	W = wtFactor*w
2155	dy = y-yc
2156	xB = np.searchsorted(x,Limits[0])
2157	xF = np.searchsorted(x,Limits[1])
2158	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2159	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2160	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2161	Dy = dy[xB:xF][np.newaxis,:]
2162	dMdvh *= Wt
2163	if dlg:
2164	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2165	if len(Hess):
2166	Hess += np.inner(dMdvh,dMdvh)
2167	dMdvh = WtDy
2168	Vec += np.sum(dMdvh,axis=1)
2169	else:
2170	Hess = np.inner(dMdvh,dMdvh)
2171	dMdvh = WtDy
2172	Vec = np.sum(dMdvh,axis=1)
2173	elif 'HKLF' in histogram[:4]:
2174	Histogram = Histograms[histogram]
2175	phase = Histogram['Reflection Lists']
2176	Phase = Phases[phase]
2177	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2178	hId = Histogram['hId']
2179	hfx = ':%d:'%(Histogram['hId'])
2180	wtFactor = calcControls[hfx+'wtFactor']
2181	# now process derivatives in constraints
2182	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2183	# print 'matrix build time: %.3f'%(time.time()-time0)
2184
2185	if dlg:
2186	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2187	if len(Hess):
2188	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
2189	Hess += wtFactor*np.inner(dMdvh,dMdvh)
2190	else:
2191	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
2192	Hess = wtFactor*np.inner(dMdvh,dMdvh)
2193	else:
2194	continue #skip non-histogram entries
2195	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2196	if np.any(pVals):
2197	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2198	Vec += np.sum(dpdvpWtpVals,axis=1)
2199	Hess += np.inner(dpdv*pWt,dpdv)
2200	return Vec,Hess
2201
2202	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2203	'Needs a doc string'
2204	Values2Dict(parmDict, varylist, values)
2205	G2mv.Dict2Map(parmDict,varylist)
2206	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2207	M = np.empty(0)
2208	SumwYo = 0
2209	Nobs = 0
2210	ApplyRBModels(parmDict,Phases,rigidbodyDict)
2211	histoList = Histograms.keys()
2212	histoList.sort()
2213	for histogram in histoList:
2214	if 'PWDR' in histogram[:4]:
2215	Histogram = Histograms[histogram]
2216	hId = Histogram['hId']
2217	hfx = ':%d:'%(hId)
2218	wtFactor = calcControls[hfx+'wtFactor']
2219	Limits = calcControls[hfx+'Limits']
2220	x,y,w,yc,yb,yd = Histogram['Data']
2221	yc *= 0.0 #zero full calcd profiles
2222	yb *= 0.0
2223	yd *= 0.0
2224	xB = np.searchsorted(x,Limits[0])
2225	xF = np.searchsorted(x,Limits[1])
2226	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2227	varylist,Histogram,Phases,calcControls,pawleyLookup)
2228	yc[xB:xF] += yb[xB:xF]
2229	if not np.any(y): #fill dummy data
2230	rv = st.poisson(yc[xB:xF])
2231	y[xB:xF] = rv.rvs()
2232	Z = np.ones_like(yc[xB:xF])
2233	Z[1::2] *= -1
2234	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2235	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2236	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2237	W = wtFactor*w
2238	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2239	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2240	Nobs += Histogram['Residuals']['Nobs']
2241	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2242	SumwYo += Histogram['Residuals']['sumwYo']
2243	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2244	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2245	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2246	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2247	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.))
2248	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.))
2249	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2250	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2251	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2252	if dlg:
2253	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2254	M = np.concatenate((M,wdy))
2255	#end of PWDR processing
2256	elif 'HKLF' in histogram[:4]:
2257	Histogram = Histograms[histogram]
2258	Histogram['Residuals'] = {}
2259	phase = Histogram['Reflection Lists']
2260	Phase = Phases[phase]
2261	hId = Histogram['hId']
2262	hfx = ':%d:'%(hId)
2263	wtFactor = calcControls[hfx+'wtFactor']
2264	pfx = '%d::'%(Phase['pId'])
2265	phfx = '%d:%d:'%(Phase['pId'],hId)
2266	SGData = Phase['General']['SGData']
2267	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2268	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2269	refDict = Histogram['Data']
2270	time0 = time.time()
2271	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2272	# print 'sf-calc time: %.3f'%(time.time()-time0)
2273	df = np.zeros(len(refDict['RefList']))
2274	sumwYo = 0
2275	sumFo = 0
2276	sumFo2 = 0
2277	sumdF = 0
2278	sumdF2 = 0
2279	nobs = 0
2280	if calcControls['F**2']:
2281	for i,ref in enumerate(refDict['RefList']):
2282	if ref[6] > 0:
2283	ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2284	w = 1.0/ref[6]
2285	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinction
2286	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2287	if w*ref[5] >= calcControls['minF/sig']:
2288	Fo = np.sqrt(ref[5])
2289	sumFo += Fo
2290	sumFo2 += ref[5]
2291	sumdF += abs(Fo-np.sqrt(ref[7]))
2292	sumdF2 += abs(ref[5]-ref[7])
2293	nobs += 1
2294	df[i] = -w*(ref[5]-ref[7])
2295	sumwYo += (wref[5])*2
2296	else:
2297	for i,ref in enumerate(refDict['RefList']):
2298	if ref[5] > 0.:
2299	ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2300	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinction
2301	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2302	Fo = np.sqrt(ref[5])
2303	Fc = np.sqrt(ref[7])
2304	w = 2.0*Fo/ref[6]
2305	if w*Fo >= calcControls['minF/sig']:
2306	sumFo += Fo
2307	sumFo2 += ref[5]
2308	sumdF += abs(Fo-Fc)
2309	sumdF2 += abs(ref[5]-ref[7])
2310	nobs += 1
2311	df[i] = -w*(Fo-Fc)
2312	sumwYo += (wFo)*2
2313	Histogram['Residuals']['Nobs'] = nobs
2314	Histogram['Residuals']['sumwYo'] = sumwYo
2315	SumwYo += sumwYo
2316	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['Residuals']['sumwYo'])100.)
2317	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2318	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2319	Histogram['Residuals'][phfx+'Nref'] = nobs
2320	Nobs += nobs
2321	if dlg:
2322	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2323	M = np.concatenate((M,wtFactor*df))
2324	# end of HKLF processing
2325	Histograms['sumwYo'] = SumwYo
2326	Histograms['Nobs'] = Nobs
2327	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
2328	if dlg:
2329	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2330	if not GoOn:
2331	parmDict['saved values'] = values
2332	dlg.Destroy()
2333	raise Exception #Abort!!
2334	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2335	if len(pVals):
2336	pSum = np.sum(pWtpVals*2)
2337	for name in pWsum:
2338	if pWsum:
2339	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
2340	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2341	Nobs += len(pVals)
2342	M = np.concatenate((M,np.sqrt(pWt)*pVals))
2343	return M
2344

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

Download in other formats: