Context Navigation

source: trunk/GSASIIstrIO.py @ 1876

Last change on this file since 1876 was 1876, checked in by vondreele, 8 years ago
Flack fix
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 134.7 KB

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

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

Download in other formats: