Context Navigation

source: trunk/GSASIIstrIO.py @ 1896

Last change on this file since 1896 was 1896, checked in by vondreele, 10 years ago
add auto constraint for TwinFr? sum = 1.0 fix TwinFr? derivative - now OK & works
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 137.9 KB

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

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

Download in other formats: