Context Navigation

source: trunk/GSASIIstrMath.py @ 1559

Last change on this file since 1559 was 1559, checked in by vondreele, 8 years ago
use deepcopy for various copy operations in G2ddataGUI use complementary colors to background for cell edges trap superlattice in cubics more fixes to cell indexing routines from possible use of superlattice add HStrainVals to G2spc use it (maybe) in G2strIO, Math make arguments for GetReflPos? & GetReflPosDeriv? the same (use A not G) add GetDij? to G2strMath topas xye file comments start with "'" sometimes
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 116.6 KB

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

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

Download in other formats: