Context Navigation

source: trunk/GSASIIstrIO.py @ 1604

Last change on this file since 1604 was 1604, checked in by vondreele, 9 years ago
make modulation wave maps fix atom index bugs begin modulated structure imput to least squares
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 119.7 KB

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

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

Download in other formats: