Context Navigation

source: trunk/GSASIIstruct.py @ 420

Last change on this file since 420 was 420, checked in by vondreele, 12 years ago
small fixes including move msg = in SetVaryFlags? outside loop arrange so plot shows Yo-Yc (was other way around) put Yo-Yc at bottom in multiplot plots more on seq refinement plots
Property svn:keywords set to `Date Author Revision URL Id`
File size: 119.4 KB

Line
1	#GSASIIstructure - structure computation routines
2	########### SVN repository information ###################
3	# $Date: 2011-11-17 19:35:34 +0000 (Thu, 17 Nov 2011) $
4	# $Author: vondreele $
5	# $Revision: 420 $
6	# $URL: trunk/GSASIIstruct.py $
7	# $Id: GSASIIstruct.py 420 2011-11-17 19:35:34Z vondreele $
8	########### SVN repository information ###################
9	import sys
10	import os
11	import os.path as ospath
12	import numpy as np
13	import numpy.linalg as nl
14	import cPickle
15	import time
16	import math
17	import GSASIIpath
18	import GSASIIElem as G2el
19	import GSASIIlattice as G2lat
20	import GSASIIspc as G2spc
21	import GSASIIpwd as G2pwd
22	import GSASIImapvars as G2mv
23	import scipy.optimize as so
24
25	sind = lambda x: np.sin(x*np.pi/180.)
26	cosd = lambda x: np.cos(x*np.pi/180.)
27	tand = lambda x: np.tan(x*np.pi/180.)
28	asind = lambda x: 180.*np.arcsin(x)/np.pi
29	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
30
31
32	def ShowBanner():
33	print 80''
34	print ' General Structure Analysis System-II Crystal Structure Refinement'
35	print ' by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
36	print ' This product includes software developed by the UChicago Argonne, LLC,'
37	print ' as Operator of Argonne National Laboratory.'
38	print 80'','\n'
39
40	def GetControls(GPXfile):
41	''' Returns dictionary of control items found in GSASII gpx file
42	input:
43	GPXfile = .gpx full file name
44	return:
45	Controls = dictionary of control items
46	'''
47	Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
48	file = open(GPXfile,'rb')
49	while True:
50	try:
51	data = cPickle.load(file)
52	except EOFError:
53	break
54	datum = data[0]
55	if datum[0] == 'Controls':
56	Controls.update(datum[1])
57	file.close()
58	return Controls
59
60	def ShowControls(Controls):
61	print ' Least squares controls:'
62	print ' Derivative type: ',Controls['deriv type']
63	print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
64	print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
65
66	def GetConstraints(GPXfile):
67	constList = []
68	file = open(GPXfile,'rb')
69	while True:
70	try:
71	data = cPickle.load(file)
72	except EOFError:
73	break
74	datum = data[0]
75	if datum[0] == 'Constraints':
76	constDict = datum[1]
77	for item in constDict:
78	constList += constDict[item]
79	file.close()
80	constDict = []
81	constFlag = []
82	fixedList = []
83	for item in constList:
84	if item[-2]:
85	fixedList.append(str(item[-2]))
86	else:
87	fixedList.append('0')
88	if item[-1]:
89	constFlag.append(['VARY'])
90	else:
91	constFlag.append([])
92	itemDict = {}
93	for term in item[:-2]:
94	itemDict[term[1]] = term[0]
95	constDict.append(itemDict)
96	return constDict,constFlag,fixedList
97
98	def GetPhaseNames(GPXfile):
99	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
100	input:
101	GPXfile = gpx full file name
102	return:
103	PhaseNames = list of phase names
104	'''
105	file = open(GPXfile,'rb')
106	PhaseNames = []
107	while True:
108	try:
109	data = cPickle.load(file)
110	except EOFError:
111	break
112	datum = data[0]
113	if 'Phases' == datum[0]:
114	for datus in data[1:]:
115	PhaseNames.append(datus[0])
116	file.close()
117	return PhaseNames
118
119	def GetAllPhaseData(GPXfile,PhaseName):
120	''' Returns the entire dictionary for PhaseName from GSASII gpx file
121	input:
122	GPXfile = gpx full file name
123	PhaseName = phase name
124	return:
125	phase dictionary
126	'''
127	file = open(GPXfile,'rb')
128	General = {}
129	Atoms = []
130	while True:
131	try:
132	data = cPickle.load(file)
133	except EOFError:
134	break
135	datum = data[0]
136	if 'Phases' == datum[0]:
137	for datus in data[1:]:
138	if datus[0] == PhaseName:
139	break
140	file.close()
141	return datus[1]
142
143	def GetHistograms(GPXfile,hNames):
144	""" Returns a dictionary of histograms found in GSASII gpx file
145	input:
146	GPXfile = .gpx full file name
147	hNames = list of histogram names
148	return:
149	Histograms = dictionary of histograms (types = PWDR & HKLF)
150	"""
151	file = open(GPXfile,'rb')
152	Histograms = {}
153	while True:
154	try:
155	data = cPickle.load(file)
156	except EOFError:
157	break
158	datum = data[0]
159	hist = datum[0]
160	if hist in hNames:
161	if 'PWDR' in hist[:4]:
162	PWDRdata = {}
163	PWDRdata['Data'] = datum[1][1] #powder data arrays
164	PWDRdata[data[2][0]] = data[2][1] #Limits
165	PWDRdata[data[3][0]] = data[3][1] #Background
166	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
167	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
168	try:
169	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
170	except IndexError:
171	PWDRdata['Reflection lists'] = {}
172
173	Histograms[hist] = PWDRdata
174	elif 'HKLF' in hist[:4]:
175	HKLFdata = []
176	datum = data[0]
177	HKLFdata = datum[1:][0]
178	Histograms[hist] = HKLFdata
179	file.close()
180	return Histograms
181
182	def GetHistogramNames(GPXfile,hType):
183	""" Returns a list of histogram names found in GSASII gpx file
184	input:
185	GPXfile = .gpx full file name
186	hType = list ['PWDR','HKLF']
187	return:
188	HistogramNames = list of histogram names (types = PWDR & HKLF)
189	"""
190	file = open(GPXfile,'rb')
191	HistogramNames = []
192	while True:
193	try:
194	data = cPickle.load(file)
195	except EOFError:
196	break
197	datum = data[0]
198	if datum[0][:4] in hType:
199	HistogramNames.append(datum[0])
200	file.close()
201	return HistogramNames
202
203	def GetUsedHistogramsAndPhases(GPXfile):
204	''' Returns all histograms that are found in any phase
205	and any phase that uses a histogram
206	input:
207	GPXfile = .gpx full file name
208	return:
209	Histograms = dictionary of histograms as {name:data,...}
210	Phases = dictionary of phases that use histograms
211	'''
212	phaseNames = GetPhaseNames(GPXfile)
213	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
214	allHistograms = GetHistograms(GPXfile,histoList)
215	phaseData = {}
216	for name in phaseNames:
217	phaseData[name] = GetAllPhaseData(GPXfile,name)
218	Histograms = {}
219	Phases = {}
220	for phase in phaseData:
221	Phase = phaseData[phase]
222	if Phase['Histograms']:
223	if phase not in Phases:
224	pId = phaseNames.index(phase)
225	Phase['pId'] = pId
226	Phases[phase] = Phase
227	for hist in Phase['Histograms']:
228	if hist not in Histograms:
229	Histograms[hist] = allHistograms[hist]
230	#future restraint, etc. histograms here
231	hId = histoList.index(hist)
232	Histograms[hist]['hId'] = hId
233	return Histograms,Phases
234
235	def GPXBackup(GPXfile,makeBack=True):
236	import distutils.file_util as dfu
237	GPXpath,GPXname = ospath.split(GPXfile)
238	if GPXpath == '': GPXpath = '.'
239	Name = ospath.splitext(GPXname)[0]
240	files = os.listdir(GPXpath)
241	last = 0
242	for name in files:
243	name = name.split('.')
244	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
245	if makeBack:
246	last = max(last,int(name[1].strip('bak'))+1)
247	else:
248	last = max(last,int(name[1].strip('bak')))
249	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
250	dfu.copy_file(GPXfile,GPXback)
251	return GPXback
252
253	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
254	''' Updates gpxfile from all histograms that are found in any phase
255	and any phase that used a histogram
256	input:
257	GPXfile = .gpx full file name
258	Histograms = dictionary of histograms as {name:data,...}
259	Phases = dictionary of phases that use histograms
260	CovData = dictionary of refined variables, varyList, & covariance matrix
261	makeBack = True if new backup of .gpx file is to be made; else use the last one made
262	'''
263
264	GPXback = GPXBackup(GPXfile,makeBack)
265	print '\n',135*'-'
266	print 'Read from file:',GPXback
267	print 'Save to file :',GPXfile
268	infile = open(GPXback,'rb')
269	outfile = open(GPXfile,'wb')
270	while True:
271	try:
272	data = cPickle.load(infile)
273	except EOFError:
274	break
275	datum = data[0]
276	# print 'read: ',datum[0]
277	if datum[0] == 'Phases':
278	for iphase in range(len(data)):
279	if data[iphase][0] in Phases:
280	phaseName = data[iphase][0]
281	data[iphase][1].update(Phases[phaseName])
282	elif datum[0] == 'Covariance':
283	data[0][1] = CovData
284	try:
285	histogram = Histograms[datum[0]]
286	# print 'found ',datum[0]
287	data[0][1][1] = histogram['Data']
288	for datus in data[1:]:
289	# print ' read: ',datus[0]
290	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
291	datus[1] = histogram[datus[0]]
292	except KeyError:
293	pass
294
295	cPickle.dump(data,outfile,1)
296	infile.close()
297	outfile.close()
298	print 'refinement save successful'
299
300	def SetSeqResult(GPXfile,Histograms,SeqResult):
301	GPXback = GPXBackup(GPXfile)
302	print '\n',135*'-'
303	print 'Read from file:',GPXback
304	print 'Save to file :',GPXfile
305	infile = open(GPXback,'rb')
306	outfile = open(GPXfile,'wb')
307	while True:
308	try:
309	data = cPickle.load(infile)
310	except EOFError:
311	break
312	datum = data[0]
313	if datum[0] == 'Sequental results':
314	data[0][1] = SeqResult
315	try:
316	histogram = Histograms[datum[0]]
317	data[0][1][1] = histogram['Data']
318	for datus in data[1:]:
319	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
320	datus[1] = histogram[datus[0]]
321	except KeyError:
322	pass
323
324	cPickle.dump(data,outfile,1)
325	infile.close()
326	outfile.close()
327	print 'refinement save successful'
328
329
330	def GetPWDRdata(GPXfile,PWDRname):
331	''' Returns powder data from GSASII gpx file
332	input:
333	GPXfile = .gpx full file name
334	PWDRname = powder histogram name as obtained from GetHistogramNames
335	return:
336	PWDRdata = powder data dictionary with:
337	Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
338
339	'''
340	file = open(GPXfile,'rb')
341	PWDRdata = {}
342	while True:
343	try:
344	data = cPickle.load(file)
345	except EOFError:
346	break
347	datum = data[0]
348	if datum[0] == PWDRname:
349	PWDRdata['Data'] = datum[1][1] #powder data arrays
350	PWDRdata[data[2][0]] = data[2][1] #Limits
351	PWDRdata[data[3][0]] = data[3][1] #Background
352	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
353	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
354	try:
355	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
356	except IndexError:
357	PWDRdata['Reflection lists'] = {}
358	file.close()
359	return PWDRdata
360
361	def GetHKLFdata(GPXfile,HKLFname):
362	''' Returns single crystal data from GSASII gpx file
363	input:
364	GPXfile = .gpx full file name
365	HKLFname = single crystal histogram name as obtained from GetHistogramNames
366	return:
367	HKLFdata = single crystal data list of reflections: for each reflection:
368	HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
369	'''
370	file = open(GPXfile,'rb')
371	HKLFdata = []
372	while True:
373	try:
374	data = cPickle.load(file)
375	except EOFError:
376	break
377	datum = data[0]
378	if datum[0] == HKLFname:
379	HKLFdata = datum[1:][0]
380	file.close()
381	return HKLFdata
382
383	def GetFFtable(General):
384	''' returns a dictionary of form factor data for atom types found in General
385	input:
386	General = dictionary of phase info.; includes AtomTypes
387	return:
388	FFtable = dictionary of form factor data; key is atom type
389	'''
390	atomTypes = General['AtomTypes']
391	FFtable = {}
392	for El in atomTypes:
393	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
394	for item in FFs:
395	if item['Symbol'] == El.upper():
396	FFtable[El] = item
397	return FFtable
398
399	def GetBLtable(General):
400	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
401	input:
402	General = dictionary of phase info.; includes AtomTypes & Isotopes
403	return:
404	BLtable = dictionary of scattering length data; key is atom type
405	'''
406	atomTypes = General['AtomTypes']
407	BLtable = {}
408	isotopes = General['Isotopes']
409	isotope = General['Isotope']
410	for El in atomTypes:
411	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
412	return BLtable
413
414	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
415	if SGLaue in ['-1','2/m','mmm']:
416	return #no Pawley symmetry required constraints
417	for i,varyI in enumerate(pawleyVary):
418	refI = int(varyI.split(':')[-1])
419	ih,ik,il = PawleyRef[refI][:3]
420	for varyJ in pawleyVary[0:i]:
421	refJ = int(varyJ.split(':')[-1])
422	jh,jk,jl = PawleyRef[refJ][:3]
423	if SGLaue in ['4/m','4/mmm']:
424	isum = ih2+ik2
425	jsum = jh2+jk2
426	if abs(il) == abs(jl) and isum == jsum:
427	G2mv.StoreEquivalence(varyJ,(varyI,))
428	elif SGLaue in ['3R','3mR']:
429	isum = ih2+ik2+il**2
430	jsum = jh2+jk2jl*2
431	isum2 = ihik+ihil+ik*il
432	jsum2 = jhjk+jhjl+jk*jl
433	if isum == jsum and isum2 == jsum2:
434	G2mv.StoreEquivalence(varyJ,(varyI,))
435	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
436	isum = ih2+ik2+ih*ik
437	jsum = jh2+jk2+jh*jk
438	if abs(il) == abs(jl) and isum == jsum:
439	G2mv.StoreEquivalence(varyJ,(varyI,))
440	elif SGLaue in ['m3','m3m']:
441	isum = ih2+ik2+il**2
442	jsum = jh2+jk2+jl**2
443	if isum == jsum:
444	G2mv.StoreEquivalence(varyJ,(varyI,))
445
446	def cellVary(pfx,SGData):
447	if SGData['SGLaue'] in ['-1',]:
448	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
449	elif SGData['SGLaue'] in ['2/m',]:
450	if SGData['SGUniq'] == 'a':
451	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
452	elif SGData['SGUniq'] == 'b':
453	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
454	else:
455	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
456	elif SGData['SGLaue'] in ['mmm',]:
457	return [pfx+'A0',pfx+'A1',pfx+'A2']
458	elif SGData['SGLaue'] in ['4/m','4/mmm']:
459	return [pfx+'A0',pfx+'A2']
460	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
461	return [pfx+'A0',pfx+'A2']
462	elif SGData['SGLaue'] in ['3R', '3mR']:
463	return [pfx+'A0',pfx+'A3']
464	elif SGData['SGLaue'] in ['m3m','m3']:
465	return [pfx+'A0',]
466
467
468	def GetPhaseData(PhaseData,Print=True):
469
470	def PrintFFtable(FFtable):
471	print '\n X-ray scattering factors:'
472	print ' Symbol fa fb fc'
473	print 99*'-'
474	for Ename in FFtable:
475	ffdata = FFtable[Ename]
476	fa = ffdata['fa']
477	fb = ffdata['fb']
478	print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
479	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
480
481	def PrintBLtable(BLtable):
482	print '\n Neutron scattering factors:'
483	print ' Symbol isotope mass b resonant terms'
484	print 99*'-'
485	for Ename in BLtable:
486	bldata = BLtable[Ename]
487	isotope = bldata[0]
488	mass = bldata[1][0]
489	blen = bldata[1][1]
490	bres = []
491	if len(bldata[1]) > 2:
492	bres = bldata[1][2:]
493	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
494	for item in bres:
495	line += '%10.5g'%(item)
496	print line
497
498	def PrintAtoms(General,Atoms):
499	print '\n Atoms:'
500	line = ' name type refine? x y z '+ \
501	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
502	if General['Type'] == 'magnetic':
503	line += ' Mx My Mz'
504	elif General['Type'] == 'macromolecular':
505	line = ' res no residue chain '+line
506	print line
507	if General['Type'] == 'nuclear':
508	print 135*'-'
509	for i,at in enumerate(Atoms):
510	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
511	'%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
512	if at[9] == 'I':
513	line += '%8.4f'%(at[10])+48*' '
514	else:
515	line += 8*' '
516	for j in range(6):
517	line += '%8.4f'%(at[11+j])
518	print line
519
520	def PrintTexture(textureData):
521	topstr = '\n Spherical harmonics texture: Order:' + \
522	str(textureData['Order'])
523	if textureData['Order']:
524	print topstr+' Refine? '+str(textureData['SH Coeff'][0])
525	else:
526	print topstr
527	return
528	names = ['omega','chi','phi']
529	line = '\n'
530	for name in names:
531	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
532	print line
533	print '\n Texture coefficients:'
534	ptlbls = ' names :'
535	ptstr = ' values:'
536	SHcoeff = textureData['SH Coeff'][1]
537	for item in SHcoeff:
538	ptlbls += '%12s'%(item)
539	ptstr += '%12.4f'%(SHcoeff[item])
540	print ptlbls
541	print ptstr
542
543	if Print: print ' Phases:'
544	phaseVary = []
545	phaseDict = {}
546	phaseConstr = {}
547	pawleyLookup = {}
548	FFtables = {} #scattering factors - xrays
549	BLtables = {} # neutrons
550	Natoms = {}
551	AtMults = {}
552	AtIA = {}
553	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
554	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
555	for name in PhaseData:
556	General = PhaseData[name]['General']
557	pId = PhaseData[name]['pId']
558	pfx = str(pId)+'::'
559	FFtable = GetFFtable(General)
560	BLtable = GetBLtable(General)
561	FFtables.update(FFtable)
562	BLtables.update(BLtable)
563	Atoms = PhaseData[name]['Atoms']
564	try:
565	PawleyRef = PhaseData[name]['Pawley ref']
566	except KeyError:
567	PawleyRef = []
568	SGData = General['SGData']
569	SGtext = G2spc.SGPrint(SGData)
570	cell = General['Cell']
571	A = G2lat.cell2A(cell[1:7])
572	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
573	if cell[0]:
574	phaseVary += cellVary(pfx,SGData)
575	Natoms[pfx] = 0
576	if Atoms:
577	if General['Type'] == 'nuclear':
578	Natoms[pfx] = len(Atoms)
579	for i,at in enumerate(Atoms):
580	phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
581	pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
582	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
583	pfx+'AI/A:'+str(i):at[9],})
584	if at[9] == 'I':
585	phaseDict[pfx+'AUiso:'+str(i)] = at[10]
586	else:
587	phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
588	pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
589	if 'F' in at[2]:
590	phaseVary.append(pfx+'Afrac:'+str(i))
591	if 'X' in at[2]:
592	xId,xCoef = G2spc.GetCSxinel(at[7])
593	delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
594	for j in range(3):
595	if xId[j] > 0:
596	phaseVary.append(delnames[j])
597	for k in range(j):
598	if xId[j] == xId[k]:
599	G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),))
600	if 'U' in at[2]:
601	if at[9] == 'I':
602	phaseVary.append(pfx+'AUiso:'+str(i))
603	else:
604	uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
605	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
606	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
607	for j in range(6):
608	if uId[j] > 0:
609	phaseVary.append(names[j])
610	for k in range(j):
611	if uId[j] == uId[k]:
612	G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
613	# elif General['Type'] == 'magnetic':
614	# elif General['Type'] == 'macromolecular':
615
616
617	if 'SH Texture' in General:
618	textureData = General['SH Texture']
619	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
620	phaseDict[pfx+'SHorder'] = textureData['Order']
621	for name in ['omega','chi','phi']:
622	phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
623	if textureData['Sample '+name][0]:
624	phaseVary.append(pfx+'SH '+name)
625	for name in textureData['SH Coeff'][1]:
626	phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
627	if textureData['SH Coeff'][0]:
628	phaseVary.append(pfx+name)
629
630	if Print:
631	print '\n Phase name: ',General['Name']
632	print 135*'-'
633	PrintFFtable(FFtable)
634	PrintBLtable(BLtable)
635	print ''
636	for line in SGtext: print line
637	PrintAtoms(General,Atoms)
638	print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
639	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
640	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
641	if 'SH Texture' in General:
642	PrintTexture(textureData)
643
644	elif PawleyRef:
645	pawleyVary = []
646	for i,refl in enumerate(PawleyRef):
647	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
648	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
649	if refl[5]:
650	pawleyVary.append(pfx+'PWLref:'+str(i))
651	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
652	phaseVary += pawleyVary
653
654	return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
655
656	def getVCov(varyNames,varyList,covMatrix):
657	vcov = np.zeros((len(varyNames),len(varyNames)))
658	for i1,name1 in enumerate(varyNames):
659	for i2,name2 in enumerate(varyNames):
660	try:
661	vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
662	except ValueError:
663	vcov[i1][i2] = 0.0
664	return vcov
665
666	def cellFill(pfx,SGData,parmDict,sigDict):
667	if SGData['SGLaue'] in ['-1',]:
668	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
669	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
670	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
671	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
672	elif SGData['SGLaue'] in ['2/m',]:
673	if SGData['SGUniq'] == 'a':
674	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
675	parmDict[pfx+'A3'],0,0]
676	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
677	sigDict[pfx+'A3'],0,0]
678	elif SGData['SGUniq'] == 'b':
679	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
680	0,parmDict[pfx+'A4'],0]
681	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
682	0,sigDict[pfx+'A4'],0]
683	else:
684	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
685	0,0,parmDict[pfx+'A5']]
686	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
687	0,0,sigDict[pfx+'A5']]
688	elif SGData['SGLaue'] in ['mmm',]:
689	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
690	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
691	elif SGData['SGLaue'] in ['4/m','4/mmm']:
692	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
693	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
694	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
695	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
696	parmDict[pfx+'A0'],0,0]
697	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
698	elif SGData['SGLaue'] in ['3R', '3mR']:
699	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
700	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
701	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
702	elif SGData['SGLaue'] in ['m3m','m3']:
703	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
704	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
705	return A,sigA
706
707	def getCellEsd(pfx,SGData,A,covData):
708	dpr = 180./np.pi
709	rVsq = G2lat.calc_rVsq(A)
710	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
711	cell = np.array(G2lat.Gmat2cell(g)) #real cell
712	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
713	scos = cosd(cellst[3:6])
714	ssin = sind(cellst[3:6])
715	scot = scos/ssin
716	rcos = cosd(cell[3:6])
717	rsin = sind(cell[3:6])
718	rcot = rcos/rsin
719	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
720	varyList = covData['varyList']
721	covMatrix = covData['covMatrix']
722	vcov = getVCov(RMnames,varyList,covMatrix)
723	Ax = np.array(A)
724	Ax[3:] /= 2.
725	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
726	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
727	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
728	Vol = 1/np.sqrt(rVsq)
729	sigVol = Vol*3np.sqrt(srcvlsq)/2.
730	R123 = Ax[0]Ax[1]Ax[2]
731	dsasdg = np.zeros((3,6))
732	dadg = np.zeros((6,6))
733	for i0 in range(3): #0 1 2
734	i1 = (i0+1)%3 #1 2 0
735	i2 = (i1+1)%3 #2 0 1
736	i3 = 5-i2 #3 5 4
737	i4 = 5-i1 #4 3 5
738	i5 = 5-i0 #5 4 3
739	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
740	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
741	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
742	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
743	denom = np.sqrt(denmsq)
744	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
745	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
746	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
747	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
748	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
749	dadg[i5][i5] = -Ax[i0]/denom
750	for i0 in range(3):
751	i1 = (i0+1)%3
752	i2 = (i1+1)%3
753	i3 = 5-i2
754	for ij in range(6):
755	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
756	if ij == i0:
757	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
758	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
759	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
760	var = np.diag(sigMat)
761	CS = np.where(var>0.,np.sqrt(var),0.)
762	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
763	return cellSig
764
765	def SetPhaseData(parmDict,sigDict,Phases,covData):
766
767	def PrintAtomsAndSig(General,Atoms,atomsSig):
768	print '\n Atoms:'
769	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
770	if General['Type'] == 'magnetic':
771	line += ' Mx My Mz'
772	elif General['Type'] == 'macromolecular':
773	line = ' res no residue chain '+line
774	print line
775	if General['Type'] == 'nuclear':
776	print 135*'-'
777	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
778	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
779	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
780	for i,at in enumerate(Atoms):
781	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
782	valstr = ' values:'
783	sigstr = ' sig :'
784	for ind in [3,4,5,6]:
785	sigind = str(i)+':'+str(ind)
786	valstr += fmt[ind]%(at[ind])
787	if sigind in atomsSig:
788	sigstr += fmt[ind]%(atomsSig[sigind])
789	else:
790	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
791	if at[9] == 'I':
792	valstr += fmt[10]%(at[10])
793	if str(i)+':10' in atomsSig:
794	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
795	else:
796	sigstr += 8*' '
797	else:
798	valstr += 8*' '
799	sigstr += 8*' '
800	for ind in [11,12,13,14,15,16]:
801	sigind = str(i)+':'+str(ind)
802	valstr += fmt[ind]%(at[ind])
803	if sigind in atomsSig:
804	sigstr += fmt[ind]%(atomsSig[sigind])
805	else:
806	sigstr += 8*' '
807	print name
808	print valstr
809	print sigstr
810
811	def PrintSHtextureAndSig(textureData,SHtextureSig):
812	print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
813	names = ['omega','chi','phi']
814	namstr = ' names :'
815	ptstr = ' values:'
816	sigstr = ' esds :'
817	for name in names:
818	namstr += '%12s'%(name)
819	ptstr += '%12.3f'%(textureData['Sample '+name][1])
820	if 'Sample '+name in SHtextureSig:
821	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
822	else:
823	sigstr += 12*' '
824	print namstr
825	print ptstr
826	print sigstr
827	print '\n Texture coefficients:'
828	namstr = ' names :'
829	ptstr = ' values:'
830	sigstr = ' esds :'
831	SHcoeff = textureData['SH Coeff'][1]
832	for name in SHcoeff:
833	namstr += '%12s'%(name)
834	ptstr += '%12.3f'%(SHcoeff[name])
835	if name in SHtextureSig:
836	sigstr += '%12.3f'%(SHtextureSig[name])
837	else:
838	sigstr += 12*' '
839	print namstr
840	print ptstr
841	print sigstr
842
843
844	print '\n Phases:'
845	for phase in Phases:
846	print ' Result for phase: ',phase
847	Phase = Phases[phase]
848	General = Phase['General']
849	SGData = General['SGData']
850	Atoms = Phase['Atoms']
851	cell = General['Cell']
852	pId = Phase['pId']
853	pfx = str(pId)+'::'
854	if cell[0]:
855	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
856	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
857	print ' Reciprocal metric tensor: '
858	ptfmt = "%15.9f"
859	names = ['A11','A22','A33','A12','A13','A23']
860	namstr = ' names :'
861	ptstr = ' values:'
862	sigstr = ' esds :'
863	for name,a,siga in zip(names,A,sigA):
864	namstr += '%15s'%(name)
865	ptstr += ptfmt%(a)
866	if siga:
867	sigstr += ptfmt%(siga)
868	else:
869	sigstr += 15*' '
870	print namstr
871	print ptstr
872	print sigstr
873	cell[1:7] = G2lat.A2cell(A)
874	cell[7] = G2lat.calc_V(A)
875	print ' New unit cell:'
876	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
877	names = ['a','b','c','alpha','beta','gamma','Volume']
878	namstr = ' names :'
879	ptstr = ' values:'
880	sigstr = ' esds :'
881	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
882	namstr += '%12s'%(name)
883	ptstr += fmt%(a)
884	if siga:
885	sigstr += fmt%(siga)
886	else:
887	sigstr += 12*' '
888	print namstr
889	print ptstr
890	print sigstr
891
892	if 'Pawley' in Phase['General']['Type']:
893	pawleyRef = Phase['Pawley ref']
894	for i,refl in enumerate(pawleyRef):
895	key = pfx+'PWLref:'+str(i)
896	refl[6] = abs(parmDict[key]) #suppress negative Fsq
897	if key in sigDict:
898	refl[7] = sigDict[key]
899	else:
900	refl[7] = 0
901	else:
902	atomsSig = {}
903	if General['Type'] == 'nuclear':
904	for i,at in enumerate(Atoms):
905	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
906	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
907	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
908	for ind in [3,4,5,6]:
909	at[ind] = parmDict[names[ind]]
910	if ind in [3,4,5]:
911	name = names[ind].replace('A','dA')
912	else:
913	name = names[ind]
914	if name in sigDict:
915	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
916	if at[9] == 'I':
917	at[10] = parmDict[names[10]]
918	if names[10] in sigDict:
919	atomsSig[str(i)+':10'] = sigDict[names[10]]
920	else:
921	for ind in [11,12,13,14,15,16]:
922	at[ind] = parmDict[names[ind]]
923	if names[ind] in sigDict:
924	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
925	PrintAtomsAndSig(General,Atoms,atomsSig)
926
927	if 'SH Texture' in General:
928	textureData = General['SH Texture']
929	if textureData['Order']:
930	SHtextureSig = {}
931	for name in ['omega','chi','phi']:
932	aname = pfx+'SH '+name
933	textureData['Sample '+name][1] = parmDict[aname]
934	if aname in sigDict:
935	SHtextureSig['Sample '+name] = sigDict[aname]
936	for name in textureData['SH Coeff'][1]:
937	aname = pfx+name
938	textureData['SH Coeff'][1][name] = parmDict[aname]
939	if aname in sigDict:
940	SHtextureSig[name] = sigDict[aname]
941	PrintSHtextureAndSig(textureData,SHtextureSig)
942
943	def GetHistogramPhaseData(Phases,Histograms,Print=True):
944
945	def PrintSize(hapData):
946	if hapData[0] in ['isotropic','uniaxial']:
947	line = '\n Size model : %9s'%(hapData[0])
948	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
949	if hapData[0] == 'uniaxial':
950	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
951	print line
952	else:
953	print '\n Size model : %s'%(hapData[0])
954	Snames = ['S11','S22','S33','S12','S13','S23']
955	ptlbls = ' names :'
956	ptstr = ' values:'
957	varstr = ' refine:'
958	for i,name in enumerate(Snames):
959	ptlbls += '%12s' % (name)
960	ptstr += '%12.6f' % (hapData[4][i])
961	varstr += '%12s' % (str(hapData[5][i]))
962	print ptlbls
963	print ptstr
964	print varstr
965
966	def PrintMuStrain(hapData,SGData):
967	if hapData[0] in ['isotropic','uniaxial']:
968	line = '\n Mustrain model: %9s'%(hapData[0])
969	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
970	if hapData[0] == 'uniaxial':
971	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
972	print line
973	else:
974	print '\n Mustrain model: %s'%(hapData[0])
975	Snames = G2spc.MustrainNames(SGData)
976	ptlbls = ' names :'
977	ptstr = ' values:'
978	varstr = ' refine:'
979	for i,name in enumerate(Snames):
980	ptlbls += '%12s' % (name)
981	ptstr += '%12.6f' % (hapData[4][i])
982	varstr += '%12s' % (str(hapData[5][i]))
983	print ptlbls
984	print ptstr
985	print varstr
986
987	def PrintHStrain(hapData,SGData):
988	print '\n Hydrostatic/elastic strain: '
989	Hsnames = G2spc.HStrainNames(SGData)
990	ptlbls = ' names :'
991	ptstr = ' values:'
992	varstr = ' refine:'
993	for i,name in enumerate(Hsnames):
994	ptlbls += '%12s' % (name)
995	ptstr += '%12.6f' % (hapData[0][i])
996	varstr += '%12s' % (str(hapData[1][i]))
997	print ptlbls
998	print ptstr
999	print varstr
1000
1001	def PrintSHPO(hapData):
1002	print '\n Spherical harmonics preferred orientation: Order:' + \
1003	str(hapData[4])+' Refine? '+str(hapData[2])
1004	ptlbls = ' names :'
1005	ptstr = ' values:'
1006	for item in hapData[5]:
1007	ptlbls += '%12s'%(item)
1008	ptstr += '%12.3f'%(hapData[5][item])
1009	print ptlbls
1010	print ptstr
1011
1012	hapDict = {}
1013	hapVary = []
1014	controlDict = {}
1015	poType = {}
1016	poAxes = {}
1017	spAxes = {}
1018	spType = {}
1019
1020	for phase in Phases:
1021	HistoPhase = Phases[phase]['Histograms']
1022	SGData = Phases[phase]['General']['SGData']
1023	cell = Phases[phase]['General']['Cell'][1:7]
1024	A = G2lat.cell2A(cell)
1025	pId = Phases[phase]['pId']
1026	for histogram in HistoPhase:
1027	try:
1028	Histogram = Histograms[histogram]
1029	except KeyError:
1030	#skip if histogram not included e.g. in a sequential refinement
1031	continue
1032	hapData = HistoPhase[histogram]
1033	hId = Histogram['hId']
1034	limits = Histogram['Limits'][1]
1035	inst = Histogram['Instrument Parameters']
1036	inst = dict(zip(inst[3],inst[1]))
1037	Zero = inst['Zero']
1038	if 'C' in inst['Type']:
1039	try:
1040	wave = inst['Lam']
1041	except KeyError:
1042	wave = inst['Lam1']
1043	dmin = wave/(2.0*sind(limits[1]/2.0))
1044	pfx = str(pId)+':'+str(hId)+':'
1045	for item in ['Scale','Extinction']:
1046	hapDict[pfx+item] = hapData[item][0]
1047	if hapData[item][1]:
1048	hapVary.append(pfx+item)
1049	names = G2spc.HStrainNames(SGData)
1050	for i,name in enumerate(names):
1051	hapDict[pfx+name] = hapData['HStrain'][0][i]
1052	if hapData['HStrain'][1][i]:
1053	hapVary.append(pfx+name)
1054	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1055	if hapData['Pref.Ori.'][0] == 'MD':
1056	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1057	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1058	if hapData['Pref.Ori.'][2]:
1059	hapVary.append(pfx+'MD')
1060	else: #'SH' spherical harmonics
1061	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1062	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1063	for item in hapData['Pref.Ori.'][5]:
1064	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1065	if hapData['Pref.Ori.'][2]:
1066	hapVary.append(pfx+item)
1067	for item in ['Mustrain','Size']:
1068	controlDict[pfx+item+'Type'] = hapData[item][0]
1069	if hapData[item][0] in ['isotropic','uniaxial']:
1070	hapDict[pfx+item+':0'] = hapData[item][1][0]
1071	if hapData[item][2][0]:
1072	hapVary.append(pfx+item+':0')
1073	if hapData[item][0] == 'uniaxial':
1074	controlDict[pfx+item+'Axis'] = hapData[item][3]
1075	hapDict[pfx+item+':1'] = hapData[item][1][1]
1076	if hapData[item][2][1]:
1077	hapVary.append(pfx+item+':1')
1078	else: #generalized for mustrain or ellipsoidal for size
1079	if item == 'Mustrain':
1080	names = G2spc.MustrainNames(SGData)
1081	pwrs = []
1082	for name in names:
1083	h,k,l = name[1:]
1084	pwrs.append([int(h),int(k),int(l)])
1085	controlDict[pfx+'MuPwrs'] = pwrs
1086	for i in range(len(hapData[item][4])):
1087	sfx = ':'+str(i)
1088	hapDict[pfx+item+sfx] = hapData[item][4][i]
1089	if hapData[item][5][i]:
1090	hapVary.append(pfx+item+sfx)
1091
1092	if Print:
1093	print '\n Phase: ',phase,' in histogram: ',histogram
1094	print 135*'-'
1095	print ' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1096	print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1097	if hapData['Pref.Ori.'][0] == 'MD':
1098	Ax = hapData['Pref.Ori.'][3]
1099	print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1100	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1101	else: #'SH' for spherical harmonics
1102	PrintSHPO(hapData['Pref.Ori.'])
1103	PrintSize(hapData['Size'])
1104	PrintMuStrain(hapData['Mustrain'],SGData)
1105	PrintHStrain(hapData['HStrain'],SGData)
1106	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1107	refList = []
1108	for h,k,l,d in HKLd:
1109	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1110	if ext:
1111	continue
1112	if 'C' in inst['Type']:
1113	pos = 2.0asind(wave/(2.0d))
1114	if limits[0] < pos < limits[1]:
1115	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
1116	else:
1117	raise ValueError
1118	Histogram['Reflection Lists'][phase] = refList
1119	return hapVary,hapDict,controlDict
1120
1121	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1122
1123	def PrintSizeAndSig(hapData,sizeSig):
1124	line = '\n Size model: %9s'%(hapData[0])
1125	if hapData[0] in ['isotropic','uniaxial']:
1126	line += ' equatorial:%12.3f'%(hapData[1][0])
1127	if sizeSig[0][0]:
1128	line += ', sig: %8.3f'%(sizeSig[0][0])
1129	if hapData[0] == 'uniaxial':
1130	line += ' axial:%12.3f'%(hapData[1][1])
1131	if sizeSig[0][1]:
1132	line += ', sig: %8.3f'%(sizeSig[0][1])
1133	print line
1134	else:
1135	print line
1136	Snames = ['S11','S22','S33','S12','S13','S23']
1137	ptlbls = ' name :'
1138	ptstr = ' value :'
1139	sigstr = ' sig :'
1140	for i,name in enumerate(Snames):
1141	ptlbls += '%12s' % (name)
1142	ptstr += '%12.6f' % (hapData[4][i])
1143	if sizeSig[1][i]:
1144	sigstr += '%12.6f' % (sizeSig[1][i])
1145	else:
1146	sigstr += 12*' '
1147	print ptlbls
1148	print ptstr
1149	print sigstr
1150
1151	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1152	line = '\n Mustrain model: %9s'%(hapData[0])
1153	if hapData[0] in ['isotropic','uniaxial']:
1154	line += ' equatorial:%12.1f'%(hapData[1][0])
1155	if mustrainSig[0][0]:
1156	line += ', sig: %8.1f'%(mustrainSig[0][0])
1157	if hapData[0] == 'uniaxial':
1158	line += ' axial:%12.1f'%(hapData[1][1])
1159	if mustrainSig[0][1]:
1160	line += ', sig: %8.1f'%(mustrainSig[0][1])
1161	print line
1162	else:
1163	print line
1164	Snames = G2spc.MustrainNames(SGData)
1165	ptlbls = ' name :'
1166	ptstr = ' value :'
1167	sigstr = ' sig :'
1168	for i,name in enumerate(Snames):
1169	ptlbls += '%12s' % (name)
1170	ptstr += '%12.6f' % (hapData[4][i])
1171	if mustrainSig[1][i]:
1172	sigstr += '%12.6f' % (mustrainSig[1][i])
1173	else:
1174	sigstr += 12*' '
1175	print ptlbls
1176	print ptstr
1177	print sigstr
1178
1179	def PrintHStrainAndSig(hapData,strainSig,SGData):
1180	print '\n Hydrostatic/elastic strain: '
1181	Hsnames = G2spc.HStrainNames(SGData)
1182	ptlbls = ' name :'
1183	ptstr = ' value :'
1184	sigstr = ' sig :'
1185	for i,name in enumerate(Hsnames):
1186	ptlbls += '%12s' % (name)
1187	ptstr += '%12.6g' % (hapData[0][i])
1188	if name in strainSig:
1189	sigstr += '%12.6g' % (strainSig[name])
1190	else:
1191	sigstr += 12*' '
1192	print ptlbls
1193	print ptstr
1194	print sigstr
1195
1196	def PrintSHPOAndSig(hapData,POsig):
1197	print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1198	ptlbls = ' names :'
1199	ptstr = ' values:'
1200	sigstr = ' sig :'
1201	for item in hapData[5]:
1202	ptlbls += '%12s'%(item)
1203	ptstr += '%12.3f'%(hapData[5][item])
1204	if item in POsig:
1205	sigstr += '%12.3f'%(POsig[item])
1206	else:
1207	sigstr += 12*' '
1208	print ptlbls
1209	print ptstr
1210	print sigstr
1211
1212	for phase in Phases:
1213	HistoPhase = Phases[phase]['Histograms']
1214	SGData = Phases[phase]['General']['SGData']
1215	pId = Phases[phase]['pId']
1216	for histogram in HistoPhase:
1217	try:
1218	Histogram = Histograms[histogram]
1219	except KeyError:
1220	#skip if histogram not included e.g. in a sequential refinement
1221	continue
1222	print '\n Phase: ',phase,' in histogram: ',histogram
1223	print 130*'-'
1224	hapData = HistoPhase[histogram]
1225	hId = Histogram['hId']
1226	pfx = str(pId)+':'+str(hId)+':'
1227	print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1228	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1229
1230	PhFrExtPOSig = {}
1231	for item in ['Scale','Extinction']:
1232	hapData[item][0] = parmDict[pfx+item]
1233	if pfx+item in sigDict:
1234	PhFrExtPOSig[item] = sigDict[pfx+item]
1235	if hapData['Pref.Ori.'][0] == 'MD':
1236	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1237	if pfx+'MD' in sigDict:
1238	PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1239	else: #'SH' spherical harmonics
1240	for item in hapData['Pref.Ori.'][5]:
1241	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1242	if pfx+item in sigDict:
1243	PhFrExtPOSig[item] = sigDict[pfx+item]
1244	if Print:
1245	if 'Scale' in PhFrExtPOSig:
1246	print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1247	if 'Extinction' in PhFrExtPOSig:
1248	print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1249	if hapData['Pref.Ori.'][0] == 'MD':
1250	if 'MD' in PhFrExtPOSig:
1251	print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1252	else:
1253	PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1254	SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1255	'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
1256	'HStrain':{}}
1257	for item in ['Mustrain','Size']:
1258	if hapData[item][0] in ['isotropic','uniaxial']:
1259	hapData[item][1][0] = parmDict[pfx+item+':0']
1260	if item == 'Size':
1261	hapData[item][1][0] = min(10.,max(0.01,hapData[item][1][0]))
1262	if pfx+item+':0' in sigDict:
1263	SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
1264	if hapData[item][0] == 'uniaxial':
1265	hapData[item][1][1] = parmDict[pfx+item+':1']
1266	if item == 'Size':
1267	hapData[item][1][1] = min(10.,max(0.01,hapData[item][1][1]))
1268	if pfx+item+':1' in sigDict:
1269	SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
1270	else: #generalized for mustrain or ellipsoidal for size
1271	for i in range(len(hapData[item][4])):
1272	sfx = ':'+str(i)
1273	hapData[item][4][i] = parmDict[pfx+item+sfx]
1274	if pfx+item+sfx in sigDict:
1275	SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1276	names = G2spc.HStrainNames(SGData)
1277	for i,name in enumerate(names):
1278	hapData['HStrain'][0][i] = parmDict[pfx+name]
1279	if pfx+name in sigDict:
1280	SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1281	if Print:
1282	PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1283	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1284	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1285
1286	def GetHistogramData(Histograms,Print=True):
1287
1288	def GetBackgroundParms(hId,Background):
1289	bakType,bakFlag = Background[:2]
1290	backVals = Background[3:]
1291	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1292	if bakFlag: #returns backNames as varyList = backNames
1293	return bakType,dict(zip(backNames,backVals)),backNames
1294	else: #no background varied; varyList = []
1295	return bakType,dict(zip(backNames,backVals)),[]
1296
1297	def GetInstParms(hId,Inst):
1298	insVals,insFlags,insNames = Inst[1:4]
1299	dataType = insVals[0]
1300	instDict = {}
1301	insVary = []
1302	pfx = ':'+str(hId)+':'
1303	for i,flag in enumerate(insFlags):
1304	insName = pfx+insNames[i]
1305	instDict[insName] = insVals[i]
1306	if flag:
1307	insVary.append(insName)
1308	instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
1309	instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
1310	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1311	return dataType,instDict,insVary
1312
1313	def GetSampleParms(hId,Sample):
1314	sampVary = []
1315	hfx = ':'+str(hId)+':'
1316	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1317	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1318	Type = Sample['Type']
1319	if 'Bragg' in Type: #Bragg-Brentano
1320	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
1321	sampDict[hfx+item] = Sample[item][0]
1322	if Sample[item][1]:
1323	sampVary.append(hfx+item)
1324	elif 'Debye' in Type: #Debye-Scherrer
1325	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1326	sampDict[hfx+item] = Sample[item][0]
1327	if Sample[item][1]:
1328	sampVary.append(hfx+item)
1329	return Type,sampDict,sampVary
1330
1331	def PrintBackground(Background):
1332	print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
1333	line = ' Coefficients: '
1334	for i,back in enumerate(Background[3:]):
1335	line += '%10.3f'%(back)
1336	if i and not i%10:
1337	line += '\n'+15*' '
1338	print line
1339
1340	def PrintInstParms(Inst):
1341	print '\n Instrument Parameters:'
1342	ptlbls = ' name :'
1343	ptstr = ' value :'
1344	varstr = ' refine:'
1345	instNames = Inst[3][1:]
1346	for i,name in enumerate(instNames):
1347	ptlbls += '%12s' % (name)
1348	ptstr += '%12.6f' % (Inst[1][i+1])
1349	if name in ['Lam1','Lam2','Azimuth']:
1350	varstr += 12*' '
1351	else:
1352	varstr += '%12s' % (str(bool(Inst[2][i+1])))
1353	print ptlbls
1354	print ptstr
1355	print varstr
1356
1357	def PrintSampleParms(Sample):
1358	print '\n Sample Parameters:'
1359	print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1360	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
1361	ptlbls = ' name :'
1362	ptstr = ' value :'
1363	varstr = ' refine:'
1364	if 'Bragg' in Sample['Type']:
1365	for item in ['Scale','Shift','Transparency']:
1366	ptlbls += '%14s'%(item)
1367	ptstr += '%14.4f'%(Sample[item][0])
1368	varstr += '%14s'%(str(bool(Sample[item][1])))
1369
1370	elif 'Debye' in Type: #Debye-Scherrer
1371	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1372	ptlbls += '%14s'%(item)
1373	ptstr += '%14.4f'%(Sample[item][0])
1374	varstr += '%14s'%(str(bool(Sample[item][1])))
1375
1376	print ptlbls
1377	print ptstr
1378	print varstr
1379
1380
1381	histDict = {}
1382	histVary = []
1383	controlDict = {}
1384	for histogram in Histograms:
1385	Histogram = Histograms[histogram]
1386	hId = Histogram['hId']
1387	pfx = ':'+str(hId)+':'
1388	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1389
1390	Background = Histogram['Background'][0]
1391	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1392	controlDict[pfx+'bakType'] = Type
1393	histDict.update(bakDict)
1394	histVary += bakVary
1395
1396	Inst = Histogram['Instrument Parameters']
1397	Type,instDict,insVary = GetInstParms(hId,Inst)
1398	controlDict[pfx+'histType'] = Type
1399	if pfx+'Lam1' in instDict:
1400	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1401	else:
1402	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
1403	histDict.update(instDict)
1404	histVary += insVary
1405
1406	Sample = Histogram['Sample Parameters']
1407	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1408	controlDict[pfx+'instType'] = Type
1409	histDict.update(sampDict)
1410	histVary += sampVary
1411
1412	if Print:
1413	print '\n Histogram: ',histogram,' histogram Id: ',hId
1414	print 135*'-'
1415	Units = {'C':' deg','T':' msec'}
1416	units = Units[controlDict[pfx+'histType'][2]]
1417	Limits = controlDict[pfx+'Limits']
1418	print ' Instrument type: ',Sample['Type']
1419	print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
1420	PrintSampleParms(Sample)
1421	PrintInstParms(Inst)
1422	PrintBackground(Background)
1423
1424	return histVary,histDict,controlDict
1425
1426	def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1427
1428	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1429	lenBack = len(Background[3:])
1430	backSig = [0 for i in range(lenBack)]
1431	for i in range(lenBack):
1432	Background[3+i] = parmDict[pfx+'Back:'+str(i)]
1433	if pfx+'Back:'+str(i) in sigDict:
1434	backSig[i] = sigDict[pfx+'Back:'+str(i)]
1435	return backSig
1436
1437	def SetInstParms(pfx,Inst,parmDict,sigDict):
1438	insVals,insFlags,insNames = Inst[1:4]
1439	instSig = [0 for i in range(len(insVals))]
1440	for i,flag in enumerate(insFlags):
1441	insName = pfx+insNames[i]
1442	insVals[i] = parmDict[insName]
1443	if insName in sigDict:
1444	instSig[i] = sigDict[insName]
1445	return instSig
1446
1447	def SetSampleParms(pfx,Sample,parmDict,sigDict):
1448	if 'Bragg' in Sample['Type']: #Bragg-Brentano
1449	sampSig = [0 for i in range(3)]
1450	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
1451	Sample[item][0] = parmDict[pfx+item]
1452	if pfx+item in sigDict:
1453	sampSig[i] = sigDict[pfx+item]
1454	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1455	sampSig = [0 for i in range(4)]
1456	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1457	Sample[item][0] = parmDict[pfx+item]
1458	if pfx+item in sigDict:
1459	sampSig[i] = sigDict[pfx+item]
1460	return sampSig
1461
1462	def PrintBackgroundSig(Background,backSig):
1463	print '\n Background function: ',Background[0]
1464	valstr = ' value : '
1465	sigstr = ' sig : '
1466	for i,back in enumerate(Background[3:]):
1467	valstr += '%10.4f'%(back)
1468	if Background[1]:
1469	sigstr += '%10.4f'%(backSig[i])
1470	else:
1471	sigstr += 10*' '
1472	print valstr
1473	print sigstr
1474
1475	def PrintInstParmsSig(Inst,instSig):
1476	print '\n Instrument Parameters:'
1477	ptlbls = ' names :'
1478	ptstr = ' value :'
1479	sigstr = ' sig :'
1480	instNames = Inst[3][1:]
1481	for i,name in enumerate(instNames):
1482	ptlbls += '%12s' % (name)
1483	ptstr += '%12.6f' % (Inst[1][i+1])
1484	if instSig[i+1]:
1485	sigstr += '%12.6f' % (instSig[i+1])
1486	else:
1487	sigstr += 12*' '
1488	print ptlbls
1489	print ptstr
1490	print sigstr
1491
1492	def PrintSampleParmsSig(Sample,sampleSig):
1493	print '\n Sample Parameters:'
1494	ptlbls = ' names :'
1495	ptstr = ' values:'
1496	sigstr = ' sig :'
1497	if 'Bragg' in Sample['Type']:
1498	for i,item in enumerate(['Scale','Shift','Transparency']):
1499	ptlbls += '%14s'%(item)
1500	ptstr += '%14.4f'%(Sample[item][0])
1501	if sampleSig[i]:
1502	sigstr += '%14.4f'%(sampleSig[i])
1503	else:
1504	sigstr += 14*' '
1505
1506	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1507	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1508	ptlbls += '%14s'%(item)
1509	ptstr += '%14.4f'%(Sample[item][0])
1510	if sampleSig[i]:
1511	sigstr += '%14.4f'%(sampleSig[i])
1512	else:
1513	sigstr += 14*' '
1514
1515	print ptlbls
1516	print ptstr
1517	print sigstr
1518
1519	for histogram in Histograms:
1520	if 'PWDR' in histogram:
1521	Histogram = Histograms[histogram]
1522	hId = Histogram['hId']
1523	pfx = ':'+str(hId)+':'
1524	Background = Histogram['Background'][0]
1525	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1526
1527	Inst = Histogram['Instrument Parameters']
1528	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1529
1530	Sample = Histogram['Sample Parameters']
1531	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1532
1533	print '\n Histogram: ',histogram,' histogram Id: ',hId
1534	print 135*'-'
1535	print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1536	if Print:
1537	print ' Instrument type: ',Sample['Type']
1538	PrintSampleParmsSig(Sample,sampSig)
1539	PrintInstParmsSig(Inst,instSig)
1540	PrintBackgroundSig(Background,backSig)
1541
1542	def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
1543	Natoms = calcControls['Natoms'][pfx]
1544	Tdata = Natoms*[' ',]
1545	Mdata = np.zeros(Natoms)
1546	IAdata = Natoms*[' ',]
1547	Fdata = np.zeros(Natoms)
1548	FFdata = []
1549	BLdata = []
1550	Xdata = np.zeros((3,Natoms))
1551	dXdata = np.zeros((3,Natoms))
1552	Uisodata = np.zeros(Natoms)
1553	Uijdata = np.zeros((6,Natoms))
1554	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1555	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1556	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1557	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1558	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1559	for iatm in range(Natoms):
1560	for key in keys:
1561	parm = pfx+key+str(iatm)
1562	if parm in parmDict:
1563	keys[key][iatm] = parmDict[parm]
1564	FFdata.append(FFtables[Tdata[iatm]])
1565	BLdata.append(BLtables[Tdata[iatm]][1])
1566	return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1567
1568	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1569	''' Compute structure factors for all h,k,l for phase
1570	input:
1571	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1572	G: reciprocal metric tensor
1573	pfx: phase id string
1574	SGData: space group info. dictionary output from SpcGroup
1575	calcControls:
1576	ParmDict:
1577	puts result F^2 in each ref[8] in refList
1578	'''
1579	twopi = 2.0*np.pi
1580	twopisq = 2.0np.pi*2
1581	ast = np.sqrt(np.diag(G))
1582	Mast = twopisq*np.multiply.outer(ast,ast)
1583	FFtables = calcControls['FFtables']
1584	BLtables = calcControls['BLtables']
1585	FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1586	if 'N' in parmDict[hfx+'Type']:
1587	FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
1588	else:
1589	FP = np.array([El[hfx+'FP'] for El in FFdata])
1590	FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1591	maxPos = len(SGData['SGOps'])
1592	Uij = np.array(G2lat.U6toUij(Uijdata))
1593	bij = Mast*Uij.T
1594	for refl in refList:
1595	fbs = np.array([0,0])
1596	H = refl[:3]
1597	SQ = 1./(2.refl[4])*2
1598	if 'N' in parmDict[hfx+'Type']:
1599	FF = np.array([El[1] for El in BLdata])
1600	else: #'X'
1601	FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1602	SQfactor = 4.0SQtwopisq
1603	Uniq = refl[11]
1604	phi = refl[12]
1605	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1606	sinp = np.sin(phase)
1607	cosp = np.cos(phase)
1608	occ = Mdata*Fdata/len(Uniq)
1609	biso = -SQfactor*Uisodata
1610	Tiso = np.where(biso<1.,np.exp(biso),1.0)
1611	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1612	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1613	Tcorr = Tiso*Tuij
1614	fa = np.array([(FF+FP)occcospTcorr,-FPPoccsinpTcorr])
1615	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
1616	if not SGData['SGInv']:
1617	fb = np.array([(FF+FP)occsinpTcorr,FPPocccospTcorr])
1618	fbs = np.sum(np.sum(fb,axis=1),axis=1)
1619	fasq = fas**2
1620	fbsq = fbs**2 #imaginary
1621	refl[9] = np.sum(fasq)+np.sum(fbsq)
1622	refl[10] = atan2d(fbs[0],fas[0])
1623	return refList
1624
1625	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1626	twopi = 2.0*np.pi
1627	twopisq = 2.0np.pi*2
1628	ast = np.sqrt(np.diag(G))
1629	Mast = twopisq*np.multiply.outer(ast,ast)
1630	FFtables = calcControls['FFtables']
1631	BLtables = calcControls['BLtables']
1632	FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1633	if 'N' in parmDict[hfx+'Type']:
1634	FP = 0.
1635	FPP = 0.
1636	else:
1637	FP = np.array([El[hfx+'FP'] for El in FFdata])
1638	FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1639	maxPos = len(SGData['SGOps'])
1640	Uij = np.array(G2lat.U6toUij(Uijdata))
1641	bij = Mast*Uij.T
1642	dFdvDict = {}
1643	dFdfr = np.zeros((len(refList),len(Mdata)))
1644	dFdx = np.zeros((len(refList),len(Mdata),3))
1645	dFdui = np.zeros((len(refList),len(Mdata)))
1646	dFdua = np.zeros((len(refList),len(Mdata),6))
1647	for iref,refl in enumerate(refList):
1648	H = np.array(refl[:3])
1649	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
1650	if 'N' in parmDict[hfx+'Type']:
1651	FF = np.array([El[1] for El in BLdata])
1652	else: #'X'
1653	FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1654	SQfactor = 8.0SQnp.pi**2
1655	Uniq = refl[11]
1656	phi = refl[12]
1657	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1658	sinp = np.sin(phase)
1659	cosp = np.cos(phase)
1660	occ = Mdata*Fdata/len(Uniq)
1661	biso = -SQfactor*Uisodata
1662	Tiso = np.where(biso<1.,np.exp(biso),1.0)
1663	# HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1664	HbH = -np.inner(H,np.inner(bij,H))
1665	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1666	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1667	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1668	Tcorr = Tiso*Tuij
1669	fot = (FF+FP)occTcorr
1670	fotp = FPPoccTcorr
1671	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
1672	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
1673
1674	fas = np.sum(np.sum(fa,axis=1),axis=1)
1675	fbs = np.sum(np.sum(fb,axis=1),axis=1)
1676	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
1677	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
1678	#sum below is over Uniq
1679	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1680	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
1681	dfadui = np.sum(-SQfactor*fa,axis=2)
1682	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1683	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
1684	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
1685	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
1686	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
1687	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
1688	if not SGData['SGInv']:
1689	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
1690	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
1691	dfbdui = np.sum(-SQfactor*fb,axis=2)
1692	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1693	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
1694	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
1695	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
1696	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
1697	#loop over atoms - each dict entry is list of derivatives for all the reflections
1698	for i in range(len(Mdata)):
1699	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1700	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1701	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1702	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1703	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1704	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1705	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1706	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1707	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1708	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1709	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1710	return dFdvDict
1711
1712	def Dict2Values(parmdict, varylist):
1713	'''Use before call to leastsq to setup list of values for the parameters
1714	in parmdict, as selected by key in varylist'''
1715	return [parmdict[key] for key in varylist]
1716
1717	def Values2Dict(parmdict, varylist, values):
1718	''' Use after call to leastsq to update the parameter dictionary with
1719	values corresponding to keys in varylist'''
1720	parmdict.update(zip(varylist,values))
1721
1722	def GetNewCellParms(parmDict,varyList):
1723	newCellDict = {}
1724	Dchoices = ['D11','D22','D33','D12','D13','D23']
1725	Achoices = ['A'+str(i) for i in range(6)]
1726
1727
1728	return newCellDict
1729
1730	def ApplyXYZshifts(parmDict,varyList):
1731	''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1732	input:
1733	parmDict - parameter dictionary
1734	varyList - list of variables
1735	returns:
1736	newAtomDict - dictionary of new atomic coordinate names & values;
1737	key is parameter shift name
1738	'''
1739	newAtomDict = {}
1740	for item in parmDict:
1741	if 'dA' in item:
1742	parm = ''.join(item.split('d'))
1743	parmDict[parm] += parmDict[item]
1744	newAtomDict[item] = [parm,parmDict[parm]]
1745	return newAtomDict
1746
1747	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1748	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1749	odfCor = 1.0
1750	H = refl[:3]
1751	cell = G2lat.Gmat2cell(g)
1752	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1753	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1754	phi,beta = G2lat.CrsAng(H,cell,SGData)
1755	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1756	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1757	for item in SHnames:
1758	L,M,N = eval(item.strip('C'))
1759	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1760	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1761	Lnorm = G2lat.Lnorm(L)
1762	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1763	return odfCor
1764
1765	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1766	FORPI = 12.5663706143592
1767	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1768	odfCor = 1.0
1769	dFdODF = {}
1770	dFdSA = [0,0,0]
1771	H = refl[:3]
1772	cell = G2lat.Gmat2cell(g)
1773	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1774	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1775	phi,beta = G2lat.CrsAng(H,cell,SGData)
1776	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1777	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1778	for item in SHnames:
1779	L,M,N = eval(item.strip('C'))
1780	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1781	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1782	Lnorm = G2lat.Lnorm(L)
1783	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1784	dFdODF[pfx+item] = LnormKclKsl
1785	for i in range(3):
1786	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
1787	return odfCor,dFdODF,dFdSA
1788
1789	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1790	odfCor = 1.0
1791	H = refl[:3]
1792	cell = G2lat.Gmat2cell(g)
1793	Sangl = [0.,0.,0.]
1794	if 'Bragg' in calcControls[hfx+'instType']:
1795	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1796	IFCoup = True
1797	else:
1798	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1799	IFCoup = False
1800	phi,beta = G2lat.CrsAng(H,cell,SGData)
1801	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1802	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1803	for item in SHnames:
1804	L,N = eval(item.strip('C'))
1805	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
1806	odfCor += parmDict[phfx+item]LnormKcsl
1807	return odfCor
1808
1809	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1810	FORPI = 12.5663706143592
1811	odfCor = 1.0
1812	dFdODF = {}
1813	H = refl[:3]
1814	cell = G2lat.Gmat2cell(g)
1815	Sangl = [0.,0.,0.]
1816	if 'Bragg' in calcControls[hfx+'instType']:
1817	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1818	IFCoup = True
1819	else:
1820	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1821	IFCoup = False
1822	phi,beta = G2lat.CrsAng(H,cell,SGData)
1823	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1824	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1825	for item in SHnames:
1826	L,N = eval(item.strip('C'))
1827	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
1828	odfCor += parmDict[phfx+item]LnormKcsl
1829	dFdODF[phfx+item] = Kcsl*Lnorm
1830	return odfCor,dFdODF
1831
1832	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1833	if calcControls[phfx+'poType'] == 'MD':
1834	MD = parmDict[phfx+'MD']
1835	MDAxis = calcControls[phfx+'MDAxis']
1836	sumMD = 0
1837	for H in refl[11]:
1838	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1839	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1840	sumMD += A**3
1841	POcorr = sumMD/len(refl[11])
1842	else: #spherical harmonics
1843	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1844	return POcorr
1845
1846	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1847	POderv = {}
1848	if calcControls[phfx+'poType'] == 'MD':
1849	MD = parmDict[phfx+'MD']
1850	MDAxis = calcControls[phfx+'MDAxis']
1851	sumMD = 0
1852	sumdMD = 0
1853	for H in refl[11]:
1854	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1855	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1856	sumMD += A**3
1857	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1858	POcorr = sumMD/len(refl[11])
1859	POderv[phfx+'MD'] = sumdMD/len(refl[11])
1860	else: #spherical harmonics
1861	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1862	return POcorr,POderv
1863
1864	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1865	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
1866	if 'X' in parmDict[hfx+'Type']:
1867	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1868	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1869	if pfx+'SHorder' in parmDict:
1870	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1871	refl[13] = Icorr
1872
1873	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1874	dIdsh = 1./parmDict[hfx+'Scale']
1875	dIdsp = 1./parmDict[phfx+'Scale']
1876	if 'X' in parmDict[hfx+'Type']:
1877	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1878	dIdPola /= pola
1879	else: #'N'
1880	dIdPola = 0.0
1881	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1882	for iPO in dIdPO:
1883	dIdPO[iPO] /= POcorr
1884	dFdODF = {}
1885	dFdSA = [0,0,0]
1886	if pfx+'SHorder' in parmDict:
1887	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1888	for iSH in dFdODF:
1889	dFdODF[iSH] /= odfCor
1890	for i in range(3):
1891	dFdSA[i] /= odfCor
1892	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
1893
1894	def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1895	costh = cosd(refl[5]/2.)
1896	#crystallite size
1897	if calcControls[phfx+'SizeType'] == 'isotropic':
1898	gam = 1.8wave/(np.piparmDict[phfx+'Size:0']*costh)
1899	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1900	H = np.array(refl[:3])
1901	P = np.array(calcControls[phfx+'SizeAxis'])
1902	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1903	gam = (1.8wave/np.pi)/(parmDict[phfx+'Size:0']parmDict[phfx+'Size:1']*costh)
1904	gam = np.sqrt((cosPparmDict[phfx+'Size:1'])*2+(sinPparmDict[phfx+'Size:0'])**2)
1905	else: #ellipsoidal crystallites - wrong not make sense
1906	H = np.array(refl[:3])
1907	gam += 1.8wave/(np.picosth*np.inner(H,np.inner(sizeEllipse,H)))
1908	#microstrain
1909	if calcControls[phfx+'MustrainType'] == 'isotropic':
1910	gam += 0.018parmDict[phfx+'Mustrain:0']tand(refl[5]/2.)/np.pi
1911	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1912	H = np.array(refl[:3])
1913	P = np.array(calcControls[phfx+'MustrainAxis'])
1914	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1915	Si = parmDict[phfx+'Mustrain:0']
1916	Sa = parmDict[phfx+'Mustrain:1']
1917	gam += 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1918	else: #generalized - P.W. Stephens model
1919	pwrs = calcControls[phfx+'MuPwrs']
1920	sum = 0
1921	for i,pwr in enumerate(pwrs):
1922	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1923	gam += 0.018refl[4]2tand(refl[5]/2.)*sum
1924	return gam
1925
1926	def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1927	gamDict = {}
1928	costh = cosd(refl[5]/2.)
1929	tanth = tand(refl[5]/2.)
1930	#crystallite size derivatives
1931	if calcControls[phfx+'SizeType'] == 'isotropic':
1932	gamDict[phfx+'Size:0'] = -1.80wave/(np.picosth)
1933	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1934	H = np.array(refl[:3])
1935	P = np.array(calcControls[phfx+'SizeAxis'])
1936	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1937	Si = parmDict[phfx+'Size:0']
1938	Sa = parmDict[phfx+'Size:1']
1939	gami = (1.80wave/np.pi)/(SiSa)
1940	sqtrm = np.sqrt((cosPSa)2+(sinPSi)**2)
1941	gam = gami*sqtrm/costh
1942	gamDict[phfx+'Size:0'] = gamiSisinP*2/(sqtrmcosth)-gam/Si
1943	gamDict[phfx+'Size:1'] = gamiSacosP*2/(sqtrmcosth)-gam/Sa
1944	else: #ellipsoidal crystallites - do numerically? - not right not make sense
1945	H = np.array(refl[:3])
1946	gam = 1.8wave/(np.picosth*np.inner(H,np.inner(sizeEllipse,H)))
1947	#microstrain derivatives
1948	if calcControls[phfx+'MustrainType'] == 'isotropic':
1949	gamDict[phfx+'Mustrain:0'] = 0.018*tanth/np.pi
1950	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1951	H = np.array(refl[:3])
1952	P = np.array(calcControls[phfx+'MustrainAxis'])
1953	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1954	Si = parmDict[phfx+'Mustrain:0']
1955	Sa = parmDict[phfx+'Mustrain:1']
1956	gami = 0.018SiSa*tanth/np.pi
1957	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1958	gam = gami/sqtrm
1959	gamDict[phfx+'Mustrain:0'] = gam/Si-gamiSicosP2/sqtrm3
1960	gamDict[phfx+'Mustrain:1'] = gam/Sa-gamiSasinP2/sqtrm3
1961	else: #generalized - P.W. Stephens model
1962	pwrs = calcControls[phfx+'MuPwrs']
1963	const = 0.018refl[4]2tanth
1964	for i,pwr in enumerate(pwrs):
1965	gamDict[phfx+'Mustrain:'+str(i)] = constrefl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1966	return gamDict
1967
1968	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1969	h,k,l = refl[:3]
1970	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1971	d = np.sqrt(dsq)
1972	refl[4] = d
1973	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1974	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1975	if 'Bragg' in calcControls[hfx+'instType']:
1976	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1977	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1978	else: #Debye-Scherrer - simple but maybe not right
1979	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1980	return pos
1981
1982	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1983	dpr = 180./np.pi
1984	h,k,l = refl[:3]
1985	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1986	dst = np.sqrt(dstsq)
1987	pos = refl[5]
1988	const = dpr/np.sqrt(1.0-wave*dst/4.0)
1989	dpdw = const*dst
1990	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1991	dpdA = constwave/(2.0*dst)
1992	dpdZ = 1.0
1993	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1994	if 'Bragg' in calcControls[hfx+'instType']:
1995	dpdSh = -4.constcosd(pos/2.0)
1996	dpdTr = -constsind(pos)100.0
1997	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1998	else: #Debye-Scherrer - simple but maybe not right
1999	dpdXd = -const*cosd(pos)
2000	dpdYd = -const*sind(pos)
2001	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2002
2003	def GetHStrainShift(refl,SGData,phfx,parmDict):
2004	laue = SGData['SGLaue']
2005	uniq = SGData['SGUniq']
2006	h,k,l = refl[:3]
2007	if laue in ['m3','m3m']:
2008	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2009	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2010	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2011	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2012	elif laue in ['3R','3mR']:
2013	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2014	elif laue in ['4/m','4/mmm']:
2015	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2016	elif laue in ['mmm']:
2017	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2018	elif laue in ['2/m']:
2019	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2020	if uniq == 'a':
2021	Dij += parmDict[phfx+'D23']kl
2022	elif uniq == 'b':
2023	Dij += parmDict[phfx+'D13']hl
2024	elif uniq == 'c':
2025	Dij += parmDict[phfx+'D12']hk
2026	else:
2027	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2028	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2029	return Dijrefl[4]2tand(refl[5]/2.0)
2030
2031	def GetHStrainShiftDerv(refl,SGData,phfx):
2032	laue = SGData['SGLaue']
2033	uniq = SGData['SGUniq']
2034	h,k,l = refl[:3]
2035	if laue in ['m3','m3m']:
2036	dDijDict = {phfx+'D11':h2+k2+l**2,
2037	phfx+'eA':((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2038	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2039	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2040	elif laue in ['3R','3mR']:
2041	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2042	elif laue in ['4/m','4/mmm']:
2043	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2044	elif laue in ['mmm']:
2045	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2046	elif laue in ['2/m']:
2047	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2048	if uniq == 'a':
2049	dDijDict[phfx+'D23'] = k*l
2050	elif uniq == 'b':
2051	dDijDict[phfx+'D13'] = h*l
2052	elif uniq == 'c':
2053	dDijDict[phfx+'D12'] = h*k
2054	names.append()
2055	else:
2056	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2057	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2058	for item in dDijDict:
2059	dDijDict[item] = refl[4]2tand(refl[5]/2.0)
2060	return dDijDict
2061
2062	def GetFprime(controlDict,Histograms):
2063	FFtables = controlDict['FFtables']
2064	if not FFtables:
2065	return
2066	for histogram in Histograms:
2067	if 'PWDR' in histogram[:4]:
2068	Histogram = Histograms[histogram]
2069	hId = Histogram['hId']
2070	hfx = ':%d:'%(hId)
2071	keV = controlDict[hfx+'keV']
2072	for El in FFtables:
2073	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2074	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2075	FFtables[El][hfx+'FP'] = FP
2076	FFtables[El][hfx+'FPP'] = FPP
2077
2078	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2079
2080	def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
2081	U = parmDict[hfx+'U']
2082	V = parmDict[hfx+'V']
2083	W = parmDict[hfx+'W']
2084	X = parmDict[hfx+'X']
2085	Y = parmDict[hfx+'Y']
2086	tanPos = tand(refl[5]/2.0)
2087	sig = UtanPos2+VtanPos+W #save peak sigma
2088	sig = max(0.001,sig)
2089	gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
2090	gam = max(0.001,gam)
2091	return sig,gam
2092
2093	hId = Histogram['hId']
2094	hfx = ':%d:'%(hId)
2095	bakType = calcControls[hfx+'bakType']
2096	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2097	yc = np.zeros_like(yb)
2098
2099	if 'C' in calcControls[hfx+'histType']:
2100	shl = max(parmDict[hfx+'SH/L'],0.002)
2101	Ka2 = False
2102	if hfx+'Lam1' in parmDict.keys():
2103	wave = parmDict[hfx+'Lam1']
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	else:
2108	wave = parmDict[hfx+'Lam']
2109	else:
2110	print 'TOF Undefined at present'
2111	raise ValueError
2112	for phase in Histogram['Reflection Lists']:
2113	refList = Histogram['Reflection Lists'][phase]
2114	Phase = Phases[phase]
2115	pId = Phase['pId']
2116	pfx = '%d::'%(pId)
2117	phfx = '%d:%d:'%(pId,hId)
2118	hfx = ':%d:'%(hId)
2119	SGData = Phase['General']['SGData']
2120	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2121	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2122	Vst = np.sqrt(nl.det(G)) #V*
2123	if 'Pawley' not in Phase['General']['Type']:
2124	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2125	sizeEllipse = []
2126	if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2127	sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2128	for refl in refList:
2129	if 'C' in calcControls[hfx+'histType']:
2130	h,k,l = refl[:3]
2131	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
2132	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
2133	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
2134	refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse) #peak sig & gam
2135	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
2136	refl[13] = VstLorenz
2137	if 'Pawley' in Phase['General']['Type']:
2138	try:
2139	refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2140	except KeyError:
2141	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2142	continue
2143	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2144	iBeg = np.searchsorted(x,refl[5]-fmin)
2145	iFin = np.searchsorted(x,refl[5]+fmax)
2146	if not iBeg+iFin: #peak below low limit - skip peak
2147	continue
2148	elif not iBeg-iFin: #peak above high limit - done
2149	break
2150	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2151	if Ka2:
2152	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2153	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2154	iBeg = np.searchsorted(x,pos2-fmin)
2155	iFin = np.searchsorted(x,pos2+fmax)
2156	if not iBeg+iFin: #peak below low limit - skip peak
2157	continue
2158	elif not iBeg-iFin: #peak above high limit - done
2159	return yc,yb
2160	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here
2161	elif 'T' in calcControls[hfx+'histType']:
2162	print 'TOF Undefined at present'
2163	raise Exception #no TOF yet
2164	return yc,yb
2165
2166	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2167	for histogram in Histograms:
2168	if 'PWDR' in histogram[:4]:
2169	Histogram = Histograms[histogram]
2170	hId = Histogram['hId']
2171	hfx = ':%d:'%(hId)
2172	Limits = calcControls[hfx+'Limits']
2173	shl = max(parmDict[hfx+'SH/L'],0.002)
2174	Ka2 = False
2175	kRatio = 0.0
2176	if hfx+'Lam1' in parmDict.keys():
2177	Ka2 = True
2178	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2179	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2180	x,y,w,yc,yb,yd = Histogram['Data']
2181	ymb = np.array(y-yb)
2182	ycmb = np.array(yc-yb)
2183	ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)
2184	xB = np.searchsorted(x,Limits[0])
2185	xF = np.searchsorted(x,Limits[1])
2186	refLists = Histogram['Reflection Lists']
2187	for phase in refLists:
2188	Phase = Phases[phase]
2189	pId = Phase['pId']
2190	phfx = '%d:%d:'%(pId,hId)
2191	refList = refLists[phase]
2192	sumFo = 0.0
2193	sumdF = 0.0
2194	sumFosq = 0.0
2195	sumdFsq = 0.0
2196	for refl in refList:
2197	if 'C' in calcControls[hfx+'histType']:
2198	yp = np.zeros_like(yb)
2199	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2200	iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
2201	iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
2202	iFin2 = iFin
2203	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2204	if Ka2:
2205	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2206	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2207	iBeg2 = np.searchsorted(x,pos2-fmin)
2208	iFin2 = np.searchsorted(x,pos2+fmax)
2209	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
2210	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
2211	elif 'T' in calcControls[hfx+'histType']:
2212	print 'TOF Undefined at present'
2213	raise Exception #no TOF yet
2214	Fo = np.sqrt(np.abs(refl[8]))
2215	Fc = np.sqrt(np.abs(refl[9]))
2216	sumFo += Fo
2217	sumFosq += refl[8]**2
2218	sumdF += np.abs(Fo-Fc)
2219	sumdFsq += (refl[8]-refl[9])**2
2220	Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2221	Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2222	Histogram[phfx+'Nref'] = len(refList)
2223
2224	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2225
2226	def cellVaryDerv(pfx,SGData,dpdA):
2227	if SGData['SGLaue'] in ['-1',]:
2228	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2229	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2230	elif SGData['SGLaue'] in ['2/m',]:
2231	if SGData['SGUniq'] == 'a':
2232	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2233	elif SGData['SGUniq'] == 'b':
2234	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2235	else:
2236	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2237	elif SGData['SGLaue'] in ['mmm',]:
2238	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2239	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2240	return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2241	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2242	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2243	elif SGData['SGLaue'] in ['3R', '3mR']:
2244	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2245	elif SGData['SGLaue'] in ['m3m','m3']:
2246	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2247	# create a list of dependent variables and set up a dictionary to hold their derivatives
2248	dependentVars = G2mv.GetDependentVars()
2249	depDerivDict = {}
2250	for j in dependentVars:
2251	depDerivDict[j] = np.zeros(shape=(len(x)))
2252	#print 'dependent vars',dependentVars
2253	lenX = len(x)
2254	hId = Histogram['hId']
2255	hfx = ':%d:'%(hId)
2256	bakType = calcControls[hfx+'bakType']
2257	dMdv = np.zeros(shape=(len(varylist),len(x)))
2258	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2259	dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2260	bBpos =varylist.index(hfx+'Back:0')
2261	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2262
2263	if 'C' in calcControls[hfx+'histType']:
2264	dx = x[1]-x[0]
2265	shl = max(parmDict[hfx+'SH/L'],0.002)
2266	Ka2 = False
2267	if hfx+'Lam1' in parmDict.keys():
2268	wave = parmDict[hfx+'Lam1']
2269	Ka2 = True
2270	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2271	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2272	else:
2273	wave = parmDict[hfx+'Lam']
2274	else:
2275	print 'TOF Undefined at present'
2276	raise ValueError
2277	for phase in Histogram['Reflection Lists']:
2278	refList = Histogram['Reflection Lists'][phase]
2279	Phase = Phases[phase]
2280	SGData = Phase['General']['SGData']
2281	pId = Phase['pId']
2282	pfx = '%d::'%(pId)
2283	phfx = '%d:%d:'%(pId,hId)
2284	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2285	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2286	if 'Pawley' not in Phase['General']['Type']:
2287	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2288	sizeEllipse = []
2289	if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2290	sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2291	for iref,refl in enumerate(refList):
2292	if 'C' in calcControls[hfx+'histType']: #CW powder
2293	h,k,l = refl[:3]
2294	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2295	if 'Pawley' in Phase['General']['Type']:
2296	try:
2297	refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2298	except KeyError:
2299	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2300	continue
2301	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2302	iBeg = np.searchsorted(x,refl[5]-fmin)
2303	iFin = np.searchsorted(x,refl[5]+fmax)
2304	if not iBeg+iFin: #peak below low limit - skip peak
2305	continue
2306	elif not iBeg-iFin: #peak above high limit - done
2307	break
2308	pos = refl[5]
2309	tanth = tand(pos/2.0)
2310	costh = cosd(pos/2.0)
2311	dMdpk = np.zeros(shape=(6,len(x)))
2312	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2313	for i in range(1,5):
2314	dMdpk[i][iBeg:iFin] += 100.dxrefl[13]refl[9]dMdipk[i]
2315	dMdpk[0][iBeg:iFin] += 100.dxrefl[13]refl[9]dMdipk[0]
2316	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
2317	if Ka2:
2318	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2319	kdelt = int((pos2-refl[5])/dx)
2320	iBeg2 = min(lenX,iBeg+kdelt)
2321	iFin2 = min(lenX,iFin+kdelt)
2322	if iBeg2-iFin2:
2323	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2324	for i in range(1,5):
2325	dMdpk[i][iBeg2:iFin2] += 100.dxrefl[13]refl[9]kRatio*dMdipk2[i]
2326	dMdpk[0][iBeg2:iFin2] += 100.dxrefl[13]refl[9]kRatio*dMdipk2[0]
2327	dMdpk[5][iBeg2:iFin2] += 100.dxrefl[13]*dMdipk2[0]
2328	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
2329	if 'Pawley' in Phase['General']['Type']:
2330	try:
2331	idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2332	dMdv[idx] = dervDict['int']/refl[9]
2333	# Assuming Pawley variables not in constraints
2334	except ValueError:
2335	pass
2336	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2337	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2338	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2339	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2340	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2341	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2342	hfx+'DisplaceY':[dpdY,'pos'],}
2343	for name in names:
2344	if name in varylist:
2345	item = names[name]
2346	dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
2347	if name in dependentVars:
2348	item = names[name]
2349	depDerivDict[name] += item[0]*dervDict[item[1]]
2350
2351	for iPO in dIdPO:
2352	if iPO in varylist:
2353	dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
2354	if iPO in dependentVars:
2355	depDerivDict[iPO] = dIdPO[iPO]*dervDict['int']
2356
2357	for i,name in enumerate(['omega','chi','phi']):
2358	aname = pfx+'SH '+name
2359	if aname in varylist:
2360	dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
2361	if aname in dependentVars:
2362	depDerivDict[aname] += dFdSA[i]*dervDict['int']
2363	for iSH in dFdODF:
2364	if iSH in varylist:
2365	dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
2366	if iSH in dependentVars:
2367	depDerivDict[iSH] += dFdODF[iSH]*dervDict['int']
2368	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2369	for name,dpdA in cellDervNames:
2370	if name in varylist:
2371	dMdv[varylist.index(name)] += dpdA*dervDict['pos']
2372	if name in dependentVars:
2373	depDerivDict[name] += dpdA*dervDict['pos']
2374	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2375	for name in dDijDict:
2376	if name in varylist:
2377	dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
2378	if name in dependentVars:
2379	depDerivDict[name] += dDijDict[name]*dervDict['pos']
2380	gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
2381	for name in gamDict:
2382	if name in varylist:
2383	dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
2384	if name in dependentVars:
2385	depDerivDict[name] += gamDict[name]*dervDict['gam']
2386
2387	elif 'T' in calcControls[hfx+'histType']:
2388	print 'TOF Undefined at present'
2389	raise Exception #no TOF yet
2390	#do atom derivatives - for F,X & U so far
2391	corr = dervDict['int']/refl[9]
2392	for name in varylist+dependentVars:
2393	try:
2394	aname = name.split(pfx)[1][:2]
2395	if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2396	except IndexError:
2397	continue
2398	if name in varylist:
2399	dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
2400	if name in dependentVars:
2401	depDerivDict[name] += dFdvDict[name][iref]*corr
2402	# now process derivatives in constraints
2403	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2404	return dMdv
2405
2406	def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2407	parmdict.update(zip(varylist,values))
2408	G2mv.Dict2Map(parmdict,varylist)
2409	Histograms,Phases = HistoPhases
2410	dMdv = np.empty(0)
2411	for histogram in Histograms:
2412	if 'PWDR' in histogram[:4]:
2413	Histogram = Histograms[histogram]
2414	hId = Histogram['hId']
2415	hfx = ':%d:'%(hId)
2416	Limits = calcControls[hfx+'Limits']
2417	x,y,w,yc,yb,yd = Histogram['Data']
2418	xB = np.searchsorted(x,Limits[0])
2419	xF = np.searchsorted(x,Limits[1])
2420	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2421	varylist,Histogram,Phases,calcControls,pawleyLookup)
2422	if len(dMdv):
2423	dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2424	else:
2425	dMdv = dMdvh
2426	return dMdv
2427
2428	def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2429	parmdict.update(zip(varylist,values))
2430	Values2Dict(parmdict, varylist, values)
2431	G2mv.Dict2Map(parmdict,varylist)
2432	Histograms,Phases = HistoPhases
2433	M = np.empty(0)
2434	sumwYo = 0
2435	Nobs = 0
2436	for histogram in Histograms:
2437	if 'PWDR' in histogram[:4]:
2438	Histogram = Histograms[histogram]
2439	hId = Histogram['hId']
2440	hfx = ':%d:'%(hId)
2441	Limits = calcControls[hfx+'Limits']
2442	x,y,w,yc,yb,yd = Histogram['Data']
2443	yc *= 0.0 #zero full calcd profiles
2444	yb *= 0.0
2445	yd *= 0.0
2446	xB = np.searchsorted(x,Limits[0])
2447	xF = np.searchsorted(x,Limits[1])
2448	Histogram['Nobs'] = xF-xB
2449	Nobs += Histogram['Nobs']
2450	Histogram['sumwYo'] = np.sum(w[xB:xF]y[xB:xF]*2)
2451	sumwYo += Histogram['sumwYo']
2452	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2453	varylist,Histogram,Phases,calcControls,pawleyLookup)
2454	yc[xB:xF] += yb[xB:xF]
2455	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2456	Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2457	wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2458	Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy*2)/Histogram['sumwYo'])100.)
2459	M = np.concatenate((M,wdy))
2460	Histograms['sumwYo'] = sumwYo
2461	Histograms['Nobs'] = Nobs
2462	Rwp = min(100.,np.sqrt(np.sum(M*2)/sumwYo)100.)
2463	if dlg:
2464	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('wRp =',Rwp,'%'))[0]
2465	if not GoOn:
2466	parmDict['saved values'] = values
2467	raise Exception #Abort!!
2468	return M
2469
2470
2471	def Refine(GPXfile,dlg):
2472	import cPickle
2473	import pytexture as ptx
2474	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
2475
2476	ShowBanner()
2477	varyList = []
2478	parmDict = {}
2479	calcControls = {}
2480	G2mv.InitVars()
2481	Controls = GetControls(GPXfile)
2482	ShowControls(Controls)
2483	constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2484	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2485	if not Phases:
2486	print ' * ERROR - you have no histograms to refine! *'
2487	print ' * Refine aborted *'
2488	raise Exception
2489	if not Histograms:
2490	print ' * ERROR - you have no data to refine with! *'
2491	print ' * Refine aborted *'
2492	raise Exception
2493	Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2494	calcControls['Natoms'] = Natoms
2495	calcControls['FFtables'] = FFtables
2496	calcControls['BLtables'] = BLtables
2497	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2498	calcControls.update(controlDict)
2499	histVary,histDict,controlDict = GetHistogramData(Histograms)
2500	calcControls.update(controlDict)
2501	varyList = phaseVary+hapVary+histVary
2502	parmDict.update(phaseDict)
2503	parmDict.update(hapDict)
2504	parmDict.update(histDict)
2505	GetFprime(calcControls,Histograms)
2506	# do constraint processing
2507	try:
2508	groups,parmlist = G2mv.GroupConstraints(constrDict)
2509	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2510	except:
2511	print ' * ERROR - your constraints are internally inconsistent *'
2512	print ' * Refine aborted *'
2513	raise Exception
2514	# check to see which generated parameters are fully varied
2515	msg = G2mv.SetVaryFlags(varyList)
2516	if msg:
2517	print ' * ERROR - you have not set the refine flags for constraints consistently! *'
2518	print msg
2519	print ' * Refine aborted *'
2520	raise Exception
2521	G2mv.Map2Dict(parmDict,varyList)
2522	# print G2mv.VarRemapShow(varyList)
2523
2524	while True:
2525	begin = time.time()
2526	values = np.array(Dict2Values(parmDict, varyList))
2527	Ftol = Controls['min dM/M']
2528	Factor = Controls['shift factor']
2529	if Controls['deriv type'] == 'analytic':
2530	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2531	ftol=Ftol,col_deriv=True,factor=Factor,
2532	args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2533	ncyc = int(result[2]['nfev']/2)
2534	else: #'numeric'
2535	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2536	args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2537	ncyc = int(result[2]['nfev']/len(varyList))
2538	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2539	# for item in table: print item,table[item] #useful debug - are things shifting?
2540	runtime = time.time()-begin
2541	chisq = np.sum(result[2]['fvec']**2)
2542	Values2Dict(parmDict, varyList, result[0])
2543	G2mv.Dict2Map(parmDict,varyList)
2544	newCellDict = GetNewCellParms(parmDict,varyList)
2545	newAtomDict = ApplyXYZshifts(parmDict,varyList)
2546
2547	Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100. #to %
2548	GOF = chisq/(Histograms['Nobs']-len(varyList))
2549	print '\n Refinement results:'
2550	print 135*'-'
2551	print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2552	print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2553	print ' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rwp,chisq,GOF)
2554	try:
2555	covMatrix = result[1]*GOF
2556	sig = np.sqrt(np.diag(covMatrix))
2557	if np.any(np.isnan(sig)):
2558	print '* Least squares aborted - some invalid esds possible *'
2559	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2560	# for item in table: print item,table[item] #useful debug - are things shifting?
2561	break #refinement succeeded - finish up!
2562	except TypeError: #result[1] is None on singular matrix
2563	print '** Refinement failed - singular matrix **'
2564	Ipvt = result[2]['ipvt']
2565	for i,ipvt in enumerate(Ipvt):
2566	if not np.sum(result[2]['fjac'],axis=1)[i]:
2567	print 'Removing parameter: ',varyList[ipvt-1]
2568	del(varyList[ipvt-1])
2569	break
2570
2571	# print 'dependentParmList: ',G2mv.dependentParmList
2572	# print 'arrayList: ',G2mv.arrayList
2573	# print 'invarrayList: ',G2mv.invarrayList
2574	# print 'indParmList: ',G2mv.indParmList
2575	# print 'fixedDict: ',G2mv.fixedDict
2576	# print 'test1'
2577	GetFobsSq(Histograms,Phases,parmDict,calcControls)
2578	# print 'test2'
2579	sigDict = dict(zip(varyList,sig))
2580	covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2581	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2582	SetPhaseData(parmDict,sigDict,Phases,covData)
2583	SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2584	SetHistogramData(parmDict,sigDict,Histograms)
2585	G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
2586	SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2587	#for testing purposes!!!
2588	# file = open('structTestdata.dat','wb')
2589	# cPickle.dump(parmDict,file,1)
2590	# cPickle.dump(varyList,file,1)
2591	# for histogram in Histograms:
2592	# if 'PWDR' in histogram[:4]:
2593	# Histogram = Histograms[histogram]
2594	# cPickle.dump(Histogram,file,1)
2595	# cPickle.dump(Phases,file,1)
2596	# cPickle.dump(calcControls,file,1)
2597	# cPickle.dump(pawleyLookup,file,1)
2598	# file.close()
2599
2600	def SeqRefine(GPXfile,dlg):
2601	import cPickle
2602	import pytexture as ptx
2603	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
2604
2605	ShowBanner()
2606	print ' Sequential Refinement'
2607	G2mv.InitVars()
2608	Controls = GetControls(GPXfile)
2609	ShowControls(Controls)
2610	constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2611	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2612	if not Phases:
2613	print ' * ERROR - you have no histograms to refine! *'
2614	print ' * Refine aborted *'
2615	raise Exception
2616	if not Histograms:
2617	print ' * ERROR - you have no data to refine with! *'
2618	print ' * Refine aborted *'
2619	raise Exception
2620	Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
2621	if 'Seq Data' in Controls:
2622	histNames = Controls['Seq Data']
2623	else:
2624	histNames = GetHistogramNames(GPXfile,['PWDR',])
2625	if 'Reverse Seq' in Controls:
2626	if Controls['Reverse Seq']:
2627	histNames.reverse()
2628	SeqResult = {'histNames':histNames}
2629	makeBack = True
2630	for ihst,histogram in enumerate(histNames):
2631	ifPrint = False
2632	if dlg:
2633	dlg.SetTitle('Residual for histogram '+str(ihst))
2634	calcControls = {}
2635	calcControls['Natoms'] = Natoms
2636	calcControls['FFtables'] = FFtables
2637	calcControls['BLtables'] = BLtables
2638	varyList = []
2639	parmDict = {}
2640	Histo = {histogram:Histograms[histogram],}
2641	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
2642	calcControls.update(controlDict)
2643	histVary,histDict,controlDict = GetHistogramData(Histo,False)
2644	calcControls.update(controlDict)
2645	varyList = phaseVary+hapVary+histVary
2646	if not ihst:
2647	saveVaryList = varyList[:]
2648	for i,item in enumerate(saveVaryList):
2649	items = item.split(':')
2650	if items[1]:
2651	items[1] = ''
2652	item = ':'.join(items)
2653	saveVaryList[i] = item
2654	SeqResult['varyList'] = saveVaryList
2655	else:
2656	newVaryList = varyList[:]
2657	for i,item in enumerate(newVaryList):
2658	items = item.split(':')
2659	if items[1]:
2660	items[1] = ''
2661	item = ':'.join(items)
2662	newVaryList[i] = item
2663	if newVaryList != SeqResult['varyList']:
2664	print newVaryList
2665	print SeqResult['varyList']
2666	print '**** ERROR - variable list for this histogram does not match previous'
2667	raise Exception
2668	parmDict.update(phaseDict)
2669	parmDict.update(hapDict)
2670	parmDict.update(histDict)
2671	GetFprime(calcControls,Histo)
2672	constrDict,constrFlag,fixedList = G2mv.InputParse([]) #constraints go here?
2673	groups,parmlist = G2mv.GroupConstraints(constrDict)
2674	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2675	G2mv.Map2Dict(parmDict,varyList)
2676
2677	while True:
2678	begin = time.time()
2679	values = np.array(Dict2Values(parmDict, varyList))
2680	Ftol = Controls['min dM/M']
2681	Factor = Controls['shift factor']
2682	if Controls['deriv type'] == 'analytic':
2683	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2684	ftol=Ftol,col_deriv=True,factor=Factor,
2685	args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2686	ncyc = int(result[2]['nfev']/2)
2687	else: #'numeric'
2688	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2689	args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2690	ncyc = int(result[2]['nfev']/len(varyList))
2691	runtime = time.time()-begin
2692	chisq = np.sum(result[2]['fvec']**2)
2693	Values2Dict(parmDict, varyList, result[0])
2694	G2mv.Dict2Map(parmDict,varyList)
2695	newCellDict = GetNewCellParms(parmDict,varyList)
2696	newAtomDict = ApplyXYZshifts(parmDict,varyList)
2697
2698	Rwp = np.sqrt(chisq/Histo['sumwYo'])*100. #to %
2699	GOF = chisq/(Histo['Nobs']-len(varyList))
2700	print '\n Refinement results for histogram: v'+histogram
2701	print 135*'-'
2702	print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
2703	print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2704	print ' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rwp,chisq,GOF)
2705	try:
2706	covMatrix = result[1]*GOF
2707	sig = np.sqrt(np.diag(covMatrix))
2708	if np.any(np.isnan(sig)):
2709	print '* Least squares aborted - some invalid esds possible *'
2710	ifPrint = True
2711	break #refinement succeeded - finish up!
2712	except TypeError: #result[1] is None on singular matrix
2713	print '** Refinement failed - singular matrix **'
2714	Ipvt = result[2]['ipvt']
2715	for i,ipvt in enumerate(Ipvt):
2716	if not np.sum(result[2]['fjac'],axis=1)[i]:
2717	print 'Removing parameter: ',varyList[ipvt-1]
2718	del(varyList[ipvt-1])
2719	break
2720
2721	GetFobsSq(Histo,Phases,parmDict,calcControls)
2722	sigDict = dict(zip(varyList,sig))
2723	covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2724	'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2725	SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
2726	SetHistogramData(parmDict,sigDict,Histo,ifPrint)
2727	SeqResult[histogram] = covData
2728	SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
2729	makeBack = False
2730	SetSeqResult(GPXfile,Histograms,SeqResult)
2731
2732
2733	def main():
2734	arg = sys.argv
2735	if len(arg) > 1:
2736	GPXfile = arg[1]
2737	if not ospath.exists(GPXfile):
2738	print 'ERROR - ',GPXfile," doesn't exist!"
2739	exit()
2740	GPXpath = ospath.dirname(arg[1])
2741	Refine(GPXfile,None)
2742	else:
2743	print 'ERROR - missing filename'
2744	exit()
2745	print "Done"
2746
2747	if __name__ == '__main__':
2748	main()

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

Download in other formats: