Context Navigation

source: trunk/GSASIIstrIO.py @ 1884

Last change on this file since 1884 was 1884, checked in by vondreele, 8 years ago
remove some unused imports add merohedral/pseudomerohedral Twin Laws to G2ddataGUI and G2strIO (not in G2strmath yet). allow ReImport? atoms to fill otherwise empty Atom List clarify HKL importers as Shelx HKL 4 & HKL 5 files.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 136.6 KB

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

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

Download in other formats: