Context Navigation

source: trunk/GSASIIstrIO.py @ 1496

Last change on this file since 1496 was 1496, checked in by vondreele, 9 years ago
Add sequential peak fitting. Change Back:n to Back;n for background parameters Also DebyeA, etc.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 113.7 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2014-09-15 00:48:39 +0000 (Mon, 15 Sep 2014) $
4	# $Author: vondreele $
5	# $Revision: 1496 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 1496 2014-09-15 00:48:39Z vondreele $
8	########### SVN repository information ###################
9	'''
10	GSASIIstrIO: structure I/O routines
11	-------------------------------------
12
13	'''
14	import sys
15	import os
16	import os.path as ospath
17	import time
18	import math
19	import copy
20	import random
21	import cPickle
22	import numpy as np
23	import numpy.ma as ma
24	import numpy.linalg as nl
25	import scipy.optimize as so
26	import GSASIIpath
27	GSASIIpath.SetVersionNumber("$Revision: 1496 $")
28	import GSASIIElem as G2el
29	import GSASIIgrid as G2gd
30	import GSASIIlattice as G2lat
31	import GSASIIspc as G2spc
32	import GSASIIobj as G2obj
33	import GSASIImapvars as G2mv
34	import GSASIImath as G2mth
35
36	sind = lambda x: np.sin(x*np.pi/180.)
37	cosd = lambda x: np.cos(x*np.pi/180.)
38	tand = lambda x: np.tan(x*np.pi/180.)
39	asind = lambda x: 180.*np.arcsin(x)/np.pi
40	acosd = lambda x: 180.*np.arccos(x)/np.pi
41	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42
43	ateln2 = 8.0*math.log(2.0)
44
45	def GetControls(GPXfile):
46	''' Returns dictionary of control items found in GSASII gpx file
47
48	:param str GPXfile: full .gpx file name
49	:return: dictionary of control items
50	'''
51	Controls = copy.copy(G2obj.DefaultControls)
52	fl = open(GPXfile,'rb')
53	while True:
54	try:
55	data = cPickle.load(fl)
56	except EOFError:
57	break
58	datum = data[0]
59	if datum[0] == 'Controls':
60	Controls.update(datum[1])
61	fl.close()
62	return Controls
63
64	def GetConstraints(GPXfile):
65	'''Read the constraints from the GPX file and interpret them
66
67	called in :func:`CheckConstraints`, :func:`GSASIIstrMain.Refine`
68	and :func:`GSASIIstrMain.SeqRefine`.
69	'''
70	fl = open(GPXfile,'rb')
71	while True:
72	try:
73	data = cPickle.load(fl)
74	except EOFError:
75	break
76	datum = data[0]
77	if datum[0] == 'Constraints':
78	constList = []
79	for item in datum[1]:
80	if item.startswith('_'): continue
81	constList += datum[1][item]
82	fl.close()
83	constDict,fixedList,ignored = ProcessConstraints(constList)
84	if ignored:
85	print ignored,'Constraints were rejected. Was a constrained phase, histogram or atom deleted?'
86	return constDict,fixedList
87	fl.close()
88	raise Exception,"No constraints in GPX file"
89
90	def ProcessConstraints(constList):
91	"""Interpret the constraints in the constList input into a dictionary, etc.
92	All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
93	phase/hist/atoms based on the object internals (random Ids). If this can't be
94	done (if a phase has been deleted, etc.), the variable is ignored.
95	If the constraint cannot be used due to too many dropped variables,
96	it is counted as ignored.
97
98	:param list constList: a list of lists where each item in the outer list
99	specifies a constraint of some form, as described in the :mod:`GSASIIobj`
100	:ref:`Constraint definition <Constraint_definitions_table>`.
101
102	:returns: a tuple of (constDict,fixedList,ignored) where:
103
104	* constDict (list of dicts) contains the constraint relationships
105	* fixedList (list) contains the fixed values for each type
106	of constraint.
107	* ignored (int) counts the number of invalid constraint items
108	(should always be zero!)
109	"""
110	constDict = []
111	fixedList = []
112	ignored = 0
113	for item in constList:
114	if item[-1] == 'h':
115	# process a hold
116	fixedList.append('0')
117	var = str(item[0][1])
118	if '?' not in var:
119	constDict.append({var:0.0})
120	else:
121	ignored += 1
122	elif item[-1] == 'f':
123	# process a new variable
124	fixedList.append(None)
125	D = {}
126	varyFlag = item[-2]
127	varname = item[-3]
128	for term in item[:-3]:
129	var = str(term[1])
130	if '?' not in var:
131	D[var] = term[0]
132	if len(D) > 1:
133	# add extra dict terms for input variable name and vary flag
134	if varname is not None:
135	if varname.startswith('::'):
136	varname = varname[2:].replace(':',';')
137	else:
138	varname = varname.replace(':',';')
139	D['_name'] = '::' + varname
140	D['_vary'] = varyFlag == True # force to bool
141	constDict.append(D)
142	else:
143	ignored += 1
144	#constFlag[-1] = ['Vary']
145	elif item[-1] == 'c':
146	# process a contraint relationship
147	D = {}
148	for term in item[:-3]:
149	var = str(term[1])
150	if '?' not in var:
151	D[var] = term[0]
152	if len(D) >= 1:
153	fixedList.append(str(item[-3]))
154	constDict.append(D)
155	else:
156	ignored += 1
157	elif item[-1] == 'e':
158	# process an equivalence
159	firstmult = None
160	eqlist = []
161	for term in item[:-3]:
162	if term[0] == 0: term[0] = 1.0
163	var = str(term[1])
164	if '?' in var: continue
165	if firstmult is None:
166	firstmult = term[0]
167	firstvar = var
168	else:
169	eqlist.append([var,firstmult/term[0]])
170	if len(eqlist) > 0:
171	G2mv.StoreEquivalence(firstvar,eqlist)
172	else:
173	ignored += 1
174	else:
175	ignored += 1
176	return constDict,fixedList,ignored
177
178	def CheckConstraints(GPXfile):
179	'''Load constraints and related info and return any error or warning messages'''
180	# init constraints
181	G2mv.InitVars()
182	# get variables
183	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
184	if not Phases:
185	return 'Error: No Phases!',''
186	if not Histograms:
187	return 'Error: no diffraction data',''
188	rigidbodyDict = GetRigidBodies(GPXfile)
189	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
190	rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
191	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
192	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
193	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
194	varyList = rbVary+phaseVary+hapVary+histVary
195	constrDict,fixedList = GetConstraints(GPXfile)
196	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
197
198	def GetRestraints(GPXfile):
199	'''Read the restraints from the GPX file.
200	Throws an exception if not found in the .GPX file
201	'''
202	fl = open(GPXfile,'rb')
203	while True:
204	try:
205	data = cPickle.load(fl)
206	except EOFError:
207	break
208	datum = data[0]
209	if datum[0] == 'Restraints':
210	restraintDict = datum[1]
211	fl.close()
212	return restraintDict
213
214	def GetRigidBodies(GPXfile):
215	'''Read the rigid body models from the GPX file
216	'''
217	fl = open(GPXfile,'rb')
218	while True:
219	try:
220	data = cPickle.load(fl)
221	except EOFError:
222	break
223	datum = data[0]
224	if datum[0] == 'Rigid bodies':
225	rigidbodyDict = datum[1]
226	fl.close()
227	return rigidbodyDict
228
229	def GetFprime(controlDict,Histograms):
230	'Needs a doc string'
231	FFtables = controlDict['FFtables']
232	if not FFtables:
233	return
234	histoList = Histograms.keys()
235	histoList.sort()
236	for histogram in histoList:
237	if histogram[:4] in ['PWDR','HKLF']:
238	Histogram = Histograms[histogram]
239	hId = Histogram['hId']
240	hfx = ':%d:'%(hId)
241	if 'X' in controlDict[hfx+'histType']:
242	keV = controlDict[hfx+'keV']
243	for El in FFtables:
244	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
245	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
246	FFtables[El][hfx+'FP'] = FP
247	FFtables[El][hfx+'FPP'] = FPP
248
249	def GetPhaseNames(GPXfile):
250	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
251
252	:param str GPXfile: full .gpx file name
253	:return: list of phase names
254	'''
255	fl = open(GPXfile,'rb')
256	PhaseNames = []
257	while True:
258	try:
259	data = cPickle.load(fl)
260	except EOFError:
261	break
262	datum = data[0]
263	if 'Phases' == datum[0]:
264	for datus in data[1:]:
265	PhaseNames.append(datus[0])
266	fl.close()
267	return PhaseNames
268
269	def GetAllPhaseData(GPXfile,PhaseName):
270	''' Returns the entire dictionary for PhaseName from GSASII gpx file
271
272	:param str GPXfile: full .gpx file name
273	:param str PhaseName: phase name
274	:return: phase dictionary
275	'''
276	fl = open(GPXfile,'rb')
277	General = {}
278	Atoms = []
279	while True:
280	try:
281	data = cPickle.load(fl)
282	except EOFError:
283	break
284	datum = data[0]
285	if 'Phases' == datum[0]:
286	for datus in data[1:]:
287	if datus[0] == PhaseName:
288	break
289	fl.close()
290	return datus[1]
291
292	def GetHistograms(GPXfile,hNames):
293	""" Returns a dictionary of histograms found in GSASII gpx file
294
295	:param str GPXfile: full .gpx file name
296	:param str hNames: list of histogram names
297	:return: dictionary of histograms (types = PWDR & HKLF)
298
299	"""
300	fl = open(GPXfile,'rb')
301	Histograms = {}
302	while True:
303	try:
304	data = cPickle.load(fl)
305	except EOFError:
306	break
307	datum = data[0]
308	hist = datum[0]
309	if hist in hNames:
310	if 'PWDR' in hist[:4]:
311	PWDRdata = {}
312	PWDRdata.update(datum[1][0]) #weight factor
313	PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1])) #masked powder data arrays/clear previous masks
314	PWDRdata[data[2][0]] = data[2][1] #Limits & excluded regions (if any)
315	PWDRdata[data[3][0]] = data[3][1] #Background
316	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
317	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
318	try:
319	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
320	except IndexError:
321	PWDRdata['Reflection Lists'] = {}
322	PWDRdata['Residuals'] = {}
323
324	Histograms[hist] = PWDRdata
325	elif 'HKLF' in hist[:4]:
326	HKLFdata = {}
327	HKLFdata.update(datum[1][0]) #weight factor
328	#patch
329	if isinstance(datum[1][1],list):
330	RefData = {'RefList':[],'FF':{}}
331	for ref in datum[1][1]:
332	RefData['RefList'].append(ref[:11]+[ref[13],])
333	RefData['RefList'] = np.array(RefData['RefList'])
334	datum[1][1] = RefData
335	#end patch
336	datum[1][1]['FF'] = {}
337	HKLFdata['Data'] = datum[1][1]
338	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
339	HKLFdata['Reflection Lists'] = None
340	HKLFdata['Residuals'] = {}
341	Histograms[hist] = HKLFdata
342	fl.close()
343	return Histograms
344
345	def GetHistogramNames(GPXfile,hType):
346	""" Returns a list of histogram names found in GSASII gpx file
347
348	:param str GPXfile: full .gpx file name
349	:param str hType: list of histogram types
350	:return: list of histogram names (types = PWDR & HKLF)
351
352	"""
353	fl = open(GPXfile,'rb')
354	HistogramNames = []
355	while True:
356	try:
357	data = cPickle.load(fl)
358	except EOFError:
359	break
360	datum = data[0]
361	if datum[0][:4] in hType:
362	HistogramNames.append(datum[0])
363	fl.close()
364	return HistogramNames
365
366	def GetUsedHistogramsAndPhases(GPXfile):
367	''' Returns all histograms that are found in any phase
368	and any phase that uses a histogram. This also
369	assigns numbers to used phases and histograms by the
370	order they appear in the file.
371
372	:param str GPXfile: full .gpx file name
373	:returns: (Histograms,Phases)
374
375	* Histograms = dictionary of histograms as {name:data,...}
376	* Phases = dictionary of phases that use histograms
377
378	'''
379	phaseNames = GetPhaseNames(GPXfile)
380	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
381	allHistograms = GetHistograms(GPXfile,histoList)
382	phaseData = {}
383	for name in phaseNames:
384	phaseData[name] = GetAllPhaseData(GPXfile,name)
385	Histograms = {}
386	Phases = {}
387	for phase in phaseData:
388	Phase = phaseData[phase]
389	if Phase['Histograms']:
390	if phase not in Phases:
391	pId = phaseNames.index(phase)
392	Phase['pId'] = pId
393	Phases[phase] = Phase
394	for hist in Phase['Histograms']:
395	if 'Use' not in Phase['Histograms'][hist]: #patch
396	Phase['Histograms'][hist]['Use'] = True
397	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
398	try:
399	Histograms[hist] = allHistograms[hist]
400	hId = histoList.index(hist)
401	Histograms[hist]['hId'] = hId
402	except KeyError: # would happen if a referenced histogram were
403	# renamed or deleted
404	print('For phase "'+str(phase)+
405	'" unresolved reference to histogram "'+str(hist)+'"')
406	G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
407	return Histograms,Phases
408
409	def getBackupName(GPXfile,makeBack):
410	'''
411	Get the name for the backup .gpx file name
412
413	:param str GPXfile: full .gpx file name
414	:param bool makeBack: if True the name of a new file is returned, if
415	False the name of the last file that exists is returned
416	:returns: the name of a backup file
417
418	'''
419	GPXpath,GPXname = ospath.split(GPXfile)
420	if GPXpath == '': GPXpath = '.'
421	Name = ospath.splitext(GPXname)[0]
422	files = os.listdir(GPXpath)
423	last = 0
424	for name in files:
425	name = name.split('.')
426	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
427	if makeBack:
428	last = max(last,int(name[1].strip('bak'))+1)
429	else:
430	last = max(last,int(name[1].strip('bak')))
431	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
432	return GPXback
433
434	def GPXBackup(GPXfile,makeBack=True):
435	'''
436	makes a backup of the current .gpx file (?)
437
438	:param str GPXfile: full .gpx file name
439	:param bool makeBack: if True (default), the backup is written to
440	a new file; if False, the last backup is overwritten
441	:returns: the name of the backup file that was written
442	'''
443	import distutils.file_util as dfu
444	GPXback = getBackupName(GPXfile,makeBack)
445	dfu.copy_file(GPXfile,GPXback)
446	return GPXback
447
448	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
449	''' Updates gpxfile from all histograms that are found in any phase
450	and any phase that used a histogram. Also updates rigid body definitions.
451
452
453	:param str GPXfile: full .gpx file name
454	:param dict Histograms: dictionary of histograms as {name:data,...}
455	:param dict Phases: dictionary of phases that use histograms
456	:param dict RigidBodies: dictionary of rigid bodies
457	:param dict CovData: dictionary of refined variables, varyList, & covariance matrix
458	:param bool makeBack: True if new backup of .gpx file is to be made; else use the last one made
459
460	'''
461
462	GPXback = GPXBackup(GPXfile,makeBack)
463	print 'Read from file:',GPXback
464	print 'Save to file :',GPXfile
465	infile = open(GPXback,'rb')
466	outfile = open(GPXfile,'wb')
467	while True:
468	try:
469	data = cPickle.load(infile)
470	except EOFError:
471	break
472	datum = data[0]
473	# print 'read: ',datum[0]
474	if datum[0] == 'Phases':
475	for iphase in range(len(data)):
476	if data[iphase][0] in Phases:
477	phaseName = data[iphase][0]
478	data[iphase][1].update(Phases[phaseName])
479	elif datum[0] == 'Covariance':
480	data[0][1] = CovData
481	elif datum[0] == 'Rigid bodies':
482	data[0][1] = RigidBodies
483	try:
484	histogram = Histograms[datum[0]]
485	# print 'found ',datum[0]
486	data[0][1][1] = histogram['Data']
487	data[0][1][0].update(histogram['Residuals'])
488	for datus in data[1:]:
489	# print ' read: ',datus[0]
490	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
491	datus[1] = histogram[datus[0]]
492	except KeyError:
493	pass
494
495	cPickle.dump(data,outfile,1)
496	infile.close()
497	outfile.close()
498	print 'GPX file save successful'
499
500	def SetSeqResult(GPXfile,Histograms,SeqResult):
501	'''
502	Needs doc string
503
504	:param str GPXfile: full .gpx file name
505	'''
506	GPXback = GPXBackup(GPXfile)
507	print 'Read from file:',GPXback
508	print 'Save to file :',GPXfile
509	infile = open(GPXback,'rb')
510	outfile = open(GPXfile,'wb')
511	while True:
512	try:
513	data = cPickle.load(infile)
514	except EOFError:
515	break
516	datum = data[0]
517	if datum[0] == 'Sequential results':
518	data[0][1] = SeqResult
519	# reset the Copy Next flag, since it should not be needed twice in a row
520	if datum[0] == 'Controls':
521	data[0][1]['Copy2Next'] = False
522	try:
523	histogram = Histograms[datum[0]]
524	data[0][1][1] = list(histogram['Data'])
525	for datus in data[1:]:
526	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
527	datus[1] = histogram[datus[0]]
528	except KeyError:
529	pass
530
531	cPickle.dump(data,outfile,1)
532	infile.close()
533	outfile.close()
534	print 'GPX file save successful'
535
536	def ShowBanner(pFile=None):
537	'Print authorship, copyright and citation notice'
538	print >>pFile,80''
539	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
540	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
541	print >>pFile,' Argonne National Laboratory(C), 2010'
542	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
543	print >>pFile,' as Operator of Argonne National Laboratory.'
544	print >>pFile,' Please cite:'
545	print >>pFile,' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
546
547	print >>pFile,80'','\n'
548
549	def ShowControls(Controls,pFile=None,SeqRef=False):
550	'Print controls information'
551	print >>pFile,' Least squares controls:'
552	print >>pFile,' Refinement type: ',Controls['deriv type']
553	if 'Hessian' in Controls['deriv type']:
554	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
555	else:
556	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
557	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
558	if SeqRef:
559	print >>pFile,' Sequential refinement controls:'
560	print >>pFile,' Copy of histogram results to next: ',Controls['Copy2Next']
561	print >>pFile,' Process histograms in reverse order: ',Controls['Reverse Seq']
562
563	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
564	'needs a doc string'
565	# if SGLaue in ['-1','2/m','mmm']:
566	# return #no Pawley symmetry required constraints
567	eqvDict = {}
568	for i,varyI in enumerate(pawleyVary):
569	eqvDict[varyI] = []
570	refI = int(varyI.split(':')[-1])
571	ih,ik,il = PawleyRef[refI][:3]
572	dspI = PawleyRef[refI][4]
573	for varyJ in pawleyVary[i+1:]:
574	refJ = int(varyJ.split(':')[-1])
575	jh,jk,jl = PawleyRef[refJ][:3]
576	dspJ = PawleyRef[refJ][4]
577	if SGLaue in ['4/m','4/mmm']:
578	isum = ih2+ik2
579	jsum = jh2+jk2
580	if abs(il) == abs(jl) and isum == jsum:
581	eqvDict[varyI].append(varyJ)
582	elif SGLaue in ['3R','3mR']:
583	isum = ih2+ik2+il**2
584	jsum = jh2+jk2+jl**2
585	isum2 = ihik+ihil+ik*il
586	jsum2 = jhjk+jhjl+jk*jl
587	if isum == jsum and isum2 == jsum2:
588	eqvDict[varyI].append(varyJ)
589	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
590	isum = ih2+ik2+ih*ik
591	jsum = jh2+jk2+jh*jk
592	if abs(il) == abs(jl) and isum == jsum:
593	eqvDict[varyI].append(varyJ)
594	elif SGLaue in ['m3','m3m']:
595	isum = ih2+ik2+il**2
596	jsum = jh2+jk2+jl**2
597	if isum == jsum:
598	eqvDict[varyI].append(varyJ)
599	elif abs(dspI-dspJ)/dspI < 1.e-4:
600	eqvDict[varyI].append(varyJ)
601	for item in pawleyVary:
602	if eqvDict[item]:
603	for item2 in pawleyVary:
604	if item2 in eqvDict[item]:
605	eqvDict[item2] = []
606	G2mv.StoreEquivalence(item,eqvDict[item])
607
608	def cellVary(pfx,SGData):
609	'needs a doc string'
610	if SGData['SGLaue'] in ['-1',]:
611	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
612	elif SGData['SGLaue'] in ['2/m',]:
613	if SGData['SGUniq'] == 'a':
614	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
615	elif SGData['SGUniq'] == 'b':
616	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
617	else:
618	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
619	elif SGData['SGLaue'] in ['mmm',]:
620	return [pfx+'A0',pfx+'A1',pfx+'A2']
621	elif SGData['SGLaue'] in ['4/m','4/mmm']:
622	return [pfx+'A0',pfx+'A2']
623	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
624	return [pfx+'A0',pfx+'A2']
625	elif SGData['SGLaue'] in ['3R', '3mR']:
626	return [pfx+'A0',pfx+'A3']
627	elif SGData['SGLaue'] in ['m3m','m3']:
628	return [pfx+'A0',]
629
630	################################################################################
631	##### Rigid Body Models and not General.get('doPawley')
632	################################################################################
633
634	def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
635	'needs a doc string'
636
637	def PrintResRBModel(RBModel):
638	atNames = RBModel['atNames']
639	rbRef = RBModel['rbRef']
640	rbSeq = RBModel['rbSeq']
641	print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
642	'No. times used: ',RBModel['useCount']
643	print >>pFile,' At name x y z'
644	for name,xyz in zip(atNames,RBModel['rbXYZ']):
645	print >>pFile,' %8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
646	print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
647	' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
648	if rbSeq:
649	for i,rbseq in enumerate(rbSeq):
650	print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
651	atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
652
653	def PrintVecRBModel(RBModel):
654	rbRef = RBModel['rbRef']
655	atTypes = RBModel['rbTypes']
656	print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
657	'No. times used: ',RBModel['useCount']
658	for i in range(len(RBModel['VectMag'])):
659	print >>pFile,'Vector no.: ',i,' Magnitude: ', \
660	'%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
661	print >>pFile,' No. Type vx vy vz'
662	for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
663	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
664	print >>pFile,' No. Type x y z'
665	for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
666	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
667	print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
668	' & atom ',rbRef[0],' -> atom ',rbRef[2]
669	rbVary = []
670	rbDict = {}
671	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
672	if len(rbIds['Vector']):
673	for irb,item in enumerate(rbIds['Vector']):
674	if rigidbodyDict['Vector'][item]['useCount']:
675	RBmags = rigidbodyDict['Vector'][item]['VectMag']
676	RBrefs = rigidbodyDict['Vector'][item]['VectRef']
677	for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
678	pid = '::RBV;'+str(i)+':'+str(irb)
679	rbDict[pid] = mag
680	if ref:
681	rbVary.append(pid)
682	if Print:
683	print >>pFile,'\nVector rigid body model:'
684	PrintVecRBModel(rigidbodyDict['Vector'][item])
685	if len(rbIds['Residue']):
686	for item in rbIds['Residue']:
687	if rigidbodyDict['Residue'][item]['useCount']:
688	if Print:
689	print >>pFile,'\nResidue rigid body model:'
690	PrintResRBModel(rigidbodyDict['Residue'][item])
691	return rbVary,rbDict
692
693	def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
694	'needs a doc string'
695
696	def PrintRBVectandSig(VectRB,VectSig):
697	print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
698	namstr = ' names :'
699	valstr = ' values:'
700	sigstr = ' esds :'
701	for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
702	namstr += '%12s'%('Vect '+str(i))
703	valstr += '%12.4f'%(val)
704	if sig:
705	sigstr += '%12.4f'%(sig)
706	else:
707	sigstr += 12*' '
708	print >>pFile,namstr
709	print >>pFile,valstr
710	print >>pFile,sigstr
711
712	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
713	if not RBIds['Vector']:
714	return
715	for irb,item in enumerate(RBIds['Vector']):
716	if rigidbodyDict['Vector'][item]['useCount']:
717	VectSig = []
718	RBmags = rigidbodyDict['Vector'][item]['VectMag']
719	for i,mag in enumerate(RBmags):
720	name = '::RBV;'+str(i)+':'+str(irb)
721	mag = parmDict[name]
722	if name in sigDict:
723	VectSig.append(sigDict[name])
724	PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)
725
726	################################################################################
727	##### Phase data
728	################################################################################
729
730	def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
731	'needs a doc string'
732
733	def PrintFFtable(FFtable):
734	print >>pFile,'\n X-ray scattering factors:'
735	print >>pFile,' Symbol fa fb fc'
736	print >>pFile,99*'-'
737	for Ename in FFtable:
738	ffdata = FFtable[Ename]
739	fa = ffdata['fa']
740	fb = ffdata['fb']
741	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
742	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
743
744	def PrintBLtable(BLtable):
745	print >>pFile,'\n Neutron scattering factors:'
746	print >>pFile,' Symbol isotope mass b resonant terms'
747	print >>pFile,99*'-'
748	for Ename in BLtable:
749	bldata = BLtable[Ename]
750	isotope = bldata[0]
751	mass = bldata[1]['Mass']
752	if 'SL' in bldata[1]:
753	blen = bldata[1]['SL'][0]
754	bres = []
755	else:
756	blen = 0
757	bres = bldata[1]['BW-LS']
758	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
759	for item in bres:
760	line += '%10.5g'%(item)
761	print >>pFile,line
762
763	def PrintRBObjects(resRBData,vecRBData):
764
765	def PrintRBThermals():
766	tlstr = ['11','22','33','12','13','23']
767	sstr = ['12','13','21','23','31','32','AA','BB']
768	TLS = RB['ThermalMotion'][1]
769	TLSvar = RB['ThermalMotion'][2]
770	if 'T' in RB['ThermalMotion'][0]:
771	print >>pFile,'TLS data'
772	text = ''
773	for i in range(6):
774	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
775	print >>pFile,text
776	if 'L' in RB['ThermalMotion'][0]:
777	text = ''
778	for i in range(6,12):
779	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
780	print >>pFile,text
781	if 'S' in RB['ThermalMotion'][0]:
782	text = ''
783	for i in range(12,20):
784	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
785	print >>pFile,text
786	if 'U' in RB['ThermalMotion'][0]:
787	print >>pFile,'Uiso data'
788	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
789
790	if len(resRBData):
791	for RB in resRBData:
792	Oxyz = RB['Orig'][0]
793	Qrijk = RB['Orient'][0]
794	Angle = 2.0*acosd(Qrijk[0])
795	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
796	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
797	print >>pFile,'Orientation angle,vector:', \
798	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
799	Torsions = RB['Torsions']
800	if len(Torsions):
801	text = 'Torsions: '
802	for torsion in Torsions:
803	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
804	print >>pFile,text
805	PrintRBThermals()
806	if len(vecRBData):
807	for RB in vecRBData:
808	Oxyz = RB['Orig'][0]
809	Qrijk = RB['Orient'][0]
810	Angle = 2.0*acosd(Qrijk[0])
811	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
812	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
813	print >>pFile,'Orientation angle,vector:', \
814	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
815	PrintRBThermals()
816
817	def PrintAtoms(General,Atoms):
818	cx,ct,cs,cia = General['AtomPtrs']
819	print >>pFile,'\n Atoms:'
820	line = ' name type refine? x y z '+ \
821	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
822	if General['Type'] == 'magnetic':
823	line += ' Mx My Mz'
824	elif General['Type'] == 'macromolecular':
825	line = ' res no residue chain'+line
826	print >>pFile,line
827	if General['Type'] == 'nuclear':
828	print >>pFile,135*'-'
829	for i,at in enumerate(Atoms):
830	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
831	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
832	if at[cia] == 'I':
833	line += '%8.5f'%(at[cia+1])+48*' '
834	else:
835	line += 8*' '
836	for j in range(6):
837	line += '%8.5f'%(at[cia+2+j])
838	print >>pFile,line
839	elif General['Type'] == 'macromolecular':
840	print >>pFile,135*'-'
841	for i,at in enumerate(Atoms):
842	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
843	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
844	if at[cia] == 'I':
845	line += '%8.4f'%(at[cia+1])+48*' '
846	else:
847	line += 8*' '
848	for j in range(6):
849	line += '%8.4f'%(at[cia+2+j])
850	print >>pFile,line
851
852	def PrintTexture(textureData):
853	topstr = '\n Spherical harmonics texture: Order:' + \
854	str(textureData['Order'])
855	if textureData['Order']:
856	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
857	else:
858	print >>pFile,topstr
859	return
860	names = ['omega','chi','phi']
861	line = '\n'
862	for name in names:
863	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
864	print >>pFile,line
865	print >>pFile,'\n Texture coefficients:'
866	ptlbls = ' names :'
867	ptstr = ' values:'
868	SHcoeff = textureData['SH Coeff'][1]
869	for item in SHcoeff:
870	ptlbls += '%12s'%(item)
871	ptstr += '%12.4f'%(SHcoeff[item])
872	print >>pFile,ptlbls
873	print >>pFile,ptstr
874
875	def MakeRBParms(rbKey,phaseVary,phaseDict):
876	rbid = str(rbids.index(RB['RBId']))
877	pfxRB = pfx+'RB'+rbKey+'P'
878	pstr = ['x','y','z']
879	ostr = ['a','i','j','k']
880	for i in range(3):
881	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
882	phaseDict[name] = RB['Orig'][0][i]
883	if RB['Orig'][1]:
884	phaseVary += [name,]
885	pfxRB = pfx+'RB'+rbKey+'O'
886	for i in range(4):
887	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
888	phaseDict[name] = RB['Orient'][0][i]
889	if RB['Orient'][1] == 'AV' and i:
890	phaseVary += [name,]
891	elif RB['Orient'][1] == 'A' and not i:
892	phaseVary += [name,]
893
894	def MakeRBThermals(rbKey,phaseVary,phaseDict):
895	rbid = str(rbids.index(RB['RBId']))
896	tlstr = ['11','22','33','12','13','23']
897	sstr = ['12','13','21','23','31','32','AA','BB']
898	if 'T' in RB['ThermalMotion'][0]:
899	pfxRB = pfx+'RB'+rbKey+'T'
900	for i in range(6):
901	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
902	phaseDict[name] = RB['ThermalMotion'][1][i]
903	if RB['ThermalMotion'][2][i]:
904	phaseVary += [name,]
905	if 'L' in RB['ThermalMotion'][0]:
906	pfxRB = pfx+'RB'+rbKey+'L'
907	for i in range(6):
908	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
909	phaseDict[name] = RB['ThermalMotion'][1][i+6]
910	if RB['ThermalMotion'][2][i+6]:
911	phaseVary += [name,]
912	if 'S' in RB['ThermalMotion'][0]:
913	pfxRB = pfx+'RB'+rbKey+'S'
914	for i in range(8):
915	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
916	phaseDict[name] = RB['ThermalMotion'][1][i+12]
917	if RB['ThermalMotion'][2][i+12]:
918	phaseVary += [name,]
919	if 'U' in RB['ThermalMotion'][0]:
920	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
921	phaseDict[name] = RB['ThermalMotion'][1][0]
922	if RB['ThermalMotion'][2][0]:
923	phaseVary += [name,]
924
925	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
926	rbid = str(rbids.index(RB['RBId']))
927	pfxRB = pfx+'RB'+rbKey+'Tr;'
928	for i,torsion in enumerate(RB['Torsions']):
929	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
930	phaseDict[name] = torsion[0]
931	if torsion[1]:
932	phaseVary += [name,]
933
934	if Print:
935	print >>pFile,'\n Phases:'
936	phaseVary = []
937	phaseDict = {}
938	phaseConstr = {}
939	pawleyLookup = {}
940	FFtables = {} #scattering factors - xrays
941	BLtables = {} # neutrons
942	Natoms = {}
943	AtMults = {}
944	AtIA = {}
945	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
946	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
947	atomIndx = {}
948	for name in PhaseData:
949	General = PhaseData[name]['General']
950	pId = PhaseData[name]['pId']
951	pfx = str(pId)+'::'
952	FFtable = G2el.GetFFtable(General['AtomTypes'])
953	BLtable = G2el.GetBLtable(General)
954	FFtables.update(FFtable)
955	BLtables.update(BLtable)
956	Atoms = PhaseData[name]['Atoms']
957	AtLookup = G2mth.FillAtomLookUp(Atoms)
958	PawleyRef = PhaseData[name].get('Pawley ref',[])
959	SGData = General['SGData']
960	SGtext = G2spc.SGPrint(SGData)
961	cell = General['Cell']
962	A = G2lat.cell2A(cell[1:7])
963	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
964	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
965	if cell[0]:
966	phaseVary += cellVary(pfx,SGData)
967	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
968	if resRBData:
969	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
970	for iRB,RB in enumerate(resRBData):
971	MakeRBParms('R',phaseVary,phaseDict)
972	MakeRBThermals('R',phaseVary,phaseDict)
973	MakeRBTorsions('R',phaseVary,phaseDict)
974
975	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
976	if vecRBData:
977	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
978	for iRB,RB in enumerate(vecRBData):
979	MakeRBParms('V',phaseVary,phaseDict)
980	MakeRBThermals('V',phaseVary,phaseDict)
981
982	Natoms[pfx] = 0
983	if Atoms and not General.get('doPawley'):
984	cx,ct,cs,cia = General['AtomPtrs']
985	if General['Type'] in ['nuclear','macromolecular']:
986	Natoms[pfx] = len(Atoms)
987	for i,at in enumerate(Atoms):
988	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
989	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
990	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
991	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
992	pfx+'AI/A:'+str(i):at[cia],})
993	if at[cia] == 'I':
994	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
995	else:
996	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
997	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
998	if 'F' in at[ct+1]:
999	phaseVary.append(pfx+'Afrac:'+str(i))
1000	if 'X' in at[ct+1]:
1001	try: #patch for sytsym name changes
1002	xId,xCoef = G2spc.GetCSxinel(at[cs])
1003	except KeyError:
1004	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1005	at[cs] = Sytsym
1006	xId,xCoef = G2spc.GetCSxinel(at[cs])
1007	xId,xCoef = G2spc.GetCSxinel(at[cs])
1008	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1009	equivs = [[],[],[]]
1010	for j in range(3):
1011	if xId[j] > 0:
1012	phaseVary.append(names[j])
1013	equivs[xId[j]-1].append([names[j],xCoef[j]])
1014	for equiv in equivs:
1015	if len(equiv) > 1:
1016	name = equiv[0][0]
1017	coef = equiv[0][1]
1018	for eqv in equiv[1:]:
1019	eqv[1] /= coef
1020	G2mv.StoreEquivalence(name,(eqv,))
1021	if 'U' in at[ct+1]:
1022	if at[cia] == 'I':
1023	phaseVary.append(pfx+'AUiso:'+str(i))
1024	else:
1025	try: #patch for sytsym name changes
1026	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1027	except KeyError:
1028	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1029	at[cs] = Sytsym
1030	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1031	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1032	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1033	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1034	equivs = [[],[],[],[],[],[]]
1035	for j in range(6):
1036	if uId[j] > 0:
1037	phaseVary.append(names[j])
1038	equivs[uId[j]-1].append([names[j],uCoef[j]])
1039	for equiv in equivs:
1040	if len(equiv) > 1:
1041	name = equiv[0][0]
1042	coef = equiv[0][1]
1043	for eqv in equiv[1:]:
1044	eqv[1] /= coef
1045	G2mv.StoreEquivalence(name,equiv[1:])
1046	# elif General['Type'] == 'magnetic':
1047	# elif General['Type'] == 'macromolecular':
1048	textureData = General['SH Texture']
1049	if textureData['Order']:
1050	phaseDict[pfx+'SHorder'] = textureData['Order']
1051	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1052	for item in ['omega','chi','phi']:
1053	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1054	if textureData['Sample '+item][0]:
1055	phaseVary.append(pfx+'SH '+item)
1056	for item in textureData['SH Coeff'][1]:
1057	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1058	if textureData['SH Coeff'][0]:
1059	phaseVary.append(pfx+item)
1060
1061	if Print:
1062	print >>pFile,'\n Phase name: ',General['Name']
1063	print >>pFile,135*'-'
1064	PrintFFtable(FFtable)
1065	PrintBLtable(BLtable)
1066	print >>pFile,''
1067	for line in SGtext: print >>pFile,line
1068	PrintRBObjects(resRBData,vecRBData)
1069	PrintAtoms(General,Atoms)
1070	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1071	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1072	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1073	PrintTexture(textureData)
1074	if name in RestraintDict:
1075	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1076	textureData,RestraintDict[name],pFile)
1077
1078	elif PawleyRef:
1079	pawleyVary = []
1080	for i,refl in enumerate(PawleyRef):
1081	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1082	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1083	if refl[5]:
1084	pawleyVary.append(pfx+'PWLref:'+str(i))
1085	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
1086	phaseVary += pawleyVary
1087
1088	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1089
1090	def cellFill(pfx,SGData,parmDict,sigDict):
1091	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1092	from the parameter and sig dictionaries.
1093
1094	:param str pfx: parameter prefix ("n::", where n is a phase number)
1095	:param dict SGdata: a symmetry object
1096	:param dict parmDict: a dictionary of parameters
1097	:param dict sigDict: a dictionary of uncertainties on parameters
1098
1099	:returns: A,sigA where each is a list of six terms with the A terms
1100	'''
1101	if SGData['SGLaue'] in ['-1',]:
1102	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1103	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1104	elif SGData['SGLaue'] in ['2/m',]:
1105	if SGData['SGUniq'] == 'a':
1106	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1107	parmDict[pfx+'A3'],0,0]
1108	elif SGData['SGUniq'] == 'b':
1109	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1110	0,parmDict[pfx+'A4'],0]
1111	else:
1112	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1113	0,0,parmDict[pfx+'A5']]
1114	elif SGData['SGLaue'] in ['mmm',]:
1115	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1116	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1117	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1118	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1119	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1120	parmDict[pfx+'A0'],0,0]
1121	elif SGData['SGLaue'] in ['3R', '3mR']:
1122	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1123	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1124	elif SGData['SGLaue'] in ['m3m','m3']:
1125	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1126
1127	try:
1128	if SGData['SGLaue'] in ['-1',]:
1129	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1130	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1131	elif SGData['SGLaue'] in ['2/m',]:
1132	if SGData['SGUniq'] == 'a':
1133	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1134	sigDict[pfx+'A3'],0,0]
1135	elif SGData['SGUniq'] == 'b':
1136	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1137	0,sigDict[pfx+'A4'],0]
1138	else:
1139	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1140	0,0,sigDict[pfx+'A5']]
1141	elif SGData['SGLaue'] in ['mmm',]:
1142	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1143	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1144	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1145	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1146	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1147	elif SGData['SGLaue'] in ['3R', '3mR']:
1148	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1149	elif SGData['SGLaue'] in ['m3m','m3']:
1150	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1151	except KeyError:
1152	sigA = [0,0,0,0,0,0]
1153
1154	return A,sigA
1155
1156	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1157	'needs a doc string'
1158	if phaseRest:
1159	Amat = G2lat.cell2AB(cell)[0]
1160	cx,ct,cs = AtPtrs[:3]
1161	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1162	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1163	['ChemComp','Sites'],['Texture','HKLs']]
1164	for name,rest in names:
1165	itemRest = phaseRest[name]
1166	if itemRest[rest] and itemRest['Use']:
1167	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1168	if name in ['Bond','Angle','Plane','Chiral']:
1169	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1170	for indx,ops,obs,esd in itemRest[rest]:
1171	try:
1172	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1173	AtName = ''
1174	for i,Aname in enumerate(AtNames):
1175	AtName += Aname
1176	if ops[i] == '1':
1177	AtName += '-'
1178	else:
1179	AtName += '+('+ops[i]+')-'
1180	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1181	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1182	if name == 'Bond':
1183	calc = G2mth.getRestDist(XYZ,Amat)
1184	elif name == 'Angle':
1185	calc = G2mth.getRestAngle(XYZ,Amat)
1186	elif name == 'Plane':
1187	calc = G2mth.getRestPlane(XYZ,Amat)
1188	elif name == 'Chiral':
1189	calc = G2mth.getRestChiral(XYZ,Amat)
1190	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1191	except KeyError:
1192	del itemRest[rest]
1193	elif name in ['Torsion','Rama']:
1194	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1195	coeffDict = itemRest['Coeff']
1196	for indx,ops,cofName,esd in itemRest[rest]:
1197	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1198	AtName = ''
1199	for i,Aname in enumerate(AtNames):
1200	AtName += Aname+'+('+ops[i]+')-'
1201	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1202	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1203	if name == 'Torsion':
1204	tor = G2mth.getRestTorsion(XYZ,Amat)
1205	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1206	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1207	else:
1208	phi,psi = G2mth.getRestRama(XYZ,Amat)
1209	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1210	print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
1211	elif name == 'ChemComp':
1212	print >>pFile,' atoms mul*frac factor prod'
1213	for indx,factors,obs,esd in itemRest[rest]:
1214	try:
1215	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1216	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1217	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1218	mulfrac = mul*frac
1219	calcs = mulfracfactors
1220	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1221	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1222	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1223	except KeyError:
1224	del itemRest[rest]
1225	elif name == 'Texture' and textureData['Order']:
1226	Start = False
1227	SHCoef = textureData['SH Coeff'][1]
1228	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1229	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1230	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1231	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1232	PH = np.array(hkl)
1233	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1234	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1235	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1236	Z = ma.masked_greater(Z,0.0)
1237	num = ma.count(Z)
1238	sum = 0
1239	if num:
1240	sum = np.sum(Z)
1241	print ' %d %d %d %d %8.3f %8.3f %8d %s %8.3f'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)
1242
1243	def getCellEsd(pfx,SGData,A,covData):
1244	'needs a doc string'
1245	dpr = 180./np.pi
1246	rVsq = G2lat.calc_rVsq(A)
1247	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1248	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1249	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1250	scos = cosd(cellst[3:6])
1251	ssin = sind(cellst[3:6])
1252	scot = scos/ssin
1253	rcos = cosd(cell[3:6])
1254	rsin = sind(cell[3:6])
1255	rcot = rcos/rsin
1256	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1257	varyList = covData['varyList']
1258	covMatrix = covData['covMatrix']
1259	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1260	Ax = np.array(A)
1261	Ax[3:] /= 2.
1262	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1263	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1264	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1265	Vol = 1/np.sqrt(rVsq)
1266	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1267	R123 = Ax[0]Ax[1]Ax[2]
1268	dsasdg = np.zeros((3,6))
1269	dadg = np.zeros((6,6))
1270	for i0 in range(3): #0 1 2
1271	i1 = (i0+1)%3 #1 2 0
1272	i2 = (i1+1)%3 #2 0 1
1273	i3 = 5-i2 #3 5 4
1274	i4 = 5-i1 #4 3 5
1275	i5 = 5-i0 #5 4 3
1276	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1277	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1278	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1279	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1280	denom = np.sqrt(denmsq)
1281	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1282	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1283	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1284	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1285	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1286	dadg[i5][i5] = -Ax[i0]/denom
1287	for i0 in range(3):
1288	i1 = (i0+1)%3
1289	i2 = (i1+1)%3
1290	i3 = 5-i2
1291	for ij in range(6):
1292	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1293	if ij == i0:
1294	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1295	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1296	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1297	var = np.diag(sigMat)
1298	CS = np.where(var>0.,np.sqrt(var),0.)
1299	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1300
1301	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1302	'needs a doc string'
1303
1304	def PrintAtomsAndSig(General,Atoms,atomsSig):
1305	print >>pFile,'\n Atoms:'
1306	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1307	if General['Type'] == 'magnetic':
1308	line += ' Mx My Mz'
1309	elif General['Type'] == 'macromolecular':
1310	line = ' res no residue chain '+line
1311	cx,ct,cs,cia = General['AtomPtrs']
1312	print >>pFile,line
1313	if General['Type'] == 'nuclear':
1314	print >>pFile,135*'-'
1315	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1316	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1317	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1318	for atyp in General['AtomTypes']: #zero composition
1319	General['NoAtoms'][atyp] = 0.
1320	for i,at in enumerate(Atoms):
1321	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1322	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1323	valstr = ' values:'
1324	sigstr = ' sig :'
1325	for ind in [3,4,5,6]:
1326	sigind = str(i)+':'+str(ind)
1327	valstr += fmt[ind]%(at[ind])
1328	if sigind in atomsSig:
1329	sigstr += fmt[ind]%(atomsSig[sigind])
1330	else:
1331	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1332	if at[9] == 'I':
1333	valstr += fmt[10]%(at[10])
1334	if str(i)+':10' in atomsSig:
1335	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1336	else:
1337	sigstr += 8*' '
1338	else:
1339	valstr += 8*' '
1340	sigstr += 8*' '
1341	for ind in [11,12,13,14,15,16]:
1342	sigind = str(i)+':'+str(ind)
1343	valstr += fmt[ind]%(at[ind])
1344	if sigind in atomsSig:
1345	sigstr += fmt[ind]%(atomsSig[sigind])
1346	else:
1347	sigstr += 8*' '
1348	print >>pFile,name
1349	print >>pFile,valstr
1350	print >>pFile,sigstr
1351
1352	def PrintRBObjPOAndSig(rbfx,rbsx):
1353	namstr = ' names :'
1354	valstr = ' values:'
1355	sigstr = ' esds :'
1356	for i,px in enumerate(['Px:','Py:','Pz:']):
1357	name = pfx+rbfx+px+rbsx
1358	namstr += '%12s'%('Pos '+px[1])
1359	valstr += '%12.5f'%(parmDict[name])
1360	if name in sigDict:
1361	sigstr += '%12.5f'%(sigDict[name])
1362	else:
1363	sigstr += 12*' '
1364	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1365	name = pfx+rbfx+po+rbsx
1366	namstr += '%12s'%('Ori '+po[1])
1367	valstr += '%12.5f'%(parmDict[name])
1368	if name in sigDict:
1369	sigstr += '%12.5f'%(sigDict[name])
1370	else:
1371	sigstr += 12*' '
1372	print >>pFile,namstr
1373	print >>pFile,valstr
1374	print >>pFile,sigstr
1375
1376	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1377	namstr = ' names :'
1378	valstr = ' values:'
1379	sigstr = ' esds :'
1380	if 'N' not in TLS:
1381	print >>pFile,' Thermal motion:'
1382	if 'T' in TLS:
1383	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1384	name = pfx+rbfx+pt+rbsx
1385	namstr += '%12s'%(pt[:3])
1386	valstr += '%12.5f'%(parmDict[name])
1387	if name in sigDict:
1388	sigstr += '%12.5f'%(sigDict[name])
1389	else:
1390	sigstr += 12*' '
1391	print >>pFile,namstr
1392	print >>pFile,valstr
1393	print >>pFile,sigstr
1394	if 'L' in TLS:
1395	namstr = ' names :'
1396	valstr = ' values:'
1397	sigstr = ' esds :'
1398	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1399	name = pfx+rbfx+pt+rbsx
1400	namstr += '%12s'%(pt[:3])
1401	valstr += '%12.3f'%(parmDict[name])
1402	if name in sigDict:
1403	sigstr += '%12.3f'%(sigDict[name])
1404	else:
1405	sigstr += 12*' '
1406	print >>pFile,namstr
1407	print >>pFile,valstr
1408	print >>pFile,sigstr
1409	if 'S' in TLS:
1410	namstr = ' names :'
1411	valstr = ' values:'
1412	sigstr = ' esds :'
1413	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1414	name = pfx+rbfx+pt+rbsx
1415	namstr += '%12s'%(pt[:3])
1416	valstr += '%12.4f'%(parmDict[name])
1417	if name in sigDict:
1418	sigstr += '%12.4f'%(sigDict[name])
1419	else:
1420	sigstr += 12*' '
1421	print >>pFile,namstr
1422	print >>pFile,valstr
1423	print >>pFile,sigstr
1424	if 'U' in TLS:
1425	name = pfx+rbfx+'U:'+rbsx
1426	namstr = ' names :'+'%12s'%('Uiso')
1427	valstr = ' values:'+'%12.5f'%(parmDict[name])
1428	if name in sigDict:
1429	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1430	else:
1431	sigstr = ' esds :'+12*' '
1432	print >>pFile,namstr
1433	print >>pFile,valstr
1434	print >>pFile,sigstr
1435
1436	def PrintRBObjTorAndSig(rbsx):
1437	namstr = ' names :'
1438	valstr = ' values:'
1439	sigstr = ' esds :'
1440	nTors = len(RBObj['Torsions'])
1441	if nTors:
1442	print >>pFile,' Torsions:'
1443	for it in range(nTors):
1444	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1445	namstr += '%12s'%('Tor'+str(it))
1446	valstr += '%12.4f'%(parmDict[name])
1447	if name in sigDict:
1448	sigstr += '%12.4f'%(sigDict[name])
1449	print >>pFile,namstr
1450	print >>pFile,valstr
1451	print >>pFile,sigstr
1452
1453	def PrintSHtextureAndSig(textureData,SHtextureSig):
1454	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1455	names = ['omega','chi','phi']
1456	namstr = ' names :'
1457	ptstr = ' values:'
1458	sigstr = ' esds :'
1459	for name in names:
1460	namstr += '%12s'%(name)
1461	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1462	if 'Sample '+name in SHtextureSig:
1463	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1464	else:
1465	sigstr += 12*' '
1466	print >>pFile,namstr
1467	print >>pFile,ptstr
1468	print >>pFile,sigstr
1469	print >>pFile,'\n Texture coefficients:'
1470	namstr = ' names :'
1471	ptstr = ' values:'
1472	sigstr = ' esds :'
1473	SHcoeff = textureData['SH Coeff'][1]
1474	for name in SHcoeff:
1475	namstr += '%12s'%(name)
1476	ptstr += '%12.3f'%(SHcoeff[name])
1477	if name in SHtextureSig:
1478	sigstr += '%12.3f'%(SHtextureSig[name])
1479	else:
1480	sigstr += 12*' '
1481	print >>pFile,namstr
1482	print >>pFile,ptstr
1483	print >>pFile,sigstr
1484
1485	print >>pFile,'\n Phases:'
1486	for phase in Phases:
1487	print >>pFile,' Result for phase: ',phase
1488	Phase = Phases[phase]
1489	General = Phase['General']
1490	SGData = General['SGData']
1491	Atoms = Phase['Atoms']
1492	AtLookup = G2mth.FillAtomLookUp(Atoms)
1493	cell = General['Cell']
1494	pId = Phase['pId']
1495	pfx = str(pId)+'::'
1496	if cell[0]:
1497	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1498	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1499	print >>pFile,' Reciprocal metric tensor: '
1500	ptfmt = "%15.9f"
1501	names = ['A11','A22','A33','A12','A13','A23']
1502	namstr = ' names :'
1503	ptstr = ' values:'
1504	sigstr = ' esds :'
1505	for name,a,siga in zip(names,A,sigA):
1506	namstr += '%15s'%(name)
1507	ptstr += ptfmt%(a)
1508	if siga:
1509	sigstr += ptfmt%(siga)
1510	else:
1511	sigstr += 15*' '
1512	print >>pFile,namstr
1513	print >>pFile,ptstr
1514	print >>pFile,sigstr
1515	cell[1:7] = G2lat.A2cell(A)
1516	cell[7] = G2lat.calc_V(A)
1517	print >>pFile,' New unit cell:'
1518	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1519	names = ['a','b','c','alpha','beta','gamma','Volume']
1520	namstr = ' names :'
1521	ptstr = ' values:'
1522	sigstr = ' esds :'
1523	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1524	namstr += '%12s'%(name)
1525	ptstr += fmt%(a)
1526	if siga:
1527	sigstr += fmt%(siga)
1528	else:
1529	sigstr += 12*' '
1530	print >>pFile,namstr
1531	print >>pFile,ptstr
1532	print >>pFile,sigstr
1533
1534	General['Mass'] = 0.
1535	if Phase['General'].get('doPawley'):
1536	pawleyRef = Phase['Pawley ref']
1537	for i,refl in enumerate(pawleyRef):
1538	key = pfx+'PWLref:'+str(i)
1539	refl[6] = parmDict[key]
1540	if key in sigDict:
1541	refl[7] = sigDict[key]
1542	else:
1543	refl[7] = 0
1544	else:
1545	VRBIds = RBIds['Vector']
1546	RRBIds = RBIds['Residue']
1547	RBModels = Phase['RBModels']
1548	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1549	jrb = VRBIds.index(RBObj['RBId'])
1550	rbsx = str(irb)+':'+str(jrb)
1551	print >>pFile,' Vector rigid body parameters:'
1552	PrintRBObjPOAndSig('RBV',rbsx)
1553	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1554	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1555	jrb = RRBIds.index(RBObj['RBId'])
1556	rbsx = str(irb)+':'+str(jrb)
1557	print >>pFile,' Residue rigid body parameters:'
1558	PrintRBObjPOAndSig('RBR',rbsx)
1559	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1560	PrintRBObjTorAndSig(rbsx)
1561	atomsSig = {}
1562	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1563	for i,at in enumerate(Atoms):
1564	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1565	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1566	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1567	for ind in [3,4,5,6]:
1568	at[ind] = parmDict[names[ind]]
1569	if ind in [3,4,5]:
1570	name = names[ind].replace('A','dA')
1571	else:
1572	name = names[ind]
1573	if name in sigDict:
1574	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1575	if at[9] == 'I':
1576	at[10] = parmDict[names[10]]
1577	if names[10] in sigDict:
1578	atomsSig[str(i)+':10'] = sigDict[names[10]]
1579	else:
1580	for ind in [11,12,13,14,15,16]:
1581	at[ind] = parmDict[names[ind]]
1582	if names[ind] in sigDict:
1583	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1584	ind = General['AtomTypes'].index(at[1])
1585	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1586	PrintAtomsAndSig(General,Atoms,atomsSig)
1587
1588	textureData = General['SH Texture']
1589	if textureData['Order']:
1590	SHtextureSig = {}
1591	for name in ['omega','chi','phi']:
1592	aname = pfx+'SH '+name
1593	textureData['Sample '+name][1] = parmDict[aname]
1594	if aname in sigDict:
1595	SHtextureSig['Sample '+name] = sigDict[aname]
1596	for name in textureData['SH Coeff'][1]:
1597	aname = pfx+name
1598	textureData['SH Coeff'][1][name] = parmDict[aname]
1599	if aname in sigDict:
1600	SHtextureSig[name] = sigDict[aname]
1601	PrintSHtextureAndSig(textureData,SHtextureSig)
1602	if phase in RestraintDict:
1603	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1604	textureData,RestraintDict[phase],pFile)
1605
1606	################################################################################
1607	##### Histogram & Phase data
1608	################################################################################
1609
1610	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1611	'''Loads the HAP histogram/phase information into dicts
1612
1613	:param dict Phases: phase information
1614	:param dict Histograms: Histogram information
1615	:param bool Print: prints information as it is read
1616	:param file pFile: file object to print to (the default, None causes printing to the console)
1617	:param bool resetRefList: Should the contents of the Reflection List be initialized
1618	on loading. The default, True, initializes the Reflection List as it is loaded.
1619
1620	:returns: (hapVary,hapDict,controlDict)
1621	* hapVary: list of refined variables
1622	* hapDict: dict with refined variables and their values
1623	* controlDict: dict with computation controls (?)
1624	'''
1625
1626	def PrintSize(hapData):
1627	if hapData[0] in ['isotropic','uniaxial']:
1628	line = '\n Size model : %9s'%(hapData[0])
1629	line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1630	if hapData[0] == 'uniaxial':
1631	line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1632	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1633	print >>pFile,line
1634	else:
1635	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1636	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1637	Snames = ['S11','S22','S33','S12','S13','S23']
1638	ptlbls = ' names :'
1639	ptstr = ' values:'
1640	varstr = ' refine:'
1641	for i,name in enumerate(Snames):
1642	ptlbls += '%12s' % (name)
1643	ptstr += '%12.6f' % (hapData[4][i])
1644	varstr += '%12s' % (str(hapData[5][i]))
1645	print >>pFile,ptlbls
1646	print >>pFile,ptstr
1647	print >>pFile,varstr
1648
1649	def PrintMuStrain(hapData,SGData):
1650	if hapData[0] in ['isotropic','uniaxial']:
1651	line = '\n Mustrain model: %9s'%(hapData[0])
1652	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1653	if hapData[0] == 'uniaxial':
1654	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1655	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1656	print >>pFile,line
1657	else:
1658	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1659	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1660	Snames = G2spc.MustrainNames(SGData)
1661	ptlbls = ' names :'
1662	ptstr = ' values:'
1663	varstr = ' refine:'
1664	for i,name in enumerate(Snames):
1665	ptlbls += '%12s' % (name)
1666	ptstr += '%12.6f' % (hapData[4][i])
1667	varstr += '%12s' % (str(hapData[5][i]))
1668	print >>pFile,ptlbls
1669	print >>pFile,ptstr
1670	print >>pFile,varstr
1671
1672	def PrintHStrain(hapData,SGData):
1673	print >>pFile,'\n Hydrostatic/elastic strain: '
1674	Hsnames = G2spc.HStrainNames(SGData)
1675	ptlbls = ' names :'
1676	ptstr = ' values:'
1677	varstr = ' refine:'
1678	for i,name in enumerate(Hsnames):
1679	ptlbls += '%12s' % (name)
1680	ptstr += '%12.4g' % (hapData[0][i])
1681	varstr += '%12s' % (str(hapData[1][i]))
1682	print >>pFile,ptlbls
1683	print >>pFile,ptstr
1684	print >>pFile,varstr
1685
1686	def PrintSHPO(hapData):
1687	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1688	str(hapData[4])+' Refine? '+str(hapData[2])
1689	ptlbls = ' names :'
1690	ptstr = ' values:'
1691	for item in hapData[5]:
1692	ptlbls += '%12s'%(item)
1693	ptstr += '%12.3f'%(hapData[5][item])
1694	print >>pFile,ptlbls
1695	print >>pFile,ptstr
1696
1697	def PrintBabinet(hapData):
1698	print >>pFile,'\n Babinet form factor modification: '
1699	ptlbls = ' names :'
1700	ptstr = ' values:'
1701	varstr = ' refine:'
1702	for name in ['BabA','BabU']:
1703	ptlbls += '%12s' % (name)
1704	ptstr += '%12.6f' % (hapData[name][0])
1705	varstr += '%12s' % (str(hapData[name][1]))
1706	print >>pFile,ptlbls
1707	print >>pFile,ptstr
1708	print >>pFile,varstr
1709
1710
1711	hapDict = {}
1712	hapVary = []
1713	controlDict = {}
1714	poType = {}
1715	poAxes = {}
1716	spAxes = {}
1717	spType = {}
1718
1719	for phase in Phases:
1720	HistoPhase = Phases[phase]['Histograms']
1721	SGData = Phases[phase]['General']['SGData']
1722	cell = Phases[phase]['General']['Cell'][1:7]
1723	A = G2lat.cell2A(cell)
1724	pId = Phases[phase]['pId']
1725	histoList = HistoPhase.keys()
1726	histoList.sort()
1727	for histogram in histoList:
1728	try:
1729	Histogram = Histograms[histogram]
1730	except KeyError:
1731	#skip if histogram not included e.g. in a sequential refinement
1732	continue
1733	hapData = HistoPhase[histogram]
1734	hId = Histogram['hId']
1735	if 'PWDR' in histogram:
1736	limits = Histogram['Limits'][1]
1737	inst = Histogram['Instrument Parameters'][0]
1738	Zero = inst['Zero'][0]
1739	if 'C' in inst['Type'][1]:
1740	try:
1741	wave = inst['Lam'][1]
1742	except KeyError:
1743	wave = inst['Lam1'][1]
1744	dmin = wave/(2.0*sind(limits[1]/2.0))
1745	elif 'T' in inst['Type'][0]:
1746	dmin = limits[0]/inst['difC'][1]
1747	pfx = str(pId)+':'+str(hId)+':'
1748	for item in ['Scale','Extinction']:
1749	hapDict[pfx+item] = hapData[item][0]
1750	if hapData[item][1]:
1751	hapVary.append(pfx+item)
1752	names = G2spc.HStrainNames(SGData)
1753	for i,name in enumerate(names):
1754	hapDict[pfx+name] = hapData['HStrain'][0][i]
1755	if hapData['HStrain'][1][i]:
1756	hapVary.append(pfx+name)
1757	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1758	if hapData['Pref.Ori.'][0] == 'MD':
1759	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1760	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1761	if hapData['Pref.Ori.'][2]:
1762	hapVary.append(pfx+'MD')
1763	else: #'SH' spherical harmonics
1764	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1765	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1766	for item in hapData['Pref.Ori.'][5]:
1767	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1768	if hapData['Pref.Ori.'][2]:
1769	hapVary.append(pfx+item)
1770	for item in ['Mustrain','Size']:
1771	controlDict[pfx+item+'Type'] = hapData[item][0]
1772	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1773	if hapData[item][2][2]:
1774	hapVary.append(pfx+item+';mx')
1775	if hapData[item][0] in ['isotropic','uniaxial']:
1776	hapDict[pfx+item+';i'] = hapData[item][1][0]
1777	if hapData[item][2][0]:
1778	hapVary.append(pfx+item+';i')
1779	if hapData[item][0] == 'uniaxial':
1780	controlDict[pfx+item+'Axis'] = hapData[item][3]
1781	hapDict[pfx+item+';a'] = hapData[item][1][1]
1782	if hapData[item][2][1]:
1783	hapVary.append(pfx+item+';a')
1784	else: #generalized for mustrain or ellipsoidal for size
1785	Nterms = len(hapData[item][4])
1786	if item == 'Mustrain':
1787	names = G2spc.MustrainNames(SGData)
1788	pwrs = []
1789	for name in names:
1790	h,k,l = name[1:]
1791	pwrs.append([int(h),int(k),int(l)])
1792	controlDict[pfx+'MuPwrs'] = pwrs
1793	for i in range(Nterms):
1794	sfx = ':'+str(i)
1795	hapDict[pfx+item+sfx] = hapData[item][4][i]
1796	if hapData[item][5][i]:
1797	hapVary.append(pfx+item+sfx)
1798	for bab in ['BabA','BabU']:
1799	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1800	if hapData['Babinet'][bab][1]:
1801	hapVary.append(pfx+bab)
1802
1803	if Print:
1804	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1805	print >>pFile,135*'-'
1806	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1807	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1808	if hapData['Pref.Ori.'][0] == 'MD':
1809	Ax = hapData['Pref.Ori.'][3]
1810	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1811	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1812	else: #'SH' for spherical harmonics
1813	PrintSHPO(hapData['Pref.Ori.'])
1814	PrintSize(hapData['Size'])
1815	PrintMuStrain(hapData['Mustrain'],SGData)
1816	PrintHStrain(hapData['HStrain'],SGData)
1817	if hapData['Babinet']['BabA'][0]:
1818	PrintBabinet(hapData['Babinet'])
1819	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1820	if resetRefList:
1821	refList = []
1822	Uniq = []
1823	Phi = []
1824	for h,k,l,d in HKLd:
1825	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1826	mul *= 2 # for powder overlap of Friedel pairs
1827	if ext:
1828	continue
1829	if 'C' in inst['Type'][0]:
1830	pos = 2.0asind(wave/(2.0d))+Zero
1831	if limits[0] < pos < limits[1]:
1832	refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
1833	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1834	Uniq.append(uniq)
1835	Phi.append(phi)
1836	elif 'T' in inst['Type'][0]:
1837	pos = inst['difC'][1]d+inst['difA'][1]d**2+inst['difB'][1]/d+Zero
1838	if limits[0] < pos < limits[1]:
1839	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
1840	refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
1841	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1842	Uniq.append(uniq)
1843	Phi.append(phi)
1844	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0]}
1845	elif 'HKLF' in histogram:
1846	inst = Histogram['Instrument Parameters'][0]
1847	hId = Histogram['hId']
1848	hfx = ':%d:'%(hId)
1849	for item in inst:
1850	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1851	hapDict[hfx+item] = inst[item][1]
1852	pfx = str(pId)+':'+str(hId)+':'
1853	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1854	if hapData['Scale'][1]:
1855	hapVary.append(pfx+'Scale')
1856
1857	extApprox,extType,extParms = hapData['Extinction']
1858	controlDict[pfx+'EType'] = extType
1859	controlDict[pfx+'EApprox'] = extApprox
1860	if 'C' in inst['Type'][0]:
1861	controlDict[pfx+'Tbar'] = extParms['Tbar']
1862	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1863	if 'Primary' in extType:
1864	Ekey = ['Ep',]
1865	elif 'I & II' in extType:
1866	Ekey = ['Eg','Es']
1867	elif 'Secondary Type II' == extType:
1868	Ekey = ['Es',]
1869	elif 'Secondary Type I' == extType:
1870	Ekey = ['Eg',]
1871	else: #'None'
1872	Ekey = []
1873	for eKey in Ekey:
1874	hapDict[pfx+eKey] = extParms[eKey][0]
1875	if extParms[eKey][1]:
1876	hapVary.append(pfx+eKey)
1877	for bab in ['BabA','BabU']:
1878	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1879	if hapData['Babinet'][bab][1]:
1880	hapVary.append(pfx+bab)
1881	if Print:
1882	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1883	print >>pFile,135*'-'
1884	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1885	if extType != 'None':
1886	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox)
1887	text = ' Parameters :'
1888	for eKey in Ekey:
1889	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1890	print >>pFile,text
1891	if hapData['Babinet']['BabA'][0]:
1892	PrintBabinet(hapData['Babinet'])
1893	Histogram['Reflection Lists'] = phase
1894
1895	return hapVary,hapDict,controlDict
1896
1897	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1898	'needs a doc string'
1899
1900	def PrintSizeAndSig(hapData,sizeSig):
1901	line = '\n Size model: %9s'%(hapData[0])
1902	refine = False
1903	if hapData[0] in ['isotropic','uniaxial']:
1904	line += ' equatorial:%12.5f'%(hapData[1][0])
1905	if sizeSig[0][0]:
1906	line += ', sig:%8.4f'%(sizeSig[0][0])
1907	refine = True
1908	if hapData[0] == 'uniaxial':
1909	line += ' axial:%12.4f'%(hapData[1][1])
1910	if sizeSig[0][1]:
1911	refine = True
1912	line += ', sig:%8.4f'%(sizeSig[0][1])
1913	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1914	if sizeSig[0][2]:
1915	refine = True
1916	line += ', sig:%8.4f'%(sizeSig[0][2])
1917	if refine:
1918	print >>pFile,line
1919	else:
1920	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1921	if sizeSig[0][2]:
1922	refine = True
1923	line += ', sig:%8.4f'%(sizeSig[0][2])
1924	Snames = ['S11','S22','S33','S12','S13','S23']
1925	ptlbls = ' name :'
1926	ptstr = ' value :'
1927	sigstr = ' sig :'
1928	for i,name in enumerate(Snames):
1929	ptlbls += '%12s' % (name)
1930	ptstr += '%12.6f' % (hapData[4][i])
1931	if sizeSig[1][i]:
1932	refine = True
1933	sigstr += '%12.6f' % (sizeSig[1][i])
1934	else:
1935	sigstr += 12*' '
1936	if refine:
1937	print >>pFile,line
1938	print >>pFile,ptlbls
1939	print >>pFile,ptstr
1940	print >>pFile,sigstr
1941
1942	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1943	line = '\n Mustrain model: %9s'%(hapData[0])
1944	refine = False
1945	if hapData[0] in ['isotropic','uniaxial']:
1946	line += ' equatorial:%12.1f'%(hapData[1][0])
1947	if mustrainSig[0][0]:
1948	line += ', sig:%8.1f'%(mustrainSig[0][0])
1949	refine = True
1950	if hapData[0] == 'uniaxial':
1951	line += ' axial:%12.1f'%(hapData[1][1])
1952	if mustrainSig[0][1]:
1953	line += ', sig:%8.1f'%(mustrainSig[0][1])
1954	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1955	if mustrainSig[0][2]:
1956	refine = True
1957	line += ', sig:%8.3f'%(mustrainSig[0][2])
1958	if refine:
1959	print >>pFile,line
1960	else:
1961	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1962	if mustrainSig[0][2]:
1963	refine = True
1964	line += ', sig:%8.3f'%(mustrainSig[0][2])
1965	Snames = G2spc.MustrainNames(SGData)
1966	ptlbls = ' name :'
1967	ptstr = ' value :'
1968	sigstr = ' sig :'
1969	for i,name in enumerate(Snames):
1970	ptlbls += '%12s' % (name)
1971	ptstr += '%12.6f' % (hapData[4][i])
1972	if mustrainSig[1][i]:
1973	refine = True
1974	sigstr += '%12.6f' % (mustrainSig[1][i])
1975	else:
1976	sigstr += 12*' '
1977	if refine:
1978	print >>pFile,line
1979	print >>pFile,ptlbls
1980	print >>pFile,ptstr
1981	print >>pFile,sigstr
1982
1983	def PrintHStrainAndSig(hapData,strainSig,SGData):
1984	Hsnames = G2spc.HStrainNames(SGData)
1985	ptlbls = ' name :'
1986	ptstr = ' value :'
1987	sigstr = ' sig :'
1988	refine = False
1989	for i,name in enumerate(Hsnames):
1990	ptlbls += '%12s' % (name)
1991	ptstr += '%12.4g' % (hapData[0][i])
1992	if name in strainSig:
1993	refine = True
1994	sigstr += '%12.4g' % (strainSig[name])
1995	else:
1996	sigstr += 12*' '
1997	if refine:
1998	print >>pFile,'\n Hydrostatic/elastic strain: '
1999	print >>pFile,ptlbls
2000	print >>pFile,ptstr
2001	print >>pFile,sigstr
2002
2003	def PrintSHPOAndSig(pfx,hapData,POsig):
2004	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2005	ptlbls = ' names :'
2006	ptstr = ' values:'
2007	sigstr = ' sig :'
2008	for item in hapData[5]:
2009	ptlbls += '%12s'%(item)
2010	ptstr += '%12.3f'%(hapData[5][item])
2011	if pfx+item in POsig:
2012	sigstr += '%12.3f'%(POsig[pfx+item])
2013	else:
2014	sigstr += 12*' '
2015	print >>pFile,ptlbls
2016	print >>pFile,ptstr
2017	print >>pFile,sigstr
2018
2019	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2020	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2021	text = ''
2022	for item in hapData[2]:
2023	if pfx+item in ScalExtSig:
2024	text += ' %s: '%(item)
2025	text += '%12.2e'%(hapData[2][item][0])
2026	if pfx+item in ScalExtSig:
2027	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2028	print >>pFile,text
2029
2030	def PrintBabinetAndSig(pfx,hapData,BabSig):
2031	print >>pFile,'\n Babinet form factor modification: '
2032	ptlbls = ' names :'
2033	ptstr = ' values:'
2034	sigstr = ' sig :'
2035	for item in hapData:
2036	ptlbls += '%12s'%(item)
2037	ptstr += '%12.3f'%(hapData[item][0])
2038	if pfx+item in BabSig:
2039	sigstr += '%12.3f'%(BabSig[pfx+item])
2040	else:
2041	sigstr += 12*' '
2042	print >>pFile,ptlbls
2043	print >>pFile,ptstr
2044	print >>pFile,sigstr
2045
2046	PhFrExtPOSig = {}
2047	SizeMuStrSig = {}
2048	ScalExtSig = {}
2049	BabSig = {}
2050	wtFrSum = {}
2051	for phase in Phases:
2052	HistoPhase = Phases[phase]['Histograms']
2053	General = Phases[phase]['General']
2054	SGData = General['SGData']
2055	pId = Phases[phase]['pId']
2056	histoList = HistoPhase.keys()
2057	histoList.sort()
2058	for histogram in histoList:
2059	try:
2060	Histogram = Histograms[histogram]
2061	except KeyError:
2062	#skip if histogram not included e.g. in a sequential refinement
2063	continue
2064	hapData = HistoPhase[histogram]
2065	hId = Histogram['hId']
2066	pfx = str(pId)+':'+str(hId)+':'
2067	if hId not in wtFrSum:
2068	wtFrSum[hId] = 0.
2069	if 'PWDR' in histogram:
2070	for item in ['Scale','Extinction']:
2071	hapData[item][0] = parmDict[pfx+item]
2072	if pfx+item in sigDict:
2073	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2074	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2075	if hapData['Pref.Ori.'][0] == 'MD':
2076	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2077	if pfx+'MD' in sigDict:
2078	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2079	else: #'SH' spherical harmonics
2080	for item in hapData['Pref.Ori.'][5]:
2081	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2082	if pfx+item in sigDict:
2083	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2084	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2085	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2086	pfx+'HStrain':{}})
2087	for item in ['Mustrain','Size']:
2088	hapData[item][1][2] = parmDict[pfx+item+';mx']
2089	hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2090	if pfx+item+';mx' in sigDict:
2091	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2092	if hapData[item][0] in ['isotropic','uniaxial']:
2093	hapData[item][1][0] = parmDict[pfx+item+';i']
2094	if item == 'Size':
2095	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2096	if pfx+item+';i' in sigDict:
2097	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2098	if hapData[item][0] == 'uniaxial':
2099	hapData[item][1][1] = parmDict[pfx+item+';a']
2100	if item == 'Size':
2101	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2102	if pfx+item+';a' in sigDict:
2103	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2104	else: #generalized for mustrain or ellipsoidal for size
2105	Nterms = len(hapData[item][4])
2106	for i in range(Nterms):
2107	sfx = ':'+str(i)
2108	hapData[item][4][i] = parmDict[pfx+item+sfx]
2109	if pfx+item+sfx in sigDict:
2110	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2111	names = G2spc.HStrainNames(SGData)
2112	for i,name in enumerate(names):
2113	hapData['HStrain'][0][i] = parmDict[pfx+name]
2114	if pfx+name in sigDict:
2115	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2116	for name in ['BabA','BabU']:
2117	hapData['Babinet'][name][0] = parmDict[pfx+name]
2118	if pfx+name in sigDict:
2119	BabSig[pfx+name] = sigDict[pfx+name]
2120
2121	elif 'HKLF' in histogram:
2122	for item in ['Scale',]:
2123	if parmDict.get(pfx+item):
2124	hapData[item][0] = parmDict[pfx+item]
2125	if pfx+item in sigDict:
2126	ScalExtSig[pfx+item] = sigDict[pfx+item]
2127	for item in ['Ep','Eg','Es']:
2128	if parmDict.get(pfx+item):
2129	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2130	if pfx+item in sigDict:
2131	ScalExtSig[pfx+item] = sigDict[pfx+item]
2132	for name in ['BabA','BabU']:
2133	hapData['Babinet'][name][0] = parmDict[pfx+name]
2134	if pfx+name in sigDict:
2135	BabSig[pfx+name] = sigDict[pfx+name]
2136
2137	if Print:
2138	for phase in Phases:
2139	HistoPhase = Phases[phase]['Histograms']
2140	General = Phases[phase]['General']
2141	SGData = General['SGData']
2142	pId = Phases[phase]['pId']
2143	histoList = HistoPhase.keys()
2144	histoList.sort()
2145	for histogram in histoList:
2146	try:
2147	Histogram = Histograms[histogram]
2148	except KeyError:
2149	#skip if histogram not included e.g. in a sequential refinement
2150	continue
2151	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2152	print >>pFile,130*'-'
2153	hapData = HistoPhase[histogram]
2154	hId = Histogram['hId']
2155	Histogram['Residuals'][str(pId)+'::Name'] = phase
2156	pfx = str(pId)+':'+str(hId)+':'
2157	if 'PWDR' in histogram:
2158	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2159	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2160
2161	if pfx+'Scale' in PhFrExtPOSig:
2162	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2163	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2164	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2165	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2166	if pfx+'Extinction' in PhFrExtPOSig:
2167	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2168	if hapData['Pref.Ori.'][0] == 'MD':
2169	if pfx+'MD' in PhFrExtPOSig:
2170	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2171	else:
2172	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2173	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2174	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2175	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2176	if len(BabSig):
2177	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2178
2179	elif 'HKLF' in histogram:
2180	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2181	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2182	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2183	if pfx+'Scale' in ScalExtSig:
2184	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2185	if hapData['Extinction'][0] != 'None':
2186	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2187	if len(BabSig):
2188	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2189
2190	################################################################################
2191	##### Histogram data
2192	################################################################################
2193
2194	def GetHistogramData(Histograms,Print=True,pFile=None):
2195	'needs a doc string'
2196
2197	def GetBackgroundParms(hId,Background):
2198	Back = Background[0]
2199	DebyePeaks = Background[1]
2200	bakType,bakFlag = Back[:2]
2201	backVals = Back[3:]
2202	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2203	backDict = dict(zip(backNames,backVals))
2204	backVary = []
2205	if bakFlag:
2206	backVary = backNames
2207	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2208	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2209	debyeDict = {}
2210	debyeList = []
2211	for i in range(DebyePeaks['nDebye']):
2212	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2213	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2214	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2215	debyeVary = []
2216	for item in debyeList:
2217	if item[1]:
2218	debyeVary.append(item[0])
2219	backDict.update(debyeDict)
2220	backVary += debyeVary
2221	peakDict = {}
2222	peakList = []
2223	for i in range(DebyePeaks['nPeaks']):
2224	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2225	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2226	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2227	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2228	peakVary = []
2229	for item in peakList:
2230	if item[1]:
2231	peakVary.append(item[0])
2232	backDict.update(peakDict)
2233	backVary += peakVary
2234	return bakType,backDict,backVary
2235
2236	def GetInstParms(hId,Inst):
2237	dataType = Inst['Type'][0]
2238	instDict = {}
2239	insVary = []
2240	pfx = ':'+str(hId)+':'
2241	insKeys = Inst.keys()
2242	insKeys.sort()
2243	for item in insKeys:
2244	insName = pfx+item
2245	instDict[insName] = Inst[item][1]
2246	if len(Inst[item]) > 2 and Inst[item][2]:
2247	insVary.append(insName)
2248	if 'C' in dataType:
2249	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2250	return dataType,instDict,insVary
2251
2252	def GetSampleParms(hId,Sample):
2253	sampVary = []
2254	hfx = ':'+str(hId)+':'
2255	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2256	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2257	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2258	if key in Sample:
2259	sampDict[hfx+key] = Sample[key]
2260	Type = Sample['Type']
2261	if 'Bragg' in Type: #Bragg-Brentano
2262	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2263	sampDict[hfx+item] = Sample[item][0]
2264	if Sample[item][1]:
2265	sampVary.append(hfx+item)
2266	elif 'Debye' in Type: #Debye-Scherrer
2267	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2268	sampDict[hfx+item] = Sample[item][0]
2269	if Sample[item][1]:
2270	sampVary.append(hfx+item)
2271	return Type,sampDict,sampVary
2272
2273	def PrintBackground(Background):
2274	Back = Background[0]
2275	DebyePeaks = Background[1]
2276	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2277	line = ' Coefficients: '
2278	for i,back in enumerate(Back[3:]):
2279	line += '%10.3f'%(back)
2280	if i and not i%10:
2281	line += '\n'+15*' '
2282	print >>pFile,line
2283	if DebyePeaks['nDebye']:
2284	print >>pFile,'\n Debye diffuse scattering coefficients'
2285	parms = ['DebyeA','DebyeR','DebyeU']
2286	line = ' names : '
2287	for parm in parms:
2288	line += '%8s refine?'%(parm)
2289	print >>pFile,line
2290	for j,term in enumerate(DebyePeaks['debyeTerms']):
2291	line = ' term'+'%2d'%(j)+':'
2292	for i in range(3):
2293	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2294	print >>pFile,line
2295	if DebyePeaks['nPeaks']:
2296	print >>pFile,'\n Single peak coefficients'
2297	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2298	line = ' names : '
2299	for parm in parms:
2300	line += '%8s refine?'%(parm)
2301	print >>pFile,line
2302	for j,term in enumerate(DebyePeaks['peaksList']):
2303	line = ' peak'+'%2d'%(j)+':'
2304	for i in range(4):
2305	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2306	print >>pFile,line
2307
2308	def PrintInstParms(Inst):
2309	print >>pFile,'\n Instrument Parameters:'
2310	insKeys = Inst.keys()
2311	insKeys.sort()
2312	iBeg = 0
2313	Ok = True
2314	while Ok:
2315	ptlbls = ' name :'
2316	ptstr = ' value :'
2317	varstr = ' refine:'
2318	iFin = min(iBeg+9,len(insKeys))
2319	for item in insKeys[iBeg:iFin]:
2320	if item not in ['Type','Source']:
2321	ptlbls += '%12s' % (item)
2322	ptstr += '%12.6f' % (Inst[item][1])
2323	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2324	varstr += 12*' '
2325	else:
2326	varstr += '%12s' % (str(bool(Inst[item][2])))
2327	print >>pFile,ptlbls
2328	print >>pFile,ptstr
2329	print >>pFile,varstr
2330	iBeg = iFin
2331	if iBeg == len(insKeys):
2332	Ok = False
2333	else:
2334	print >>pFile,'\n'
2335
2336	def PrintSampleParms(Sample):
2337	print >>pFile,'\n Sample Parameters:'
2338	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2339	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2340	ptlbls = ' name :'
2341	ptstr = ' value :'
2342	varstr = ' refine:'
2343	if 'Bragg' in Sample['Type']:
2344	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2345	ptlbls += '%14s'%(item)
2346	ptstr += '%14.4f'%(Sample[item][0])
2347	varstr += '%14s'%(str(bool(Sample[item][1])))
2348
2349	elif 'Debye' in Type: #Debye-Scherrer
2350	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2351	ptlbls += '%14s'%(item)
2352	ptstr += '%14.4f'%(Sample[item][0])
2353	varstr += '%14s'%(str(bool(Sample[item][1])))
2354
2355	print >>pFile,ptlbls
2356	print >>pFile,ptstr
2357	print >>pFile,varstr
2358
2359	histDict = {}
2360	histVary = []
2361	controlDict = {}
2362	histoList = Histograms.keys()
2363	histoList.sort()
2364	for histogram in histoList:
2365	if 'PWDR' in histogram:
2366	Histogram = Histograms[histogram]
2367	hId = Histogram['hId']
2368	pfx = ':'+str(hId)+':'
2369	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2370	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2371	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2372	for excl in controlDict[pfx+'Exclude']:
2373	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2374	if controlDict[pfx+'Exclude']:
2375	ma.mask_rows(Histogram['Data'])
2376	Background = Histogram['Background']
2377	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2378	controlDict[pfx+'bakType'] = Type
2379	histDict.update(bakDict)
2380	histVary += bakVary
2381
2382	Inst = Histogram['Instrument Parameters'][0]
2383	Type,instDict,insVary = GetInstParms(hId,Inst)
2384	controlDict[pfx+'histType'] = Type
2385	if 'XC' in Type:
2386	if pfx+'Lam1' in instDict:
2387	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2388	else:
2389	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2390	histDict.update(instDict)
2391	histVary += insVary
2392
2393	Sample = Histogram['Sample Parameters']
2394	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2395	controlDict[pfx+'instType'] = Type
2396	histDict.update(sampDict)
2397	histVary += sampVary
2398
2399
2400	if Print:
2401	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2402	print >>pFile,135*'-'
2403	Units = {'C':' deg','T':' msec'}
2404	units = Units[controlDict[pfx+'histType'][2]]
2405	Limits = controlDict[pfx+'Limits']
2406	print >>pFile,' Instrument type: ',Sample['Type']
2407	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2408	if len(controlDict[pfx+'Exclude']):
2409	excls = controlDict[pfx+'Exclude']
2410	for excl in excls:
2411	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2412	PrintSampleParms(Sample)
2413	PrintInstParms(Inst)
2414	PrintBackground(Background)
2415	elif 'HKLF' in histogram:
2416	Histogram = Histograms[histogram]
2417	hId = Histogram['hId']
2418	pfx = ':'+str(hId)+':'
2419	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2420	Inst = Histogram['Instrument Parameters'][0]
2421	controlDict[pfx+'histType'] = Inst['Type'][0]
2422	if 'X' in Inst['Type'][0]:
2423	histDict[pfx+'Lam'] = Inst['Lam'][1]
2424	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2425	return histVary,histDict,controlDict
2426
2427	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2428	'needs a doc string'
2429
2430	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2431	Back = Background[0]
2432	DebyePeaks = Background[1]
2433	lenBack = len(Back[3:])
2434	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2435	for i in range(lenBack):
2436	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2437	if pfx+'Back;'+str(i) in sigDict:
2438	backSig[i] = sigDict[pfx+'Back;'+str(i)]
2439	if DebyePeaks['nDebye']:
2440	for i in range(DebyePeaks['nDebye']):
2441	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2442	for j,name in enumerate(names):
2443	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2444	if name in sigDict:
2445	backSig[lenBack+3*i+j] = sigDict[name]
2446	if DebyePeaks['nPeaks']:
2447	for i in range(DebyePeaks['nPeaks']):
2448	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2449	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2450	for j,name in enumerate(names):
2451	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2452	if name in sigDict:
2453	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2454	return backSig
2455
2456	def SetInstParms(pfx,Inst,parmDict,sigDict):
2457	instSig = {}
2458	insKeys = Inst.keys()
2459	insKeys.sort()
2460	for item in insKeys:
2461	insName = pfx+item
2462	Inst[item][1] = parmDict[insName]
2463	if insName in sigDict:
2464	instSig[item] = sigDict[insName]
2465	else:
2466	instSig[item] = 0
2467	return instSig
2468
2469	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2470	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2471	sampSig = [0 for i in range(5)]
2472	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2473	Sample[item][0] = parmDict[pfx+item]
2474	if pfx+item in sigDict:
2475	sampSig[i] = sigDict[pfx+item]
2476	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2477	sampSig = [0 for i in range(4)]
2478	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2479	Sample[item][0] = parmDict[pfx+item]
2480	if pfx+item in sigDict:
2481	sampSig[i] = sigDict[pfx+item]
2482	return sampSig
2483
2484	def PrintBackgroundSig(Background,backSig):
2485	Back = Background[0]
2486	DebyePeaks = Background[1]
2487	lenBack = len(Back[3:])
2488	valstr = ' value : '
2489	sigstr = ' sig : '
2490	refine = False
2491	for i,back in enumerate(Back[3:]):
2492	valstr += '%10.4g'%(back)
2493	if Back[1]:
2494	refine = True
2495	sigstr += '%10.4g'%(backSig[i])
2496	else:
2497	sigstr += 10*' '
2498	if refine:
2499	print >>pFile,'\n Background function: ',Back[0]
2500	print >>pFile,valstr
2501	print >>pFile,sigstr
2502	if DebyePeaks['nDebye']:
2503	ifAny = False
2504	ptfmt = "%12.3f"
2505	names = ' names :'
2506	ptstr = ' values:'
2507	sigstr = ' esds :'
2508	for item in sigDict:
2509	if 'Debye' in item:
2510	ifAny = True
2511	names += '%12s'%(item)
2512	ptstr += ptfmt%(parmDict[item])
2513	sigstr += ptfmt%(sigDict[item])
2514	if ifAny:
2515	print >>pFile,'\n Debye diffuse scattering coefficients'
2516	print >>pFile,names
2517	print >>pFile,ptstr
2518	print >>pFile,sigstr
2519	if DebyePeaks['nPeaks']:
2520	ifAny = False
2521	ptfmt = "%14.3f"
2522	names = ' names :'
2523	ptstr = ' values:'
2524	sigstr = ' esds :'
2525	for item in sigDict:
2526	if 'BkPk' in item:
2527	ifAny = True
2528	names += '%14s'%(item)
2529	ptstr += ptfmt%(parmDict[item])
2530	sigstr += ptfmt%(sigDict[item])
2531	if ifAny:
2532	print >>pFile,'\n Single peak coefficients'
2533	print >>pFile,names
2534	print >>pFile,ptstr
2535	print >>pFile,sigstr
2536
2537	def PrintInstParmsSig(Inst,instSig):
2538	refine = False
2539	insKeys = instSig.keys()
2540	insKeys.sort()
2541	iBeg = 0
2542	Ok = True
2543	while Ok:
2544	ptlbls = ' names :'
2545	ptstr = ' value :'
2546	sigstr = ' sig :'
2547	iFin = min(iBeg+9,len(insKeys))
2548	for name in insKeys[iBeg:iFin]:
2549	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2550	ptlbls += '%12s' % (name)
2551	ptstr += '%12.6f' % (Inst[name][1])
2552	if instSig[name]:
2553	refine = True
2554	sigstr += '%12.6f' % (instSig[name])
2555	else:
2556	sigstr += 12*' '
2557	if refine:
2558	print >>pFile,'\n Instrument Parameters:'
2559	print >>pFile,ptlbls
2560	print >>pFile,ptstr
2561	print >>pFile,sigstr
2562	iBeg = iFin
2563	if iBeg == len(insKeys):
2564	Ok = False
2565
2566	def PrintSampleParmsSig(Sample,sampleSig):
2567	ptlbls = ' names :'
2568	ptstr = ' values:'
2569	sigstr = ' sig :'
2570	refine = False
2571	if 'Bragg' in Sample['Type']:
2572	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2573	ptlbls += '%14s'%(item)
2574	ptstr += '%14.4f'%(Sample[item][0])
2575	if sampleSig[i]:
2576	refine = True
2577	sigstr += '%14.4f'%(sampleSig[i])
2578	else:
2579	sigstr += 14*' '
2580
2581	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2582	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2583	ptlbls += '%14s'%(item)
2584	ptstr += '%14.4f'%(Sample[item][0])
2585	if sampleSig[i]:
2586	refine = True
2587	sigstr += '%14.4f'%(sampleSig[i])
2588	else:
2589	sigstr += 14*' '
2590
2591	if refine:
2592	print >>pFile,'\n Sample Parameters:'
2593	print >>pFile,ptlbls
2594	print >>pFile,ptstr
2595	print >>pFile,sigstr
2596
2597	histoList = Histograms.keys()
2598	histoList.sort()
2599	for histogram in histoList:
2600	if 'PWDR' in histogram:
2601	Histogram = Histograms[histogram]
2602	hId = Histogram['hId']
2603	pfx = ':'+str(hId)+':'
2604	Background = Histogram['Background']
2605	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2606
2607	Inst = Histogram['Instrument Parameters'][0]
2608	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2609
2610	Sample = Histogram['Sample Parameters']
2611	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2612
2613	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2614	print >>pFile,135*'-'
2615	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2616	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2617	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2618	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2619	if Print:
2620	print >>pFile,' Instrument type: ',Sample['Type']
2621	PrintSampleParmsSig(Sample,sampSig)
2622	PrintInstParmsSig(Inst,instSig)
2623	PrintBackgroundSig(Background,backSig)
2624

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

Download in other formats: