Context Navigation

source: trunk/GSASIIstruct.py @ 609

Last change on this file since 609 was 609, checked in by vondreele, 11 years ago
adding more tutorials & fix some uncovered bugs
Property svn:keywords set to `Date Author Revision URL Id`
File size: 148.7 KB

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

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

Download in other formats: