Context Navigation

source: trunk/GSASIIstrIO.py @ 1944

Last change on this file since 1944 was 1944, checked in by vondreele, 8 years ago
add G2HistDataDialog for editing e.g. orientation angles over all images - could be used elsewhere fix problem of bad RefList? key on plotting new histograms fix failure (file busy!) in GPXBackup
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 139.1 KB

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

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

Download in other formats: