Context Navigation

source: trunk/GSASIIstrMath.py @ 1474

Last change on this file since 1474 was 1474, checked in by vondreele, 9 years ago
1st MC/SA tutorial various MC/SA fixes fix to background peak fitting for CW & TOF
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 115.4 KB

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

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

Download in other formats: