Context Navigation

source: trunk/GSASIIstrMath.py @ 1465

Last change on this file since 1465 was 1465, checked in by vondreele, 9 years ago
added mod to ImageIntegrate? requested by Steven Weigand - return NST array if requested. Powder extinction now correct for TOF - not same as for GSAS; that code has an error!
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 114.9 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2014-08-18 15:06:22 +0000 (Mon, 18 Aug 2014) $
8	# $Author: vondreele $
9	# $Revision: 1465 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 1465 2014-08-18 15:06:22Z 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: 1465 $")
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[7]parmDict[hfx+'Lam']**2
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]
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,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	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*refl[14]/refl[16] #wave/abs corr
1152	dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[17] #/ext corr
1153	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1154
1155	def GetSampleSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1156	'Needs a doc string'
1157	if 'C' in calcControls[hfx+'histType']:
1158	costh = cosd(refl[5]/2.)
1159	#crystallite size
1160	if calcControls[phfx+'SizeType'] == 'isotropic':
1161	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1162	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1163	H = np.array(refl[:3])
1164	P = np.array(calcControls[phfx+'SizeAxis'])
1165	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1166	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1167	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1168	else: #ellipsoidal crystallites
1169	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1170	H = np.array(refl[:3])
1171	lenR = G2pwd.ellipseSize(H,Sij,GB)
1172	Sgam = 1.8wave/(np.picosth*lenR)
1173	#microstrain
1174	if calcControls[phfx+'MustrainType'] == 'isotropic':
1175	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1176	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1177	H = np.array(refl[:3])
1178	P = np.array(calcControls[phfx+'MustrainAxis'])
1179	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1180	Si = parmDict[phfx+'Mustrain;i']
1181	Sa = parmDict[phfx+'Mustrain;a']
1182	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1183	else: #generalized - P.W. Stephens model
1184	pwrs = calcControls[phfx+'MuPwrs']
1185	sum = 0
1186	for i,pwr in enumerate(pwrs):
1187	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1188	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1189	elif 'T' in calcControls[hfx+'histType']:
1190	#crystallite size
1191	if calcControls[phfx+'SizeType'] == 'isotropic':
1192	Sgam = 1.e-4parmDict[hfx+'difC']parmDict[phfx+'Size;i']
1193	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1194	H = np.array(refl[:3])
1195	P = np.array(calcControls[phfx+'SizeAxis'])
1196	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1197	Sgam = 1.e-4parmDict[hfx+'difC'](parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1198	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1199	else: #ellipsoidal crystallites
1200	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1201	H = np.array(refl[:3])
1202	lenR = G2pwd.ellipseSize(H,Sij,GB)
1203	Sgam = 1.e-4parmDict[hfx+'difC'](refl[4]*2lenR)
1204	#microstrain
1205	if calcControls[phfx+'MustrainType'] == 'isotropic':
1206	Mgam = 1.e-6parmDict[hfx+'difC']parmDict[phfx+'Mustrain;i']
1207	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1208	H = np.array(refl[:3])
1209	P = np.array(calcControls[phfx+'MustrainAxis'])
1210	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1211	Si = parmDict[phfx+'Mustrain;i']
1212	Sa = parmDict[phfx+'Mustrain;a']
1213	Mgam = 1.e-6parmDict[hfx+'difC']SiSa/np.sqrt((SicosP)*2+(SasinP)**2)
1214	else: #generalized - P.W. Stephens model
1215	pwrs = calcControls[phfx+'MuPwrs']
1216	sum = 0
1217	for i,pwr in enumerate(pwrs):
1218	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1219	Mgam = 1.e-6parmDict[hfx+'difC']refl[4]*2sum
1220
1221	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1222	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1223	sig /= ateln2
1224	return sig,gam
1225
1226	def GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1227	'Needs a doc string'
1228	gamDict = {}
1229	sigDict = {}
1230	if 'C' in calcControls[hfx+'histType']:
1231	costh = cosd(refl[5]/2.)
1232	tanth = tand(refl[5]/2.)
1233	#crystallite size derivatives
1234	if calcControls[phfx+'SizeType'] == 'isotropic':
1235	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1236	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1237	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1238	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1239	H = np.array(refl[:3])
1240	P = np.array(calcControls[phfx+'SizeAxis'])
1241	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1242	Si = parmDict[phfx+'Size;i']
1243	Sa = parmDict[phfx+'Size;a']
1244	gami = (1.8wave/np.pi)/(SiSa)
1245	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1246	Sgam = gami*sqtrm
1247	gam = Sgam/costh
1248	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
1249	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
1250	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1251	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1252	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1253	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1254	else: #ellipsoidal crystallites
1255	const = 1.8wave/(np.picosth)
1256	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1257	H = np.array(refl[:3])
1258	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1259	Sgam = const/lenR
1260	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1261	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1262	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1263	gamDict[phfx+'Size;mx'] = Sgam
1264	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1265
1266	#microstrain derivatives
1267	if calcControls[phfx+'MustrainType'] == 'isotropic':
1268	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1269	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1270	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1271	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1272	H = np.array(refl[:3])
1273	P = np.array(calcControls[phfx+'MustrainAxis'])
1274	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1275	Si = parmDict[phfx+'Mustrain;i']
1276	Sa = parmDict[phfx+'Mustrain;a']
1277	gami = 0.018SiSa*tanth/np.pi
1278	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1279	Mgam = gami/sqtrm
1280	dsi = -gamiSicosP2/sqtrm3
1281	dsa = -gamiSasinP2/sqtrm3
1282	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1283	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1284	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1285	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1286	else: #generalized - P.W. Stephens model
1287	pwrs = calcControls[phfx+'MuPwrs']
1288	const = 0.018refl[4]2tanth
1289	sum = 0
1290	for i,pwr in enumerate(pwrs):
1291	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1292	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1293	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
1294	sigDict[phfx+'Mustrain:'+str(i)] = \
1295	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1296	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1297	for i in range(len(pwrs)):
1298	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1299	gamDict[phfx+'Mustrain;mx'] = Mgam
1300	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1301	else: #'T'OF
1302	if calcControls[phfx+'SizeType'] == 'isotropic':
1303	Sgam = 1.e-4parmDict[hfx+'difC']parmDict[phfx+'Size;i']
1304	gamDict[phfx+'Size;i'] = 1.e-4parmDict[hfx+'difC']parmDict[phfx+'Size;mx']
1305	sigDict[phfx+'Size;i'] = 1.e-4parmDict[hfx+'difC']Sgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1306	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1307	const = 1.e-4*parmDict[hfx+'difC']
1308	H = np.array(refl[:3])
1309	P = np.array(calcControls[phfx+'SizeAxis'])
1310	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1311	Si = parmDict[phfx+'Size;i']
1312	Sa = parmDict[phfx+'Size;a']
1313	gami = const(SiSa)
1314	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1315	Sgam = gami*sqtrm
1316	dsi = gamiSicosP**2/sqtrm-gam/Si
1317	dsa = gamiSasinP**2/sqtrm-gam/Sa
1318	gamDict[phfx+'Size;i'] = constparmDict[phfx+'Size;mx']Sa
1319	gamDict[phfx+'Size;a'] = constparmDict[phfx+'Size;mx']Si
1320	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1321	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1322	else: #ellipsoidal crystallites
1323	const = 1.e-4*parmDict[hfx+'difC']
1324	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1325	H = np.array(refl[:3])
1326	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1327	Sgam = const/lenR
1328	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1329	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1330	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1331	gamDict[phfx+'Size;mx'] = Sgam
1332	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1333
1334	#microstrain derivatives
1335	if calcControls[phfx+'MustrainType'] == 'isotropic':
1336	Mgam = 1.e-6parmDict[hfx+'difC']parmDict[phfx+'Mustrain;i']
1337	gamDict[phfx+'Mustrain;i'] = 1.e-6parmDict[hfx+'difC']parmDict[phfx+'Mustrain;mx']
1338	sigDict[phfx+'Mustrain;i'] = 2.e-6parmDict[hfx+'difC']Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1339	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1340	H = np.array(refl[:3])
1341	P = np.array(calcControls[phfx+'MustrainAxis'])
1342	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1343	Si = parmDict[phfx+'Mustrain;i']
1344	Sa = parmDict[phfx+'Mustrain;a']
1345	gami = 1.e-6parmDict[hfx+'difC']Si*Sa
1346	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1347	Mgam = gami/sqtrm
1348	dsi = -gamiSicosP2/sqtrm3
1349	dsa = -gamiSasinP2/sqtrm3
1350	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1351	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1352	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1353	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1354	else: #generalized - P.W. Stephens model
1355	pwrs = calcControls[phfx+'MuPwrs']
1356	const = 1.e-6parmDict[hfx+'difC']refl[4]**2
1357	sum = 0
1358	for i,pwr in enumerate(pwrs):
1359	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1360	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1361	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
1362	sigDict[phfx+'Mustrain:'+str(i)] = \
1363	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1364	Mgam = const*sum
1365	for i in range(len(pwrs)):
1366	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1367	gamDict[phfx+'Mustrain;mx'] = Mgam
1368	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1369
1370	return sigDict,gamDict
1371
1372	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1373	'Needs a doc string'
1374	h,k,l = refl[:3]
1375	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1376	d = np.sqrt(dsq)
1377
1378	refl[4] = d
1379	if 'C' in calcControls[hfx+'histType']:
1380	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1381	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1382	if 'Bragg' in calcControls[hfx+'instType']:
1383	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1384	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1385	else: #Debye-Scherrer - simple but maybe not right
1386	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1387	elif 'T' in calcControls[hfx+'histType']:
1388	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d*2+parmDict[hfx+'difB']d**3+parmDict[hfx+'Zero']
1389	#do I need sample position effects - maybe?
1390	return pos
1391
1392	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1393	'Needs a doc string'
1394	dpr = 180./np.pi
1395	h,k,l = refl[:3]
1396	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1397	dst = np.sqrt(dstsq)
1398	dsp = 1./dst
1399	if 'C' in calcControls[hfx+'histType']:
1400	pos = refl[5]-parmDict[hfx+'Zero']
1401	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1402	dpdw = const*dst
1403	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1404	dpdA = constwave/(2.0*dst)
1405	dpdZ = 1.0
1406	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1407	if 'Bragg' in calcControls[hfx+'instType']:
1408	dpdSh = -4.constcosd(pos/2.0)
1409	dpdTr = -constsind(pos)100.0
1410	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1411	else: #Debye-Scherrer - simple but maybe not right
1412	dpdXd = -const*cosd(pos)
1413	dpdYd = -const*sind(pos)
1414	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1415	elif 'T' in calcControls[hfx+'histType']:
1416	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1417	dpdZ = 1.0
1418	dpdDC = dsp
1419	dpdDA = dsp**2
1420	dpdDB = dsp**3
1421	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
1422
1423	def GetHStrainShift(refl,SGData,phfx,parmDict):
1424	'Needs a doc string'
1425	laue = SGData['SGLaue']
1426	uniq = SGData['SGUniq']
1427	h,k,l = refl[:3]
1428	if laue in ['m3','m3m']:
1429	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
1430	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
1431	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1432	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
1433	elif laue in ['3R','3mR']:
1434	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
1435	elif laue in ['4/m','4/mmm']:
1436	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
1437	elif laue in ['mmm']:
1438	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1439	elif laue in ['2/m']:
1440	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1441	if uniq == 'a':
1442	Dij += parmDict[phfx+'D23']kl
1443	elif uniq == 'b':
1444	Dij += parmDict[phfx+'D13']hl
1445	elif uniq == 'c':
1446	Dij += parmDict[phfx+'D12']hk
1447	else:
1448	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
1449	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
1450	return -Dijrefl[4]2tand(refl[5]/2.0)
1451
1452	def GetHStrainShiftDerv(refl,SGData,phfx):
1453	'Needs a doc string'
1454	laue = SGData['SGLaue']
1455	uniq = SGData['SGUniq']
1456	h,k,l = refl[:3]
1457	if laue in ['m3','m3m']:
1458	dDijDict = {phfx+'D11':h2+k2+l**2,
1459	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
1460	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1461	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
1462	elif laue in ['3R','3mR']:
1463	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
1464	elif laue in ['4/m','4/mmm']:
1465	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
1466	elif laue in ['mmm']:
1467	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1468	elif laue in ['2/m']:
1469	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1470	if uniq == 'a':
1471	dDijDict[phfx+'D23'] = k*l
1472	elif uniq == 'b':
1473	dDijDict[phfx+'D13'] = h*l
1474	elif uniq == 'c':
1475	dDijDict[phfx+'D12'] = h*k
1476	names.append()
1477	else:
1478	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
1479	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
1480	for item in dDijDict:
1481	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
1482	return dDijDict
1483
1484	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1485	'Needs a doc string'
1486	histoList = Histograms.keys()
1487	histoList.sort()
1488	for histogram in histoList:
1489	if 'PWDR' in histogram[:4]:
1490	Histogram = Histograms[histogram]
1491	hId = Histogram['hId']
1492	hfx = ':%d:'%(hId)
1493	Limits = calcControls[hfx+'Limits']
1494	if 'C' in calcControls[hfx+'histType']:
1495	shl = max(parmDict[hfx+'SH/L'],0.0005)
1496	Ka2 = False
1497	kRatio = 0.0
1498	if hfx+'Lam1' in parmDict.keys():
1499	Ka2 = True
1500	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1501	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1502	x,y,w,yc,yb,yd = Histogram['Data']
1503	xB = np.searchsorted(x,Limits[0])
1504	xF = np.searchsorted(x,Limits[1])
1505	ymb = np.array(y-yb)
1506	ymb = np.where(ymb,ymb,1.0)
1507	ycmb = np.array(yc-yb)
1508	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
1509	refLists = Histogram['Reflection Lists']
1510	for phase in refLists:
1511	Phase = Phases[phase]
1512	pId = Phase['pId']
1513	phfx = '%d:%d:'%(pId,hId)
1514	refDict = refLists[phase]
1515	sumFo = 0.0
1516	sumdF = 0.0
1517	sumFosq = 0.0
1518	sumdFsq = 0.0
1519	for refl in refDict['RefList']:
1520	if 'C' in calcControls[hfx+'histType']:
1521	yp = np.zeros_like(yb)
1522	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1523	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1524	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1525	iFin2 = iFin
1526	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
1527	if Ka2:
1528	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1529	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1530	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1531	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1532	if not iBeg2+iFin2: #peak below low limit - skip peak
1533	continue
1534	elif not iBeg2-iFin2: #peak above high limit - done
1535	break
1536	yp[iBeg2:iFin2] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2])) #and here
1537	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11](1.+kRatio)),0.0))
1538	elif 'T' in calcControls[hfx+'histType']:
1539	yp = np.zeros_like(yb)
1540	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1541	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1542	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1543	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
1544	refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
1545	Fo = np.sqrt(np.abs(refl[8]))
1546	Fc = np.sqrt(np.abs(refl[9]))
1547	sumFo += Fo
1548	sumFosq += refl[8]**2
1549	sumdF += np.abs(Fo-Fc)
1550	sumdFsq += (refl[8]-refl[9])**2
1551	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1552	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1553	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1554	Histogram['Residuals']['hId'] = hId
1555	elif 'HKLF' in histogram[:4]:
1556	Histogram = Histograms[histogram]
1557	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1558
1559	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1560	'Needs a doc string'
1561
1562	def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
1563	U = parmDict[hfx+'U']
1564	V = parmDict[hfx+'V']
1565	W = parmDict[hfx+'W']
1566	X = parmDict[hfx+'X']
1567	Y = parmDict[hfx+'Y']
1568	tanPos = tand(refl[5]/2.0)
1569	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)
1570	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
1571	sig = max(0.001,sig)
1572	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
1573	gam = max(0.001,gam)
1574	return sig,gam
1575
1576	def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
1577	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4]*2+ \
1578	parmDict[hfx+'sig-2']refl[4]4+parmDict[hfx+'sig-q']/refl[4]*2
1579	gam = parmDict[hfx+'X']refl[4]+parmDict[hfx+'Y']refl[4]**2
1580	Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict)
1581	sig += Ssig #save peak sigma
1582	sig = max(0.001,sig)
1583	gam += Sgam #save peak gamma
1584	gam = max(0.001,gam)
1585	return sig,gam
1586
1587	def GetReflAlpBet(refl,hfx,parmDict):
1588	alp = parmDict[hfx+'alpha']/refl[4]
1589	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]4+parmDict[hfx+'beta-q']/refl[4]2
1590	return alp,bet
1591
1592	hId = Histogram['hId']
1593	hfx = ':%d:'%(hId)
1594	bakType = calcControls[hfx+'bakType']
1595	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1596	yc = np.zeros_like(yb)
1597	cw = np.diff(x)
1598	cw = np.append(cw,cw[-1])
1599
1600	if 'C' in calcControls[hfx+'histType']:
1601	shl = max(parmDict[hfx+'SH/L'],0.002)
1602	Ka2 = False
1603	if hfx+'Lam1' in parmDict.keys():
1604	wave = parmDict[hfx+'Lam1']
1605	Ka2 = True
1606	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1607	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1608	else:
1609	wave = parmDict[hfx+'Lam']
1610	for phase in Histogram['Reflection Lists']:
1611	refDict = Histogram['Reflection Lists'][phase]
1612	Phase = Phases[phase]
1613	pId = Phase['pId']
1614	pfx = '%d::'%(pId)
1615	phfx = '%d:%d:'%(pId,hId)
1616	hfx = ':%d:'%(hId)
1617	SGData = Phase['General']['SGData']
1618	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1619	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #Do I want to modify by Dij?
1620	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1621	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1622	Vst = np.sqrt(nl.det(G)) #V*
1623	if not Phase['General'].get('doPawley'):
1624	time0 = time.time()
1625	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1626	# print 'sf calc time: %.3fs'%(time.time()-time0)
1627	time0 = time.time()
1628	for iref,refl in enumerate(refDict['RefList']):
1629	if 'C' in calcControls[hfx+'histType']:
1630	h,k,l = refl[:3]
1631	Uniq = np.inner(refl[:3],SGMT)
1632	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
1633	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
1634	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
1635	refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
1636	refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1637	refl[11] = VstLorenz
1638
1639	if Phase['General'].get('doPawley'):
1640	try:
1641	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1642	refl[9] = parmDict[pInd]
1643	except KeyError:
1644	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1645	continue
1646	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1647	iBeg = np.searchsorted(x,refl[5]-fmin)
1648	iFin = np.searchsorted(x,refl[5]+fmax)
1649	if not iBeg+iFin: #peak below low limit - skip peak
1650	continue
1651	elif not iBeg-iFin: #peak above high limit - done
1652	break
1653	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
1654	if Ka2:
1655	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1656	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1657	iBeg = np.searchsorted(x,pos2-fmin)
1658	iFin = np.searchsorted(x,pos2+fmax)
1659	if not iBeg+iFin: #peak below low limit - skip peak
1660	continue
1661	elif not iBeg-iFin: #peak above high limit - done
1662	return yc,yb
1663	yc[iBeg:iFin] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here
1664	elif 'T' in calcControls[hfx+'histType']:
1665	h,k,l = refl[:3]
1666	Uniq = np.inner(refl[:3],SGMT)
1667	refl[5] = GetReflPos(refl,0.0,G,hfx,calcControls,parmDict) #corrected reflection position
1668	Lorenz = sind(parmDict[hfx+'2-theta']/2)refl[4]*4 #TOF Lorentz correction
1669	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
1670	refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict) #peak sig & gam
1671	refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
1672	refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1673	refl[11] = VstLorenz
1674	if Phase['General'].get('doPawley'):
1675	try:
1676	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1677	refl[9] = parmDict[pInd]
1678	except KeyError:
1679	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1680	continue
1681	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1682	iBeg = np.searchsorted(x,refl[5]-fmin)
1683	iFin = np.searchsorted(x,refl[5]+fmax)
1684	if not iBeg+iFin: #peak below low limit - skip peak
1685	continue
1686	elif not iBeg-iFin: #peak above high limit - done
1687	break
1688	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]
1689	# print 'profile calc time: %.3fs'%(time.time()-time0)
1690	return yc,yb
1691
1692	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1693	'Needs a doc string'
1694
1695	def cellVaryDerv(pfx,SGData,dpdA):
1696	if SGData['SGLaue'] in ['-1',]:
1697	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1698	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1699	elif SGData['SGLaue'] in ['2/m',]:
1700	if SGData['SGUniq'] == 'a':
1701	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1702	elif SGData['SGUniq'] == 'b':
1703	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1704	else:
1705	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1706	elif SGData['SGLaue'] in ['mmm',]:
1707	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1708	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1709	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1710	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1711	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1712	elif SGData['SGLaue'] in ['3R', '3mR']:
1713	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
1714	elif SGData['SGLaue'] in ['m3m','m3']:
1715	return [[pfx+'A0',dpdA[0]]]
1716
1717	# create a list of dependent variables and set up a dictionary to hold their derivatives
1718	dependentVars = G2mv.GetDependentVars()
1719	depDerivDict = {}
1720	for j in dependentVars:
1721	depDerivDict[j] = np.zeros(shape=(len(x)))
1722	#print 'dependent vars',dependentVars
1723	lenX = len(x)
1724	hId = Histogram['hId']
1725	hfx = ':%d:'%(hId)
1726	bakType = calcControls[hfx+'bakType']
1727	dMdv = np.zeros(shape=(len(varylist),len(x)))
1728	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1729	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1730	bBpos =varylist.index(hfx+'Back:0')
1731	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1732	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1733	for name in varylist:
1734	if 'Debye' in name:
1735	id = int(name.split(':')[-1])
1736	parm = name[:int(name.rindex(':'))]
1737	ip = names.index(parm)
1738	dMdv[varylist.index(name)] = dMddb[3*id+ip]
1739	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1740	for name in varylist:
1741	if 'BkPk' in name:
1742	parm,id = name.split(';')
1743	id = int(id)
1744	if parm in names:
1745	ip = names.index(parm)
1746	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1747	cw = np.diff(x)
1748	cw = np.append(cw,cw[-1])
1749	Ka2 = False #also for TOF!
1750	if 'C' in calcControls[hfx+'histType']:
1751	shl = max(parmDict[hfx+'SH/L'],0.002)
1752	if hfx+'Lam1' in parmDict.keys():
1753	wave = parmDict[hfx+'Lam1']
1754	Ka2 = True
1755	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1756	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1757	else:
1758	wave = parmDict[hfx+'Lam']
1759	for phase in Histogram['Reflection Lists']:
1760	refDict = Histogram['Reflection Lists'][phase]
1761	Phase = Phases[phase]
1762	SGData = Phase['General']['SGData']
1763	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1764	pId = Phase['pId']
1765	pfx = '%d::'%(pId)
1766	phfx = '%d:%d:'%(pId,hId)
1767	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #And modify here by Dij?
1768	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1769	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1770	if not Phase['General'].get('doPawley'):
1771	time0 = time.time()
1772	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1773	# print 'sf-derv time %.3fs'%(time.time()-time0)
1774	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1775	time0 = time.time()
1776	for iref,refl in enumerate(refDict['RefList']):
1777	h,k,l = refl[:3]
1778	Uniq = np.inner(refl[:3],SGMT)
1779	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1780	if 'C' in calcControls[hfx+'histType']: #CW powder
1781	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1782	else: #'T'OF
1783	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
1784	iBeg = np.searchsorted(x,refl[5]-fmin)
1785	iFin = np.searchsorted(x,refl[5]+fmax)
1786	if not iBeg+iFin: #peak below low limit - skip peak
1787	continue
1788	elif not iBeg-iFin: #peak above high limit - done
1789	break
1790	pos = refl[5]
1791	if 'C' in calcControls[hfx+'histType']:
1792	tanth = tand(pos/2.0)
1793	costh = cosd(pos/2.0)
1794	lenBF = iFin-iBeg
1795	dMdpk = np.zeros(shape=(6,lenBF))
1796	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1797	for i in range(5):
1798	dMdpk[i] += 100.cw[iBeg:iFin]refl[11]refl[9]dMdipk[i]
1799	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1800	if Ka2:
1801	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
1802	iBeg2 = np.searchsorted(x,pos2-fmin)
1803	iFin2 = np.searchsorted(x,pos2+fmax)
1804	if iBeg2-iFin2:
1805	lenBF2 = iFin2-iBeg2
1806	dMdpk2 = np.zeros(shape=(6,lenBF2))
1807	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1808	for i in range(5):
1809	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11]refl[9]kRatio*dMdipk2[i]
1810	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11]*dMdipk2[0]
1811	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1812	else: #'T'OF
1813	lenBF = iFin-iBeg
1814	dMdpk = np.zeros(shape=(6,lenBF))
1815	dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
1816	for i in range(6):
1817	dMdpk[i] += refl[11]refl[9]dMdipk[i] #cw[iBeg:iFin]*
1818	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1819	if Phase['General'].get('doPawley'):
1820	dMdpw = np.zeros(len(x))
1821	try:
1822	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1823	idx = varylist.index(pIdx)
1824	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1825	if Ka2: #not for TOF either
1826	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1827	dMdv[idx] = dMdpw
1828	except: # ValueError:
1829	pass
1830	if 'C' in calcControls[hfx+'histType']:
1831	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1832	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1833	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1834	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1835	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1836	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1837	hfx+'DisplaceY':[dpdY,'pos'],}
1838	if 'Bragg' in calcControls[hfx+'instType']:
1839	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1840	hfx+'SurfRoughB':[dFdAb[1],'int'],})
1841	else:
1842	names.update({hfx+'Absorption':[dFdAb,'int'],})
1843	else: #'T'OF
1844	dpdA,dpdZ,dpdDC,dpdDA,dpdDB = GetReflPosDerv(refl,0.0,A,hfx,calcControls,parmDict)
1845	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
1846	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
1847	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
1848	hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
1849	hfx+'beta-q':[1./refl[4]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]2,'sig'],
1850	hfx+'sig-2':[refl[4]4,'sig'],hfx+'sig-q':[1./refl[4]2,'sig'],
1851	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
1852	for name in names:
1853	item = names[name]
1854	if name in varylist:
1855	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1856	if Ka2:
1857	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1858	elif name in dependentVars:
1859	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1860	if Ka2:
1861	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1862	for iPO in dIdPO:
1863	if iPO in varylist:
1864	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1865	if Ka2:
1866	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1867	elif iPO in dependentVars:
1868	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1869	if Ka2:
1870	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1871	for i,name in enumerate(['omega','chi','phi']):
1872	aname = pfx+'SH '+name
1873	if aname in varylist:
1874	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1875	if Ka2:
1876	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1877	elif aname in dependentVars:
1878	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1879	if Ka2:
1880	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1881	for iSH in dFdODF:
1882	if iSH in varylist:
1883	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1884	if Ka2:
1885	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1886	elif iSH in dependentVars:
1887	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1888	if Ka2:
1889	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1890	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1891	for name,dpdA in cellDervNames:
1892	if name in varylist:
1893	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1894	if Ka2:
1895	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1896	elif name in dependentVars:
1897	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1898	if Ka2:
1899	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1900	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1901	for name in dDijDict:
1902	if name in varylist:
1903	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1904	if Ka2:
1905	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1906	elif name in dependentVars:
1907	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1908	if Ka2:
1909	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1910	if 'C' in calcControls[hfx+'histType']:
1911	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)
1912	else: #'T'OF
1913	sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,hfx,phfx,calcControls,parmDict)
1914	for name in gamDict:
1915	if name in varylist:
1916	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1917	if Ka2:
1918	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1919	elif name in dependentVars:
1920	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1921	if Ka2:
1922	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1923	for name in sigDict:
1924	if name in varylist:
1925	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1926	if Ka2:
1927	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1928	elif name in dependentVars:
1929	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1930	if Ka2:
1931	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1932	for name in ['BabA','BabU']:
1933	if refl[9]:
1934	if phfx+name in varylist:
1935	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1936	if Ka2:
1937	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1938	elif phfx+name in dependentVars:
1939	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1940	if Ka2:
1941	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1942	if not Phase['General'].get('doPawley'):
1943	#do atom derivatives - for RB,F,X & U so far
1944	corr = dervDict['int']/refl[9]
1945	if Ka2:
1946	corr2 = dervDict2['int']/refl[9]
1947	for name in varylist+dependentVars:
1948	if '::RBV;' in name:
1949	pass
1950	else:
1951	try:
1952	aname = name.split(pfx)[1][:2]
1953	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1954	except IndexError:
1955	continue
1956	if name in varylist:
1957	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1958	if Ka2:
1959	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1960	elif name in dependentVars:
1961	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1962	if Ka2:
1963	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1964	# print 'profile derv time: %.3fs'%(time.time()-time0)
1965	# now process derivatives in constraints
1966	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1967	return dMdv
1968
1969	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1970	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
1971	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
1972	are returned to either dervRefine or HessRefine.
1973
1974	:returns:
1975	'''
1976	nobs = Histogram['Residuals']['Nobs']
1977	hId = Histogram['hId']
1978	hfx = ':%d:'%(hId)
1979	pfx = '%d::'%(Phase['pId'])
1980	phfx = '%d:%d:'%(Phase['pId'],hId)
1981	SGData = Phase['General']['SGData']
1982	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1983	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1984	refDict = Histogram['Data']
1985	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1986	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1987	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1988	dependentVars = G2mv.GetDependentVars()
1989	depDerivDict = {}
1990	for j in dependentVars:
1991	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1992	wdf = np.zeros(len(refDict['RefList']))
1993	if calcControls['F**2']:
1994	for iref,ref in enumerate(refDict['RefList']):
1995	if ref[6] > 0:
1996	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:]
1997	w = 1.0/ref[6]
1998	if w*ref[5] >= calcControls['minF/sig']:
1999	wdf[iref] = w*(ref[5]-ref[7])
2000	for j,var in enumerate(varylist):
2001	if var in dFdvDict:
2002	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']dervCorref[11]
2003	for var in dependentVars:
2004	if var in dFdvDict:
2005	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']dervCorref[11]
2006	if phfx+'Scale' in varylist:
2007	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
2008	elif phfx+'Scale' in dependentVars:
2009	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
2010	for item in ['Ep','Es','Eg']:
2011	if phfx+item in varylist and dervDict:
2012	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor
2013	elif phfx+item in dependentVars and dervDict:
2014	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
2015	for item in ['BabA','BabU']:
2016	if phfx+item in varylist:
2017	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
2018	elif phfx+item in dependentVars:
2019	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
2020	else:
2021	for iref,ref in enumerate(refDict['RefList']):
2022	if ref[5] > 0.:
2023	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[1:]
2024	Fo = np.sqrt(ref[5])
2025	Fc = np.sqrt(ref[7])
2026	w = 1.0/ref[6]
2027	if 2.0Fow*Fo >= calcControls['minF/sig']:
2028	wdf[iref] = 2.0Fow*(Fo-Fc)
2029	for j,var in enumerate(varylist):
2030	if var in dFdvDict:
2031	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']ref[11] #dervCor
2032	for var in dependentVars:
2033	if var in dFdvDict:
2034	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']ref[11] #dervCor
2035	if phfx+'Scale' in varylist:
2036	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]ref[11] #*dervCor
2037	elif phfx+'Scale' in dependentVars:
2038	depDerivDict[phfx+'Scale'][iref] = wref[9]ref[11] #*dervCor
2039	for item in ['Ep','Es','Eg']:
2040	if phfx+item in varylist and dervDict:
2041	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor #correct
2042	elif phfx+item in dependentVars and dervDict:
2043	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
2044	for item in ['BabA','BabU']:
2045	if phfx+item in varylist:
2046	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
2047	elif phfx+item in dependentVars:
2048	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
2049	return dMdvh,depDerivDict,wdf
2050
2051
2052	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2053	'''Loop over histograms and compute derivatives of the fitting
2054	model (M) with respect to all parameters. Results are returned in
2055	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2056	n is the number of parameters and m is the number of data
2057	points. This can exceed memory when m gets large. This routine is
2058	used when refinement derivatives are selected as "analtytic
2059	Jacobian" in Controls.
2060
2061	:returns: Jacobian numpy.array dMdv for all histograms concatinated
2062	'''
2063	parmDict.update(zip(varylist,values))
2064	G2mv.Dict2Map(parmDict,varylist)
2065	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2066	nvar = len(varylist)
2067	dMdv = np.empty(0)
2068	histoList = Histograms.keys()
2069	histoList.sort()
2070	for histogram in histoList:
2071	if 'PWDR' in histogram[:4]:
2072	Histogram = Histograms[histogram]
2073	hId = Histogram['hId']
2074	hfx = ':%d:'%(hId)
2075	wtFactor = calcControls[hfx+'wtFactor']
2076	Limits = calcControls[hfx+'Limits']
2077	x,y,w,yc,yb,yd = Histogram['Data']
2078	W = wtFactor*w
2079	xB = np.searchsorted(x,Limits[0])
2080	xF = np.searchsorted(x,Limits[1])
2081	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2082	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2083	elif 'HKLF' in histogram[:4]:
2084	Histogram = Histograms[histogram]
2085	phase = Histogram['Reflection Lists']
2086	Phase = Phases[phase]
2087	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2088	hfx = ':%d:'%(Histogram['hId'])
2089	wtFactor = calcControls[hfx+'wtFactor']
2090	# now process derivatives in constraints
2091	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2092	else:
2093	continue #skip non-histogram entries
2094	if len(dMdv):
2095	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2096	else:
2097	dMdv = np.sqrt(wtFactor)*dMdvh
2098
2099	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2100	if np.any(pVals):
2101	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2102	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2103
2104	return dMdv
2105
2106	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2107	'''Loop over histograms and compute derivatives of the fitting
2108	model (M) with respect to all parameters. For each histogram, the
2109	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2110	number of parameters and m is the number of data points *in the
2111	histogram*. The (n by n) Hessian is computed from each Jacobian
2112	and it is returned. This routine is used when refinement
2113	derivatives are selected as "analtytic Hessian" in Controls.
2114
2115	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2116	'''
2117	parmDict.update(zip(varylist,values))
2118	G2mv.Dict2Map(parmDict,varylist)
2119	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2120	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
2121	nvar = len(varylist)
2122	Hess = np.empty(0)
2123	histoList = Histograms.keys()
2124	histoList.sort()
2125	for histogram in histoList:
2126	if 'PWDR' in histogram[:4]:
2127	Histogram = Histograms[histogram]
2128	hId = Histogram['hId']
2129	hfx = ':%d:'%(hId)
2130	wtFactor = calcControls[hfx+'wtFactor']
2131	Limits = calcControls[hfx+'Limits']
2132	x,y,w,yc,yb,yd = Histogram['Data']
2133	W = wtFactor*w
2134	dy = y-yc
2135	xB = np.searchsorted(x,Limits[0])
2136	xF = np.searchsorted(x,Limits[1])
2137	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2138	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2139	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2140	Dy = dy[xB:xF][np.newaxis,:]
2141	dMdvh *= Wt
2142	if dlg:
2143	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2144	if len(Hess):
2145	Hess += np.inner(dMdvh,dMdvh)
2146	dMdvh = WtDy
2147	Vec += np.sum(dMdvh,axis=1)
2148	else:
2149	Hess = np.inner(dMdvh,dMdvh)
2150	dMdvh = WtDy
2151	Vec = np.sum(dMdvh,axis=1)
2152	elif 'HKLF' in histogram[:4]:
2153	Histogram = Histograms[histogram]
2154	phase = Histogram['Reflection Lists']
2155	Phase = Phases[phase]
2156	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2157	hId = Histogram['hId']
2158	hfx = ':%d:'%(Histogram['hId'])
2159	wtFactor = calcControls[hfx+'wtFactor']
2160	# now process derivatives in constraints
2161	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2162	# print 'matrix build time: %.3f'%(time.time()-time0)
2163
2164	if dlg:
2165	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2166	if len(Hess):
2167	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
2168	Hess += wtFactor*np.inner(dMdvh,dMdvh)
2169	else:
2170	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
2171	Hess = wtFactor*np.inner(dMdvh,dMdvh)
2172	else:
2173	continue #skip non-histogram entries
2174	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2175	if np.any(pVals):
2176	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
2177	Vec += np.sum(dpdvpWtpVals,axis=1)
2178	Hess += np.inner(dpdv*pWt,dpdv)
2179	return Vec,Hess
2180
2181	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2182	'Needs a doc string'
2183	Values2Dict(parmDict, varylist, values)
2184	G2mv.Dict2Map(parmDict,varylist)
2185	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2186	M = np.empty(0)
2187	SumwYo = 0
2188	Nobs = 0
2189	ApplyRBModels(parmDict,Phases,rigidbodyDict)
2190	histoList = Histograms.keys()
2191	histoList.sort()
2192	for histogram in histoList:
2193	if 'PWDR' in histogram[:4]:
2194	Histogram = Histograms[histogram]
2195	hId = Histogram['hId']
2196	hfx = ':%d:'%(hId)
2197	wtFactor = calcControls[hfx+'wtFactor']
2198	Limits = calcControls[hfx+'Limits']
2199	x,y,w,yc,yb,yd = Histogram['Data']
2200	yc *= 0.0 #zero full calcd profiles
2201	yb *= 0.0
2202	yd *= 0.0
2203	xB = np.searchsorted(x,Limits[0])
2204	xF = np.searchsorted(x,Limits[1])
2205	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2206	varylist,Histogram,Phases,calcControls,pawleyLookup)
2207	yc[xB:xF] += yb[xB:xF]
2208	if not np.any(y): #fill dummy data
2209	rv = st.poisson(yc[xB:xF])
2210	y[xB:xF] = rv.rvs()
2211	Z = np.ones_like(yc[xB:xF])
2212	Z[1::2] *= -1
2213	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2214	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2215	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2216	W = wtFactor*w
2217	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2218	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2219	Nobs += Histogram['Residuals']['Nobs']
2220	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2221	SumwYo += Histogram['Residuals']['sumwYo']
2222	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2223	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2224	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2225	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2226	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.))
2227	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.))
2228	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2229	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2230	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2231	if dlg:
2232	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2233	M = np.concatenate((M,wdy))
2234	#end of PWDR processing
2235	elif 'HKLF' in histogram[:4]:
2236	Histogram = Histograms[histogram]
2237	Histogram['Residuals'] = {}
2238	phase = Histogram['Reflection Lists']
2239	Phase = Phases[phase]
2240	hId = Histogram['hId']
2241	hfx = ':%d:'%(hId)
2242	wtFactor = calcControls[hfx+'wtFactor']
2243	pfx = '%d::'%(Phase['pId'])
2244	phfx = '%d:%d:'%(Phase['pId'],hId)
2245	SGData = Phase['General']['SGData']
2246	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2247	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2248	refDict = Histogram['Data']
2249	time0 = time.time()
2250	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2251	# print 'sf-calc time: %.3f'%(time.time()-time0)
2252	df = np.zeros(len(refDict['RefList']))
2253	sumwYo = 0
2254	sumFo = 0
2255	sumFo2 = 0
2256	sumdF = 0
2257	sumdF2 = 0
2258	nobs = 0
2259	if calcControls['F**2']:
2260	for i,ref in enumerate(refDict['RefList']):
2261	if ref[6] > 0:
2262	ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2263	w = 1.0/ref[6]
2264	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinction
2265	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2266	if w*ref[5] >= calcControls['minF/sig']:
2267	Fo = np.sqrt(ref[5])
2268	sumFo += Fo
2269	sumFo2 += ref[5]
2270	sumdF += abs(Fo-np.sqrt(ref[7]))
2271	sumdF2 += abs(ref[5]-ref[7])
2272	nobs += 1
2273	df[i] = -w*(ref[5]-ref[7])
2274	sumwYo += (wref[5])*2
2275	else:
2276	for i,ref in enumerate(refDict['RefList']):
2277	if ref[5] > 0.:
2278	ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2279	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinction
2280	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2281	Fo = np.sqrt(ref[5])
2282	Fc = np.sqrt(ref[7])
2283	w = 2.0*Fo/ref[6]
2284	if w*Fo >= calcControls['minF/sig']:
2285	sumFo += Fo
2286	sumFo2 += ref[5]
2287	sumdF += abs(Fo-Fc)
2288	sumdF2 += abs(ref[5]-ref[7])
2289	nobs += 1
2290	df[i] = -w*(Fo-Fc)
2291	sumwYo += (wFo)*2
2292	Histogram['Residuals']['Nobs'] = nobs
2293	Histogram['Residuals']['sumwYo'] = sumwYo
2294	SumwYo += sumwYo
2295	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['Residuals']['sumwYo'])100.)
2296	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2297	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2298	Histogram['Residuals'][phfx+'Nref'] = nobs
2299	Nobs += nobs
2300	if dlg:
2301	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2302	M = np.concatenate((M,wtFactor*df))
2303	# end of HKLF processing
2304	Histograms['sumwYo'] = SumwYo
2305	Histograms['Nobs'] = Nobs
2306	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
2307	if dlg:
2308	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2309	if not GoOn:
2310	parmDict['saved values'] = values
2311	dlg.Destroy()
2312	raise Exception #Abort!!
2313	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2314	if len(pVals):
2315	pSum = np.sum(pWtpVals*2)
2316	for name in pWsum:
2317	if pWsum:
2318	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
2319	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2320	Nobs += len(pVals)
2321	M = np.concatenate((M,np.sqrt(pWt)*pVals))
2322	return M
2323

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

Download in other formats: