Context Navigation

source: trunk/GSASIIstrIO.py @ 2478

Last change on this file since 2478 was 2478, checked in by vondreele, 7 years ago
work up to magnetic structure factor calcs.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 143.8 KB

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

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

Download in other formats: