Context Navigation

source: trunk/GSASIIstrIO.py @ 1887

Last change on this file since 1887 was 1887, checked in by vondreele, 10 years ago
default new twin to inversion twin - usual case fix twin results output twin functions now work - derivatives work but incorrect.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 137.0 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2015-06-12 20:13:05 +0000 (Fri, 12 Jun 2015) $
4	# $Author: vondreele $
5	# $Revision: 1887 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 1887 2015-06-12 20:13:05Z 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: 1887 $")
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	hapDict = {}
2030	hapVary = []
2031	controlDict = {}
2032	poType = {}
2033	poAxes = {}
2034	spAxes = {}
2035	spType = {}
2036
2037	for phase in Phases:
2038	HistoPhase = Phases[phase]['Histograms']
2039	SGData = Phases[phase]['General']['SGData']
2040	cell = Phases[phase]['General']['Cell'][1:7]
2041	A = G2lat.cell2A(cell)
2042	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
2043	SSGData = Phases[phase]['General']['SSGData']
2044	Vec,x,maxH = Phases[phase]['General']['SuperVec']
2045	pId = Phases[phase]['pId']
2046	histoList = HistoPhase.keys()
2047	histoList.sort()
2048	for histogram in histoList:
2049	try:
2050	Histogram = Histograms[histogram]
2051	except KeyError:
2052	#skip if histogram not included e.g. in a sequential refinement
2053	continue
2054	hapData = HistoPhase[histogram]
2055	hId = Histogram['hId']
2056	if 'PWDR' in histogram:
2057	limits = Histogram['Limits'][1]
2058	inst = Histogram['Instrument Parameters'][0]
2059	Zero = inst['Zero'][0]
2060	if 'C' in inst['Type'][1]:
2061	try:
2062	wave = inst['Lam'][1]
2063	except KeyError:
2064	wave = inst['Lam1'][1]
2065	dmin = wave/(2.0*sind(limits[1]/2.0))
2066	elif 'T' in inst['Type'][0]:
2067	dmin = limits[0]/inst['difC'][1]
2068	pfx = str(pId)+':'+str(hId)+':'
2069	for item in ['Scale','Extinction']:
2070	hapDict[pfx+item] = hapData[item][0]
2071	if hapData[item][1]:
2072	hapVary.append(pfx+item)
2073	names = G2spc.HStrainNames(SGData)
2074	HSvals = []
2075	for i,name in enumerate(names):
2076	hapDict[pfx+name] = hapData['HStrain'][0][i]
2077	HSvals.append(hapDict[pfx+name])
2078	if hapData['HStrain'][1][i]:
2079	hapVary.append(pfx+name)
2080	DIJS = G2spc.HStrainVals(HSvals,SGData)
2081	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2082	if hapData['Pref.Ori.'][0] == 'MD':
2083	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2084	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2085	if hapData['Pref.Ori.'][2]:
2086	hapVary.append(pfx+'MD')
2087	else: #'SH' spherical harmonics
2088	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2089	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2090	controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2091	controlDict[pfx+'SHhkl'] = []
2092	try: #patch for old Pref.Ori. items
2093	controlDict[pfx+'SHtoler'] = 0.1
2094	if hapData['Pref.Ori.'][6][0] != '':
2095	controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2096	controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2097	except IndexError:
2098	pass
2099	for item in hapData['Pref.Ori.'][5]:
2100	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2101	if hapData['Pref.Ori.'][2]:
2102	hapVary.append(pfx+item)
2103	for item in ['Mustrain','Size']:
2104	controlDict[pfx+item+'Type'] = hapData[item][0]
2105	hapDict[pfx+item+';mx'] = hapData[item][1][2]
2106	if hapData[item][2][2]:
2107	hapVary.append(pfx+item+';mx')
2108	if hapData[item][0] in ['isotropic','uniaxial']:
2109	hapDict[pfx+item+';i'] = hapData[item][1][0]
2110	if hapData[item][2][0]:
2111	hapVary.append(pfx+item+';i')
2112	if hapData[item][0] == 'uniaxial':
2113	controlDict[pfx+item+'Axis'] = hapData[item][3]
2114	hapDict[pfx+item+';a'] = hapData[item][1][1]
2115	if hapData[item][2][1]:
2116	hapVary.append(pfx+item+';a')
2117	else: #generalized for mustrain or ellipsoidal for size
2118	Nterms = len(hapData[item][4])
2119	if item == 'Mustrain':
2120	names = G2spc.MustrainNames(SGData)
2121	pwrs = []
2122	for name in names:
2123	h,k,l = name[1:]
2124	pwrs.append([int(h),int(k),int(l)])
2125	controlDict[pfx+'MuPwrs'] = pwrs
2126	for i in range(Nterms):
2127	sfx = ':'+str(i)
2128	hapDict[pfx+item+sfx] = hapData[item][4][i]
2129	if hapData[item][5][i]:
2130	hapVary.append(pfx+item+sfx)
2131	for bab in ['BabA','BabU']:
2132	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2133	if hapData['Babinet'][bab][1]:
2134	hapVary.append(pfx+bab)
2135
2136	if Print:
2137	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2138	print >>pFile,135*'-'
2139	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2140	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2141	if hapData['Pref.Ori.'][0] == 'MD':
2142	Ax = hapData['Pref.Ori.'][3]
2143	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2144	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2145	else: #'SH' for spherical harmonics
2146	PrintSHPO(hapData['Pref.Ori.'])
2147	print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2148	PrintSize(hapData['Size'])
2149	PrintMuStrain(hapData['Mustrain'],SGData)
2150	PrintHStrain(hapData['HStrain'],SGData)
2151	if hapData['Babinet']['BabA'][0]:
2152	PrintBabinet(hapData['Babinet'])
2153	if resetRefList:
2154	refList = []
2155	Uniq = []
2156	Phi = []
2157	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
2158	ifSuper = True
2159	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2160	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2161	for h,k,l,m,d in HKLd:
2162	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData) #is this right for SS refl.??
2163	mul *= 2 # for powder overlap of Friedel pairs
2164	if m or not ext:
2165	if 'C' in inst['Type'][0]:
2166	pos = G2lat.Dsp2pos(inst,d)
2167	if limits[0] < pos < limits[1]:
2168	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])
2169	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2170	Uniq.append(uniq)
2171	Phi.append(phi)
2172	elif 'T' in inst['Type'][0]:
2173	pos = G2lat.Dsp2pos(inst,d)
2174	if limits[0] < pos < limits[1]:
2175	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2176	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])
2177	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2178	Uniq.append(uniq)
2179	Phi.append(phi)
2180	else:
2181	ifSuper = False
2182	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2183	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2184	for h,k,l,d in HKLd:
2185	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2186	mul *= 2 # for powder overlap of Friedel pairs
2187	if ext:
2188	continue
2189	if 'C' in inst['Type'][0]:
2190	pos = G2lat.Dsp2pos(inst,d)
2191	if limits[0] < pos < limits[1]:
2192	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])
2193	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2194	Uniq.append(uniq)
2195	Phi.append(phi)
2196	elif 'T' in inst['Type'][0]:
2197	pos = G2lat.Dsp2pos(inst,d)
2198	if limits[0] < pos < limits[1]:
2199	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2200	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])
2201	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2202	Uniq.append(uniq)
2203	Phi.append(phi)
2204	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2205	elif 'HKLF' in histogram:
2206	inst = Histogram['Instrument Parameters'][0]
2207	hId = Histogram['hId']
2208	hfx = ':%d:'%(hId)
2209	for item in inst:
2210	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2211	hapDict[hfx+item] = inst[item][1]
2212	pfx = str(pId)+':'+str(hId)+':'
2213	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2214	if hapData['Scale'][1]:
2215	hapVary.append(pfx+'Scale')
2216
2217	extApprox,extType,extParms = hapData['Extinction']
2218	controlDict[pfx+'EType'] = extType
2219	controlDict[pfx+'EApprox'] = extApprox
2220	if 'C' in inst['Type'][0]:
2221	controlDict[pfx+'Tbar'] = extParms['Tbar']
2222	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2223	if 'Primary' in extType:
2224	Ekey = ['Ep',]
2225	elif 'I & II' in extType:
2226	Ekey = ['Eg','Es']
2227	elif 'Secondary Type II' == extType:
2228	Ekey = ['Es',]
2229	elif 'Secondary Type I' == extType:
2230	Ekey = ['Eg',]
2231	else: #'None'
2232	Ekey = []
2233	for eKey in Ekey:
2234	hapDict[pfx+eKey] = extParms[eKey][0]
2235	if extParms[eKey][1]:
2236	hapVary.append(pfx+eKey)
2237	for bab in ['BabA','BabU']:
2238	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2239	if hapData['Babinet'][bab][1]:
2240	hapVary.append(pfx+bab)
2241	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2242	if hapData.get('Flack',[0,False])[1]:
2243	hapVary.append(pfx+'Flack')
2244	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False]],])
2245	sumTwFr = 0.
2246	controlDict[pfx+'TwinLaw'] = []
2247	for it,twin in enumerate(Twins):
2248	controlDict[pfx+'TwinLaw'].append(twin[0])
2249	hapDict[pfx+'TwinFr;'+str(it)] = twin[1][0]
2250	if it:
2251	sumTwFr += twin[1][0]
2252	if twin[1][1]:
2253	hapVary.append(pfx+'TwinFr;'+str(it))
2254	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2255	if len(Twins) > 1: #force sum to unity
2256	hapDict[pfx+'TwinFr;0'] = 1.-sumTwFr
2257	if Print:
2258	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2259	print >>pFile,135*'-'
2260	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2261	if extType != 'None':
2262	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox)
2263	text = ' Parameters :'
2264	for eKey in Ekey:
2265	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2266	print >>pFile,text
2267	if hapData['Babinet']['BabA'][0]:
2268	PrintBabinet(hapData['Babinet'])
2269	if not SGData['SGInv']:
2270	print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2271	if len(Twins) > 1:
2272	for it,twin in enumerate(Twins):
2273	print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr;'+str(it)]),twin[1][1]
2274
2275	Histogram['Reflection Lists'] = phase
2276
2277	return hapVary,hapDict,controlDict
2278
2279	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2280	'needs a doc string'
2281
2282	def PrintSizeAndSig(hapData,sizeSig):
2283	line = '\n Size model: %9s'%(hapData[0])
2284	refine = False
2285	if hapData[0] in ['isotropic','uniaxial']:
2286	line += ' equatorial:%12.5f'%(hapData[1][0])
2287	if sizeSig[0][0]:
2288	line += ', sig:%8.4f'%(sizeSig[0][0])
2289	refine = True
2290	if hapData[0] == 'uniaxial':
2291	line += ' axial:%12.4f'%(hapData[1][1])
2292	if sizeSig[0][1]:
2293	refine = True
2294	line += ', sig:%8.4f'%(sizeSig[0][1])
2295	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2296	if sizeSig[0][2]:
2297	refine = True
2298	line += ', sig:%8.4f'%(sizeSig[0][2])
2299	if refine:
2300	print >>pFile,line
2301	else:
2302	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2303	if sizeSig[0][2]:
2304	refine = True
2305	line += ', sig:%8.4f'%(sizeSig[0][2])
2306	Snames = ['S11','S22','S33','S12','S13','S23']
2307	ptlbls = ' name :'
2308	ptstr = ' value :'
2309	sigstr = ' sig :'
2310	for i,name in enumerate(Snames):
2311	ptlbls += '%12s' % (name)
2312	ptstr += '%12.6f' % (hapData[4][i])
2313	if sizeSig[1][i]:
2314	refine = True
2315	sigstr += '%12.6f' % (sizeSig[1][i])
2316	else:
2317	sigstr += 12*' '
2318	if refine:
2319	print >>pFile,line
2320	print >>pFile,ptlbls
2321	print >>pFile,ptstr
2322	print >>pFile,sigstr
2323
2324	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2325	line = '\n Mustrain model: %9s'%(hapData[0])
2326	refine = False
2327	if hapData[0] in ['isotropic','uniaxial']:
2328	line += ' equatorial:%12.1f'%(hapData[1][0])
2329	if mustrainSig[0][0]:
2330	line += ', sig:%8.1f'%(mustrainSig[0][0])
2331	refine = True
2332	if hapData[0] == 'uniaxial':
2333	line += ' axial:%12.1f'%(hapData[1][1])
2334	if mustrainSig[0][1]:
2335	line += ', sig:%8.1f'%(mustrainSig[0][1])
2336	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2337	if mustrainSig[0][2]:
2338	refine = True
2339	line += ', sig:%8.3f'%(mustrainSig[0][2])
2340	if refine:
2341	print >>pFile,line
2342	else:
2343	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2344	if mustrainSig[0][2]:
2345	refine = True
2346	line += ', sig:%8.3f'%(mustrainSig[0][2])
2347	Snames = G2spc.MustrainNames(SGData)
2348	ptlbls = ' name :'
2349	ptstr = ' value :'
2350	sigstr = ' sig :'
2351	for i,name in enumerate(Snames):
2352	ptlbls += '%12s' % (name)
2353	ptstr += '%12.6f' % (hapData[4][i])
2354	if mustrainSig[1][i]:
2355	refine = True
2356	sigstr += '%12.6f' % (mustrainSig[1][i])
2357	else:
2358	sigstr += 12*' '
2359	if refine:
2360	print >>pFile,line
2361	print >>pFile,ptlbls
2362	print >>pFile,ptstr
2363	print >>pFile,sigstr
2364
2365	def PrintHStrainAndSig(hapData,strainSig,SGData):
2366	Hsnames = G2spc.HStrainNames(SGData)
2367	ptlbls = ' name :'
2368	ptstr = ' value :'
2369	sigstr = ' sig :'
2370	refine = False
2371	for i,name in enumerate(Hsnames):
2372	ptlbls += '%12s' % (name)
2373	ptstr += '%12.4g' % (hapData[0][i])
2374	if name in strainSig:
2375	refine = True
2376	sigstr += '%12.4g' % (strainSig[name])
2377	else:
2378	sigstr += 12*' '
2379	if refine:
2380	print >>pFile,'\n Hydrostatic/elastic strain: '
2381	print >>pFile,ptlbls
2382	print >>pFile,ptstr
2383	print >>pFile,sigstr
2384
2385	def PrintSHPOAndSig(pfx,hapData,POsig):
2386	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2387	ptlbls = ' names :'
2388	ptstr = ' values:'
2389	sigstr = ' sig :'
2390	for item in hapData[5]:
2391	ptlbls += '%12s'%(item)
2392	ptstr += '%12.3f'%(hapData[5][item])
2393	if pfx+item in POsig:
2394	sigstr += '%12.3f'%(POsig[pfx+item])
2395	else:
2396	sigstr += 12*' '
2397	print >>pFile,ptlbls
2398	print >>pFile,ptstr
2399	print >>pFile,sigstr
2400
2401	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2402	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2403	text = ''
2404	for item in hapData[2]:
2405	if pfx+item in ScalExtSig:
2406	text += ' %s: '%(item)
2407	text += '%12.2e'%(hapData[2][item][0])
2408	if pfx+item in ScalExtSig:
2409	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2410	print >>pFile,text
2411
2412	def PrintBabinetAndSig(pfx,hapData,BabSig):
2413	print >>pFile,'\n Babinet form factor modification: '
2414	ptlbls = ' names :'
2415	ptstr = ' values:'
2416	sigstr = ' sig :'
2417	for item in hapData:
2418	ptlbls += '%12s'%(item)
2419	ptstr += '%12.3f'%(hapData[item][0])
2420	if pfx+item in BabSig:
2421	sigstr += '%12.3f'%(BabSig[pfx+item])
2422	else:
2423	sigstr += 12*' '
2424	print >>pFile,ptlbls
2425	print >>pFile,ptstr
2426	print >>pFile,sigstr
2427
2428	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2429	print >>pFile,'\n Twin Law fractions : '
2430	ptlbls = ' names :'
2431	ptstr = ' values:'
2432	sigstr = ' sig :'
2433	for it,item in enumerate(twinData):
2434	ptlbls += ' twin #%d'%(it)
2435	ptstr += '%12.3f'%(item[1][0])
2436	if pfx+'TwinFr;'+str(it) in TwinSig:
2437	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr;'+str(it)])
2438	else:
2439	sigstr += 12*' '
2440	print >>pFile,ptlbls
2441	print >>pFile,ptstr
2442	print >>pFile,sigstr
2443
2444
2445	PhFrExtPOSig = {}
2446	SizeMuStrSig = {}
2447	ScalExtSig = {}
2448	BabSig = {}
2449	TwinFrSig = {}
2450	wtFrSum = {}
2451	for phase in Phases:
2452	HistoPhase = Phases[phase]['Histograms']
2453	General = Phases[phase]['General']
2454	SGData = General['SGData']
2455	pId = Phases[phase]['pId']
2456	histoList = HistoPhase.keys()
2457	histoList.sort()
2458	for histogram in histoList:
2459	try:
2460	Histogram = Histograms[histogram]
2461	except KeyError:
2462	#skip if histogram not included e.g. in a sequential refinement
2463	continue
2464	hapData = HistoPhase[histogram]
2465	hId = Histogram['hId']
2466	pfx = str(pId)+':'+str(hId)+':'
2467	if hId not in wtFrSum:
2468	wtFrSum[hId] = 0.
2469	if 'PWDR' in histogram:
2470	for item in ['Scale','Extinction']:
2471	hapData[item][0] = parmDict[pfx+item]
2472	if pfx+item in sigDict:
2473	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2474	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2475	if hapData['Pref.Ori.'][0] == 'MD':
2476	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2477	if pfx+'MD' in sigDict:
2478	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2479	else: #'SH' spherical harmonics
2480	for item in hapData['Pref.Ori.'][5]:
2481	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2482	if pfx+item in sigDict:
2483	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2484	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2485	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2486	pfx+'HStrain':{}})
2487	for item in ['Mustrain','Size']:
2488	hapData[item][1][2] = parmDict[pfx+item+';mx']
2489	hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2490	if pfx+item+';mx' in sigDict:
2491	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2492	if hapData[item][0] in ['isotropic','uniaxial']:
2493	hapData[item][1][0] = parmDict[pfx+item+';i']
2494	if item == 'Size':
2495	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2496	if pfx+item+';i' in sigDict:
2497	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2498	if hapData[item][0] == 'uniaxial':
2499	hapData[item][1][1] = parmDict[pfx+item+';a']
2500	if item == 'Size':
2501	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2502	if pfx+item+';a' in sigDict:
2503	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2504	else: #generalized for mustrain or ellipsoidal for size
2505	Nterms = len(hapData[item][4])
2506	for i in range(Nterms):
2507	sfx = ':'+str(i)
2508	hapData[item][4][i] = parmDict[pfx+item+sfx]
2509	if pfx+item+sfx in sigDict:
2510	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2511	names = G2spc.HStrainNames(SGData)
2512	for i,name in enumerate(names):
2513	hapData['HStrain'][0][i] = parmDict[pfx+name]
2514	if pfx+name in sigDict:
2515	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2516	for name in ['BabA','BabU']:
2517	hapData['Babinet'][name][0] = parmDict[pfx+name]
2518	if pfx+name in sigDict:
2519	BabSig[pfx+name] = sigDict[pfx+name]
2520
2521	elif 'HKLF' in histogram:
2522	for item in ['Scale','Flack']:
2523	if parmDict.get(pfx+item):
2524	hapData[item][0] = parmDict[pfx+item]
2525	if pfx+item in sigDict:
2526	ScalExtSig[pfx+item] = sigDict[pfx+item]
2527	for item in ['Ep','Eg','Es']:
2528	if parmDict.get(pfx+item):
2529	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2530	if pfx+item in sigDict:
2531	ScalExtSig[pfx+item] = sigDict[pfx+item]
2532	for item in ['BabA','BabU']:
2533	hapData['Babinet'][item][0] = parmDict[pfx+item]
2534	if pfx+item in sigDict:
2535	BabSig[pfx+item] = sigDict[pfx+item]
2536	if 'Twins' in hapData:
2537	it = 0
2538	sumTwFr = 0.
2539	while True:
2540	try:
2541	hapData['Twins'][it][1][0] = parmDict[pfx+'TwinFr;'+str(it)]
2542	if pfx+'TwinFr;'+str(it) in sigDict:
2543	TwinFrSig[pfx+'TwinFr;'+str(it)] = sigDict[pfx+'TwinFr;'+str(it)]
2544	if it:
2545	sumTwFr += hapData['Twins'][it][1][0]
2546	it += 1
2547	except KeyError:
2548	break
2549	hapData['Twins'][0][1][0] = 1.-sumTwFr
2550
2551	if Print:
2552	for phase in Phases:
2553	HistoPhase = Phases[phase]['Histograms']
2554	General = Phases[phase]['General']
2555	SGData = General['SGData']
2556	pId = Phases[phase]['pId']
2557	histoList = HistoPhase.keys()
2558	histoList.sort()
2559	for histogram in histoList:
2560	try:
2561	Histogram = Histograms[histogram]
2562	except KeyError:
2563	#skip if histogram not included e.g. in a sequential refinement
2564	continue
2565	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2566	print >>pFile,130*'-'
2567	hapData = HistoPhase[histogram]
2568	hId = Histogram['hId']
2569	Histogram['Residuals'][str(pId)+'::Name'] = phase
2570	pfx = str(pId)+':'+str(hId)+':'
2571	hfx = ':%s:'%(hId)
2572	if 'PWDR' in histogram:
2573	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2574	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2575	print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2576
2577	if pfx+'Scale' in PhFrExtPOSig:
2578	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2579	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2580	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2581	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2582	if pfx+'Extinction' in PhFrExtPOSig:
2583	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2584	if hapData['Pref.Ori.'][0] == 'MD':
2585	if pfx+'MD' in PhFrExtPOSig:
2586	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2587	else:
2588	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2589	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2590	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2591	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2592	if len(BabSig):
2593	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2594
2595	elif 'HKLF' in histogram:
2596	Inst = Histogram['Instrument Parameters'][0]
2597	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)' \
2598	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2599	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2600	if FFtables != None and 'N' not in Inst['Type'][0]:
2601	PrintFprime(FFtables,hfx,pFile)
2602	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2603	if pfx+'Scale' in ScalExtSig:
2604	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2605	if hapData['Extinction'][0] != 'None':
2606	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2607	if len(BabSig):
2608	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2609	if pfx+'Flack' in ScalExtSig:
2610	print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2611	if len(TwinFrSig):
2612	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2613
2614	################################################################################
2615	##### Histogram data
2616	################################################################################
2617
2618	def GetHistogramData(Histograms,Print=True,pFile=None):
2619	'needs a doc string'
2620
2621	def GetBackgroundParms(hId,Background):
2622	Back = Background[0]
2623	DebyePeaks = Background[1]
2624	bakType,bakFlag = Back[:2]
2625	backVals = Back[3:]
2626	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2627	backDict = dict(zip(backNames,backVals))
2628	backVary = []
2629	if bakFlag:
2630	backVary = backNames
2631	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2632	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2633	debyeDict = {}
2634	debyeList = []
2635	for i in range(DebyePeaks['nDebye']):
2636	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2637	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2638	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2639	debyeVary = []
2640	for item in debyeList:
2641	if item[1]:
2642	debyeVary.append(item[0])
2643	backDict.update(debyeDict)
2644	backVary += debyeVary
2645	peakDict = {}
2646	peakList = []
2647	for i in range(DebyePeaks['nPeaks']):
2648	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2649	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2650	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2651	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2652	peakVary = []
2653	for item in peakList:
2654	if item[1]:
2655	peakVary.append(item[0])
2656	backDict.update(peakDict)
2657	backVary += peakVary
2658	return bakType,backDict,backVary
2659
2660	def GetInstParms(hId,Inst):
2661	dataType = Inst['Type'][0]
2662	instDict = {}
2663	insVary = []
2664	pfx = ':'+str(hId)+':'
2665	insKeys = Inst.keys()
2666	insKeys.sort()
2667	for item in insKeys:
2668	insName = pfx+item
2669	instDict[insName] = Inst[item][1]
2670	if len(Inst[item]) > 2 and Inst[item][2]:
2671	insVary.append(insName)
2672	if 'C' in dataType:
2673	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2674	return dataType,instDict,insVary
2675
2676	def GetSampleParms(hId,Sample):
2677	sampVary = []
2678	hfx = ':'+str(hId)+':'
2679	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2680	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2681	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2682	if key in Sample:
2683	sampDict[hfx+key] = Sample[key]
2684	Type = Sample['Type']
2685	if 'Bragg' in Type: #Bragg-Brentano
2686	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2687	sampDict[hfx+item] = Sample[item][0]
2688	if Sample[item][1]:
2689	sampVary.append(hfx+item)
2690	elif 'Debye' in Type: #Debye-Scherrer
2691	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2692	sampDict[hfx+item] = Sample[item][0]
2693	if Sample[item][1]:
2694	sampVary.append(hfx+item)
2695	return Type,sampDict,sampVary
2696
2697	def PrintBackground(Background):
2698	Back = Background[0]
2699	DebyePeaks = Background[1]
2700	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2701	line = ' Coefficients: '
2702	for i,back in enumerate(Back[3:]):
2703	line += '%10.3f'%(back)
2704	if i and not i%10:
2705	line += '\n'+15*' '
2706	print >>pFile,line
2707	if DebyePeaks['nDebye']:
2708	print >>pFile,'\n Debye diffuse scattering coefficients'
2709	parms = ['DebyeA','DebyeR','DebyeU']
2710	line = ' names : '
2711	for parm in parms:
2712	line += '%8s refine?'%(parm)
2713	print >>pFile,line
2714	for j,term in enumerate(DebyePeaks['debyeTerms']):
2715	line = ' term'+'%2d'%(j)+':'
2716	for i in range(3):
2717	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2718	print >>pFile,line
2719	if DebyePeaks['nPeaks']:
2720	print >>pFile,'\n Single peak coefficients'
2721	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2722	line = ' names : '
2723	for parm in parms:
2724	line += '%8s refine?'%(parm)
2725	print >>pFile,line
2726	for j,term in enumerate(DebyePeaks['peaksList']):
2727	line = ' peak'+'%2d'%(j)+':'
2728	for i in range(4):
2729	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2730	print >>pFile,line
2731
2732	def PrintInstParms(Inst):
2733	print >>pFile,'\n Instrument Parameters:'
2734	insKeys = Inst.keys()
2735	insKeys.sort()
2736	iBeg = 0
2737	Ok = True
2738	while Ok:
2739	ptlbls = ' name :'
2740	ptstr = ' value :'
2741	varstr = ' refine:'
2742	iFin = min(iBeg+9,len(insKeys))
2743	for item in insKeys[iBeg:iFin]:
2744	if item not in ['Type','Source']:
2745	ptlbls += '%12s' % (item)
2746	ptstr += '%12.6f' % (Inst[item][1])
2747	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2748	varstr += 12*' '
2749	else:
2750	varstr += '%12s' % (str(bool(Inst[item][2])))
2751	print >>pFile,ptlbls
2752	print >>pFile,ptstr
2753	print >>pFile,varstr
2754	iBeg = iFin
2755	if iBeg == len(insKeys):
2756	Ok = False
2757	else:
2758	print >>pFile,'\n'
2759
2760	def PrintSampleParms(Sample):
2761	print >>pFile,'\n Sample Parameters:'
2762	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2763	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2764	ptlbls = ' name :'
2765	ptstr = ' value :'
2766	varstr = ' refine:'
2767	if 'Bragg' in Sample['Type']:
2768	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2769	ptlbls += '%14s'%(item)
2770	ptstr += '%14.4f'%(Sample[item][0])
2771	varstr += '%14s'%(str(bool(Sample[item][1])))
2772
2773	elif 'Debye' in Type: #Debye-Scherrer
2774	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2775	ptlbls += '%14s'%(item)
2776	ptstr += '%14.4f'%(Sample[item][0])
2777	varstr += '%14s'%(str(bool(Sample[item][1])))
2778
2779	print >>pFile,ptlbls
2780	print >>pFile,ptstr
2781	print >>pFile,varstr
2782
2783	histDict = {}
2784	histVary = []
2785	controlDict = {}
2786	histoList = Histograms.keys()
2787	histoList.sort()
2788	for histogram in histoList:
2789	if 'PWDR' in histogram:
2790	Histogram = Histograms[histogram]
2791	hId = Histogram['hId']
2792	pfx = ':'+str(hId)+':'
2793	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2794	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2795	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2796	for excl in controlDict[pfx+'Exclude']:
2797	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2798	if controlDict[pfx+'Exclude']:
2799	ma.mask_rows(Histogram['Data'])
2800	Background = Histogram['Background']
2801	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2802	controlDict[pfx+'bakType'] = Type
2803	histDict.update(bakDict)
2804	histVary += bakVary
2805
2806	Inst = Histogram['Instrument Parameters'][0]
2807	Type,instDict,insVary = GetInstParms(hId,Inst)
2808	controlDict[pfx+'histType'] = Type
2809	if 'XC' in Type:
2810	if pfx+'Lam1' in instDict:
2811	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2812	else:
2813	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2814	histDict.update(instDict)
2815	histVary += insVary
2816
2817	Sample = Histogram['Sample Parameters']
2818	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2819	controlDict[pfx+'instType'] = Type
2820	histDict.update(sampDict)
2821	histVary += sampVary
2822
2823
2824	if Print:
2825	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2826	print >>pFile,135*'-'
2827	Units = {'C':' deg','T':' msec'}
2828	units = Units[controlDict[pfx+'histType'][2]]
2829	Limits = controlDict[pfx+'Limits']
2830	print >>pFile,' Instrument type: ',Sample['Type']
2831	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2832	if len(controlDict[pfx+'Exclude']):
2833	excls = controlDict[pfx+'Exclude']
2834	for excl in excls:
2835	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2836	PrintSampleParms(Sample)
2837	PrintInstParms(Inst)
2838	PrintBackground(Background)
2839	elif 'HKLF' in histogram:
2840	Histogram = Histograms[histogram]
2841	hId = Histogram['hId']
2842	pfx = ':'+str(hId)+':'
2843	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2844	Inst = Histogram['Instrument Parameters'][0]
2845	controlDict[pfx+'histType'] = Inst['Type'][0]
2846	if 'X' in Inst['Type'][0]:
2847	histDict[pfx+'Lam'] = Inst['Lam'][1]
2848	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2849	return histVary,histDict,controlDict
2850
2851	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
2852	'needs a doc string'
2853
2854	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2855	Back = Background[0]
2856	DebyePeaks = Background[1]
2857	lenBack = len(Back[3:])
2858	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2859	for i in range(lenBack):
2860	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2861	if pfx+'Back;'+str(i) in sigDict:
2862	backSig[i] = sigDict[pfx+'Back;'+str(i)]
2863	if DebyePeaks['nDebye']:
2864	for i in range(DebyePeaks['nDebye']):
2865	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
2866	for j,name in enumerate(names):
2867	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2868	if name in sigDict:
2869	backSig[lenBack+3*i+j] = sigDict[name]
2870	if DebyePeaks['nPeaks']:
2871	for i in range(DebyePeaks['nPeaks']):
2872	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2873	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2874	for j,name in enumerate(names):
2875	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2876	if name in sigDict:
2877	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2878	return backSig
2879
2880	def SetInstParms(pfx,Inst,parmDict,sigDict):
2881	instSig = {}
2882	insKeys = Inst.keys()
2883	insKeys.sort()
2884	for item in insKeys:
2885	insName = pfx+item
2886	Inst[item][1] = parmDict[insName]
2887	if insName in sigDict:
2888	instSig[item] = sigDict[insName]
2889	else:
2890	instSig[item] = 0
2891	return instSig
2892
2893	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2894	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2895	sampSig = [0 for i in range(5)]
2896	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2897	Sample[item][0] = parmDict[pfx+item]
2898	if pfx+item in sigDict:
2899	sampSig[i] = sigDict[pfx+item]
2900	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2901	sampSig = [0 for i in range(4)]
2902	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2903	Sample[item][0] = parmDict[pfx+item]
2904	if pfx+item in sigDict:
2905	sampSig[i] = sigDict[pfx+item]
2906	return sampSig
2907
2908	def PrintBackgroundSig(Background,backSig):
2909	Back = Background[0]
2910	DebyePeaks = Background[1]
2911	lenBack = len(Back[3:])
2912	valstr = ' value : '
2913	sigstr = ' sig : '
2914	refine = False
2915	for i,back in enumerate(Back[3:]):
2916	valstr += '%10.4g'%(back)
2917	if Back[1]:
2918	refine = True
2919	sigstr += '%10.4g'%(backSig[i])
2920	else:
2921	sigstr += 10*' '
2922	if refine:
2923	print >>pFile,'\n Background function: ',Back[0]
2924	print >>pFile,valstr
2925	print >>pFile,sigstr
2926	if DebyePeaks['nDebye']:
2927	ifAny = False
2928	ptfmt = "%12.3f"
2929	names = ' names :'
2930	ptstr = ' values:'
2931	sigstr = ' esds :'
2932	for item in sigDict:
2933	if 'Debye' in item:
2934	ifAny = True
2935	names += '%12s'%(item)
2936	ptstr += ptfmt%(parmDict[item])
2937	sigstr += ptfmt%(sigDict[item])
2938	if ifAny:
2939	print >>pFile,'\n Debye diffuse scattering coefficients'
2940	print >>pFile,names
2941	print >>pFile,ptstr
2942	print >>pFile,sigstr
2943	if DebyePeaks['nPeaks']:
2944	print >>pFile,'\n Single peak coefficients:'
2945	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2946	line = ' peak no. '
2947	for parm in parms:
2948	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
2949	print >>pFile,line
2950	for ip in range(DebyePeaks['nPeaks']):
2951	ptstr = ' %4d '%(ip)
2952	for parm in parms:
2953	fmt = '%14.3f'
2954	efmt = '%12.3f'
2955	if parm == 'BkPkpos':
2956	fmt = '%14.4f'
2957	efmt = '%12.4f'
2958	name = pfx+parm+';%d'%(ip)
2959	ptstr += fmt%(parmDict[name])
2960	if name in sigDict:
2961	ptstr += efmt%(sigDict[name])
2962	else:
2963	ptstr += 12*' '
2964	print >>pFile,ptstr
2965	sumBk = np.array(Histogram['sumBk'])
2966	print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g' \
2967	%(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
2968
2969	def PrintInstParmsSig(Inst,instSig):
2970	refine = False
2971	insKeys = instSig.keys()
2972	insKeys.sort()
2973	iBeg = 0
2974	Ok = True
2975	while Ok:
2976	ptlbls = ' names :'
2977	ptstr = ' value :'
2978	sigstr = ' sig :'
2979	iFin = min(iBeg+9,len(insKeys))
2980	for name in insKeys[iBeg:iFin]:
2981	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2982	ptlbls += '%12s' % (name)
2983	ptstr += '%12.6f' % (Inst[name][1])
2984	if instSig[name]:
2985	refine = True
2986	sigstr += '%12.6f' % (instSig[name])
2987	else:
2988	sigstr += 12*' '
2989	if refine:
2990	print >>pFile,'\n Instrument Parameters:'
2991	print >>pFile,ptlbls
2992	print >>pFile,ptstr
2993	print >>pFile,sigstr
2994	iBeg = iFin
2995	if iBeg == len(insKeys):
2996	Ok = False
2997
2998	def PrintSampleParmsSig(Sample,sampleSig):
2999	ptlbls = ' names :'
3000	ptstr = ' values:'
3001	sigstr = ' sig :'
3002	refine = False
3003	if 'Bragg' in Sample['Type']:
3004	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3005	ptlbls += '%14s'%(item)
3006	ptstr += '%14.4f'%(Sample[item][0])
3007	if sampleSig[i]:
3008	refine = True
3009	sigstr += '%14.4f'%(sampleSig[i])
3010	else:
3011	sigstr += 14*' '
3012
3013	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3014	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3015	ptlbls += '%14s'%(item)
3016	ptstr += '%14.4f'%(Sample[item][0])
3017	if sampleSig[i]:
3018	refine = True
3019	sigstr += '%14.4f'%(sampleSig[i])
3020	else:
3021	sigstr += 14*' '
3022
3023	if refine:
3024	print >>pFile,'\n Sample Parameters:'
3025	print >>pFile,ptlbls
3026	print >>pFile,ptstr
3027	print >>pFile,sigstr
3028
3029	histoList = Histograms.keys()
3030	histoList.sort()
3031	for histogram in histoList:
3032	if 'PWDR' in histogram:
3033	Histogram = Histograms[histogram]
3034	hId = Histogram['hId']
3035	pfx = ':'+str(hId)+':'
3036	Background = Histogram['Background']
3037	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3038
3039	Inst = Histogram['Instrument Parameters'][0]
3040	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3041
3042	Sample = Histogram['Sample Parameters']
3043	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3044
3045	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3046	print >>pFile,135*'-'
3047	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3048	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3049	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3050	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3051	if Print:
3052	print >>pFile,' Instrument type: ',Sample['Type']
3053	if FFtables != None and 'N' not in Inst['Type'][0]:
3054	PrintFprime(FFtables,pfx,pFile)
3055	PrintSampleParmsSig(Sample,sampSig)
3056	PrintInstParmsSig(Inst,instSig)
3057	PrintBackgroundSig(Background,backSig)
3058

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

Download in other formats: