Context Navigation

source: trunk/GSASIIstrIO.py @ 1625

Last change on this file since 1625 was 1625, checked in by vondreele, 9 years ago
add a parameter to result from G2stIO.GetPhaseData? add modulation functions to G2Math add modulation names to G2obj implement various wave types for modulations plot position modulation wave function on contour plots implement shift of modulation plot by +/-/0 keys temporarily default G2spc.GetSSfxuinel to 1,-1 site symm. move maxSSwave dict out of parmDict - now in controlDict implement import of Sawtooth parms from J2K files.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 124.2 KB

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

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

Download in other formats: