Context Navigation

source: trunk/GSASIIstrMath.py @ 1484

Last change on this file since 1484 was 1484, checked in by vondreele, 9 years ago
change lower limit of LGmix to 0 from 0.1 - recompile Win-32 & Win-64 fortran sources fix ellipseSizeDeriv - didn't do it correctly before fix mustrain & size fxns & derivatives for CW (minor) & TOF major; all checked against numerical derivatives for both lorentzian & gaussian contributions as well as mixing coeff.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 115.9 KB

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

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

Download in other formats: