Context Navigation

source: trunk/GSASIIstrIO.py @ 1597

Last change on this file since 1597 was 1597, checked in by vondreele, 9 years ago
Incommensurate Pawley refinement works, including refinement of modulation vector.
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-06 21:12:32 +0000 (Sat, 06 Dec 2014) $
4	# $Author: vondreele $
5	# $Revision: 1597 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 1597 2014-12-06 21:12:32Z 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: 1597 $")
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	AtLookup = G2mth.FillAtomLookUp(Atoms)
966	PawleyRef = PhaseData[name].get('Pawley ref',[])
967	SGData = General['SGData']
968	SGtext,SGtable = G2spc.SGPrint(SGData)
969	cell = General['Cell']
970	A = G2lat.cell2A(cell[1:7])
971	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
972	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
973	if cell[0]:
974	phaseVary += cellVary(pfx,SGData)
975	SSGtext = [] #no superstructure
976	im = 0
977	if General['Type'] in ['modulated','magnetic']:
978	im = 1 #refl offset
979	Vec,vRef,maxH = General['SuperVec']
980	phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
981	SSGData = General['SSGData']
982	SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
983	if vRef:
984	phaseVary += modVary(pfx,SSGData)
985	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
986	if resRBData:
987	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
988	for iRB,RB in enumerate(resRBData):
989	MakeRBParms('R',phaseVary,phaseDict)
990	MakeRBThermals('R',phaseVary,phaseDict)
991	MakeRBTorsions('R',phaseVary,phaseDict)
992
993	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
994	if vecRBData:
995	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
996	for iRB,RB in enumerate(vecRBData):
997	MakeRBParms('V',phaseVary,phaseDict)
998	MakeRBThermals('V',phaseVary,phaseDict)
999
1000	Natoms[pfx] = 0
1001	if Atoms and not General.get('doPawley'):
1002	cx,ct,cs,cia = General['AtomPtrs']
1003	if General['Type'] in ['nuclear','macromolecular']:
1004	Natoms[pfx] = len(Atoms)
1005	for i,at in enumerate(Atoms):
1006	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
1007	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1008	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1009	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
1010	pfx+'AI/A:'+str(i):at[cia],})
1011	if at[cia] == 'I':
1012	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1013	else:
1014	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1015	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1016	if 'F' in at[ct+1]:
1017	phaseVary.append(pfx+'Afrac:'+str(i))
1018	if 'X' in at[ct+1]:
1019	try: #patch for sytsym name changes
1020	xId,xCoef = G2spc.GetCSxinel(at[cs])
1021	except KeyError:
1022	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1023	at[cs] = Sytsym
1024	xId,xCoef = G2spc.GetCSxinel(at[cs])
1025	xId,xCoef = G2spc.GetCSxinel(at[cs])
1026	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1027	equivs = [[],[],[]]
1028	for j in range(3):
1029	if xId[j] > 0:
1030	phaseVary.append(names[j])
1031	equivs[xId[j]-1].append([names[j],xCoef[j]])
1032	for equiv in equivs:
1033	if len(equiv) > 1:
1034	name = equiv[0][0]
1035	coef = equiv[0][1]
1036	for eqv in equiv[1:]:
1037	eqv[1] /= coef
1038	G2mv.StoreEquivalence(name,(eqv,))
1039	if 'U' in at[ct+1]:
1040	if at[cia] == 'I':
1041	phaseVary.append(pfx+'AUiso:'+str(i))
1042	else:
1043	try: #patch for sytsym name changes
1044	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1045	except KeyError:
1046	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1047	at[cs] = Sytsym
1048	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1049	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1050	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1051	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1052	equivs = [[],[],[],[],[],[]]
1053	for j in range(6):
1054	if uId[j] > 0:
1055	phaseVary.append(names[j])
1056	equivs[uId[j]-1].append([names[j],uCoef[j]])
1057	for equiv in equivs:
1058	if len(equiv) > 1:
1059	name = equiv[0][0]
1060	coef = equiv[0][1]
1061	for eqv in equiv[1:]:
1062	eqv[1] /= coef
1063	G2mv.StoreEquivalence(name,equiv[1:])
1064	# elif General['Type'] == 'magnetic':
1065	# elif General['Type'] == 'macromolecular':
1066	textureData = General['SH Texture']
1067	if textureData['Order']:
1068	phaseDict[pfx+'SHorder'] = textureData['Order']
1069	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1070	for item in ['omega','chi','phi']:
1071	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1072	if textureData['Sample '+item][0]:
1073	phaseVary.append(pfx+'SH '+item)
1074	for item in textureData['SH Coeff'][1]:
1075	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1076	if textureData['SH Coeff'][0]:
1077	phaseVary.append(pfx+item)
1078
1079	if Print:
1080	print >>pFile,'\n Phase name: ',General['Name']
1081	print >>pFile,135*'-'
1082	PrintFFtable(FFtable)
1083	PrintBLtable(BLtable)
1084	print >>pFile,''
1085	if len(SSGtext): #if superstructure
1086	for line in SSGtext: print >>pFile,line
1087	if len(SSGtable):
1088	for item in SSGtable:
1089	line = ' %s '%(item)
1090	print >>pFile,line
1091	else:
1092	print >>pFile,' ( 1) %s'%(SSGtable[0])
1093	else:
1094	for line in SGtext: print >>pFile,line
1095	if len(SGtable):
1096	for item in SGtable:
1097	line = ' %s '%(item)
1098	print >>pFile,line
1099	else:
1100	print >>pFile,' ( 1) %s'%(SGtable[0])
1101	PrintRBObjects(resRBData,vecRBData)
1102	PrintAtoms(General,Atoms)
1103	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1104	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1105	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1106	if len(SSGtext): #if superstructure
1107	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1108	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1109	PrintTexture(textureData)
1110	if name in RestraintDict:
1111	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1112	textureData,RestraintDict[name],pFile)
1113
1114	elif PawleyRef:
1115	if Print:
1116	print >>pFile,'\n Phase name: ',General['Name']
1117	print >>pFile,135*'-'
1118	print >>pFile,''
1119	if len(SSGtext): #if superstructure
1120	for line in SSGtext: print >>pFile,line
1121	if len(SSGtable):
1122	for item in SSGtable:
1123	line = ' %s '%(item)
1124	print >>pFile,line
1125	else:
1126	print >>pFile,' ( 1) %s'%(SSGtable[0])
1127	else:
1128	for line in SGtext: print >>pFile,line
1129	if len(SGtable):
1130	for item in SGtable:
1131	line = ' %s '%(item)
1132	print >>pFile,line
1133	else:
1134	print >>pFile,' ( 1) %s'%(SGtable[0])
1135	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1136	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1137	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1138	if len(SSGtext): #if superstructure
1139	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1140	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1141	pawleyVary = []
1142	for i,refl in enumerate(PawleyRef):
1143	phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1144	if im:
1145	pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1146	else:
1147	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1148	if refl[5+im]:
1149	pawleyVary.append(pfx+'PWLref:'+str(i))
1150	GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary) #does G2mv.StoreEquivalence
1151	phaseVary += pawleyVary
1152
1153	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1154
1155	def cellFill(pfx,SGData,parmDict,sigDict):
1156	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1157	from the parameter and sig dictionaries.
1158
1159	:param str pfx: parameter prefix ("n::", where n is a phase number)
1160	:param dict SGdata: a symmetry object
1161	:param dict parmDict: a dictionary of parameters
1162	:param dict sigDict: a dictionary of uncertainties on parameters
1163
1164	:returns: A,sigA where each is a list of six terms with the A terms
1165	'''
1166	if SGData['SGLaue'] in ['-1',]:
1167	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1168	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1169	elif SGData['SGLaue'] in ['2/m',]:
1170	if SGData['SGUniq'] == 'a':
1171	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1172	parmDict[pfx+'A3'],0,0]
1173	elif SGData['SGUniq'] == 'b':
1174	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1175	0,parmDict[pfx+'A4'],0]
1176	else:
1177	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1178	0,0,parmDict[pfx+'A5']]
1179	elif SGData['SGLaue'] in ['mmm',]:
1180	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1181	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1182	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1183	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1184	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1185	parmDict[pfx+'A0'],0,0]
1186	elif SGData['SGLaue'] in ['3R', '3mR']:
1187	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1188	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1189	elif SGData['SGLaue'] in ['m3m','m3']:
1190	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1191
1192	try:
1193	if SGData['SGLaue'] in ['-1',]:
1194	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1195	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1196	elif SGData['SGLaue'] in ['2/m',]:
1197	if SGData['SGUniq'] == 'a':
1198	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1199	sigDict[pfx+'A3'],0,0]
1200	elif SGData['SGUniq'] == 'b':
1201	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1202	0,sigDict[pfx+'A4'],0]
1203	else:
1204	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1205	0,0,sigDict[pfx+'A5']]
1206	elif SGData['SGLaue'] in ['mmm',]:
1207	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1208	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1209	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1210	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1211	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1212	elif SGData['SGLaue'] in ['3R', '3mR']:
1213	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1214	elif SGData['SGLaue'] in ['m3m','m3']:
1215	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1216	except KeyError:
1217	sigA = [0,0,0,0,0,0]
1218
1219	return A,sigA
1220
1221	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1222	'needs a doc string'
1223	if phaseRest:
1224	Amat = G2lat.cell2AB(cell)[0]
1225	cx,ct,cs = AtPtrs[:3]
1226	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1227	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1228	['ChemComp','Sites'],['Texture','HKLs']]
1229	for name,rest in names:
1230	itemRest = phaseRest[name]
1231	if itemRest[rest] and itemRest['Use']:
1232	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1233	if name in ['Bond','Angle','Plane','Chiral']:
1234	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1235	for indx,ops,obs,esd in itemRest[rest]:
1236	try:
1237	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1238	AtName = ''
1239	for i,Aname in enumerate(AtNames):
1240	AtName += Aname
1241	if ops[i] == '1':
1242	AtName += '-'
1243	else:
1244	AtName += '+('+ops[i]+')-'
1245	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1246	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1247	if name == 'Bond':
1248	calc = G2mth.getRestDist(XYZ,Amat)
1249	elif name == 'Angle':
1250	calc = G2mth.getRestAngle(XYZ,Amat)
1251	elif name == 'Plane':
1252	calc = G2mth.getRestPlane(XYZ,Amat)
1253	elif name == 'Chiral':
1254	calc = G2mth.getRestChiral(XYZ,Amat)
1255	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1256	except KeyError:
1257	del itemRest[rest]
1258	elif name in ['Torsion','Rama']:
1259	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1260	coeffDict = itemRest['Coeff']
1261	for indx,ops,cofName,esd in itemRest[rest]:
1262	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1263	AtName = ''
1264	for i,Aname in enumerate(AtNames):
1265	AtName += Aname+'+('+ops[i]+')-'
1266	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1267	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1268	if name == 'Torsion':
1269	tor = G2mth.getRestTorsion(XYZ,Amat)
1270	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1271	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1272	else:
1273	phi,psi = G2mth.getRestRama(XYZ,Amat)
1274	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1275	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])
1276	elif name == 'ChemComp':
1277	print >>pFile,' atoms mul*frac factor prod'
1278	for indx,factors,obs,esd in itemRest[rest]:
1279	try:
1280	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1281	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1282	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1283	mulfrac = mul*frac
1284	calcs = mulfracfactors
1285	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1286	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1287	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1288	except KeyError:
1289	del itemRest[rest]
1290	elif name == 'Texture' and textureData['Order']:
1291	Start = False
1292	SHCoef = textureData['SH Coeff'][1]
1293	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1294	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1295	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1296	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1297	PH = np.array(hkl)
1298	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1299	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1300	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1301	Z = ma.masked_greater(Z,0.0)
1302	num = ma.count(Z)
1303	sum = 0
1304	if num:
1305	sum = np.sum(Z)
1306	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)
1307
1308	def getCellEsd(pfx,SGData,A,covData):
1309	'needs a doc string'
1310	dpr = 180./np.pi
1311	rVsq = G2lat.calc_rVsq(A)
1312	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1313	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1314	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1315	scos = cosd(cellst[3:6])
1316	ssin = sind(cellst[3:6])
1317	scot = scos/ssin
1318	rcos = cosd(cell[3:6])
1319	rsin = sind(cell[3:6])
1320	rcot = rcos/rsin
1321	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1322	varyList = covData['varyList']
1323	covMatrix = covData['covMatrix']
1324	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1325	Ax = np.array(A)
1326	Ax[3:] /= 2.
1327	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1328	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1329	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1330	Vol = 1/np.sqrt(rVsq)
1331	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1332	R123 = Ax[0]Ax[1]Ax[2]
1333	dsasdg = np.zeros((3,6))
1334	dadg = np.zeros((6,6))
1335	for i0 in range(3): #0 1 2
1336	i1 = (i0+1)%3 #1 2 0
1337	i2 = (i1+1)%3 #2 0 1
1338	i3 = 5-i2 #3 5 4
1339	i4 = 5-i1 #4 3 5
1340	i5 = 5-i0 #5 4 3
1341	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1342	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1343	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1344	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1345	denom = np.sqrt(denmsq)
1346	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1347	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1348	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1349	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1350	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1351	dadg[i5][i5] = -Ax[i0]/denom
1352	for i0 in range(3):
1353	i1 = (i0+1)%3
1354	i2 = (i1+1)%3
1355	i3 = 5-i2
1356	for ij in range(6):
1357	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1358	if ij == i0:
1359	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1360	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1361	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1362	var = np.diag(sigMat)
1363	CS = np.where(var>0.,np.sqrt(var),0.)
1364	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1365
1366	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1367	'needs a doc string'
1368
1369	def PrintAtomsAndSig(General,Atoms,atomsSig):
1370	print >>pFile,'\n Atoms:'
1371	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1372	if General['Type'] == 'magnetic':
1373	line += ' Mx My Mz'
1374	elif General['Type'] == 'macromolecular':
1375	line = ' res no residue chain '+line
1376	cx,ct,cs,cia = General['AtomPtrs']
1377	print >>pFile,line
1378	if General['Type'] == 'nuclear':
1379	print >>pFile,135*'-'
1380	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1381	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1382	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1383	for atyp in General['AtomTypes']: #zero composition
1384	General['NoAtoms'][atyp] = 0.
1385	for i,at in enumerate(Atoms):
1386	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1387	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1388	valstr = ' values:'
1389	sigstr = ' sig :'
1390	for ind in [3,4,5,6]:
1391	sigind = str(i)+':'+str(ind)
1392	valstr += fmt[ind]%(at[ind])
1393	if sigind in atomsSig:
1394	sigstr += fmt[ind]%(atomsSig[sigind])
1395	else:
1396	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1397	if at[9] == 'I':
1398	valstr += fmt[10]%(at[10])
1399	if str(i)+':10' in atomsSig:
1400	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1401	else:
1402	sigstr += 8*' '
1403	else:
1404	valstr += 8*' '
1405	sigstr += 8*' '
1406	for ind in [11,12,13,14,15,16]:
1407	sigind = str(i)+':'+str(ind)
1408	valstr += fmt[ind]%(at[ind])
1409	if sigind in atomsSig:
1410	sigstr += fmt[ind]%(atomsSig[sigind])
1411	else:
1412	sigstr += 8*' '
1413	print >>pFile,name
1414	print >>pFile,valstr
1415	print >>pFile,sigstr
1416
1417	def PrintRBObjPOAndSig(rbfx,rbsx):
1418	namstr = ' names :'
1419	valstr = ' values:'
1420	sigstr = ' esds :'
1421	for i,px in enumerate(['Px:','Py:','Pz:']):
1422	name = pfx+rbfx+px+rbsx
1423	namstr += '%12s'%('Pos '+px[1])
1424	valstr += '%12.5f'%(parmDict[name])
1425	if name in sigDict:
1426	sigstr += '%12.5f'%(sigDict[name])
1427	else:
1428	sigstr += 12*' '
1429	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1430	name = pfx+rbfx+po+rbsx
1431	namstr += '%12s'%('Ori '+po[1])
1432	valstr += '%12.5f'%(parmDict[name])
1433	if name in sigDict:
1434	sigstr += '%12.5f'%(sigDict[name])
1435	else:
1436	sigstr += 12*' '
1437	print >>pFile,namstr
1438	print >>pFile,valstr
1439	print >>pFile,sigstr
1440
1441	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1442	namstr = ' names :'
1443	valstr = ' values:'
1444	sigstr = ' esds :'
1445	if 'N' not in TLS:
1446	print >>pFile,' Thermal motion:'
1447	if 'T' in TLS:
1448	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1449	name = pfx+rbfx+pt+rbsx
1450	namstr += '%12s'%(pt[:3])
1451	valstr += '%12.5f'%(parmDict[name])
1452	if name in sigDict:
1453	sigstr += '%12.5f'%(sigDict[name])
1454	else:
1455	sigstr += 12*' '
1456	print >>pFile,namstr
1457	print >>pFile,valstr
1458	print >>pFile,sigstr
1459	if 'L' in TLS:
1460	namstr = ' names :'
1461	valstr = ' values:'
1462	sigstr = ' esds :'
1463	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1464	name = pfx+rbfx+pt+rbsx
1465	namstr += '%12s'%(pt[:3])
1466	valstr += '%12.3f'%(parmDict[name])
1467	if name in sigDict:
1468	sigstr += '%12.3f'%(sigDict[name])
1469	else:
1470	sigstr += 12*' '
1471	print >>pFile,namstr
1472	print >>pFile,valstr
1473	print >>pFile,sigstr
1474	if 'S' in TLS:
1475	namstr = ' names :'
1476	valstr = ' values:'
1477	sigstr = ' esds :'
1478	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1479	name = pfx+rbfx+pt+rbsx
1480	namstr += '%12s'%(pt[:3])
1481	valstr += '%12.4f'%(parmDict[name])
1482	if name in sigDict:
1483	sigstr += '%12.4f'%(sigDict[name])
1484	else:
1485	sigstr += 12*' '
1486	print >>pFile,namstr
1487	print >>pFile,valstr
1488	print >>pFile,sigstr
1489	if 'U' in TLS:
1490	name = pfx+rbfx+'U:'+rbsx
1491	namstr = ' names :'+'%12s'%('Uiso')
1492	valstr = ' values:'+'%12.5f'%(parmDict[name])
1493	if name in sigDict:
1494	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1495	else:
1496	sigstr = ' esds :'+12*' '
1497	print >>pFile,namstr
1498	print >>pFile,valstr
1499	print >>pFile,sigstr
1500
1501	def PrintRBObjTorAndSig(rbsx):
1502	namstr = ' names :'
1503	valstr = ' values:'
1504	sigstr = ' esds :'
1505	nTors = len(RBObj['Torsions'])
1506	if nTors:
1507	print >>pFile,' Torsions:'
1508	for it in range(nTors):
1509	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1510	namstr += '%12s'%('Tor'+str(it))
1511	valstr += '%12.4f'%(parmDict[name])
1512	if name in sigDict:
1513	sigstr += '%12.4f'%(sigDict[name])
1514	print >>pFile,namstr
1515	print >>pFile,valstr
1516	print >>pFile,sigstr
1517
1518	def PrintSHtextureAndSig(textureData,SHtextureSig):
1519	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1520	names = ['omega','chi','phi']
1521	namstr = ' names :'
1522	ptstr = ' values:'
1523	sigstr = ' esds :'
1524	for name in names:
1525	namstr += '%12s'%(name)
1526	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1527	if 'Sample '+name in SHtextureSig:
1528	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1529	else:
1530	sigstr += 12*' '
1531	print >>pFile,namstr
1532	print >>pFile,ptstr
1533	print >>pFile,sigstr
1534	print >>pFile,'\n Texture coefficients:'
1535	namstr = ' names :'
1536	ptstr = ' values:'
1537	sigstr = ' esds :'
1538	SHcoeff = textureData['SH Coeff'][1]
1539	for name in SHcoeff:
1540	namstr += '%12s'%(name)
1541	ptstr += '%12.3f'%(SHcoeff[name])
1542	if name in SHtextureSig:
1543	sigstr += '%12.3f'%(SHtextureSig[name])
1544	else:
1545	sigstr += 12*' '
1546	print >>pFile,namstr
1547	print >>pFile,ptstr
1548	print >>pFile,sigstr
1549
1550	print >>pFile,'\n Phases:'
1551	for phase in Phases:
1552	print >>pFile,' Result for phase: ',phase
1553	Phase = Phases[phase]
1554	General = Phase['General']
1555	SGData = General['SGData']
1556	Atoms = Phase['Atoms']
1557	AtLookup = G2mth.FillAtomLookUp(Atoms)
1558	cell = General['Cell']
1559	pId = Phase['pId']
1560	pfx = str(pId)+'::'
1561	if cell[0]:
1562	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1563	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1564	print >>pFile,' Reciprocal metric tensor: '
1565	ptfmt = "%15.9f"
1566	names = ['A11','A22','A33','A12','A13','A23']
1567	namstr = ' names :'
1568	ptstr = ' values:'
1569	sigstr = ' esds :'
1570	for name,a,siga in zip(names,A,sigA):
1571	namstr += '%15s'%(name)
1572	ptstr += ptfmt%(a)
1573	if siga:
1574	sigstr += ptfmt%(siga)
1575	else:
1576	sigstr += 15*' '
1577	print >>pFile,namstr
1578	print >>pFile,ptstr
1579	print >>pFile,sigstr
1580	cell[1:7] = G2lat.A2cell(A)
1581	cell[7] = G2lat.calc_V(A)
1582	print >>pFile,' New unit cell:'
1583	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1584	names = ['a','b','c','alpha','beta','gamma','Volume']
1585	namstr = ' names :'
1586	ptstr = ' values:'
1587	sigstr = ' esds :'
1588	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1589	namstr += '%12s'%(name)
1590	ptstr += fmt%(a)
1591	if siga:
1592	sigstr += fmt%(siga)
1593	else:
1594	sigstr += 12*' '
1595	print >>pFile,namstr
1596	print >>pFile,ptstr
1597	print >>pFile,sigstr
1598	ik = 6 #for Pawley stuff below
1599	if General['Type'] in ['modulated','magnetic']:
1600	ik = 7
1601	Vec,vRef,maxH = General['SuperVec']
1602	if vRef:
1603	print >>pFile,' New modulation vector:'
1604	namstr = ' names :'
1605	ptstr = ' values:'
1606	sigstr = ' esds :'
1607	for var in ['mV0','mV1','mV2']:
1608	namstr += '%12s'%(pfx+var)
1609	ptstr += '%12.6f'%(parmDict[pfx+var])
1610	if pfx+var in sigDict:
1611	sigstr += '%12.6f'%(sigDict[pfx+var])
1612	else:
1613	sigstr += 12*' '
1614	print >>pFile,namstr
1615	print >>pFile,ptstr
1616	print >>pFile,sigstr
1617
1618	General['Mass'] = 0.
1619	if Phase['General'].get('doPawley'):
1620	pawleyRef = Phase['Pawley ref']
1621	for i,refl in enumerate(pawleyRef):
1622	key = pfx+'PWLref:'+str(i)
1623	refl[ik] = parmDict[key]
1624	if key in sigDict:
1625	refl[ik+1] = sigDict[key]
1626	else:
1627	refl[ik+1] = 0
1628	else:
1629	VRBIds = RBIds['Vector']
1630	RRBIds = RBIds['Residue']
1631	RBModels = Phase['RBModels']
1632	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1633	jrb = VRBIds.index(RBObj['RBId'])
1634	rbsx = str(irb)+':'+str(jrb)
1635	print >>pFile,' Vector rigid body parameters:'
1636	PrintRBObjPOAndSig('RBV',rbsx)
1637	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1638	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1639	jrb = RRBIds.index(RBObj['RBId'])
1640	rbsx = str(irb)+':'+str(jrb)
1641	print >>pFile,' Residue rigid body parameters:'
1642	PrintRBObjPOAndSig('RBR',rbsx)
1643	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1644	PrintRBObjTorAndSig(rbsx)
1645	atomsSig = {}
1646	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1647	for i,at in enumerate(Atoms):
1648	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1649	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1650	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1651	for ind in [3,4,5,6]:
1652	at[ind] = parmDict[names[ind]]
1653	if ind in [3,4,5]:
1654	name = names[ind].replace('A','dA')
1655	else:
1656	name = names[ind]
1657	if name in sigDict:
1658	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1659	if at[9] == 'I':
1660	at[10] = parmDict[names[10]]
1661	if names[10] in sigDict:
1662	atomsSig[str(i)+':10'] = sigDict[names[10]]
1663	else:
1664	for ind in [11,12,13,14,15,16]:
1665	at[ind] = parmDict[names[ind]]
1666	if names[ind] in sigDict:
1667	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1668	ind = General['AtomTypes'].index(at[1])
1669	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1670	PrintAtomsAndSig(General,Atoms,atomsSig)
1671
1672	textureData = General['SH Texture']
1673	if textureData['Order']:
1674	SHtextureSig = {}
1675	for name in ['omega','chi','phi']:
1676	aname = pfx+'SH '+name
1677	textureData['Sample '+name][1] = parmDict[aname]
1678	if aname in sigDict:
1679	SHtextureSig['Sample '+name] = sigDict[aname]
1680	for name in textureData['SH Coeff'][1]:
1681	aname = pfx+name
1682	textureData['SH Coeff'][1][name] = parmDict[aname]
1683	if aname in sigDict:
1684	SHtextureSig[name] = sigDict[aname]
1685	PrintSHtextureAndSig(textureData,SHtextureSig)
1686	if phase in RestraintDict:
1687	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1688	textureData,RestraintDict[phase],pFile)
1689
1690	################################################################################
1691	##### Histogram & Phase data
1692	################################################################################
1693
1694	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1695	'''Loads the HAP histogram/phase information into dicts
1696
1697	:param dict Phases: phase information
1698	:param dict Histograms: Histogram information
1699	:param bool Print: prints information as it is read
1700	:param file pFile: file object to print to (the default, None causes printing to the console)
1701	:param bool resetRefList: Should the contents of the Reflection List be initialized
1702	on loading. The default, True, initializes the Reflection List as it is loaded.
1703
1704	:returns: (hapVary,hapDict,controlDict)
1705	* hapVary: list of refined variables
1706	* hapDict: dict with refined variables and their values
1707	* controlDict: dict with computation controls (?)
1708	'''
1709
1710	def PrintSize(hapData):
1711	if hapData[0] in ['isotropic','uniaxial']:
1712	line = '\n Size model : %9s'%(hapData[0])
1713	line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1714	if hapData[0] == 'uniaxial':
1715	line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1716	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1717	print >>pFile,line
1718	else:
1719	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1720	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1721	Snames = ['S11','S22','S33','S12','S13','S23']
1722	ptlbls = ' names :'
1723	ptstr = ' values:'
1724	varstr = ' refine:'
1725	for i,name in enumerate(Snames):
1726	ptlbls += '%12s' % (name)
1727	ptstr += '%12.6f' % (hapData[4][i])
1728	varstr += '%12s' % (str(hapData[5][i]))
1729	print >>pFile,ptlbls
1730	print >>pFile,ptstr
1731	print >>pFile,varstr
1732
1733	def PrintMuStrain(hapData,SGData):
1734	if hapData[0] in ['isotropic','uniaxial']:
1735	line = '\n Mustrain model: %9s'%(hapData[0])
1736	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1737	if hapData[0] == 'uniaxial':
1738	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1739	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1740	print >>pFile,line
1741	else:
1742	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1743	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1744	Snames = G2spc.MustrainNames(SGData)
1745	ptlbls = ' names :'
1746	ptstr = ' values:'
1747	varstr = ' refine:'
1748	for i,name in enumerate(Snames):
1749	ptlbls += '%12s' % (name)
1750	ptstr += '%12.6f' % (hapData[4][i])
1751	varstr += '%12s' % (str(hapData[5][i]))
1752	print >>pFile,ptlbls
1753	print >>pFile,ptstr
1754	print >>pFile,varstr
1755
1756	def PrintHStrain(hapData,SGData):
1757	print >>pFile,'\n Hydrostatic/elastic strain: '
1758	Hsnames = G2spc.HStrainNames(SGData)
1759	ptlbls = ' names :'
1760	ptstr = ' values:'
1761	varstr = ' refine:'
1762	for i,name in enumerate(Hsnames):
1763	ptlbls += '%12s' % (name)
1764	ptstr += '%12.4g' % (hapData[0][i])
1765	varstr += '%12s' % (str(hapData[1][i]))
1766	print >>pFile,ptlbls
1767	print >>pFile,ptstr
1768	print >>pFile,varstr
1769
1770	def PrintSHPO(hapData):
1771	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1772	str(hapData[4])+' Refine? '+str(hapData[2])
1773	ptlbls = ' names :'
1774	ptstr = ' values:'
1775	for item in hapData[5]:
1776	ptlbls += '%12s'%(item)
1777	ptstr += '%12.3f'%(hapData[5][item])
1778	print >>pFile,ptlbls
1779	print >>pFile,ptstr
1780
1781	def PrintBabinet(hapData):
1782	print >>pFile,'\n Babinet form factor modification: '
1783	ptlbls = ' names :'
1784	ptstr = ' values:'
1785	varstr = ' refine:'
1786	for name in ['BabA','BabU']:
1787	ptlbls += '%12s' % (name)
1788	ptstr += '%12.6f' % (hapData[name][0])
1789	varstr += '%12s' % (str(hapData[name][1]))
1790	print >>pFile,ptlbls
1791	print >>pFile,ptstr
1792	print >>pFile,varstr
1793
1794
1795	hapDict = {}
1796	hapVary = []
1797	controlDict = {}
1798	poType = {}
1799	poAxes = {}
1800	spAxes = {}
1801	spType = {}
1802
1803	for phase in Phases:
1804	HistoPhase = Phases[phase]['Histograms']
1805	SGData = Phases[phase]['General']['SGData']
1806	cell = Phases[phase]['General']['Cell'][1:7]
1807	A = G2lat.cell2A(cell)
1808	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
1809	SSGData = Phases[phase]['General']['SSGData']
1810	Vec,x,maxH = Phases[phase]['General']['SuperVec']
1811	pId = Phases[phase]['pId']
1812	histoList = HistoPhase.keys()
1813	histoList.sort()
1814	for histogram in histoList:
1815	try:
1816	Histogram = Histograms[histogram]
1817	except KeyError:
1818	#skip if histogram not included e.g. in a sequential refinement
1819	continue
1820	hapData = HistoPhase[histogram]
1821	hId = Histogram['hId']
1822	if 'PWDR' in histogram:
1823	limits = Histogram['Limits'][1]
1824	inst = Histogram['Instrument Parameters'][0]
1825	Zero = inst['Zero'][0]
1826	if 'C' in inst['Type'][1]:
1827	try:
1828	wave = inst['Lam'][1]
1829	except KeyError:
1830	wave = inst['Lam1'][1]
1831	dmin = wave/(2.0*sind(limits[1]/2.0))
1832	elif 'T' in inst['Type'][0]:
1833	dmin = limits[0]/inst['difC'][1]
1834	pfx = str(pId)+':'+str(hId)+':'
1835	for item in ['Scale','Extinction']:
1836	hapDict[pfx+item] = hapData[item][0]
1837	if hapData[item][1]:
1838	hapVary.append(pfx+item)
1839	names = G2spc.HStrainNames(SGData)
1840	HSvals = []
1841	for i,name in enumerate(names):
1842	hapDict[pfx+name] = hapData['HStrain'][0][i]
1843	HSvals.append(hapDict[pfx+name])
1844	if hapData['HStrain'][1][i]:
1845	hapVary.append(pfx+name)
1846	DIJS = G2spc.HStrainVals(HSvals,SGData)
1847	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1848	if hapData['Pref.Ori.'][0] == 'MD':
1849	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1850	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1851	if hapData['Pref.Ori.'][2]:
1852	hapVary.append(pfx+'MD')
1853	else: #'SH' spherical harmonics
1854	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1855	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1856	for item in hapData['Pref.Ori.'][5]:
1857	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1858	if hapData['Pref.Ori.'][2]:
1859	hapVary.append(pfx+item)
1860	for item in ['Mustrain','Size']:
1861	controlDict[pfx+item+'Type'] = hapData[item][0]
1862	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1863	if hapData[item][2][2]:
1864	hapVary.append(pfx+item+';mx')
1865	if hapData[item][0] in ['isotropic','uniaxial']:
1866	hapDict[pfx+item+';i'] = hapData[item][1][0]
1867	if hapData[item][2][0]:
1868	hapVary.append(pfx+item+';i')
1869	if hapData[item][0] == 'uniaxial':
1870	controlDict[pfx+item+'Axis'] = hapData[item][3]
1871	hapDict[pfx+item+';a'] = hapData[item][1][1]
1872	if hapData[item][2][1]:
1873	hapVary.append(pfx+item+';a')
1874	else: #generalized for mustrain or ellipsoidal for size
1875	Nterms = len(hapData[item][4])
1876	if item == 'Mustrain':
1877	names = G2spc.MustrainNames(SGData)
1878	pwrs = []
1879	for name in names:
1880	h,k,l = name[1:]
1881	pwrs.append([int(h),int(k),int(l)])
1882	controlDict[pfx+'MuPwrs'] = pwrs
1883	for i in range(Nterms):
1884	sfx = ':'+str(i)
1885	hapDict[pfx+item+sfx] = hapData[item][4][i]
1886	if hapData[item][5][i]:
1887	hapVary.append(pfx+item+sfx)
1888	for bab in ['BabA','BabU']:
1889	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1890	if hapData['Babinet'][bab][1]:
1891	hapVary.append(pfx+bab)
1892
1893	if Print:
1894	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1895	print >>pFile,135*'-'
1896	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1897	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1898	if hapData['Pref.Ori.'][0] == 'MD':
1899	Ax = hapData['Pref.Ori.'][3]
1900	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1901	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1902	else: #'SH' for spherical harmonics
1903	PrintSHPO(hapData['Pref.Ori.'])
1904	PrintSize(hapData['Size'])
1905	PrintMuStrain(hapData['Mustrain'],SGData)
1906	PrintHStrain(hapData['HStrain'],SGData)
1907	if hapData['Babinet']['BabA'][0]:
1908	PrintBabinet(hapData['Babinet'])
1909	if resetRefList:
1910	refList = []
1911	Uniq = []
1912	Phi = []
1913	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
1914	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
1915	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
1916	for h,k,l,m,d in HKLd:
1917	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData) #is this right for SS refl.??
1918	mul *= 2 # for powder overlap of Friedel pairs
1919	if m or not ext:
1920	if 'C' in inst['Type'][0]:
1921	pos = G2lat.Dsp2pos(inst,d)
1922	if limits[0] < pos < limits[1]:
1923	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])
1924	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1925	Uniq.append(uniq)
1926	Phi.append(phi)
1927	elif 'T' in inst['Type'][0]:
1928	pos = G2lat.Dsp2pos(inst,d)
1929	if limits[0] < pos < limits[1]:
1930	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
1931	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])
1932	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1933	Uniq.append(uniq)
1934	Phi.append(phi)
1935	else:
1936	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1937	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
1938	for h,k,l,d in HKLd:
1939	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1940	mul *= 2 # for powder overlap of Friedel pairs
1941	if ext:
1942	continue
1943	if 'C' in inst['Type'][0]:
1944	pos = G2lat.Dsp2pos(inst,d)
1945	if limits[0] < pos < limits[1]:
1946	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])
1947	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1948	Uniq.append(uniq)
1949	Phi.append(phi)
1950	elif 'T' in inst['Type'][0]:
1951	pos = G2lat.Dsp2pos(inst,d)
1952	if limits[0] < pos < limits[1]:
1953	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
1954	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])
1955	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1956	Uniq.append(uniq)
1957	Phi.append(phi)
1958	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0]}
1959	elif 'HKLF' in histogram:
1960	inst = Histogram['Instrument Parameters'][0]
1961	hId = Histogram['hId']
1962	hfx = ':%d:'%(hId)
1963	for item in inst:
1964	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1965	hapDict[hfx+item] = inst[item][1]
1966	pfx = str(pId)+':'+str(hId)+':'
1967	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1968	if hapData['Scale'][1]:
1969	hapVary.append(pfx+'Scale')
1970
1971	extApprox,extType,extParms = hapData['Extinction']
1972	controlDict[pfx+'EType'] = extType
1973	controlDict[pfx+'EApprox'] = extApprox
1974	if 'C' in inst['Type'][0]:
1975	controlDict[pfx+'Tbar'] = extParms['Tbar']
1976	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1977	if 'Primary' in extType:
1978	Ekey = ['Ep',]
1979	elif 'I & II' in extType:
1980	Ekey = ['Eg','Es']
1981	elif 'Secondary Type II' == extType:
1982	Ekey = ['Es',]
1983	elif 'Secondary Type I' == extType:
1984	Ekey = ['Eg',]
1985	else: #'None'
1986	Ekey = []
1987	for eKey in Ekey:
1988	hapDict[pfx+eKey] = extParms[eKey][0]
1989	if extParms[eKey][1]:
1990	hapVary.append(pfx+eKey)
1991	for bab in ['BabA','BabU']:
1992	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1993	if hapData['Babinet'][bab][1]:
1994	hapVary.append(pfx+bab)
1995	if Print:
1996	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1997	print >>pFile,135*'-'
1998	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1999	if extType != 'None':
2000	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox)
2001	text = ' Parameters :'
2002	for eKey in Ekey:
2003	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2004	print >>pFile,text
2005	if hapData['Babinet']['BabA'][0]:
2006	PrintBabinet(hapData['Babinet'])
2007	Histogram['Reflection Lists'] = phase
2008
2009	return hapVary,hapDict,controlDict
2010
2011	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
2012	'needs a doc string'
2013
2014	def PrintSizeAndSig(hapData,sizeSig):
2015	line = '\n Size model: %9s'%(hapData[0])
2016	refine = False
2017	if hapData[0] in ['isotropic','uniaxial']:
2018	line += ' equatorial:%12.5f'%(hapData[1][0])
2019	if sizeSig[0][0]:
2020	line += ', sig:%8.4f'%(sizeSig[0][0])
2021	refine = True
2022	if hapData[0] == 'uniaxial':
2023	line += ' axial:%12.4f'%(hapData[1][1])
2024	if sizeSig[0][1]:
2025	refine = True
2026	line += ', sig:%8.4f'%(sizeSig[0][1])
2027	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2028	if sizeSig[0][2]:
2029	refine = True
2030	line += ', sig:%8.4f'%(sizeSig[0][2])
2031	if refine:
2032	print >>pFile,line
2033	else:
2034	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2035	if sizeSig[0][2]:
2036	refine = True
2037	line += ', sig:%8.4f'%(sizeSig[0][2])
2038	Snames = ['S11','S22','S33','S12','S13','S23']
2039	ptlbls = ' name :'
2040	ptstr = ' value :'
2041	sigstr = ' sig :'
2042	for i,name in enumerate(Snames):
2043	ptlbls += '%12s' % (name)
2044	ptstr += '%12.6f' % (hapData[4][i])
2045	if sizeSig[1][i]:
2046	refine = True
2047	sigstr += '%12.6f' % (sizeSig[1][i])
2048	else:
2049	sigstr += 12*' '
2050	if refine:
2051	print >>pFile,line
2052	print >>pFile,ptlbls
2053	print >>pFile,ptstr
2054	print >>pFile,sigstr
2055
2056	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2057	line = '\n Mustrain model: %9s'%(hapData[0])
2058	refine = False
2059	if hapData[0] in ['isotropic','uniaxial']:
2060	line += ' equatorial:%12.1f'%(hapData[1][0])
2061	if mustrainSig[0][0]:
2062	line += ', sig:%8.1f'%(mustrainSig[0][0])
2063	refine = True
2064	if hapData[0] == 'uniaxial':
2065	line += ' axial:%12.1f'%(hapData[1][1])
2066	if mustrainSig[0][1]:
2067	line += ', sig:%8.1f'%(mustrainSig[0][1])
2068	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2069	if mustrainSig[0][2]:
2070	refine = True
2071	line += ', sig:%8.3f'%(mustrainSig[0][2])
2072	if refine:
2073	print >>pFile,line
2074	else:
2075	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2076	if mustrainSig[0][2]:
2077	refine = True
2078	line += ', sig:%8.3f'%(mustrainSig[0][2])
2079	Snames = G2spc.MustrainNames(SGData)
2080	ptlbls = ' name :'
2081	ptstr = ' value :'
2082	sigstr = ' sig :'
2083	for i,name in enumerate(Snames):
2084	ptlbls += '%12s' % (name)
2085	ptstr += '%12.6f' % (hapData[4][i])
2086	if mustrainSig[1][i]:
2087	refine = True
2088	sigstr += '%12.6f' % (mustrainSig[1][i])
2089	else:
2090	sigstr += 12*' '
2091	if refine:
2092	print >>pFile,line
2093	print >>pFile,ptlbls
2094	print >>pFile,ptstr
2095	print >>pFile,sigstr
2096
2097	def PrintHStrainAndSig(hapData,strainSig,SGData):
2098	Hsnames = G2spc.HStrainNames(SGData)
2099	ptlbls = ' name :'
2100	ptstr = ' value :'
2101	sigstr = ' sig :'
2102	refine = False
2103	for i,name in enumerate(Hsnames):
2104	ptlbls += '%12s' % (name)
2105	ptstr += '%12.4g' % (hapData[0][i])
2106	if name in strainSig:
2107	refine = True
2108	sigstr += '%12.4g' % (strainSig[name])
2109	else:
2110	sigstr += 12*' '
2111	if refine:
2112	print >>pFile,'\n Hydrostatic/elastic strain: '
2113	print >>pFile,ptlbls
2114	print >>pFile,ptstr
2115	print >>pFile,sigstr
2116
2117	def PrintSHPOAndSig(pfx,hapData,POsig):
2118	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2119	ptlbls = ' names :'
2120	ptstr = ' values:'
2121	sigstr = ' sig :'
2122	for item in hapData[5]:
2123	ptlbls += '%12s'%(item)
2124	ptstr += '%12.3f'%(hapData[5][item])
2125	if pfx+item in POsig:
2126	sigstr += '%12.3f'%(POsig[pfx+item])
2127	else:
2128	sigstr += 12*' '
2129	print >>pFile,ptlbls
2130	print >>pFile,ptstr
2131	print >>pFile,sigstr
2132
2133	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2134	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2135	text = ''
2136	for item in hapData[2]:
2137	if pfx+item in ScalExtSig:
2138	text += ' %s: '%(item)
2139	text += '%12.2e'%(hapData[2][item][0])
2140	if pfx+item in ScalExtSig:
2141	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2142	print >>pFile,text
2143
2144	def PrintBabinetAndSig(pfx,hapData,BabSig):
2145	print >>pFile,'\n Babinet form factor modification: '
2146	ptlbls = ' names :'
2147	ptstr = ' values:'
2148	sigstr = ' sig :'
2149	for item in hapData:
2150	ptlbls += '%12s'%(item)
2151	ptstr += '%12.3f'%(hapData[item][0])
2152	if pfx+item in BabSig:
2153	sigstr += '%12.3f'%(BabSig[pfx+item])
2154	else:
2155	sigstr += 12*' '
2156	print >>pFile,ptlbls
2157	print >>pFile,ptstr
2158	print >>pFile,sigstr
2159
2160	PhFrExtPOSig = {}
2161	SizeMuStrSig = {}
2162	ScalExtSig = {}
2163	BabSig = {}
2164	wtFrSum = {}
2165	for phase in Phases:
2166	HistoPhase = Phases[phase]['Histograms']
2167	General = Phases[phase]['General']
2168	SGData = General['SGData']
2169	pId = Phases[phase]['pId']
2170	histoList = HistoPhase.keys()
2171	histoList.sort()
2172	for histogram in histoList:
2173	try:
2174	Histogram = Histograms[histogram]
2175	except KeyError:
2176	#skip if histogram not included e.g. in a sequential refinement
2177	continue
2178	hapData = HistoPhase[histogram]
2179	hId = Histogram['hId']
2180	pfx = str(pId)+':'+str(hId)+':'
2181	if hId not in wtFrSum:
2182	wtFrSum[hId] = 0.
2183	if 'PWDR' in histogram:
2184	for item in ['Scale','Extinction']:
2185	hapData[item][0] = parmDict[pfx+item]
2186	if pfx+item in sigDict:
2187	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2188	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2189	if hapData['Pref.Ori.'][0] == 'MD':
2190	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2191	if pfx+'MD' in sigDict:
2192	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2193	else: #'SH' spherical harmonics
2194	for item in hapData['Pref.Ori.'][5]:
2195	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2196	if pfx+item in sigDict:
2197	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2198	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2199	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2200	pfx+'HStrain':{}})
2201	for item in ['Mustrain','Size']:
2202	hapData[item][1][2] = parmDict[pfx+item+';mx']
2203	hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2204	if pfx+item+';mx' in sigDict:
2205	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2206	if hapData[item][0] in ['isotropic','uniaxial']:
2207	hapData[item][1][0] = parmDict[pfx+item+';i']
2208	if item == 'Size':
2209	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2210	if pfx+item+';i' in sigDict:
2211	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2212	if hapData[item][0] == 'uniaxial':
2213	hapData[item][1][1] = parmDict[pfx+item+';a']
2214	if item == 'Size':
2215	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2216	if pfx+item+';a' in sigDict:
2217	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2218	else: #generalized for mustrain or ellipsoidal for size
2219	Nterms = len(hapData[item][4])
2220	for i in range(Nterms):
2221	sfx = ':'+str(i)
2222	hapData[item][4][i] = parmDict[pfx+item+sfx]
2223	if pfx+item+sfx in sigDict:
2224	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2225	names = G2spc.HStrainNames(SGData)
2226	for i,name in enumerate(names):
2227	hapData['HStrain'][0][i] = parmDict[pfx+name]
2228	if pfx+name in sigDict:
2229	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2230	for name in ['BabA','BabU']:
2231	hapData['Babinet'][name][0] = parmDict[pfx+name]
2232	if pfx+name in sigDict:
2233	BabSig[pfx+name] = sigDict[pfx+name]
2234
2235	elif 'HKLF' in histogram:
2236	for item in ['Scale',]:
2237	if parmDict.get(pfx+item):
2238	hapData[item][0] = parmDict[pfx+item]
2239	if pfx+item in sigDict:
2240	ScalExtSig[pfx+item] = sigDict[pfx+item]
2241	for item in ['Ep','Eg','Es']:
2242	if parmDict.get(pfx+item):
2243	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2244	if pfx+item in sigDict:
2245	ScalExtSig[pfx+item] = sigDict[pfx+item]
2246	for name in ['BabA','BabU']:
2247	hapData['Babinet'][name][0] = parmDict[pfx+name]
2248	if pfx+name in sigDict:
2249	BabSig[pfx+name] = sigDict[pfx+name]
2250
2251	if Print:
2252	for phase in Phases:
2253	HistoPhase = Phases[phase]['Histograms']
2254	General = Phases[phase]['General']
2255	SGData = General['SGData']
2256	pId = Phases[phase]['pId']
2257	histoList = HistoPhase.keys()
2258	histoList.sort()
2259	for histogram in histoList:
2260	try:
2261	Histogram = Histograms[histogram]
2262	except KeyError:
2263	#skip if histogram not included e.g. in a sequential refinement
2264	continue
2265	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2266	print >>pFile,130*'-'
2267	hapData = HistoPhase[histogram]
2268	hId = Histogram['hId']
2269	Histogram['Residuals'][str(pId)+'::Name'] = phase
2270	pfx = str(pId)+':'+str(hId)+':'
2271	if 'PWDR' in histogram:
2272	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2273	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2274
2275	if pfx+'Scale' in PhFrExtPOSig:
2276	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2277	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2278	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2279	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2280	if pfx+'Extinction' in PhFrExtPOSig:
2281	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2282	if hapData['Pref.Ori.'][0] == 'MD':
2283	if pfx+'MD' in PhFrExtPOSig:
2284	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2285	else:
2286	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2287	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2288	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2289	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2290	if len(BabSig):
2291	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2292
2293	elif 'HKLF' in histogram:
2294	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2295	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2296	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2297	if pfx+'Scale' in ScalExtSig:
2298	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2299	if hapData['Extinction'][0] != 'None':
2300	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2301	if len(BabSig):
2302	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2303
2304	################################################################################
2305	##### Histogram data
2306	################################################################################
2307
2308	def GetHistogramData(Histograms,Print=True,pFile=None):
2309	'needs a doc string'
2310
2311	def GetBackgroundParms(hId,Background):
2312	Back = Background[0]
2313	DebyePeaks = Background[1]
2314	bakType,bakFlag = Back[:2]
2315	backVals = Back[3:]
2316	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2317	backDict = dict(zip(backNames,backVals))
2318	backVary = []
2319	if bakFlag:
2320	backVary = backNames
2321	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2322	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2323	debyeDict = {}
2324	debyeList = []
2325	for i in range(DebyePeaks['nDebye']):
2326	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2327	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2328	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2329	debyeVary = []
2330	for item in debyeList:
2331	if item[1]:
2332	debyeVary.append(item[0])
2333	backDict.update(debyeDict)
2334	backVary += debyeVary
2335	peakDict = {}
2336	peakList = []
2337	for i in range(DebyePeaks['nPeaks']):
2338	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2339	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2340	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2341	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2342	peakVary = []
2343	for item in peakList:
2344	if item[1]:
2345	peakVary.append(item[0])
2346	backDict.update(peakDict)
2347	backVary += peakVary
2348	return bakType,backDict,backVary
2349
2350	def GetInstParms(hId,Inst):
2351	dataType = Inst['Type'][0]
2352	instDict = {}
2353	insVary = []
2354	pfx = ':'+str(hId)+':'
2355	insKeys = Inst.keys()
2356	insKeys.sort()
2357	for item in insKeys:
2358	insName = pfx+item
2359	instDict[insName] = Inst[item][1]
2360	if len(Inst[item]) > 2 and Inst[item][2]:
2361	insVary.append(insName)
2362	if 'C' in dataType:
2363	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2364	return dataType,instDict,insVary
2365
2366	def GetSampleParms(hId,Sample):
2367	sampVary = []
2368	hfx = ':'+str(hId)+':'
2369	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2370	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2371	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2372	if key in Sample:
2373	sampDict[hfx+key] = Sample[key]
2374	Type = Sample['Type']
2375	if 'Bragg' in Type: #Bragg-Brentano
2376	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2377	sampDict[hfx+item] = Sample[item][0]
2378	if Sample[item][1]:
2379	sampVary.append(hfx+item)
2380	elif 'Debye' in Type: #Debye-Scherrer
2381	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2382	sampDict[hfx+item] = Sample[item][0]
2383	if Sample[item][1]:
2384	sampVary.append(hfx+item)
2385	return Type,sampDict,sampVary
2386
2387	def PrintBackground(Background):
2388	Back = Background[0]
2389	DebyePeaks = Background[1]
2390	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2391	line = ' Coefficients: '
2392	for i,back in enumerate(Back[3:]):
2393	line += '%10.3f'%(back)
2394	if i and not i%10:
2395	line += '\n'+15*' '
2396	print >>pFile,line
2397	if DebyePeaks['nDebye']:
2398	print >>pFile,'\n Debye diffuse scattering coefficients'
2399	parms = ['DebyeA','DebyeR','DebyeU']
2400	line = ' names : '
2401	for parm in parms:
2402	line += '%8s refine?'%(parm)
2403	print >>pFile,line
2404	for j,term in enumerate(DebyePeaks['debyeTerms']):
2405	line = ' term'+'%2d'%(j)+':'
2406	for i in range(3):
2407	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2408	print >>pFile,line
2409	if DebyePeaks['nPeaks']:
2410	print >>pFile,'\n Single peak coefficients'
2411	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2412	line = ' names : '
2413	for parm in parms:
2414	line += '%8s refine?'%(parm)
2415	print >>pFile,line
2416	for j,term in enumerate(DebyePeaks['peaksList']):
2417	line = ' peak'+'%2d'%(j)+':'
2418	for i in range(4):
2419	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2420	print >>pFile,line
2421
2422	def PrintInstParms(Inst):
2423	print >>pFile,'\n Instrument Parameters:'
2424	insKeys = Inst.keys()
2425	insKeys.sort()
2426	iBeg = 0
2427	Ok = True
2428	while Ok:
2429	ptlbls = ' name :'
2430	ptstr = ' value :'
2431	varstr = ' refine:'
2432	iFin = min(iBeg+9,len(insKeys))
2433	for item in insKeys[iBeg:iFin]:
2434	if item not in ['Type','Source']:
2435	ptlbls += '%12s' % (item)
2436	ptstr += '%12.6f' % (Inst[item][1])
2437	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2438	varstr += 12*' '
2439	else:
2440	varstr += '%12s' % (str(bool(Inst[item][2])))
2441	print >>pFile,ptlbls
2442	print >>pFile,ptstr
2443	print >>pFile,varstr
2444	iBeg = iFin
2445	if iBeg == len(insKeys):
2446	Ok = False
2447	else:
2448	print >>pFile,'\n'
2449
2450	def PrintSampleParms(Sample):
2451	print >>pFile,'\n Sample Parameters:'
2452	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2453	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2454	ptlbls = ' name :'
2455	ptstr = ' value :'
2456	varstr = ' refine:'
2457	if 'Bragg' in Sample['Type']:
2458	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2459	ptlbls += '%14s'%(item)
2460	ptstr += '%14.4f'%(Sample[item][0])
2461	varstr += '%14s'%(str(bool(Sample[item][1])))
2462
2463	elif 'Debye' in Type: #Debye-Scherrer
2464	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2465	ptlbls += '%14s'%(item)
2466	ptstr += '%14.4f'%(Sample[item][0])
2467	varstr += '%14s'%(str(bool(Sample[item][1])))
2468
2469	print >>pFile,ptlbls
2470	print >>pFile,ptstr
2471	print >>pFile,varstr
2472
2473	histDict = {}
2474	histVary = []
2475	controlDict = {}
2476	histoList = Histograms.keys()
2477	histoList.sort()
2478	for histogram in histoList:
2479	if 'PWDR' in histogram:
2480	Histogram = Histograms[histogram]
2481	hId = Histogram['hId']
2482	pfx = ':'+str(hId)+':'
2483	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2484	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2485	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2486	for excl in controlDict[pfx+'Exclude']:
2487	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2488	if controlDict[pfx+'Exclude']:
2489	ma.mask_rows(Histogram['Data'])
2490	Background = Histogram['Background']
2491	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2492	controlDict[pfx+'bakType'] = Type
2493	histDict.update(bakDict)
2494	histVary += bakVary
2495
2496	Inst = Histogram['Instrument Parameters'][0]
2497	Type,instDict,insVary = GetInstParms(hId,Inst)
2498	controlDict[pfx+'histType'] = Type
2499	if 'XC' in Type:
2500	if pfx+'Lam1' in instDict:
2501	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2502	else:
2503	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2504	histDict.update(instDict)
2505	histVary += insVary
2506
2507	Sample = Histogram['Sample Parameters']
2508	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2509	controlDict[pfx+'instType'] = Type
2510	histDict.update(sampDict)
2511	histVary += sampVary
2512
2513
2514	if Print:
2515	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2516	print >>pFile,135*'-'
2517	Units = {'C':' deg','T':' msec'}
2518	units = Units[controlDict[pfx+'histType'][2]]
2519	Limits = controlDict[pfx+'Limits']
2520	print >>pFile,' Instrument type: ',Sample['Type']
2521	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2522	if len(controlDict[pfx+'Exclude']):
2523	excls = controlDict[pfx+'Exclude']
2524	for excl in excls:
2525	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2526	PrintSampleParms(Sample)
2527	PrintInstParms(Inst)
2528	PrintBackground(Background)
2529	elif 'HKLF' in histogram:
2530	Histogram = Histograms[histogram]
2531	hId = Histogram['hId']
2532	pfx = ':'+str(hId)+':'
2533	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2534	Inst = Histogram['Instrument Parameters'][0]
2535	controlDict[pfx+'histType'] = Inst['Type'][0]
2536	if 'X' in Inst['Type'][0]:
2537	histDict[pfx+'Lam'] = Inst['Lam'][1]
2538	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2539	return histVary,histDict,controlDict
2540
2541	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2542	'needs a doc string'
2543
2544	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2545	Back = Background[0]
2546	DebyePeaks = Background[1]
2547	lenBack = len(Back[3:])
2548	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2549	for i in range(lenBack):
2550	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2551	if pfx+'Back;'+str(i) in sigDict:
2552	backSig[i] = sigDict[pfx+'Back;'+str(i)]
2553	if DebyePeaks['nDebye']:
2554	for i in range(DebyePeaks['nDebye']):
2555	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2556	for j,name in enumerate(names):
2557	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2558	if name in sigDict:
2559	backSig[lenBack+3*i+j] = sigDict[name]
2560	if DebyePeaks['nPeaks']:
2561	for i in range(DebyePeaks['nPeaks']):
2562	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2563	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2564	for j,name in enumerate(names):
2565	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2566	if name in sigDict:
2567	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2568	return backSig
2569
2570	def SetInstParms(pfx,Inst,parmDict,sigDict):
2571	instSig = {}
2572	insKeys = Inst.keys()
2573	insKeys.sort()
2574	for item in insKeys:
2575	insName = pfx+item
2576	Inst[item][1] = parmDict[insName]
2577	if insName in sigDict:
2578	instSig[item] = sigDict[insName]
2579	else:
2580	instSig[item] = 0
2581	return instSig
2582
2583	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2584	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2585	sampSig = [0 for i in range(5)]
2586	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2587	Sample[item][0] = parmDict[pfx+item]
2588	if pfx+item in sigDict:
2589	sampSig[i] = sigDict[pfx+item]
2590	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2591	sampSig = [0 for i in range(4)]
2592	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2593	Sample[item][0] = parmDict[pfx+item]
2594	if pfx+item in sigDict:
2595	sampSig[i] = sigDict[pfx+item]
2596	return sampSig
2597
2598	def PrintBackgroundSig(Background,backSig):
2599	Back = Background[0]
2600	DebyePeaks = Background[1]
2601	lenBack = len(Back[3:])
2602	valstr = ' value : '
2603	sigstr = ' sig : '
2604	refine = False
2605	for i,back in enumerate(Back[3:]):
2606	valstr += '%10.4g'%(back)
2607	if Back[1]:
2608	refine = True
2609	sigstr += '%10.4g'%(backSig[i])
2610	else:
2611	sigstr += 10*' '
2612	if refine:
2613	print >>pFile,'\n Background function: ',Back[0]
2614	print >>pFile,valstr
2615	print >>pFile,sigstr
2616	if DebyePeaks['nDebye']:
2617	ifAny = False
2618	ptfmt = "%12.3f"
2619	names = ' names :'
2620	ptstr = ' values:'
2621	sigstr = ' esds :'
2622	for item in sigDict:
2623	if 'Debye' in item:
2624	ifAny = True
2625	names += '%12s'%(item)
2626	ptstr += ptfmt%(parmDict[item])
2627	sigstr += ptfmt%(sigDict[item])
2628	if ifAny:
2629	print >>pFile,'\n Debye diffuse scattering coefficients'
2630	print >>pFile,names
2631	print >>pFile,ptstr
2632	print >>pFile,sigstr
2633	if DebyePeaks['nPeaks']:
2634	ifAny = False
2635	ptfmt = "%14.3f"
2636	names = ' names :'
2637	ptstr = ' values:'
2638	sigstr = ' esds :'
2639	for item in sigDict:
2640	if 'BkPk' in item:
2641	ifAny = True
2642	names += '%14s'%(item)
2643	ptstr += ptfmt%(parmDict[item])
2644	sigstr += ptfmt%(sigDict[item])
2645	if ifAny:
2646	print >>pFile,'\n Single peak coefficients'
2647	print >>pFile,names
2648	print >>pFile,ptstr
2649	print >>pFile,sigstr
2650
2651	def PrintInstParmsSig(Inst,instSig):
2652	refine = False
2653	insKeys = instSig.keys()
2654	insKeys.sort()
2655	iBeg = 0
2656	Ok = True
2657	while Ok:
2658	ptlbls = ' names :'
2659	ptstr = ' value :'
2660	sigstr = ' sig :'
2661	iFin = min(iBeg+9,len(insKeys))
2662	for name in insKeys[iBeg:iFin]:
2663	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2664	ptlbls += '%12s' % (name)
2665	ptstr += '%12.6f' % (Inst[name][1])
2666	if instSig[name]:
2667	refine = True
2668	sigstr += '%12.6f' % (instSig[name])
2669	else:
2670	sigstr += 12*' '
2671	if refine:
2672	print >>pFile,'\n Instrument Parameters:'
2673	print >>pFile,ptlbls
2674	print >>pFile,ptstr
2675	print >>pFile,sigstr
2676	iBeg = iFin
2677	if iBeg == len(insKeys):
2678	Ok = False
2679
2680	def PrintSampleParmsSig(Sample,sampleSig):
2681	ptlbls = ' names :'
2682	ptstr = ' values:'
2683	sigstr = ' sig :'
2684	refine = False
2685	if 'Bragg' in Sample['Type']:
2686	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2687	ptlbls += '%14s'%(item)
2688	ptstr += '%14.4f'%(Sample[item][0])
2689	if sampleSig[i]:
2690	refine = True
2691	sigstr += '%14.4f'%(sampleSig[i])
2692	else:
2693	sigstr += 14*' '
2694
2695	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2696	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2697	ptlbls += '%14s'%(item)
2698	ptstr += '%14.4f'%(Sample[item][0])
2699	if sampleSig[i]:
2700	refine = True
2701	sigstr += '%14.4f'%(sampleSig[i])
2702	else:
2703	sigstr += 14*' '
2704
2705	if refine:
2706	print >>pFile,'\n Sample Parameters:'
2707	print >>pFile,ptlbls
2708	print >>pFile,ptstr
2709	print >>pFile,sigstr
2710
2711	histoList = Histograms.keys()
2712	histoList.sort()
2713	for histogram in histoList:
2714	if 'PWDR' in histogram:
2715	Histogram = Histograms[histogram]
2716	hId = Histogram['hId']
2717	pfx = ':'+str(hId)+':'
2718	Background = Histogram['Background']
2719	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2720
2721	Inst = Histogram['Instrument Parameters'][0]
2722	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2723
2724	Sample = Histogram['Sample Parameters']
2725	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2726
2727	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2728	print >>pFile,135*'-'
2729	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2730	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2731	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2732	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2733	if Print:
2734	print >>pFile,' Instrument type: ',Sample['Type']
2735	PrintSampleParmsSig(Sample,sampSig)
2736	PrintInstParmsSig(Inst,instSig)
2737	PrintBackgroundSig(Background,backSig)
2738

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

Download in other formats: