Context Navigation

source: trunk/GSASIIstrMath.py @ 2052

Last change on this file since 2052 was 2052, checked in by vondreele, 7 years ago
relax cell order rules for C- ortho & C-mono cell indexing change default view of modulation plot work on powder incommensurate fxn & derivs.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 161.6 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2015-11-14 22:13:43 +0000 (Sat, 14 Nov 2015) $
8	# $Author: vondreele $
9	# $Revision: 2052 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 2052 2015-11-14 22:13:43Z vondreele $
12	########### SVN repository information ###################
13	import time
14	import copy
15	import numpy as np
16	import numpy.ma as ma
17	import numpy.linalg as nl
18	import scipy.optimize as so
19	import scipy.stats as st
20	import GSASIIpath
21	GSASIIpath.SetVersionNumber("$Revision: 2052 $")
22	import GSASIIElem as G2el
23	import GSASIIlattice as G2lat
24	import GSASIIspc as G2spc
25	import GSASIIpwd as G2pwd
26	import GSASIImapvars as G2mv
27	import GSASIImath as G2mth
28	import GSASIIobj as G2obj
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*np.log(2.0)
38	twopi = 2.0*np.pi
39	twopisq = 2.0np.pi*2
40
41	################################################################################
42	##### Rigid Body Models
43	################################################################################
44
45	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
46	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
47	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
48	'''
49	atxIds = ['Ax:','Ay:','Az:']
50	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
51	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
52	if not RBIds['Vector'] and not RBIds['Residue']:
53	return
54	VRBIds = RBIds['Vector']
55	RRBIds = RBIds['Residue']
56	if Update:
57	RBData = rigidbodyDict
58	else:
59	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
60	if RBIds['Vector']: # first update the vector magnitudes
61	VRBData = RBData['Vector']
62	for i,rbId in enumerate(VRBIds):
63	if VRBData[rbId]['useCount']:
64	for j in range(len(VRBData[rbId]['VectMag'])):
65	name = '::RBV;'+str(j)+':'+str(i)
66	VRBData[rbId]['VectMag'][j] = parmDict[name]
67	for phase in Phases:
68	Phase = Phases[phase]
69	General = Phase['General']
70	cx,ct,cs,cia = General['AtomPtrs']
71	cell = General['Cell'][1:7]
72	Amat,Bmat = G2lat.cell2AB(cell)
73	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
74	pfx = str(Phase['pId'])+'::'
75	if Update:
76	RBModels = Phase['RBModels']
77	else:
78	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
79	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
80	jrb = VRBIds.index(RBObj['RBId'])
81	rbsx = str(irb)+':'+str(jrb)
82	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
83	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
84	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
85	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
86	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
87	TLS = RBObj['ThermalMotion']
88	if 'T' in TLS[0]:
89	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
90	TLS[1][i] = parmDict[pfx+pt+rbsx]
91	if 'L' in TLS[0]:
92	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
93	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
94	if 'S' in TLS[0]:
95	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
96	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
97	if 'U' in TLS[0]:
98	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
99	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
100	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
101	for i,x in enumerate(XYZ):
102	atId = RBObj['Ids'][i]
103	for j in [0,1,2]:
104	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
105	if UIJ[i][0] == 'A':
106	for j in range(6):
107	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
108	elif UIJ[i][0] == 'I':
109	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
110
111	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
112	jrb = RRBIds.index(RBObj['RBId'])
113	rbsx = str(irb)+':'+str(jrb)
114	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
115	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
116	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
117	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
118	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
119	TLS = RBObj['ThermalMotion']
120	if 'T' in TLS[0]:
121	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
122	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
123	if 'L' in TLS[0]:
124	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
125	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
126	if 'S' in TLS[0]:
127	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
128	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
129	if 'U' in TLS[0]:
130	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
131	for itors,tors in enumerate(RBObj['Torsions']):
132	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
133	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
134	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
135	for i,x in enumerate(XYZ):
136	atId = RBObj['Ids'][i]
137	for j in [0,1,2]:
138	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
139	if UIJ[i][0] == 'A':
140	for j in range(6):
141	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
142	elif UIJ[i][0] == 'I':
143	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
144
145	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
146	'Needs a doc string'
147	atxIds = ['dAx:','dAy:','dAz:']
148	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
149	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
150	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
151	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
152	PIds = ['Px:','Py:','Pz:']
153	OIds = ['Oa:','Oi:','Oj:','Ok:']
154	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
155	if not RBIds['Vector'] and not RBIds['Residue']:
156	return
157	VRBIds = RBIds['Vector']
158	RRBIds = RBIds['Residue']
159	RBData = rigidbodyDict
160	for item in parmDict:
161	if 'RB' in item:
162	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
163	General = Phase['General']
164	cx,ct,cs,cia = General['AtomPtrs']
165	cell = General['Cell'][1:7]
166	Amat,Bmat = G2lat.cell2AB(cell)
167	rpd = np.pi/180.
168	rpd2 = rpd**2
169	g = nl.inv(np.inner(Bmat,Bmat))
170	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
171	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
172	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
173	pfx = str(Phase['pId'])+'::'
174	RBModels = Phase['RBModels']
175	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
176	VModel = RBData['Vector'][RBObj['RBId']]
177	Q = RBObj['Orient'][0]
178	Pos = RBObj['Orig'][0]
179	jrb = VRBIds.index(RBObj['RBId'])
180	rbsx = str(irb)+':'+str(jrb)
181	dXdv = []
182	for iv in range(len(VModel['VectMag'])):
183	dCdv = []
184	for vec in VModel['rbVect'][iv]:
185	dCdv.append(G2mth.prodQVQ(Q,vec))
186	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
187	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
188	for ia,atId in enumerate(RBObj['Ids']):
189	atNum = AtLookup[atId]
190	dx = 0.00001
191	for iv in range(len(VModel['VectMag'])):
192	for ix in [0,1,2]:
193	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
194	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
195	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
196	for iv in range(4):
197	Q[iv] -= dx
198	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
199	Q[iv] += 2.*dx
200	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201	Q[iv] -= dx
202	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
203	for ix in [0,1,2]:
204	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
205	X = G2mth.prodQVQ(Q,Cart[ia])
206	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
207	dFdu = G2lat.U6toUij(dFdu.T)
208	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
209	dFdu = G2lat.UijtoU6(dFdu)
210	atNum = AtLookup[atId]
211	if 'T' in RBObj['ThermalMotion'][0]:
212	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
213	dFdvDict[pfx+name+rbsx] += dFdu[i]
214	if 'L' in RBObj['ThermalMotion'][0]:
215	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
216	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
217	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
218	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
219	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
220	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
221	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
222	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
223	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
224	if 'S' in RBObj['ThermalMotion'][0]:
225	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
226	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
227	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
228	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
229	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
230	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
231	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
232	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
233	if 'U' in RBObj['ThermalMotion'][0]:
234	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
235
236
237	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
238	Q = RBObj['Orient'][0]
239	Pos = RBObj['Orig'][0]
240	jrb = RRBIds.index(RBObj['RBId'])
241	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
242	rbsx = str(irb)+':'+str(jrb)
243	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
244	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
245	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
246	orId,pvId = torData[itors][:2]
247	pivotVec = Cart[orId]-Cart[pvId]
248	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
249	QB = G2mth.AVdeg2Q(0.001,pivotVec)
250	for ir in torData[itors][3]:
251	atNum = AtLookup[RBObj['Ids'][ir]]
252	rVec = Cart[ir]-Cart[pvId]
253	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
254	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
255	for ix in [0,1,2]:
256	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
257	for ia,atId in enumerate(RBObj['Ids']):
258	atNum = AtLookup[atId]
259	dx = 0.00001
260	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
261	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
262	for iv in range(4):
263	Q[iv] -= dx
264	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
265	Q[iv] += 2.*dx
266	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
267	Q[iv] -= dx
268	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
269	for ix in [0,1,2]:
270	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
271	X = G2mth.prodQVQ(Q,Cart[ia])
272	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
273	dFdu = G2lat.U6toUij(dFdu.T)
274	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
275	dFdu = G2lat.UijtoU6(dFdu)
276	atNum = AtLookup[atId]
277	if 'T' in RBObj['ThermalMotion'][0]:
278	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
279	dFdvDict[pfx+name+rbsx] += dFdu[i]
280	if 'L' in RBObj['ThermalMotion'][0]:
281	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
282	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
283	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
284	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
285	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
286	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
287	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
288	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
289	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
290	if 'S' in RBObj['ThermalMotion'][0]:
291	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
292	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
293	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
294	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
295	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
296	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
297	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
298	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
299	if 'U' in RBObj['ThermalMotion'][0]:
300	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
301
302	################################################################################
303	##### Penalty & restraint functions
304	################################################################################
305
306	def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
307	'Needs a doc string'
308	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
309	pNames = []
310	pVals = []
311	pWt = []
312	negWt = {}
313	pWsum = {}
314	for phase in Phases:
315	pId = Phases[phase]['pId']
316	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
317	General = Phases[phase]['General']
318	cx,ct,cs,cia = General['AtomPtrs']
319	textureData = General['SH Texture']
320	SGData = General['SGData']
321	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
322	cell = General['Cell'][1:7]
323	Amat,Bmat = G2lat.cell2AB(cell)
324	if phase not in restraintDict:
325	continue
326	phaseRest = restraintDict[phase]
327	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
328	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
329	['ChemComp','Sites'],['Texture','HKLs'],]
330	for name,rest in names:
331	pWsum[name] = 0.
332	itemRest = phaseRest[name]
333	if itemRest[rest] and itemRest['Use']:
334	wt = itemRest['wtFactor']
335	if name in ['Bond','Angle','Plane','Chiral']:
336	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
337	pNames.append(str(pId)+':'+name+':'+str(i))
338	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
339	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
340	if name == 'Bond':
341	calc = G2mth.getRestDist(XYZ,Amat)
342	elif name == 'Angle':
343	calc = G2mth.getRestAngle(XYZ,Amat)
344	elif name == 'Plane':
345	calc = G2mth.getRestPlane(XYZ,Amat)
346	elif name == 'Chiral':
347	calc = G2mth.getRestChiral(XYZ,Amat)
348	pVals.append(obs-calc)
349	pWt.append(wt/esd**2)
350	pWsum[name] += wt((obs-calc)/esd)*2
351	elif name in ['Torsion','Rama']:
352	coeffDict = itemRest['Coeff']
353	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
354	pNames.append(str(pId)+':'+name+':'+str(i))
355	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
356	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
357	if name == 'Torsion':
358	tor = G2mth.getRestTorsion(XYZ,Amat)
359	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
360	else:
361	phi,psi = G2mth.getRestRama(XYZ,Amat)
362	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
363	pVals.append(restr)
364	pWt.append(wt/esd**2)
365	pWsum[name] += wt(restr/esd)*2
366	elif name == 'ChemComp':
367	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
368	pNames.append(str(pId)+':'+name+':'+str(i))
369	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
370	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
371	calc = np.sum(mulfracfactors)
372	pVals.append(obs-calc)
373	pWt.append(wt/esd**2)
374	pWsum[name] += wt((obs-calc)/esd)*2
375	elif name == 'Texture':
376	SHkeys = textureData['SH Coeff'][1].keys()
377	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
378	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
379	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
380	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
381	PH = np.array(hkl)
382	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
383	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
384	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
385	Z1 = -ma.masked_greater(Z,0.0)
386	IndZ1 = np.array(ma.nonzero(Z1))
387	for ind in IndZ1.T:
388	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
389	pVals.append(Z1[ind[0]][ind[1]])
390	pWt.append(wt/esd1**2)
391	pWsum[name] += wt(-Z1[ind[0]][ind[1]]/esd1)*2
392	if ifesd2:
393	Z2 = 1.-Z
394	for ind in np.ndindex(grid,grid):
395	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
396	pVals.append(Z1[ind[0]][ind[1]])
397	pWt.append(wt/esd2**2)
398	pWsum[name] += wt(Z2/esd2)*2
399
400	for phase in Phases:
401	name = 'SH-Pref.Ori.'
402	pId = Phases[phase]['pId']
403	General = Phases[phase]['General']
404	SGData = General['SGData']
405	cell = General['Cell'][1:7]
406	pWsum[name] = 0.0
407	for hist in Phases[phase]['Histograms']:
408	if hist in Histograms and 'PWDR' in hist:
409	hId = Histograms[hist]['hId']
410	phfx = '%d:%d:'%(pId,hId)
411	if calcControls[phfx+'poType'] == 'SH':
412	toler = calcControls[phfx+'SHtoler']
413	wt = 1./toler**2
414	HKLs = np.array(calcControls[phfx+'SHhkl'])
415	SHnames = calcControls[phfx+'SHnames']
416	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
417	for i,PH in enumerate(HKLs):
418	phi,beta = G2lat.CrsAng(PH,cell,SGData)
419	SH3Coef = {}
420	for item in SHcof:
421	L,N = eval(item.strip('C'))
422	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
423	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
424	X = np.linspace(0,90.0,26)
425	Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
426	IndY = ma.nonzero(Y)
427	for ind in IndY[0]:
428	pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
429	pVals.append(Y[ind])
430	pWt.append(wt)
431	pWsum[name] += wt(Y[ind])*2
432	pWsum['PWLref'] = 0.
433	for item in varyList:
434	if 'PWLref' in item and parmDict[item] < 0.:
435	pId = int(item.split(':')[0])
436	if negWt[pId]:
437	pNames.append(item)
438	pVals.append(-parmDict[item])
439	pWt.append(negWt[pId])
440	pWsum['PWLref'] += negWt[pId](-parmDict[item])*2
441	pVals = np.array(pVals)
442	pWt = np.array(pWt) #should this be np.sqrt?
443	return pNames,pVals,pWt,pWsum
444
445	def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
446	'Needs a doc string'
447	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
448	pDerv = np.zeros((len(varyList),len(pVal)))
449	for phase in Phases:
450	# if phase not in restraintDict:
451	# continue
452	pId = Phases[phase]['pId']
453	General = Phases[phase]['General']
454	cx,ct,cs,cia = General['AtomPtrs']
455	SGData = General['SGData']
456	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
457	cell = General['Cell'][1:7]
458	Amat,Bmat = G2lat.cell2AB(cell)
459	textureData = General['SH Texture']
460
461	SHkeys = textureData['SH Coeff'][1].keys()
462	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
463	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
464	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
465	sam = SamSym[textureData['Model']]
466	phaseRest = restraintDict.get(phase,{})
467	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
468	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
469	'ChemComp':'Sites','Texture':'HKLs'}
470	lasthkl = np.array([0,0,0])
471	for ip,pName in enumerate(pNames):
472	pnames = pName.split(':')
473	if pId == int(pnames[0]):
474	name = pnames[1]
475	if 'PWL' in pName:
476	pDerv[varyList.index(pName)][ip] += 1.
477	continue
478	elif 'SH-' in pName:
479	continue
480	id = int(pnames[2])
481	itemRest = phaseRest[name]
482	if name in ['Bond','Angle','Plane','Chiral']:
483	indx,ops,obs,esd = itemRest[names[name]][id]
484	dNames = []
485	for ind in indx:
486	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
487	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
488	if name == 'Bond':
489	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
490	elif name == 'Angle':
491	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
492	elif name == 'Plane':
493	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
494	elif name == 'Chiral':
495	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
496	elif name in ['Torsion','Rama']:
497	coffDict = itemRest['Coeff']
498	indx,ops,cofName,esd = itemRest[names[name]][id]
499	dNames = []
500	for ind in indx:
501	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
502	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
503	if name == 'Torsion':
504	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
505	else:
506	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
507	elif name == 'ChemComp':
508	indx,factors,obs,esd = itemRest[names[name]][id]
509	dNames = []
510	for ind in indx:
511	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
512	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
513	deriv = mul*factors
514	elif 'Texture' in name:
515	deriv = []
516	dNames = []
517	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
518	hkl = np.array(hkl)
519	if np.any(lasthkl-hkl):
520	PH = np.array(hkl)
521	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
522	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
523	lasthkl = copy.copy(hkl)
524	if 'unit' in name:
525	pass
526	else:
527	gam = float(pnames[3])
528	psi = float(pnames[4])
529	for SHname in ODFln:
530	l,m,n = eval(SHname[1:])
531	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
532	dNames += [str(pId)+'::'+SHname]
533	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
534	for dName,drv in zip(dNames,deriv):
535	try:
536	ind = varyList.index(dName)
537	pDerv[ind][ip] += drv
538	except ValueError:
539	pass
540
541	lasthkl = np.array([0,0,0])
542	for ip,pName in enumerate(pNames):
543	deriv = []
544	dNames = []
545	pnames = pName.split(':')
546	if 'SH-' in pName and pId == int(pnames[0]):
547	hId = int(pnames[1])
548	phfx = '%d:%d:'%(pId,hId)
549	psi = float(pnames[4])
550	HKLs = calcControls[phfx+'SHhkl']
551	SHnames = calcControls[phfx+'SHnames']
552	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
553	hkl = np.array(HKLs[int(pnames[3])])
554	if np.any(lasthkl-hkl):
555	PH = np.array(hkl)
556	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
557	SH3Coef = {}
558	for item in SHcof:
559	L,N = eval(item.strip('C'))
560	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
561	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
562	lasthkl = copy.copy(hkl)
563	for SHname in SHnames:
564	l,n = eval(SHname[1:])
565	SH3name = 'C%d,0,%d'%(l,n)
566	Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
567	dNames += [phfx+SHname]
568	deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
569	for dName,drv in zip(dNames,deriv):
570	try:
571	ind = varyList.index(dName)
572	pDerv[ind][ip] += drv
573	except ValueError:
574	pass
575	return pDerv
576
577	################################################################################
578	##### Function & derivative calculations
579	################################################################################
580
581	def GetAtomFXU(pfx,calcControls,parmDict):
582	'Needs a doc string'
583	Natoms = calcControls['Natoms'][pfx]
584	Tdata = Natoms*[' ',]
585	Mdata = np.zeros(Natoms)
586	IAdata = Natoms*[' ',]
587	Fdata = np.zeros(Natoms)
588	FFdata = []
589	BLdata = []
590	Xdata = np.zeros((3,Natoms))
591	dXdata = np.zeros((3,Natoms))
592	Uisodata = np.zeros(Natoms)
593	Uijdata = np.zeros((6,Natoms))
594	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
595	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
596	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
597	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
598	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
599	for iatm in range(Natoms):
600	for key in keys:
601	parm = pfx+key+str(iatm)
602	if parm in parmDict:
603	keys[key][iatm] = parmDict[parm]
604	Fdata = np.where(Fdata,Fdata,1.e-8) #avoid divide by zero in derivative calc.
605	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
606
607	def GetAtomSSFXU(pfx,calcControls,parmDict):
608	'Needs a doc string'
609	Natoms = calcControls['Natoms'][pfx]
610	maxSSwave = calcControls['maxSSwave'][pfx]
611	Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
612	'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
613	XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
614	FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
615	USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
616	MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
617	waveTypes = []
618	keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
619	'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
620	'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
621	'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
622	'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
623	'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
624	for iatm in range(Natoms):
625	waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
626	for key in keys:
627	for m in range(Nwave[key[0]]):
628	parm = pfx+key+str(iatm)+':%d'%(m)
629	if parm in parmDict:
630	keys[key][m][iatm] = parmDict[parm]
631	return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
632
633	def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
634
635	Natoms = calcControls['Natoms'][pfx]
636	maxSSwave = calcControls['maxSSwave'][pfx]
637	Smult = np.zeros((Natoms,len(SGOps)))
638	TauT = np.zeros((Natoms,len(SGOps)))
639	for ix,xyz in enumerate(XData.T):
640	for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
641	sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
642	Smult[ix][isym] = sdet
643	TauT[ix][isym] = tauT
644	return Smult,TauT
645
646	def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
647	''' Compute structure factors for all h,k,l for phase
648	puts the result, F^2, in each ref[8] in refList
649	input:
650
651	:param dict refDict: where
652	'RefList' list where each ref = h,k,l,m,d,...
653	'FF' dict of form factors - filed in below
654	:param np.array G: reciprocal metric tensor
655	:param str pfx: phase id string
656	:param dict SGData: space group info. dictionary output from SpcGroup
657	:param dict calcControls:
658	:param dict ParmDict:
659
660	'''
661	phfx = pfx.split(':')[0]+hfx
662	ast = np.sqrt(np.diag(G))
663	Mast = twopisq*np.multiply.outer(ast,ast)
664	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
665	SGT = np.array([ops[1] for ops in SGData['SGOps']])
666	FFtables = calcControls['FFtables']
667	BLtables = calcControls['BLtables']
668	Flack = 1.0
669	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
670	Flack = 1.-2.*parmDict[phfx+'Flack']
671	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
672	TwDict = refDict.get('TwDict',{})
673	if 'S' in calcControls[hfx+'histType']:
674	NTL = calcControls[phfx+'NTL']
675	NM = calcControls[phfx+'TwinNMN']+1
676	TwinLaw = calcControls[phfx+'TwinLaw']
677	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
678	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
679	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
680	FF = np.zeros(len(Tdata))
681	if 'NC' in calcControls[hfx+'histType']:
682	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
683	elif 'X' in calcControls[hfx+'histType']:
684	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
685	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
686	Uij = np.array(G2lat.U6toUij(Uijdata))
687	bij = Mast*Uij.T
688	blkSize = 100 #no. of reflections in a block
689	nRef = refDict['RefList'].shape[0]
690	if not len(refDict['FF']): #no form factors - 1st time thru StructureFactor
691	if 'N' in calcControls[hfx+'histType']:
692	dat = G2el.getBLvalues(BLtables)
693	refDict['FF']['El'] = dat.keys()
694	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
695	else: #'X'
696	dat = G2el.getFFvalues(FFtables,0.)
697	refDict['FF']['El'] = dat.keys()
698	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
699	for iref,ref in enumerate(refDict['RefList']):
700	SQ = 1./(2.ref[4])*2
701	dat = G2el.getFFvalues(FFtables,SQ)
702	refDict['FF']['FF'][iref] *= dat.values()
703	#reflection processing begins here - big arrays!
704	iBeg = 0
705	while iBeg < nRef:
706	iFin = min(iBeg+blkSize,nRef)
707	refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems)
708	H = refl.T[:3] #array(blkSize,3)
709	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,3) or (blkSize,3)
710	TwMask = np.any(H,axis=-1)
711	if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
712	for ir in range(blkSize):
713	iref = ir+iBeg
714	if iref in TwDict:
715	for i in TwDict[iref]:
716	for n in range(NTL):
717	H[ir][i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
718	TwMask = np.any(H,axis=-1)
719	SQ = 1./(2.refl.T[4])*2 #array(blkSize)
720	SQfactor = 4.0SQtwopisq #ditto prev.
721	if 'T' in calcControls[hfx+'histType']:
722	if 'P' in calcControls[hfx+'histType']:
723	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
724	else:
725	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
726	FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
727	FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
728	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),len(SGT)*len(TwinLaw))
729	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
730	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
731	Uniq = np.inner(H,SGMT)
732	Phi = np.inner(H,SGT)
733	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
734	sinp = np.sin(phase)
735	cosp = np.cos(phase)
736	biso = -SQfactor*Uisodata[:,np.newaxis]
737	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
738	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
739	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
740	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/len(SGMT)
741	if 'T' in calcControls[hfx+'histType']:
742	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-np.reshape(FlackFPP,sinp.shape)sinp*Tcorr])
743	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,np.reshape(FlackFPP,cosp.shape)cosp*Tcorr])
744	else:
745	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-FlackFPPsinp*Tcorr])
746	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,FlackFPPcosp*Tcorr])
747	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl
748	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
749	if SGData['SGInv']: #centrosymmetric; B=0
750	fbs[0] *= 0.
751	if 'P' in calcControls[hfx+'histType']:
752	refl.T[9] = np.sum(fas2,axis=0)+np.sum(fbs2,axis=0)
753	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
754	else:
755	if len(TwinLaw) > 1:
756	refl.T[9] = np.sum(fas[:,:,0],axis=0)2+np.sum(fbs[:,:,0],axis=0)2 #FcT from primary twin element
757	refl.T[7] = np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fas,axis=0)**2,axis=-1)+ \
758	np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fbs,axis=0)**2,axis=-1) #Fc sum over twins
759	refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f"
760	else:
761	refl.T[9] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2
762	refl.T[7] = np.copy(refl.T[9])
763	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
764	# refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
765	iBeg += blkSize
766
767	def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
768	'Needs a doc string'
769	phfx = pfx.split(':')[0]+hfx
770	ast = np.sqrt(np.diag(G))
771	Mast = twopisq*np.multiply.outer(ast,ast)
772	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
773	SGT = np.array([ops[1] for ops in SGData['SGOps']])
774	FFtables = calcControls['FFtables']
775	BLtables = calcControls['BLtables']
776	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
777	TwDict = refDict.get('TwDict',{})
778	if 'S' in calcControls[hfx+'histType']:
779	NTL = calcControls[phfx+'NTL']
780	NM = calcControls[phfx+'TwinNMN']+1
781	TwinLaw = calcControls[phfx+'TwinLaw']
782	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
783	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
784	nTwin = len(TwinLaw)
785	nRef = len(refDict['RefList'])
786	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
787	mSize = len(Mdata)
788	FF = np.zeros(len(Tdata))
789	if 'NC' in calcControls[hfx+'histType']:
790	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
791	elif 'X' in calcControls[hfx+'histType']:
792	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
793	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
794	Uij = np.array(G2lat.U6toUij(Uijdata))
795	bij = Mast*Uij.T
796	dFdvDict = {}
797	if nTwin > 1:
798	dFdfr = np.zeros((nRef,nTwin,mSize))
799	dFdx = np.zeros((nRef,nTwin,mSize,3))
800	dFdui = np.zeros((nRef,nTwin,mSize))
801	dFdua = np.zeros((nRef,nTwin,mSize,6))
802	dFdbab = np.zeros((nRef,nTwin,2))
803	dFdfl = np.zeros((nRef,nTwin))
804	dFdtw = np.zeros((nRef,nTwin))
805	else:
806	dFdfr = np.zeros((nRef,mSize))
807	dFdx = np.zeros((nRef,mSize,3))
808	dFdui = np.zeros((nRef,mSize))
809	dFdua = np.zeros((nRef,mSize,6))
810	dFdbab = np.zeros((nRef,2))
811	dFdfl = np.zeros((nRef))
812	dFdtw = np.zeros((nRef))
813	Flack = 1.0
814	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
815	Flack = 1.-2.*parmDict[phfx+'Flack']
816	time0 = time.time()
817	nref = len(refDict['RefList'])/100
818	for iref,refl in enumerate(refDict['RefList']):
819	if 'T' in calcControls[hfx+'histType']:
820	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
821	H = np.array(refl[:3])
822	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)
823	TwMask = np.any(H,axis=-1)
824	if TwinLaw.shape[0] > 1 and TwDict:
825	if iref in TwDict:
826	for i in TwDict[iref]:
827	for n in range(NTL):
828	H[i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
829	TwMask = np.any(H,axis=-1)
830	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
831	SQfactor = 8.0SQnp.pi**2
832	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
833	Bab = parmDict[phfx+'BabA']*dBabdA
834	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
835	FF = refDict['FF']['FF'][iref].T[Tindx].T
836	Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3)
837	Phi = np.inner(H,SGT)
838	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
839	sinp = np.sin(phase)
840	cosp = np.cos(phase)
841	occ = Mdata*Fdata/len(SGT)
842	biso = -SQfactor*Uisodata[:,np.newaxis]
843	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
844	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
845	Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
846	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
847	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
848	Tcorr = (np.reshape(Tiso,Tuij.shape)Tuij).Tocc
849	fot = (FF+FP-Bab)*Tcorr
850	fotp = FPP*Tcorr
851	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr])
852	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr])
853	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins)
854	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
855	fax = np.array([-fotsinp,-fotpcosp]) #positions array(2,ntwi,nEqv,nAtoms)
856	fbx = np.array([fotcosp,-fotpsinp])
857	#sum below is over Uniq
858	dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
859	dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
860	dfadui = np.sum(-SQfactor*fa,axis=-2)
861	if nTwin > 1:
862	dfadx = np.array([np.sum(twopiUniq[it]np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
863	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
864	# array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
865	else:
866	dfadx = np.sum(twopiUniqnp.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
867	dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
868	# array(2,nAtom,3) & array(2,nAtom,6)
869	if not SGData['SGInv']:
870	dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem
871	dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
872	dfbdui = np.sum(-SQfactor*fb,axis=-2)
873	if len(TwinLaw) > 1:
874	dfbdx = np.array([np.sum(twopiUniq[it]np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
875	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
876	else:
877	dfadfl = np.sum(-FPPTcorrsinp)
878	dfbdfl = np.sum(FPPTcorrcosp)
879	dfbdx = np.sum(twopiUniqnp.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)
880	dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
881	else:
882	dfbdfr = np.zeros_like(dfadfr)
883	dfbdx = np.zeros_like(dfadx)
884	dfbdui = np.zeros_like(dfadui)
885	dfbdua = np.zeros_like(dfadua)
886	dfbdba = np.zeros_like(dfadba)
887	dfadfl = 0.0
888	dfbdfl = 0.0
889	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
890	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
891	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
892	2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
893	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
894	2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
895	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
896	2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
897	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
898	2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
899	else:
900	SA = fas[0]+fas[1]
901	SB = fbs[0]+fbs[1]
902	if nTwin > 1:
903	dFdfr[iref] = [2.TwMask[it](SA[it]dfadfr[0][it]+SA[it]dfadfr[1][it]+SB[it]dfbdfr[0][it]+SB[it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
904	dFdx[iref] = [2.TwMask[it](SA[it]dfadx[it][0]+SA[it]dfadx[it][1]+SB[it]dfbdx[it][0]+SB[it]dfbdx[it][1]) for it in range(nTwin)]
905	dFdui[iref] = [2.TwMask[it](SA[it]dfadui[it][0]+SA[it]dfadui[it][1]+SB[it]dfbdui[it][0]+SB[it]dfbdui[it][1]) for it in range(nTwin)]
906	dFdua[iref] = [2.TwMask[it](SA[it]dfadua[it][0]+SA[it]dfadua[it][1]+SB[it]dfbdua[it][0]+SB[it]dfbdua[it][1]) for it in range(nTwin)]
907	dFdtw[iref] = np.sum(TwMaskfas,axis=0)2+np.sum(TwMaskfbs,axis=0)**2
908
909	else: #these are good for no twin single crystals
910	dFdfr[iref] = (2.SA(dfadfr[0]+dfadfr[1])+2.SB(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
911	dFdx[iref] = 2.SA(dfadx[0]+dfadx[1])+2.SB(dfbdx[0]+dfbdx[1])
912	dFdui[iref] = 2.SA(dfadui[0]+dfadui[1])+2.SB(dfbdui[0]+dfbdui[1])
913	dFdua[iref] = 2.SA(dfadua[0]+dfadua[1])+2.SB(dfbdua[0]+dfbdua[1])
914	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
915	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
916	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
917	# GSASIIpath.IPyBreak()
918	if not iref%100 :
919	print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
920	#loop over atoms - each dict entry is list of derivatives for all the reflections
921	if nTwin > 1:
922	for i in range(len(Mdata)): #these all OK?
923	dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
924	dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
925	dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
926	dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
927	dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
928	dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
929	dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
930	dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
931	dFdvDict[pfx+'AU12:'+str(i)] = 2.np.sum(dFdua.T[3][i]TwinFr[:,np.newaxis],axis=0)
932	dFdvDict[pfx+'AU13:'+str(i)] = 2.np.sum(dFdua.T[4][i]TwinFr[:,np.newaxis],axis=0)
933	dFdvDict[pfx+'AU23:'+str(i)] = 2.np.sum(dFdua.T[5][i]TwinFr[:,np.newaxis],axis=0)
934	else:
935	for i in range(len(Mdata)):
936	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
937	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
938	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
939	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
940	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
941	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
942	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
943	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
944	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
945	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
946	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
947	dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
948	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
949	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
950	if nTwin > 1:
951	for i in range(nTwin):
952	dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
953	return dFdvDict
954
955	def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
956	'''
957	Compute super structure factors for all h,k,l,m for phase
958	puts the result, F^2, in each ref[8+im] in refList
959	input:
960
961	:param dict refDict: where
962	'RefList' list where each ref = h,k,l,m,d,...
963	'FF' dict of form factors - filed in below
964	:param np.array G: reciprocal metric tensor
965	:param str pfx: phase id string
966	:param dict SGData: space group info. dictionary output from SpcGroup
967	:param dict calcControls:
968	:param dict ParmDict:
969
970	'''
971	nxs = np.newaxis
972	phfx = pfx.split(':')[0]+hfx
973	ast = np.sqrt(np.diag(G))
974	Mast = twopisq*np.multiply.outer(ast,ast)
975	SGInv = SGData['SGInv']
976	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
977	SGT = np.array([ops[1] for ops in SGData['SGOps']])
978	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
979	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
980	# eps = SSGMT[:,3,3]
981	FFtables = calcControls['FFtables']
982	BLtables = calcControls['BLtables']
983	Flack = 1.0
984	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
985	Flack = 1.-2.*parmDict[phfx+'Flack']
986	TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],]) #4D?
987	TwDict = refDict.get('TwDict',{})
988	if 'S' in calcControls[hfx+'histType']:
989	NTL = calcControls[phfx+'NTL']
990	NM = calcControls[phfx+'TwinNMN']+1
991	TwinLaw = calcControls[phfx+'TwinLaw'] #this'll have to be 4D also...
992	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
993	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
994	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
995	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
996	# Smult,TauT = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
997	# GSASIIpath.IPyBreak()
998	# if SGInv:
999	# TauT = np.hstack((TauT,-TauT))
1000	# eps = np.concatenate((eps,-eps))
1001	modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1002	FF = np.zeros(len(Tdata))
1003	if 'NC' in calcControls[hfx+'histType']:
1004	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1005	else:
1006	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1007	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1008	Uij = np.array(G2lat.U6toUij(Uijdata)).T
1009	bij = Mast*Uij
1010	blkSize = 100 #no. of reflections in a block
1011	nRef = refDict['RefList'].shape[0]
1012	if not len(refDict['FF']):
1013	if 'N' in calcControls[hfx+'histType']:
1014	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
1015	refDict['FF']['El'] = dat.keys()
1016	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1017	else:
1018	dat = G2el.getFFvalues(FFtables,0.)
1019	refDict['FF']['El'] = dat.keys()
1020	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
1021	for iref,ref in enumerate(refDict['RefList']):
1022	SQ = 1./(2.ref[5])*2
1023	dat = G2el.getFFvalues(FFtables,SQ)
1024	refDict['FF']['FF'][iref] *= dat.values()
1025	time0 = time.time()
1026	#reflection processing begins here - big arrays!
1027	iBeg = 0
1028	while iBeg < nRef:
1029	iFin = min(iBeg+blkSize,nRef)
1030	refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems)
1031	H = refl.T[:4] #array(blkSize,4)
1032	HP = H[:3]+modQ[:,nxs]*H[3:] #projected hklm to hkl
1033	# H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4)
1034	# TwMask = np.any(H,axis=-1)
1035	# if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1036	# for ir in range(blkSize):
1037	# iref = ir+iBeg
1038	# if iref in TwDict:
1039	# for i in TwDict[iref]:
1040	# for n in range(NTL):
1041	# H[ir][i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
1042	# TwMask = np.any(H,axis=-1)
1043	SQ = 1./(2.refl.T[5])*2 #array(blkSize)
1044	SQfactor = 4.0SQtwopisq #ditto prev.
1045	Uniq = np.inner(H.T,SSGMT)
1046	UniqP = np.inner(HP.T,SGMT)
1047	Phi = np.inner(H.T,SSGT)
1048	# PhiP = np.inner(HP.T,SGT)
1049	if SGInv: #if centro - expand HKL sets
1050	Uniq = np.hstack((Uniq,-Uniq))
1051	Phi = np.hstack((Phi,-Phi))
1052	UniqP = np.hstack((UniqP,-UniqP))
1053	# PhiP = np.hstack((PhiP,-PhiP))
1054	if 'T' in calcControls[hfx+'histType']:
1055	if 'P' in calcControls[hfx+'histType']:
1056	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1057	else:
1058	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1059	FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1060	FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1061	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),Uniq.shape[1]*len(TwinLaw))
1062	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1063	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1064	phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1065	sinp = np.sin(phase)
1066	cosp = np.cos(phase)
1067	biso = -SQfactor*Uisodata[:,nxs]
1068	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1069	HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1) #use hklt proj to hkl
1070	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1071	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/Uniq.shape[1] #refBlk x ops x atoms
1072	if 'T' in calcControls[hfx+'histType']:
1073	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-np.reshape(FlackFPP,sinp.shape)sinp*Tcorr])
1074	fb = np.array([np.reshape(FlackFPP,cosp.shape)cospTcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinp*Tcorr])
1075	else:
1076	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-FlackFPPsinp*Tcorr])
1077	fb = np.array([FlackFPPcospTcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinp*Tcorr])
1078	GfpuA = G2mth.Modulation(waveTypes,Uniq,UniqP,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
1079	fag = faGfpuA[0]-fbGfpuA[1] #real; 2 x refBlk x sym x atoms
1080	fbg = fbGfpuA[0]+faGfpuA[1]
1081	fas = np.sum(np.sum(fag,axis=-1),axis=-1) #2 x refBlk; sum sym & atoms
1082	fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1083	# GSASIIpath.IPyBreak()
1084	if 'P' in calcControls[hfx+'histType']:
1085	refl.T[10] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2 #square of sums
1086	# refl.T[10] = np.sum(fas2,axis=0)+np.sum(fbs2,axis=0)
1087	refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f"
1088	else:
1089	if len(TwinLaw) > 1:
1090	refl.T[10] = np.sum(fas[:,:,0],axis=0)2+np.sum(fbs[:,:,0],axis=0)2 #FcT from primary twin element
1091	refl.T[8] = np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fas,axis=0)**2,axis=-1)+ \
1092	np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fbs,axis=0)**2,axis=-1) #Fc sum over twins
1093	refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f"
1094	else:
1095	refl.T[10] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2 #square of sums
1096	refl.T[8] = np.copy(refl.T[10])
1097	refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f"
1098	iBeg += blkSize
1099	print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1100
1101	def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1102	'Needs a doc string'
1103	nxs = np.newaxis
1104	phfx = pfx.split(':')[0]+hfx
1105	ast = np.sqrt(np.diag(G))
1106	Mast = twopisq*np.multiply.outer(ast,ast)
1107	SGInv = SGData['SGInv']
1108	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1109	SGT = np.array([ops[1] for ops in SGData['SGOps']])
1110	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1111	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1112	FFtables = calcControls['FFtables']
1113	BLtables = calcControls['BLtables']
1114	TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
1115	TwDict = refDict.get('TwDict',{})
1116	if 'S' in calcControls[hfx+'histType']:
1117	NTL = calcControls[phfx+'NTL']
1118	NM = calcControls[phfx+'TwinNMN']+1
1119	TwinLaw = calcControls[phfx+'TwinLaw']
1120	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1121	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1122	nTwin = len(TwinLaw)
1123	nRef = len(refDict['RefList'])
1124	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1125	mSize = len(Mdata) #no. atoms
1126	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1127	# Smult,TauT = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
1128	# if SGInv:
1129	# TauT = np.hstack((TauT,-TauT))
1130	modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1131	FF = np.zeros(len(Tdata))
1132	if 'NC' in calcControls[hfx+'histType']:
1133	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1134	elif 'X' in calcControls[hfx+'histType']:
1135	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1136	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1137	Uij = np.array(G2lat.U6toUij(Uijdata)).T
1138	bij = Mast*Uij
1139	if not len(refDict['FF']):
1140	if 'N' in calcControls[hfx+'histType']:
1141	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
1142	else:
1143	dat = G2el.getFFvalues(FFtables,0.)
1144	refDict['FF']['El'] = dat.keys()
1145	refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1146	dFdvDict = {}
1147	if nTwin > 1:
1148	dFdfr = np.zeros((nRef,nTwin,mSize))
1149	dFdx = np.zeros((nRef,nTwin,mSize,3))
1150	dFdui = np.zeros((nRef,nTwin,mSize))
1151	dFdua = np.zeros((nRef,nTwin,mSize,6))
1152	dFdbab = np.zeros((nRef,nTwin,2))
1153	dFdfl = np.zeros((nRef,nTwin))
1154	dFdtw = np.zeros((nRef,nTwin))
1155	dFdGfA = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
1156	dFdGfB = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
1157	dFdGxA = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
1158	dFdGxB = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
1159	dFdGuA = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
1160	dFdGuB = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
1161	else:
1162	dFdfr = np.zeros((nRef,mSize))
1163	dFdx = np.zeros((nRef,mSize,3))
1164	dFdui = np.zeros((nRef,mSize))
1165	dFdua = np.zeros((nRef,mSize,6))
1166	dFdbab = np.zeros((nRef,2))
1167	dFdfl = np.zeros((nRef))
1168	dFdtw = np.zeros((nRef))
1169	dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1170	dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1171	dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1172	Flack = 1.0
1173	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1174	Flack = 1.-2.*parmDict[phfx+'Flack']
1175	time0 = time.time()
1176	nRef = len(refDict['RefList'])/100
1177	for iref,refl in enumerate(refDict['RefList']):
1178	if 'T' in calcControls[hfx+'histType']:
1179	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1180	H = np.array(refl[:4])
1181	HP = H[:3]+modQ*H[3:] #projected hklm to hkl
1182	# H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(4,nTwins) or (4)
1183	# TwMask = np.any(H,axis=-1)
1184	# if TwinLaw.shape[0] > 1 and TwDict:
1185	# if iref in TwDict:
1186	# for i in TwDict[iref]:
1187	# for n in range(NTL):
1188	# H[i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
1189	# TwMask = np.any(H,axis=-1)
1190	SQ = 1./(2.refl[4+im])2 # or (sin(theta)/lambda)*2
1191	SQfactor = 8.0SQnp.pi**2
1192	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1193	Bab = parmDict[phfx+'BabA']*dBabdA
1194	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1195	FF = refDict['FF']['FF'][iref].T[Tindx]
1196	Uniq = np.inner(H,SSGMT)
1197	Phi = np.inner(H,SSGT)
1198	UniqP = np.inner(HP,SGMT)
1199	if SGInv: #if centro - expand HKL sets
1200	Uniq = np.vstack((Uniq,-Uniq))
1201	Phi = np.hstack((Phi,-Phi))
1202	UniqP = np.vstack((UniqP,-UniqP))
1203	phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1204	sinp = np.sin(phase)
1205	cosp = np.cos(phase)
1206	occ = Mdata*Fdata/Uniq.shape[0]
1207	biso = -SQfactor*Uisodata[:,nxs]
1208	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T #ops x atoms
1209	HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms
1210	Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1211	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
1212	Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms
1213	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/Uniq.shape[0] #ops x atoms
1214	fot = (FF+FP-Bab)*Tcorr #ops x atoms
1215	fotp = FPP*Tcorr #ops x atoms
1216	GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,UniqP,Hij,FSSdata,XSSdata,USSdata,Mast)
1217	# GfpuA is 2 x ops x atoms
1218	# derivs are: ops x atoms x waves x 2,6,or 12 parms as [real,imag] parts
1219	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
1220	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr]) #or array(2,nEqv,nAtoms)
1221	fag = faGfpuA[0]-fbGfpuA[1]
1222	fbg = fbGfpuA[0]+faGfpuA[1]
1223
1224	fas = np.sum(np.sum(fag,axis=1),axis=1) # 2 x twin
1225	fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1226	fax = np.array([-fotsinp,-fotpcosp]) #positions; 2 x twin x ops x atoms
1227	fbx = np.array([fotcosp,-fotpsinp])
1228	fax = faxGfpuA[0]-fbxGfpuA[1]
1229	fbx = fbxGfpuA[0]+faxGfpuA[1]
1230	#sum below is over Uniq
1231	dfadfr = np.sum(fag/occ,axis=1) #Fdata != 0 ever avoids /0. problem
1232	dfbdfr = np.sum(fbg/occ,axis=1) #Fdata != 0 avoids /0. problem
1233	dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1234	dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1235	dfadui = np.sum(-SQfactor*fag,axis=1)
1236	dfbdui = np.sum(-SQfactor*fbg,axis=1)
1237	if nTwin > 1:
1238	dfadx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1239	dfbdx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1240	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1241	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1242	# array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
1243	dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1244	dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1245	dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1246	dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1247	dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1248	dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1249	else:
1250	dfadx = np.sum(twopiUniq[:,:3]np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 3xyz; sum nOps
1251	dfbdx = np.sum(twopiUniq[:,:3]np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)
1252	dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 6Uij; sum nOps
1253	dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2) #these are correct also for twins above
1254	# array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1255	dfadGf = np.sum(fa[:,:,:,nxs,nxs]dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdf[1][nxs,:,:,:,:],axis=1)
1256	dfbdGf = np.sum(fb[:,:,:,nxs,nxs]dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdf[1][nxs,:,:,:,:],axis=1)
1257	dfadGx = np.sum(fa[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1258	dfbdGx = np.sum(fb[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1259	dfadGu = np.sum(fa[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1260	dfbdGu = np.sum(fb[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1261	if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative
1262	dfadfl = np.sum(-FPPTcorrsinp)
1263	dfbdfl = np.sum(FPPTcorrcosp)
1264	else:
1265	dfadfl = 1.0
1266	dfbdfl = 1.0
1267	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
1268	SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin])
1269	SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin])
1270	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1271	dFdfr[iref] = 2.(SAdfadfr[0]+SAdfadfr[1]+SBdfbdfr[0]+SBdfbdfr[1])Mdata/len(Uniq) #array(nRef,nAtom)
1272	dFdx[iref] = 2.(SAdfadx[0]+SAdfadx[1]+SBdfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3)
1273	dFdui[iref] = 2.(SAdfadui[0]+SAdfadui[1]+SBdfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom)
1274	dFdua[iref] = 2.(SAdfadua[0]+SAdfadua[1]+SBdfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6)
1275	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
1276
1277	dFdGf[iref] = 2.(SAdfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2)
1278	dFdGx[iref] = 2.(SAdfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6)
1279	dFdGu[iref] = 2.(SAdfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12)
1280	# dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
1281	# 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
1282	# dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
1283	# 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
1284	# dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
1285	# 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
1286	# dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
1287	# 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
1288	# dFdGf[iref] = 2.(fas[0]dfadGf[0]+fas[1]*dfadGf[1])+ \
1289	# 2.(fbs[0]dfbdGf[0]+fbs[1]*dfbdGf[1])
1290	# dFdGx[iref] = 2.(fas[0]dfadGx[0]+fas[1]*dfadGx[1])+ \
1291	# 2.(fbs[0]dfbdGx[0]+fbs[1]*dfbdGx[1])
1292	# dFdGu[iref] = 2.(fas[0]dfadGu[0]+fas[1]*dfadGu[1])+ \
1293	# 2.(fbs[0]dfbdGu[0]+fbs[1]*dfbdGu[1])
1294	else:
1295	if nTwin > 1:
1296	dFdfr[iref] = [2.TwMask[it](SA[it]dfadfr[0][it]+SA[it]dfadfr[1][it]+SB[it]dfbdfr[0][it]+SB[it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
1297	dFdx[iref] = [2.TwMask[it](SA[it]dfadx[it][0]+SA[it]dfadx[it][1]+SB[it]dfbdx[it][0]+SB[it]dfbdx[it][1]) for it in range(nTwin)]
1298	dFdui[iref] = [2.TwMask[it](SA[it]dfadui[it][0]+SA[it]dfadui[it][1]+SB[it]dfbdui[it][0]+SB[it]dfbdui[it][1]) for it in range(nTwin)]
1299	dFdua[iref] = [2.TwMask[it](SA[it]dfadua[it][0]+SA[it]dfadua[it][1]+SB[it]dfbdua[it][0]+SB[it]dfbdua[it][1]) for it in range(nTwin)]
1300	dFdtw[iref] = np.sum(TwMaskfas,axis=0)2+np.sum(TwMaskfbs,axis=0)**2
1301
1302	dFdGf[iref] = [2.TwMask[it](SA[it]dfadGf[1]+SB[it]dfbdGf[1]) for it in range(nTwin)]
1303	dFdGx[iref] = [2.TwMask[it](SA[it]dfadGx[1]+SB[it]dfbdGx[1]) for it in range(nTwin)]
1304	dFdGu[iref] = [2.TwMask[it](SA[it]dfadGu[1]+SB[it]dfbdGu[1]) for it in range(nTwin)]
1305	else: #these are good for no twin single crystals
1306	dFdfr[iref] = 2.(SAdfadfr[0]+SAdfadfr[1]+SBdfbdfr[0]+SBdfbdfr[1])Mdata/len(Uniq) #array(nRef,nAtom)
1307	dFdx[iref] = 2.(SAdfadx[0]+SAdfadx[1]+SBdfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3)
1308	dFdui[iref] = 2.(SAdfadui[0]+SAdfadui[1]+SBdfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom)
1309	dFdua[iref] = 2.(SAdfadua[0]+SAdfadua[1]+SBdfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6)
1310	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
1311
1312	dFdGf[iref] = 2.(SAdfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2)
1313	dFdGx[iref] = 2.(SAdfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6)
1314	dFdGu[iref] = 2.(SAdfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12)
1315	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
1316	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
1317	#loop over atoms - each dict entry is list of derivatives for all the reflections
1318	# GSASIIpath.IPyBreak()
1319	if not iref%100 :
1320	print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1321	for i in range(len(Mdata)): #loop over atoms
1322	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1323	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1324	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1325	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1326	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1327	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1328	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1329	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1330	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1331	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1332	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1333	for j in range(FSSdata.shape[1]): #loop over waves Fzero & Fwid?
1334	dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1335	dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1336	for j in range(XSSdata.shape[1]): #loop over waves 'Tzero','Xslope','Yslope','Zslope'?
1337	dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i]
1338	dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i]
1339	dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j)] = dFdGx.T[2][j][i]
1340	dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j)] = dFdGx.T[3][j][i]
1341	dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
1342	dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
1343	for j in range(USSdata.shape[1]): #loop over waves
1344	dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1345	dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1346	dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1347	dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1348	dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1349	dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1350	dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1351	dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1352	dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1353	dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1354	dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1355	dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1356
1357	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1358	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1359	return dFdvDict
1360
1361	def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1362	''' Single crystal extinction function; returns extinction & derivative
1363	'''
1364	extCor = 1.0
1365	dervDict = {}
1366	dervCor = 1.0
1367	if calcControls[phfx+'EType'] != 'None':
1368	SQ = 1/(4.ref[4+im]*2)
1369	if 'C' in parmDict[hfx+'Type']:
1370	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
1371	else: #'T'
1372	cos2T = 1.0-2.SQref[12+im]**2 #cos(2theta)
1373	if 'SXC' in parmDict[hfx+'Type']:
1374	AV = 7.9406e5/parmDict[pfx+'Vol']**2
1375	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1376	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
1377	PLZ = AVP12ref[9+im]parmDict[hfx+'Lam']*2
1378	elif 'SNT' in parmDict[hfx+'Type']:
1379	AV = 1.e7/parmDict[pfx+'Vol']**2
1380	PL = SQ
1381	PLZ = AVref[9+im]ref[12+im]**2
1382	elif 'SNC' in parmDict[hfx+'Type']:
1383	AV = 1.e7/parmDict[pfx+'Vol']**2
1384	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1385	PLZ = AVref[9+im]parmDict[hfx+'Lam']**2
1386
1387	if 'Primary' in calcControls[phfx+'EType']:
1388	PLZ *= 1.5
1389	else:
1390	if 'C' in parmDict[hfx+'Type']:
1391	PLZ *= calcControls[phfx+'Tbar']
1392	else: #'T'
1393	PLZ *= ref[13+im] #t-bar
1394	if 'Primary' in calcControls[phfx+'EType']:
1395	PLZ *= 1.5
1396	PSIG = parmDict[phfx+'Ep']
1397	elif 'I & II' in calcControls[phfx+'EType']:
1398	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
1399	elif 'Type II' in calcControls[phfx+'EType']:
1400	PSIG = parmDict[phfx+'Es']
1401	else: # 'Secondary Type I'
1402	PSIG = parmDict[phfx+'Eg']/PL
1403
1404	AG = 0.58+0.48cos2T+0.24cos2T**2
1405	AL = 0.025+0.285*cos2T
1406	BG = 0.02-0.025*cos2T
1407	BL = 0.15-0.2(0.75-cos2T)*2
1408	if cos2T < 0.:
1409	BL = -0.45*cos2T
1410	CG = 2.
1411	CL = 2.
1412	PF = PLZ*PSIG
1413
1414	if 'Gaussian' in calcControls[phfx+'EApprox']:
1415	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
1416	extCor = np.sqrt(PF4)
1417	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
1418	else:
1419	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
1420	extCor = np.sqrt(PF4)
1421	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
1422
1423	dervCor = (1.+PF)*PF3 #extinction corr for other derivatives
1424	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1425	dervDict[phfx+'Ep'] = -ref[7+im]PLZPF3
1426	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1427	dervDict[phfx+'Es'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
1428	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1429	dervDict[phfx+'Eg'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
1430
1431	return 1./extCor,dervDict,dervCor
1432
1433	def Dict2Values(parmdict, varylist):
1434	'''Use before call to leastsq to setup list of values for the parameters
1435	in parmdict, as selected by key in varylist'''
1436	return [parmdict[key] for key in varylist]
1437
1438	def Values2Dict(parmdict, varylist, values):
1439	''' Use after call to leastsq to update the parameter dictionary with
1440	values corresponding to keys in varylist'''
1441	parmdict.update(zip(varylist,values))
1442
1443	def GetNewCellParms(parmDict,varyList):
1444	'Needs a doc string'
1445	newCellDict = {}
1446	Anames = ['A'+str(i) for i in range(6)]
1447	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1448	for item in varyList:
1449	keys = item.split(':')
1450	if keys[2] in Ddict:
1451	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
1452	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
1453	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1454	return newCellDict # is e.g. {'0::D11':A0-D11}
1455
1456	def ApplyXYZshifts(parmDict,varyList):
1457	'''
1458	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1459
1460	:param dict parmDict: parameter dictionary
1461	:param list varyList: list of variables (not used!)
1462	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1463
1464	'''
1465	newAtomDict = {}
1466	for item in parmDict:
1467	if 'dA' in item:
1468	parm = ''.join(item.split('d'))
1469	parmDict[parm] += parmDict[item]
1470	newAtomDict[item] = [parm,parmDict[parm]]
1471	return newAtomDict
1472
1473	def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1474	'Spherical harmonics texture'
1475	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1476	if 'T' in calcControls[hfx+'histType']:
1477	tth = parmDict[hfx+'2-theta']
1478	else:
1479	tth = refl[5+im]
1480	odfCor = 1.0
1481	H = refl[:3]
1482	cell = G2lat.Gmat2cell(g)
1483	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1484	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1485	phi,beta = G2lat.CrsAng(H,cell,SGData)
1486	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1487	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1488	for item in SHnames:
1489	L,M,N = eval(item.strip('C'))
1490	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1491	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1492	Lnorm = G2lat.Lnorm(L)
1493	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1494	return odfCor
1495
1496	def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1497	'Spherical harmonics texture derivatives'
1498	if 'T' in calcControls[hfx+'histType']:
1499	tth = parmDict[hfx+'2-theta']
1500	else:
1501	tth = refl[5+im]
1502	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1503	odfCor = 1.0
1504	dFdODF = {}
1505	dFdSA = [0,0,0]
1506	H = refl[:3]
1507	cell = G2lat.Gmat2cell(g)
1508	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1509	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1510	phi,beta = G2lat.CrsAng(H,cell,SGData)
1511	psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1512	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1513	for item in SHnames:
1514	L,M,N = eval(item.strip('C'))
1515	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1516	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1517	Lnorm = G2lat.Lnorm(L)
1518	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1519	dFdODF[pfx+item] = LnormKclKsl
1520	for i in range(3):
1521	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
1522	return odfCor,dFdODF,dFdSA
1523
1524	def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1525	'spherical harmonics preferred orientation (cylindrical symmetry only)'
1526	if 'T' in calcControls[hfx+'histType']:
1527	tth = parmDict[hfx+'2-theta']
1528	else:
1529	tth = refl[5+im]
1530	odfCor = 1.0
1531	H = refl[:3]
1532	cell = G2lat.Gmat2cell(g)
1533	Sangls = [0.,0.,0.]
1534	if 'Bragg' in calcControls[hfx+'instType']:
1535	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1536	IFCoup = True
1537	else:
1538	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1539	IFCoup = False
1540	phi,beta = G2lat.CrsAng(H,cell,SGData)
1541	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1542	SHnames = calcControls[phfx+'SHnames']
1543	for item in SHnames:
1544	L,N = eval(item.strip('C'))
1545	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1546	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1547	Lnorm = G2lat.Lnorm(L)
1548	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1549	return np.squeeze(odfCor)
1550
1551	def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1552	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1553	if 'T' in calcControls[hfx+'histType']:
1554	tth = parmDict[hfx+'2-theta']
1555	else:
1556	tth = refl[5+im]
1557	odfCor = 1.0
1558	dFdODF = {}
1559	H = refl[:3]
1560	cell = G2lat.Gmat2cell(g)
1561	Sangls = [0.,0.,0.]
1562	if 'Bragg' in calcControls[hfx+'instType']:
1563	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1564	IFCoup = True
1565	else:
1566	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1567	IFCoup = False
1568	phi,beta = G2lat.CrsAng(H,cell,SGData)
1569	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1570	SHnames = calcControls[phfx+'SHnames']
1571	for item in SHnames:
1572	L,N = eval(item.strip('C'))
1573	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1574	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1575	Lnorm = G2lat.Lnorm(L)
1576	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1577	dFdODF[phfx+item] = KclKslLnorm
1578	return odfCor,dFdODF
1579
1580	def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1581	'March-Dollase preferred orientation correction'
1582	POcorr = 1.0
1583	MD = parmDict[phfx+'MD']
1584	if MD != 1.0:
1585	MDAxis = calcControls[phfx+'MDAxis']
1586	sumMD = 0
1587	for H in uniq:
1588	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1589	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1590	sumMD += A**3
1591	POcorr = sumMD/len(uniq)
1592	return POcorr
1593
1594	def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1595	'Needs a doc string'
1596	POcorr = 1.0
1597	POderv = {}
1598	if calcControls[phfx+'poType'] == 'MD':
1599	MD = parmDict[phfx+'MD']
1600	MDAxis = calcControls[phfx+'MDAxis']
1601	sumMD = 0
1602	sumdMD = 0
1603	for H in uniq:
1604	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1605	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1606	sumMD += A**3
1607	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1608	POcorr = sumMD/len(uniq)
1609	POderv[phfx+'MD'] = sumdMD/len(uniq)
1610	else: #spherical harmonics
1611	if calcControls[phfx+'SHord']:
1612	POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1613	return POcorr,POderv
1614
1615	def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1616	'Needs a doc string'
1617	if 'Debye' in calcControls[hfx+'instType']:
1618	if 'T' in calcControls[hfx+'histType']:
1619	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1620	else:
1621	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1622	else:
1623	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1624
1625	def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1626	'Needs a doc string'
1627	if 'Debye' in calcControls[hfx+'instType']:
1628	if 'T' in calcControls[hfx+'histType']:
1629	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1630	else:
1631	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1632	else:
1633	return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1634
1635	def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1636	'Needs a doc string'
1637	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1638	pi2 = np.sqrt(2./np.pi)
1639	if 'T' in calcControls[hfx+'histType']:
1640	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1641	wave = refl[14+im]
1642	else: #'C'W
1643	sth2 = sind(refl[5+im]/2.)**2
1644	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1645	c2th = 1.-2.0*sth2
1646	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1647	if 'X' in calcControls[hfx+'histType']:
1648	flv2 = 0.079411(1.0+c2th**2)/2.0
1649	xfac = flv2*parmDict[phfx+'Extinction']
1650	exb = 1.0
1651	if xfac > -1.:
1652	exb = 1./np.sqrt(1.+xfac)
1653	exl = 1.0
1654	if 0 < xfac <= 1.:
1655	xn = np.array([xfac**(i+1) for i in range(6)])
1656	exl += np.sum(xn*coef)
1657	elif xfac > 1.:
1658	xfac2 = 1./np.sqrt(xfac)
1659	exl = pi2(1.-0.125/xfac)xfac2
1660	return exbsth2+exl(1.-sth2)
1661
1662	def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1663	'Needs a doc string'
1664	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1665	pi2 = np.sqrt(2./np.pi)
1666	if 'T' in calcControls[hfx+'histType']:
1667	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1668	wave = refl[14+im]
1669	else: #'C'W
1670	sth2 = sind(refl[5+im]/2.)**2
1671	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1672	c2th = 1.-2.0*sth2
1673	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1674	if 'X' in calcControls[hfx+'histType']:
1675	flv2 = 0.079411(1.0+c2th**2)/2.0
1676	xfac = flv2*parmDict[phfx+'Extinction']
1677	dbde = -500.*flv2
1678	if xfac > -1.:
1679	dbde = -0.5flv2/np.sqrt(1.+xfac)*3
1680	dlde = 0.
1681	if 0 < xfac <= 1.:
1682	xn = np.array([iflv2xfac**i for i in [1,2,3,4,5,6]])
1683	dlde = np.sum(xn*coef)
1684	elif xfac > 1.:
1685	xfac2 = 1./np.sqrt(xfac)
1686	dlde = flv2pi2xfac2(-1./xfac+0.375/xfac*2)
1687
1688	return dbdesth2+dlde(1.-sth2)
1689
1690	def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1691	'Needs a doc string' #need powder extinction!
1692	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3+im] #scale*multiplicity
1693	if 'X' in parmDict[hfx+'Type']:
1694	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1695	POcorr = 1.0
1696	if pfx+'SHorder' in parmDict: #generalized spherical harmonics texture - takes precidence
1697	POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1698	elif calcControls[phfx+'poType'] == 'MD': #March-Dollase
1699	POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1700	elif calcControls[phfx+'SHord']: #cylindrical spherical harmonics
1701	POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1702	Icorr *= POcorr
1703	AbsCorr = 1.0
1704	AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1705	Icorr *= AbsCorr
1706	ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1707	Icorr *= ExtCorr
1708	return Icorr,POcorr,AbsCorr,ExtCorr
1709
1710	def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1711	'Needs a doc string' #need powder extinction derivs!
1712	dIdsh = 1./parmDict[hfx+'Scale']
1713	dIdsp = 1./parmDict[phfx+'Scale']
1714	if 'X' in parmDict[hfx+'Type']:
1715	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1716	dIdPola /= pola
1717	else: #'N'
1718	dIdPola = 0.0
1719	dFdODF = {}
1720	dFdSA = [0,0,0]
1721	dIdPO = {}
1722	if pfx+'SHorder' in parmDict:
1723	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1724	for iSH in dFdODF:
1725	dFdODF[iSH] /= odfCor
1726	for i in range(3):
1727	dFdSA[i] /= odfCor
1728	elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1729	POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1730	for iPO in dIdPO:
1731	dIdPO[iPO] /= POcorr
1732	if 'T' in parmDict[hfx+'Type']:
1733	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1734	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im] #/ext corr
1735	else:
1736	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1737	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im] #/ext corr
1738	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1739
1740	def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1741	'Needs a doc string'
1742	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1743	costh = cosd(refl[5+im]/2.)
1744	#crystallite size
1745	if calcControls[phfx+'SizeType'] == 'isotropic':
1746	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1747	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1748	H = np.array(refl[:3])
1749	P = np.array(calcControls[phfx+'SizeAxis'])
1750	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1751	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1752	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1753	else: #ellipsoidal crystallites
1754	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1755	H = np.array(refl[:3])
1756	lenR = G2pwd.ellipseSize(H,Sij,GB)
1757	Sgam = 1.8wave/(np.picosth*lenR)
1758	#microstrain
1759	if calcControls[phfx+'MustrainType'] == 'isotropic':
1760	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1761	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1762	H = np.array(refl[:3])
1763	P = np.array(calcControls[phfx+'MustrainAxis'])
1764	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1765	Si = parmDict[phfx+'Mustrain;i']
1766	Sa = parmDict[phfx+'Mustrain;a']
1767	Mgam = 0.018SiSatand(refl[5+im]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1768	else: #generalized - P.W. Stephens model
1769	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1770	Sum = 0
1771	for i,strm in enumerate(Strms):
1772	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1773	Mgam = 0.018refl[4+im]2tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1774	elif 'T' in calcControls[hfx+'histType']: #All checked & OK
1775	#crystallite size
1776	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1777	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1778	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1779	H = np.array(refl[:3])
1780	P = np.array(calcControls[phfx+'SizeAxis'])
1781	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1782	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]*2/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a'])
1783	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1784	else: #ellipsoidal crystallites #OK
1785	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1786	H = np.array(refl[:3])
1787	lenR = G2pwd.ellipseSize(H,Sij,GB)
1788	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/lenR
1789	#microstrain
1790	if calcControls[phfx+'MustrainType'] == 'isotropic': #OK
1791	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1792	elif calcControls[phfx+'MustrainType'] == 'uniaxial': #OK
1793	H = np.array(refl[:3])
1794	P = np.array(calcControls[phfx+'MustrainAxis'])
1795	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1796	Si = parmDict[phfx+'Mustrain;i']
1797	Sa = parmDict[phfx+'Mustrain;a']
1798	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa/np.sqrt((SicosP)2+(SasinP)**2)
1799	else: #generalized - P.W. Stephens model OK
1800	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1801	Sum = 0
1802	for i,strm in enumerate(Strms):
1803	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1804	Mgam = 1.e-6parmDict[hfx+'difC']np.sqrt(Sum)refl[4+im]*3
1805
1806	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1807	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1808	sig /= ateln2
1809	return sig,gam
1810
1811	def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1812	'Needs a doc string'
1813	gamDict = {}
1814	sigDict = {}
1815	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1816	costh = cosd(refl[5+im]/2.)
1817	tanth = tand(refl[5+im]/2.)
1818	#crystallite size derivatives
1819	if calcControls[phfx+'SizeType'] == 'isotropic':
1820	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1821	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1822	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1823	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1824	H = np.array(refl[:3])
1825	P = np.array(calcControls[phfx+'SizeAxis'])
1826	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1827	Si = parmDict[phfx+'Size;i']
1828	Sa = parmDict[phfx+'Size;a']
1829	gami = 1.8wave/(costhnp.piSiSa)
1830	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1831	Sgam = gami*sqtrm
1832	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1833	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1834	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1835	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1836	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1837	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1838	else: #ellipsoidal crystallites
1839	const = 1.8wave/(np.picosth)
1840	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1841	H = np.array(refl[:3])
1842	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1843	Sgam = const/lenR
1844	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1845	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1846	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1847	gamDict[phfx+'Size;mx'] = Sgam
1848	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1849
1850	#microstrain derivatives
1851	if calcControls[phfx+'MustrainType'] == 'isotropic':
1852	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1853	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1854	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1855	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1856	H = np.array(refl[:3])
1857	P = np.array(calcControls[phfx+'MustrainAxis'])
1858	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1859	Si = parmDict[phfx+'Mustrain;i']
1860	Sa = parmDict[phfx+'Mustrain;a']
1861	gami = 0.018SiSa*tanth/np.pi
1862	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1863	Mgam = gami/sqtrm
1864	dsi = -gamiSicosP2/sqtrm3
1865	dsa = -gamiSasinP2/sqtrm3
1866	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1867	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1868	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1869	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1870	else: #generalized - P.W. Stephens model
1871	const = 0.018refl[4+im]2tanth/np.pi
1872	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1873	Sum = 0
1874	for i,strm in enumerate(Strms):
1875	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1876	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1877	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1878	Mgam = const*np.sqrt(Sum)
1879	for i in range(len(Strms)):
1880	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1881	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1882	gamDict[phfx+'Mustrain;mx'] = Mgam
1883	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1884	else: #'T'OF - All checked & OK
1885	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1886	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1887	gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1888	sigDict[phfx+'Size;i'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])*2/(ateln2parmDict[phfx+'Size;i'])
1889	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1890	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1891	H = np.array(refl[:3])
1892	P = np.array(calcControls[phfx+'SizeAxis'])
1893	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1894	Si = parmDict[phfx+'Size;i']
1895	Sa = parmDict[phfx+'Size;a']
1896	gami = const/(Si*Sa)
1897	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1898	Sgam = gami*sqtrm
1899	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1900	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1901	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1902	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1903	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1904	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1905	else: #OK ellipsoidal crystallites
1906	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1907	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1908	H = np.array(refl[:3])
1909	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1910	Sgam = const/lenR
1911	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1912	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1913	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1914	gamDict[phfx+'Size;mx'] = Sgam #OK
1915	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2 #OK
1916
1917	#microstrain derivatives
1918	if calcControls[phfx+'MustrainType'] == 'isotropic':
1919	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1920	gamDict[phfx+'Mustrain;i'] = 1.e-6refl[4+im]parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] #OK
1921	sigDict[phfx+'Mustrain;i'] = 2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])*2/(ateln2parmDict[phfx+'Mustrain;i'])
1922	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1923	H = np.array(refl[:3])
1924	P = np.array(calcControls[phfx+'MustrainAxis'])
1925	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1926	Si = parmDict[phfx+'Mustrain;i']
1927	Sa = parmDict[phfx+'Mustrain;a']
1928	gami = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa
1929	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1930	Mgam = gami/sqtrm
1931	dsi = -gamiSicosP2/sqtrm3
1932	dsa = -gamiSasinP2/sqtrm3
1933	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1934	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1935	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1936	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1937	else: #generalized - P.W. Stephens model OK
1938	pwrs = calcControls[phfx+'MuPwrs']
1939	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1940	const = 1.e-6parmDict[hfx+'difC']refl[4+im]**3
1941	Sum = 0
1942	for i,strm in enumerate(Strms):
1943	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1944	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1945	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1946	Mgam = const*np.sqrt(Sum)
1947	for i in range(len(Strms)):
1948	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1949	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1950	gamDict[phfx+'Mustrain;mx'] = Mgam
1951	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1952
1953	return sigDict,gamDict
1954
1955	def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1956	'Needs a doc string'
1957	if im:
1958	h,k,l,m = refl[:4]
1959	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1960	d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1961	else:
1962	h,k,l = refl[:3]
1963	d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1964	refl[4+im] = d
1965	if 'C' in calcControls[hfx+'histType']:
1966	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1967	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1968	if 'Bragg' in calcControls[hfx+'instType']:
1969	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1970	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1971	else: #Debye-Scherrer - simple but maybe not right
1972	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1973	elif 'T' in calcControls[hfx+'histType']:
1974	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1975	#do I need sample position effects - maybe?
1976	return pos
1977
1978	def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1979	'Needs a doc string'
1980	dpr = 180./np.pi
1981	if im:
1982	h,k,l,m = refl[:4]
1983	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1984	dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1985	else:
1986	m = 0
1987	h,k,l = refl[:3]
1988	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1989	dst = np.sqrt(dstsq)
1990	dsp = 1./dst
1991	if 'C' in calcControls[hfx+'histType']:
1992	pos = refl[5+im]-parmDict[hfx+'Zero']
1993	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1994	dpdw = const*dst
1995	dpdA = np.array([h2,k2,l*2,hk,hl,kl])constwave/(2.0*dst)
1996	dpdZ = 1.0
1997	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
1998	2lA[2]+hA[4]+kA[5]])mconstwave/(2.0dst)
1999	shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2000	if 'Bragg' in calcControls[hfx+'instType']:
2001	dpdSh = -4.shftcosd(pos/2.0)
2002	dpdTr = -shftsind(pos)100.0
2003	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2004	else: #Debye-Scherrer - simple but maybe not right
2005	dpdXd = -shft*cosd(pos)
2006	dpdYd = -shft*sind(pos)
2007	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2008	elif 'T' in calcControls[hfx+'histType']:
2009	dpdA = -np.array([h2,k2,l*2,hk,hl,kl])parmDict[hfx+'difC']dsp**3/2.
2010	dpdZ = 1.0
2011	dpdDC = dsp
2012	dpdDA = dsp**2
2013	dpdDB = 1./dsp
2014	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
2015	2lA[2]+hA[4]+kA[5]])mparmDict[hfx+'difC']dsp**3/2.
2016	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2017
2018	def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2019	'Needs a doc string'
2020	laue = SGData['SGLaue']
2021	uniq = SGData['SGUniq']
2022	h,k,l = refl[:3]
2023	if laue in ['m3','m3m']:
2024	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2025	refl[4+im]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2026	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2027	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2028	elif laue in ['3R','3mR']:
2029	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2030	elif laue in ['4/m','4/mmm']:
2031	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2032	elif laue in ['mmm']:
2033	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2034	elif laue in ['2/m']:
2035	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2036	if uniq == 'a':
2037	Dij += parmDict[phfx+'D23']kl
2038	elif uniq == 'b':
2039	Dij += parmDict[phfx+'D13']hl
2040	elif uniq == 'c':
2041	Dij += parmDict[phfx+'D12']hk
2042	else:
2043	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2044	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2045	if 'C' in calcControls[hfx+'histType']:
2046	return -180.Dijrefl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2047	else:
2048	return -DijparmDict[hfx+'difC']refl[4+im]**2/2.
2049
2050	def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2051	'Needs a doc string'
2052	laue = SGData['SGLaue']
2053	uniq = SGData['SGUniq']
2054	h,k,l = refl[:3]
2055	if laue in ['m3','m3m']:
2056	dDijDict = {phfx+'D11':h2+k2+l**2,
2057	phfx+'eA':refl[4+im]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2058	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2059	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2060	elif laue in ['3R','3mR']:
2061	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2062	elif laue in ['4/m','4/mmm']:
2063	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2064	elif laue in ['mmm']:
2065	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2066	elif laue in ['2/m']:
2067	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2068	if uniq == 'a':
2069	dDijDict[phfx+'D23'] = k*l
2070	elif uniq == 'b':
2071	dDijDict[phfx+'D13'] = h*l
2072	elif uniq == 'c':
2073	dDijDict[phfx+'D12'] = h*k
2074	else:
2075	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2076	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2077	if 'C' in calcControls[hfx+'histType']:
2078	for item in dDijDict:
2079	dDijDict[item] = 180.0refl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2080	else:
2081	for item in dDijDict:
2082	dDijDict[item] = -parmDict[hfx+'difC']refl[4+im]**3/2.
2083	return dDijDict
2084
2085	def GetDij(phfx,SGData,parmDict):
2086	HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2087	return G2spc.HStrainVals(HSvals,SGData)
2088
2089	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2090	'Needs a doc string'
2091	histoList = Histograms.keys()
2092	histoList.sort()
2093	for histogram in histoList:
2094	if 'PWDR' in histogram[:4]:
2095	Histogram = Histograms[histogram]
2096	hId = Histogram['hId']
2097	hfx = ':%d:'%(hId)
2098	Limits = calcControls[hfx+'Limits']
2099	if 'C' in calcControls[hfx+'histType']:
2100	shl = max(parmDict[hfx+'SH/L'],0.0005)
2101	Ka2 = False
2102	kRatio = 0.0
2103	if hfx+'Lam1' in parmDict.keys():
2104	Ka2 = True
2105	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2106	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2107	x,y,w,yc,yb,yd = Histogram['Data']
2108	xB = np.searchsorted(x,Limits[0])
2109	xF = np.searchsorted(x,Limits[1])
2110	ymb = np.array(y-yb)
2111	ymb = np.where(ymb,ymb,1.0)
2112	ycmb = np.array(yc-yb)
2113	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2114	refLists = Histogram['Reflection Lists']
2115	for phase in refLists:
2116	if phase not in Phases: #skips deleted or renamed phases silently!
2117	continue
2118	Phase = Phases[phase]
2119	im = 0
2120	if Phase['General']['Type'] in ['modulated','magnetic']:
2121	im = 1
2122	pId = Phase['pId']
2123	phfx = '%d:%d:'%(pId,hId)
2124	refDict = refLists[phase]
2125	sumFo = 0.0
2126	sumdF = 0.0
2127	sumFosq = 0.0
2128	sumdFsq = 0.0
2129	sumInt = 0.0
2130	for refl in refDict['RefList']:
2131	if 'C' in calcControls[hfx+'histType']:
2132	yp = np.zeros_like(yb)
2133	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2134	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2135	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2136	iFin2 = iFin
2137	if not iBeg+iFin: #peak below low limit - skip peak
2138	continue
2139	elif not iBeg-iFin: #peak above high limit - done
2140	break
2141	elif iBeg < iFin:
2142	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2143	sumInt += refl[11+im]*refl[9+im]
2144	if Ka2:
2145	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2146	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2147	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2148	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2149	if iFin2 > iBeg2:
2150	yp[iBeg2:iFin2] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2])) #and here
2151	sumInt += refl[11+im]refl[9+im]kRatio
2152	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11+im](1.+kRatio)),0.0))
2153
2154	elif 'T' in calcControls[hfx+'histType']:
2155	yp = np.zeros_like(yb)
2156	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2157	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2158	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2159	if iBeg < iFin:
2160	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2161	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
2162	sumInt += refl[11+im]*refl[9+im]
2163	Fo = np.sqrt(np.abs(refl[8+im]))
2164	Fc = np.sqrt(np.abs(refl[9]+im))
2165	sumFo += Fo
2166	sumFosq += refl[8+im]**2
2167	sumdF += np.abs(Fo-Fc)
2168	sumdFsq += (refl[8+im]-refl[9+im])**2
2169	if sumFo:
2170	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2171	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2172	else:
2173	Histogram['Residuals'][phfx+'Rf'] = 100.
2174	Histogram['Residuals'][phfx+'Rf^2'] = 100.
2175	Histogram['Residuals'][phfx+'sumInt'] = sumInt
2176	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
2177	Histogram['Residuals']['hId'] = hId
2178	elif 'HKLF' in histogram[:4]:
2179	Histogram = Histograms[histogram]
2180	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
2181
2182	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2183	'Needs a doc string'
2184
2185	def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
2186	U = parmDict[hfx+'U']
2187	V = parmDict[hfx+'V']
2188	W = parmDict[hfx+'W']
2189	X = parmDict[hfx+'X']
2190	Y = parmDict[hfx+'Y']
2191	tanPos = tand(refl[5+im]/2.0)
2192	Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2193	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
2194	sig = max(0.001,sig)
2195	gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam #save peak gamma
2196	gam = max(0.001,gam)
2197	return sig,gam
2198
2199	def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
2200	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4+im]*2+ \
2201	parmDict[hfx+'sig-2']refl[4+im]4+parmDict[hfx+'sig-q']/refl[4+im]*2
2202	gam = parmDict[hfx+'X']refl[4+im]+parmDict[hfx+'Y']refl[4+im]**2
2203	Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2204	sig += Ssig
2205	gam += Sgam
2206	return sig,gam
2207
2208	def GetReflAlpBet(refl,im,hfx,parmDict):
2209	alp = parmDict[hfx+'alpha']/refl[4+im]
2210	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]4+parmDict[hfx+'beta-q']/refl[4+im]2
2211	return alp,bet
2212
2213	hId = Histogram['hId']
2214	hfx = ':%d:'%(hId)
2215	bakType = calcControls[hfx+'bakType']
2216	yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2217	yc = np.zeros_like(yb)
2218	cw = np.diff(x)
2219	cw = np.append(cw,cw[-1])
2220
2221	if 'C' in calcControls[hfx+'histType']:
2222	shl = max(parmDict[hfx+'SH/L'],0.002)
2223	Ka2 = False
2224	if hfx+'Lam1' in parmDict.keys():
2225	wave = parmDict[hfx+'Lam1']
2226	Ka2 = True
2227	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2228	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2229	else:
2230	wave = parmDict[hfx+'Lam']
2231	for phase in Histogram['Reflection Lists']:
2232	refDict = Histogram['Reflection Lists'][phase]
2233	if phase not in Phases: #skips deleted or renamed phases silently!
2234	continue
2235	Phase = Phases[phase]
2236	pId = Phase['pId']
2237	pfx = '%d::'%(pId)
2238	phfx = '%d:%d:'%(pId,hId)
2239	hfx = ':%d:'%(hId)
2240	SGData = Phase['General']['SGData']
2241	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2242	im = 0
2243	if Phase['General']['Type'] in ['modulated','magnetic']:
2244	SSGData = Phase['General']['SSGData']
2245	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2246	im = 1 #offset in SS reflection list
2247	#??
2248	Dij = GetDij(phfx,SGData,parmDict)
2249	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2250	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2251	if np.any(np.diag(G)<0.):
2252	raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2253	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2254	Vst = np.sqrt(nl.det(G)) #V*
2255	if not Phase['General'].get('doPawley'):
2256	time0 = time.time()
2257	if im:
2258	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2259	else:
2260	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2261	# print 'sf calc time: %.3fs'%(time.time()-time0)
2262	time0 = time.time()
2263	badPeak = False
2264	for iref,refl in enumerate(refDict['RefList']):
2265	if 'C' in calcControls[hfx+'histType']:
2266	if im:
2267	h,k,l,m = refl[:4]
2268	else:
2269	h,k,l = refl[:3]
2270	Uniq = np.inner(refl[:3],SGMT)
2271	refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2272	Lorenz = 1./(2.sind(refl[5+im]/2.)2cosd(refl[5+im]/2.)) #Lorentz correction
2273	refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2274	refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2275	refl[11+im] = VstLorenz
2276
2277	if Phase['General'].get('doPawley'):
2278	try:
2279	if im:
2280	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2281	else:
2282	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2283	refl[9+im] = parmDict[pInd]
2284	except KeyError:
2285	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2286	continue
2287	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2288	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2289	iFin = np.searchsorted(x,refl[5+im]+fmax)
2290	if not iBeg+iFin: #peak below low limit - skip peak
2291	continue
2292	elif not iBeg-iFin: #peak above high limit - done
2293	break
2294	elif iBeg > iFin: #bad peak coeff - skip
2295	badPeak = True
2296	continue
2297	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2298	if Ka2:
2299	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2300	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2301	iBeg = np.searchsorted(x,pos2-fmin)
2302	iFin = np.searchsorted(x,pos2+fmax)
2303	if not iBeg+iFin: #peak below low limit - skip peak
2304	continue
2305	elif not iBeg-iFin: #peak above high limit - done
2306	return yc,yb
2307	elif iBeg > iFin: #bad peak coeff - skip
2308	continue
2309	yc[iBeg:iFin] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #and here
2310	elif 'T' in calcControls[hfx+'histType']:
2311	h,k,l = refl[:3]
2312	Uniq = np.inner(refl[:3],SGMT)
2313	refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2314	Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)refl[4+im]*4 #TOF Lorentz correction
2315	# refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
2316	refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2317	refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2318	refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2319	refl[11+im] = VstLorenz
2320	if Phase['General'].get('doPawley'):
2321	try:
2322	if im:
2323	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2324	else:
2325	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2326	refl[9+im] = parmDict[pInd]
2327	except KeyError:
2328	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2329	continue
2330	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2331	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2332	iFin = np.searchsorted(x,refl[5+im]+fmax)
2333	if not iBeg+iFin: #peak below low limit - skip peak
2334	continue
2335	elif not iBeg-iFin: #peak above high limit - done
2336	break
2337	elif iBeg > iFin: #bad peak coeff - skip
2338	badPeak = True
2339	continue
2340	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
2341	# print 'profile calc time: %.3fs'%(time.time()-time0)
2342	if badPeak:
2343	print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
2344	return yc,yb
2345
2346	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2347	'Needs a doc string'
2348
2349	def cellVaryDerv(pfx,SGData,dpdA):
2350	if SGData['SGLaue'] in ['-1',]:
2351	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2352	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2353	elif SGData['SGLaue'] in ['2/m',]:
2354	if SGData['SGUniq'] == 'a':
2355	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2356	elif SGData['SGUniq'] == 'b':
2357	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2358	else:
2359	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2360	elif SGData['SGLaue'] in ['mmm',]:
2361	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2362	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2363	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2364	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2365	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2366	elif SGData['SGLaue'] in ['3R', '3mR']:
2367	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2368	elif SGData['SGLaue'] in ['m3m','m3']:
2369	return [[pfx+'A0',dpdA[0]]]
2370
2371	# create a list of dependent variables and set up a dictionary to hold their derivatives
2372	dependentVars = G2mv.GetDependentVars()
2373	depDerivDict = {}
2374	for j in dependentVars:
2375	depDerivDict[j] = np.zeros(shape=(len(x)))
2376	#print 'dependent vars',dependentVars
2377	lenX = len(x)
2378	hId = Histogram['hId']
2379	hfx = ':%d:'%(hId)
2380	bakType = calcControls[hfx+'bakType']
2381	dMdv = np.zeros(shape=(len(varylist),len(x)))
2382	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2383	if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2384	bBpos =varylist.index(hfx+'Back;0')
2385	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2386	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2387	for name in varylist:
2388	if 'Debye' in name:
2389	id = int(name.split(';')[-1])
2390	parm = name[:int(name.rindex(';'))]
2391	ip = names.index(parm)
2392	dMdv[varylist.index(name)] = dMddb[3*id+ip]
2393	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2394	for name in varylist:
2395	if 'BkPk' in name:
2396	parm,id = name.split(';')
2397	id = int(id)
2398	if parm in names:
2399	ip = names.index(parm)
2400	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2401	cw = np.diff(x)
2402	cw = np.append(cw,cw[-1])
2403	Ka2 = False #also for TOF!
2404	if 'C' in calcControls[hfx+'histType']:
2405	shl = max(parmDict[hfx+'SH/L'],0.002)
2406	if hfx+'Lam1' in parmDict.keys():
2407	wave = parmDict[hfx+'Lam1']
2408	Ka2 = True
2409	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2410	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2411	else:
2412	wave = parmDict[hfx+'Lam']
2413	for phase in Histogram['Reflection Lists']:
2414	refDict = Histogram['Reflection Lists'][phase]
2415	if phase not in Phases: #skips deleted or renamed phases silently!
2416	continue
2417	Phase = Phases[phase]
2418	SGData = Phase['General']['SGData']
2419	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2420	im = 0
2421	if Phase['General']['Type'] in ['modulated','magnetic']:
2422	SSGData = Phase['General']['SSGData']
2423	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2424	im = 1 #offset in SS reflection list
2425	#??
2426	pId = Phase['pId']
2427	pfx = '%d::'%(pId)
2428	phfx = '%d:%d:'%(pId,hId)
2429	Dij = GetDij(phfx,SGData,parmDict)
2430	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2431	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2432	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2433	if not Phase['General'].get('doPawley'):
2434	time0 = time.time()
2435	if im:
2436	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2437	else:
2438	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2439	# print 'sf-derv time %.3fs'%(time.time()-time0)
2440	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2441	time0 = time.time()
2442	for iref,refl in enumerate(refDict['RefList']):
2443	if im:
2444	h,k,l,m = refl[:4]
2445	else:
2446	h,k,l = refl[:3]
2447	Uniq = np.inner(refl[:3],SGMT)
2448	if 'T' in calcControls[hfx+'histType']:
2449	wave = refl[14+im]
2450	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2451	if 'C' in calcControls[hfx+'histType']: #CW powder
2452	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2453	else: #'T'OF
2454	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2455	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2456	iFin = np.searchsorted(x,refl[5+im]+fmax)
2457	if not iBeg+iFin: #peak below low limit - skip peak
2458	continue
2459	elif not iBeg-iFin: #peak above high limit - done
2460	break
2461	pos = refl[5+im]
2462	if 'C' in calcControls[hfx+'histType']:
2463	tanth = tand(pos/2.0)
2464	costh = cosd(pos/2.0)
2465	lenBF = iFin-iBeg
2466	dMdpk = np.zeros(shape=(6,lenBF))
2467	dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2468	for i in range(5):
2469	dMdpk[i] += 100.cw[iBeg:iFin]refl[11+im]refl[9+im]dMdipk[i]
2470	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2471	if Ka2:
2472	pos2 = refl[5+im]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2473	iBeg2 = np.searchsorted(x,pos2-fmin)
2474	iFin2 = np.searchsorted(x,pos2+fmax)
2475	if iBeg2-iFin2:
2476	lenBF2 = iFin2-iBeg2
2477	dMdpk2 = np.zeros(shape=(6,lenBF2))
2478	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2479	for i in range(5):
2480	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11+im]refl[9+im]kRatio*dMdipk2[i]
2481	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11+im]*dMdipk2[0]
2482	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2483	else: #'T'OF
2484	lenBF = iFin-iBeg
2485	if lenBF < 0: #bad peak coeff
2486	break
2487	dMdpk = np.zeros(shape=(6,lenBF))
2488	dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2489	for i in range(6):
2490	dMdpk[i] += refl[11+im]refl[9+im]dMdipk[i] #cw[iBeg:iFin]*
2491	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
2492	if Phase['General'].get('doPawley'):
2493	dMdpw = np.zeros(len(x))
2494	try:
2495	if im:
2496	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2497	else:
2498	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2499	idx = varylist.index(pIdx)
2500	dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2501	if Ka2: #not for TOF either
2502	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2503	dMdv[idx] = dMdpw
2504	except: # ValueError:
2505	pass
2506	if 'C' in calcControls[hfx+'histType']:
2507	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2508	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2509	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2510	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2511	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2512	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2513	hfx+'DisplaceY':[dpdY,'pos'],}
2514	if 'Bragg' in calcControls[hfx+'instType']:
2515	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2516	hfx+'SurfRoughB':[dFdAb[1],'int'],})
2517	else:
2518	names.update({hfx+'Absorption':[dFdAb,'int'],})
2519	else: #'T'OF
2520	dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2521	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2522	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2523	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2524	hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2525	hfx+'beta-q':[1./refl[4+im]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]2,'sig'],
2526	hfx+'sig-2':[refl[4+im]4,'sig'],hfx+'sig-q':[1./refl[4+im]2,'sig'],
2527	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2528	for name in names:
2529	item = names[name]
2530	if name in varylist:
2531	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2532	if Ka2 and iFin2-iBeg2:
2533	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2534	elif name in dependentVars:
2535	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2536	if Ka2 and iFin2-iBeg2:
2537	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2538	for iPO in dIdPO:
2539	if iPO in varylist:
2540	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2541	if Ka2 and iFin2-iBeg2:
2542	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2543	elif iPO in dependentVars:
2544	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2545	if Ka2 and iFin2-iBeg2:
2546	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2547	for i,name in enumerate(['omega','chi','phi']):
2548	aname = pfx+'SH '+name
2549	if aname in varylist:
2550	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2551	if Ka2 and iFin2-iBeg2:
2552	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2553	elif aname in dependentVars:
2554	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2555	if Ka2 and iFin2-iBeg2:
2556	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2557	for iSH in dFdODF:
2558	if iSH in varylist:
2559	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2560	if Ka2 and iFin2-iBeg2:
2561	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2562	elif iSH in dependentVars:
2563	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2564	if Ka2 and iFin2-iBeg2:
2565	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2566	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2567	for name,dpdA in cellDervNames:
2568	if name in varylist:
2569	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2570	if Ka2 and iFin2-iBeg2:
2571	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2572	elif name in dependentVars:
2573	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2574	if Ka2 and iFin2-iBeg2:
2575	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2576	dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2577	for name in dDijDict:
2578	if name in varylist:
2579	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2580	if Ka2 and iFin2-iBeg2:
2581	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2582	elif name in dependentVars:
2583	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2584	if Ka2 and iFin2-iBeg2:
2585	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2586	for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2587	if name in varylist:
2588	dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2589	if Ka2 and iFin2-iBeg2:
2590	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2591	elif name in dependentVars:
2592	depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2593	if Ka2 and iFin2-iBeg2:
2594	depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2595	if 'C' in calcControls[hfx+'histType']:
2596	sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2597	else: #'T'OF
2598	sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2599	for name in gamDict:
2600	if name in varylist:
2601	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2602	if Ka2 and iFin2-iBeg2:
2603	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2604	elif name in dependentVars:
2605	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2606	if Ka2 and iFin2-iBeg2:
2607	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2608	for name in sigDict:
2609	if name in varylist:
2610	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2611	if Ka2 and iFin2-iBeg2:
2612	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2613	elif name in dependentVars:
2614	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2615	if Ka2 and iFin2-iBeg2:
2616	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2617	for name in ['BabA','BabU']:
2618	if refl[9+im]:
2619	if phfx+name in varylist:
2620	dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2621	if Ka2 and iFin2-iBeg2:
2622	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2623	elif phfx+name in dependentVars:
2624	depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2625	if Ka2 and iFin2-iBeg2:
2626	depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2627	if not Phase['General'].get('doPawley'):
2628	#do atom derivatives - for RB,F,X & U so far
2629	corr = dervDict['int']/refl[9+im]
2630	if Ka2 and iFin2-iBeg2:
2631	corr2 = dervDict2['int']/refl[9+im]
2632	for name in varylist+dependentVars:
2633	if '::RBV;' in name:
2634	pass
2635	else:
2636	try:
2637	aname = name.split(pfx)[1][:2]
2638	if aname not in ['Af','dA','AU','RB','Xs','Xc','Ys','Yc','Zs','Zc','U1','U2','U3']: continue # skip anything not an atom or rigid body param
2639	except IndexError:
2640	continue
2641	if name in varylist:
2642	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2643	if Ka2 and iFin2-iBeg2:
2644	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2645	elif name in dependentVars:
2646	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2647	if Ka2 and iFin2-iBeg2:
2648	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2649	# print 'profile derv time: %.3fs'%(time.time()-time0)
2650	# now process derivatives in constraints
2651	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2652	return dMdv
2653
2654	def UserRejectHKL(ref,im,userReject):
2655	if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2656	return False
2657	elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2658	return False
2659	elif ref[11+im] < userReject['MinExt']:
2660	return False
2661	elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2662	return False
2663	return True
2664
2665	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2666	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2667	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
2668	are returned to either dervRefine or HessRefine.
2669
2670	:returns:
2671	'''
2672	nobs = Histogram['Residuals']['Nobs']
2673	hId = Histogram['hId']
2674	hfx = ':%d:'%(hId)
2675	pfx = '%d::'%(Phase['pId'])
2676	phfx = '%d:%d:'%(Phase['pId'],hId)
2677	SGData = Phase['General']['SGData']
2678	im = 0
2679	if Phase['General']['Type'] in ['modulated','magnetic']:
2680	SSGData = Phase['General']['SSGData']
2681	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2682	im = 1 #offset in SS reflection list
2683	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2684	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2685	refDict = Histogram['Data']
2686	if parmDict[phfx+'Scale'] < 0.:
2687	parmDict[phfx+'Scale'] = .001
2688	if im:
2689	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2690	else:
2691	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2692	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2693	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2694	dependentVars = G2mv.GetDependentVars()
2695	depDerivDict = {}
2696	for j in dependentVars:
2697	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2698	wdf = np.zeros(len(refDict['RefList']))
2699	if calcControls['F**2']:
2700	for iref,ref in enumerate(refDict['RefList']):
2701	if ref[6+im] > 0:
2702	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2703	w = 1.0/ref[6+im]
2704	if ref[3+im] > 0:
2705	wdf[iref] = w*(ref[5+im]-ref[7+im])
2706	for j,var in enumerate(varylist):
2707	if var in dFdvDict:
2708	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2709	for var in dependentVars:
2710	if var in dFdvDict:
2711	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2712	if phfx+'Scale' in varylist:
2713	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2714	elif phfx+'Scale' in dependentVars:
2715	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2716	for item in ['Ep','Es','Eg']:
2717	if phfx+item in varylist and phfx+item in dervDict:
2718	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2719	elif phfx+item in dependentVars and phfx+item in dervDict:
2720	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2721	for item in ['BabA','BabU']:
2722	if phfx+item in varylist:
2723	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2724	elif phfx+item in dependentVars:
2725	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2726	else: #F refinement
2727	for iref,ref in enumerate(refDict['RefList']):
2728	if ref[5+im] > 0.:
2729	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2730	Fo = np.sqrt(ref[5+im])
2731	Fc = np.sqrt(ref[7+im])
2732	w = 1.0/ref[6+im]
2733	if ref[3+im] > 0:
2734	wdf[iref] = 2.0Fcw*(Fo-Fc)
2735	for j,var in enumerate(varylist):
2736	if var in dFdvDict:
2737	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2738	for var in dependentVars:
2739	if var in dFdvDict:
2740	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2741	if phfx+'Scale' in varylist:
2742	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2743	elif phfx+'Scale' in dependentVars:
2744	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2745	for item in ['Ep','Es','Eg']: #OK!
2746	if phfx+item in varylist and phfx+item in dervDict:
2747	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]
2748	elif phfx+item in dependentVars and phfx+item in dervDict:
2749	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2750	for item in ['BabA','BabU']:
2751	if phfx+item in varylist:
2752	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2753	elif phfx+item in dependentVars:
2754	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2755	return dMdvh,depDerivDict,wdf
2756
2757
2758	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2759	'''Loop over histograms and compute derivatives of the fitting
2760	model (M) with respect to all parameters. Results are returned in
2761	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2762	n is the number of parameters and m is the number of data
2763	points. This can exceed memory when m gets large. This routine is
2764	used when refinement derivatives are selected as "analtytic
2765	Jacobian" in Controls.
2766
2767	:returns: Jacobian numpy.array dMdv for all histograms concatinated
2768	'''
2769	parmDict.update(zip(varylist,values))
2770	G2mv.Dict2Map(parmDict,varylist)
2771	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2772	nvar = len(varylist)
2773	dMdv = np.empty(0)
2774	histoList = Histograms.keys()
2775	histoList.sort()
2776	for histogram in histoList:
2777	if 'PWDR' in histogram[:4]:
2778	Histogram = Histograms[histogram]
2779	hId = Histogram['hId']
2780	hfx = ':%d:'%(hId)
2781	wtFactor = calcControls[hfx+'wtFactor']
2782	Limits = calcControls[hfx+'Limits']
2783	x,y,w,yc,yb,yd = Histogram['Data']
2784	xB = np.searchsorted(x,Limits[0])
2785	xF = np.searchsorted(x,Limits[1])
2786	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2787	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2788	elif 'HKLF' in histogram[:4]:
2789	Histogram = Histograms[histogram]
2790	phase = Histogram['Reflection Lists']
2791	Phase = Phases[phase]
2792	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2793	hfx = ':%d:'%(Histogram['hId'])
2794	wtFactor = calcControls[hfx+'wtFactor']
2795	# now process derivatives in constraints
2796	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2797	else:
2798	continue #skip non-histogram entries
2799	if len(dMdv):
2800	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2801	else:
2802	dMdv = np.sqrt(wtFactor)*dMdvh
2803
2804	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2805	if np.any(pVals):
2806	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2807	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2808
2809	return dMdv
2810
2811	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2812	'''Loop over histograms and compute derivatives of the fitting
2813	model (M) with respect to all parameters. For each histogram, the
2814	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2815	number of parameters and m is the number of data points *in the
2816	histogram*. The (n by n) Hessian is computed from each Jacobian
2817	and it is returned. This routine is used when refinement
2818	derivatives are selected as "analtytic Hessian" in Controls.
2819
2820	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2821	'''
2822	parmDict.update(zip(varylist,values))
2823	G2mv.Dict2Map(parmDict,varylist)
2824	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2825	#fixup H atom positions here?
2826	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
2827	nvar = len(varylist)
2828	Hess = np.empty(0)
2829	histoList = Histograms.keys()
2830	histoList.sort()
2831	for histogram in histoList:
2832	if 'PWDR' in histogram[:4]:
2833	Histogram = Histograms[histogram]
2834	hId = Histogram['hId']
2835	hfx = ':%d:'%(hId)
2836	wtFactor = calcControls[hfx+'wtFactor']
2837	Limits = calcControls[hfx+'Limits']
2838	x,y,w,yc,yb,yd = Histogram['Data']
2839	W = wtFactor*w
2840	dy = y-yc
2841	xB = np.searchsorted(x,Limits[0])
2842	xF = np.searchsorted(x,Limits[1])
2843	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2844	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2845	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2846	Dy = dy[xB:xF][np.newaxis,:]
2847	dMdvh *= Wt
2848	if dlg:
2849	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2850	if len(Hess):
2851	Hess += np.inner(dMdvh,dMdvh)
2852	dMdvh = WtDy
2853	Vec += np.sum(dMdvh,axis=1)
2854	else:
2855	Hess = np.inner(dMdvh,dMdvh)
2856	dMdvh = WtDy
2857	Vec = np.sum(dMdvh,axis=1)
2858	elif 'HKLF' in histogram[:4]:
2859	Histogram = Histograms[histogram]
2860	phase = Histogram['Reflection Lists']
2861	Phase = Phases[phase]
2862	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2863	hId = Histogram['hId']
2864	hfx = ':%d:'%(Histogram['hId'])
2865	wtFactor = calcControls[hfx+'wtFactor']
2866	# now process derivatives in constraints
2867	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2868	# print 'matrix build time: %.3f'%(time.time()-time0)
2869
2870	if dlg:
2871	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2872	if len(Hess):
2873	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
2874	Hess += wtFactor*np.inner(dMdvh,dMdvh)
2875	else:
2876	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
2877	Hess = wtFactor*np.inner(dMdvh,dMdvh)
2878	else:
2879	continue #skip non-histogram entries
2880	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2881	if np.any(pVals):
2882	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2883	Vec += np.sum(dpdvpWtpVals,axis=1)
2884	Hess += np.inner(dpdv*pWt,dpdv)
2885	return Vec,Hess
2886
2887	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):
2888	'''Computes the point-by-point discrepancies between every data point in every histogram
2889	and the observed value
2890
2891	:returns: an np array of differences between observed and computed diffraction values.
2892	'''
2893	Values2Dict(parmDict, varylist, values)
2894	G2mv.Dict2Map(parmDict,varylist)
2895	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2896	M = np.empty(0)
2897	SumwYo = 0
2898	Nobs = 0
2899	Nrej = 0
2900	Next = 0
2901	ApplyRBModels(parmDict,Phases,rigidbodyDict)
2902	#fixup Hatom positions here....
2903	histoList = Histograms.keys()
2904	histoList.sort()
2905	for histogram in histoList:
2906	if 'PWDR' in histogram[:4]:
2907	Histogram = Histograms[histogram]
2908	hId = Histogram['hId']
2909	hfx = ':%d:'%(hId)
2910	wtFactor = calcControls[hfx+'wtFactor']
2911	Limits = calcControls[hfx+'Limits']
2912	x,y,w,yc,yb,yd = Histogram['Data']
2913	yc *= 0.0 #zero full calcd profiles
2914	yb *= 0.0
2915	yd *= 0.0
2916	xB = np.searchsorted(x,Limits[0])
2917	xF = np.searchsorted(x,Limits[1])
2918	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2919	varylist,Histogram,Phases,calcControls,pawleyLookup)
2920	yc[xB:xF] += yb[xB:xF]
2921	if not np.any(y): #fill dummy data
2922	rv = st.poisson(yc[xB:xF])
2923	y[xB:xF] = rv.rvs()
2924	Z = np.ones_like(yc[xB:xF])
2925	Z[1::2] *= -1
2926	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2927	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2928	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2929	W = wtFactor*w
2930	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2931	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2932	Nobs += Histogram['Residuals']['Nobs']
2933	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2934	SumwYo += Histogram['Residuals']['sumwYo']
2935	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2936	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2937	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2938	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2939	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.))
2940	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.))
2941	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2942	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2943	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2944	if dlg:
2945	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2946	M = np.concatenate((M,wdy))
2947	#end of PWDR processing
2948	elif 'HKLF' in histogram[:4]:
2949	Histogram = Histograms[histogram]
2950	Histogram['Residuals'] = {}
2951	phase = Histogram['Reflection Lists']
2952	Phase = Phases[phase]
2953	hId = Histogram['hId']
2954	hfx = ':%d:'%(hId)
2955	wtFactor = calcControls[hfx+'wtFactor']
2956	pfx = '%d::'%(Phase['pId'])
2957	phfx = '%d:%d:'%(Phase['pId'],hId)
2958	SGData = Phase['General']['SGData']
2959	im = 0
2960	if parmDict[phfx+'Scale'] < 0.:
2961	parmDict[phfx+'Scale'] = .001
2962	if Phase['General']['Type'] in ['modulated','magnetic']:
2963	SSGData = Phase['General']['SSGData']
2964	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2965	im = 1 #offset in SS reflection list
2966	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2967	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2968	refDict = Histogram['Data']
2969	time0 = time.time()
2970	if im:
2971	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2972	else:
2973	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2974	# print 'sf-calc time: %.3f'%(time.time()-time0)
2975	df = np.zeros(len(refDict['RefList']))
2976	sumwYo = 0
2977	sumFo = 0
2978	sumFo2 = 0
2979	sumdF = 0
2980	sumdF2 = 0
2981	if im:
2982	sumSSFo = np.zeros(10)
2983	sumSSFo2 = np.zeros(10)
2984	sumSSdF = np.zeros(10)
2985	sumSSdF2 = np.zeros(10)
2986	sumSSwYo = np.zeros(10)
2987	sumSSwdf2 = np.zeros(10)
2988	SSnobs = np.zeros(10)
2989	nobs = 0
2990	nrej = 0
2991	next = 0
2992	maxH = 0
2993	if calcControls['F**2']:
2994	for i,ref in enumerate(refDict['RefList']):
2995	if ref[6+im] > 0:
2996	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2997	w = 1.0/ref[6+im] # 1/sig(F^2)
2998	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
2999	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
3000	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
3001	ref[3+im] = abs(ref[3+im]) #mark as allowed
3002	Fo = np.sqrt(ref[5+im])
3003	sumFo += Fo
3004	sumFo2 += ref[5+im]
3005	sumdF += abs(Fo-np.sqrt(ref[7+im]))
3006	sumdF2 += abs(ref[5+im]-ref[7+im])
3007	nobs += 1
3008	df[i] = -w*(ref[5+im]-ref[7+im])
3009	sumwYo += (wref[5+im])2 #wFo^2
3010	if im: #accumulate super lattice sums
3011	ind = int(abs(ref[3]))
3012	sumSSFo[ind] += Fo
3013	sumSSFo2[ind] += ref[5+im]
3014	sumSSdF[ind] += abs(Fo-np.sqrt(ref[7+im]))
3015	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
3016	sumSSwYo[ind] += (wref[5+im])2 #wFo^2
3017	sumSSwdf2[ind] += df[i]**2
3018	SSnobs[ind] += 1
3019	maxH = max(maxH,ind)
3020	else:
3021	if ref[3+im]:
3022	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3023	nrej += 1
3024	else: #sp.gp.extinct
3025	next += 1
3026	else:
3027	for i,ref in enumerate(refDict['RefList']):
3028	if ref[5+im] > 0.:
3029	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
3030	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
3031	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
3032	Fo = np.sqrt(ref[5+im])
3033	Fc = np.sqrt(ref[7+im])
3034	w = 2.0*Fo/ref[6+im] # 1/sig(F)?
3035	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
3036	ref[3+im] = abs(ref[3+im]) #mark as allowed
3037	sumFo += Fo
3038	sumFo2 += ref[5+im]
3039	sumdF += abs(Fo-Fc)
3040	sumdF2 += abs(ref[5+im]-ref[7+im])
3041	nobs += 1
3042	df[i] = -w*(Fo-Fc)
3043	sumwYo += (wFo)*2
3044	if im:
3045	ind = int(abs(ref[3]))
3046	sumSSFo[ind] += Fo
3047	sumSSFo2[ind] += ref[5+im]
3048	sumSSdF[ind] += abs(Fo-Fc)
3049	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
3050	sumSSwYo[ind] += (wFo)*2
3051	sumSSwdf2[ind] += df[i]**2
3052	SSnobs[ind] += 1
3053	maxH = max(maxH,ind)
3054	else:
3055	if ref[3+im]:
3056	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3057	nrej += 1
3058	else: #sp.gp.extinct
3059	next += 1
3060	Histogram['Residuals']['Nobs'] = nobs
3061	Histogram['Residuals']['sumwYo'] = sumwYo
3062	SumwYo += sumwYo
3063	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/sumwYo)100.)
3064	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
3065	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3066	Histogram['Residuals'][phfx+'Nref'] = nobs
3067	Histogram['Residuals'][phfx+'Nrej'] = nrej
3068	Histogram['Residuals'][phfx+'Next'] = next
3069	if im:
3070	Histogram['Residuals'][phfx+'SSRf'] = 100.*sumSSdF[:maxH+1]/sumSSFo[:maxH+1]
3071	Histogram['Residuals'][phfx+'SSRf^2'] = 100.*sumSSdF2[:maxH+1]/sumSSFo2[:maxH+1]
3072	Histogram['Residuals'][phfx+'SSNref'] = SSnobs[:maxH+1]
3073	Histogram['Residuals']['SSwR'] = np.sqrt(sumSSwdf2[:maxH+1]/sumSSwYo[:maxH+1])*100.
3074	Nobs += nobs
3075	Nrej += nrej
3076	Next += next
3077	if dlg:
3078	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
3079	M = np.concatenate((M,wtFactor*df))
3080	# end of HKLF processing
3081	Histograms['sumwYo'] = SumwYo
3082	Histograms['Nobs'] = Nobs
3083	Histograms['Nrej'] = Nrej
3084	Histograms['Next'] = Next
3085	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
3086	if dlg:
3087	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3088	if not GoOn:
3089	parmDict['saved values'] = values
3090	dlg.Destroy()
3091	raise G2obj.G2Exception('User abort') #Abort!!
3092	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3093	if len(pVals):
3094	pSum = np.sum(pWtpVals*2)
3095	for name in pWsum:
3096	if pWsum[name]:
3097	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
3098	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
3099	Nobs += len(pVals)
3100	M = np.concatenate((M,np.sqrt(pWt)*pVals))
3101	return M
3102

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

Download in other formats: