Context Navigation

source: trunk/GSASIIstrIO.py @ 2155

Last change on this file since 2155 was 2155, checked in by vondreele, 9 years ago
fix bug - Mustrain: should be Mustrain; everywhere
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 139.7 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2016-02-23 16:20:56 +0000 (Tue, 23 Feb 2016) $
4	# $Author: vondreele $
5	# $Revision: 2155 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 2155 2016-02-23 16:20:56Z 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: 2155 $")
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','Tmin','Tmax','Xmax','Ymax','Zmax'],
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	for iw,wave in enumerate(Waves):
931	line = ''
932	if AtomSS['waveType'] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
933	for item in names[Stype][6:]:
934	line += '%8s '%(item)
935	else:
936	if Stype == 'Spos':
937	for item in names[Stype][:6]:
938	line += '%8s '%(item)
939	else:
940	for item in names[Stype]:
941	line += '%8s '%(item)
942	print >>pFile,line
943	line = ''
944	for item in wave[0]:
945	line += '%8.4f '%(item)
946	line += ' Refine? '+str(wave[1])
947	print >>pFile,line
948
949	def PrintTexture(textureData):
950	topstr = '\n Spherical harmonics texture: Order:' + \
951	str(textureData['Order'])
952	if textureData['Order']:
953	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
954	else:
955	print >>pFile,topstr
956	return
957	names = ['omega','chi','phi']
958	line = '\n'
959	for name in names:
960	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
961	print >>pFile,line
962	print >>pFile,'\n Texture coefficients:'
963	SHcoeff = textureData['SH Coeff'][1]
964	SHkeys = SHcoeff.keys()
965	nCoeff = len(SHcoeff)
966	nBlock = nCoeff/10+1
967	iBeg = 0
968	iFin = min(iBeg+10,nCoeff)
969	for block in range(nBlock):
970	ptlbls = ' names :'
971	ptstr = ' values:'
972	for item in SHkeys[iBeg:iFin]:
973	ptlbls += '%12s'%(item)
974	ptstr += '%12.4f'%(SHcoeff[item])
975	print >>pFile,ptlbls
976	print >>pFile,ptstr
977	iBeg += 10
978	iFin = min(iBeg+10,nCoeff)
979
980	def MakeRBParms(rbKey,phaseVary,phaseDict):
981	rbid = str(rbids.index(RB['RBId']))
982	pfxRB = pfx+'RB'+rbKey+'P'
983	pstr = ['x','y','z']
984	ostr = ['a','i','j','k']
985	for i in range(3):
986	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
987	phaseDict[name] = RB['Orig'][0][i]
988	if RB['Orig'][1]:
989	phaseVary += [name,]
990	pfxRB = pfx+'RB'+rbKey+'O'
991	for i in range(4):
992	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
993	phaseDict[name] = RB['Orient'][0][i]
994	if RB['Orient'][1] == 'AV' and i:
995	phaseVary += [name,]
996	elif RB['Orient'][1] == 'A' and not i:
997	phaseVary += [name,]
998
999	def MakeRBThermals(rbKey,phaseVary,phaseDict):
1000	rbid = str(rbids.index(RB['RBId']))
1001	tlstr = ['11','22','33','12','13','23']
1002	sstr = ['12','13','21','23','31','32','AA','BB']
1003	if 'T' in RB['ThermalMotion'][0]:
1004	pfxRB = pfx+'RB'+rbKey+'T'
1005	for i in range(6):
1006	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1007	phaseDict[name] = RB['ThermalMotion'][1][i]
1008	if RB['ThermalMotion'][2][i]:
1009	phaseVary += [name,]
1010	if 'L' in RB['ThermalMotion'][0]:
1011	pfxRB = pfx+'RB'+rbKey+'L'
1012	for i in range(6):
1013	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1014	phaseDict[name] = RB['ThermalMotion'][1][i+6]
1015	if RB['ThermalMotion'][2][i+6]:
1016	phaseVary += [name,]
1017	if 'S' in RB['ThermalMotion'][0]:
1018	pfxRB = pfx+'RB'+rbKey+'S'
1019	for i in range(8):
1020	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1021	phaseDict[name] = RB['ThermalMotion'][1][i+12]
1022	if RB['ThermalMotion'][2][i+12]:
1023	phaseVary += [name,]
1024	if 'U' in RB['ThermalMotion'][0]:
1025	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1026	phaseDict[name] = RB['ThermalMotion'][1][0]
1027	if RB['ThermalMotion'][2][0]:
1028	phaseVary += [name,]
1029
1030	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1031	rbid = str(rbids.index(RB['RBId']))
1032	pfxRB = pfx+'RB'+rbKey+'Tr;'
1033	for i,torsion in enumerate(RB['Torsions']):
1034	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1035	phaseDict[name] = torsion[0]
1036	if torsion[1]:
1037	phaseVary += [name,]
1038
1039	if Print:
1040	print >>pFile,'\n Phases:'
1041	phaseVary = []
1042	phaseDict = {}
1043	phaseConstr = {}
1044	pawleyLookup = {}
1045	FFtables = {} #scattering factors - xrays
1046	BLtables = {} # neutrons
1047	Natoms = {}
1048	AtMults = {}
1049	AtIA = {}
1050	maxSSwave = {}
1051	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1052	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1053	atomIndx = {}
1054	for name in PhaseData:
1055	General = PhaseData[name]['General']
1056	pId = PhaseData[name]['pId']
1057	pfx = str(pId)+'::'
1058	FFtable = G2el.GetFFtable(General['AtomTypes'])
1059	BLtable = G2el.GetBLtable(General)
1060	FFtables.update(FFtable)
1061	BLtables.update(BLtable)
1062	Atoms = PhaseData[name]['Atoms']
1063	if Atoms and not General.get('doPawley'):
1064	cx,ct,cs,cia = General['AtomPtrs']
1065	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1066	PawleyRef = PhaseData[name].get('Pawley ref',[])
1067	SGData = General['SGData']
1068	SGtext,SGtable = G2spc.SGPrint(SGData)
1069	cell = General['Cell']
1070	A = G2lat.cell2A(cell[1:7])
1071	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1072	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1073	if cell[0]:
1074	phaseVary += cellVary(pfx,SGData) #also fills in symmetry required constraints
1075	SSGtext = [] #no superstructure
1076	im = 0
1077	if General['Type'] in ['modulated','magnetic']:
1078	im = 1 #refl offset
1079	Vec,vRef,maxH = General['SuperVec']
1080	phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1081	SSGData = General['SSGData']
1082	SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1083	if vRef:
1084	phaseVary += modVary(pfx,SSGData)
1085	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1086	if resRBData:
1087	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
1088	for iRB,RB in enumerate(resRBData):
1089	MakeRBParms('R',phaseVary,phaseDict)
1090	MakeRBThermals('R',phaseVary,phaseDict)
1091	MakeRBTorsions('R',phaseVary,phaseDict)
1092
1093	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1094	if vecRBData:
1095	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
1096	for iRB,RB in enumerate(vecRBData):
1097	MakeRBParms('V',phaseVary,phaseDict)
1098	MakeRBThermals('V',phaseVary,phaseDict)
1099
1100	Natoms[pfx] = 0
1101	maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1102	if Atoms and not General.get('doPawley'):
1103	cx,ct,cs,cia = General['AtomPtrs']
1104	Natoms[pfx] = len(Atoms)
1105	for i,at in enumerate(Atoms):
1106	atomIndx[at[cia+8]] = [pfx,i] #lookup table for restraints
1107	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1108	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1109	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
1110	pfx+'AI/A:'+str(i):at[cia],})
1111	if at[cia] == 'I':
1112	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1113	else:
1114	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1115	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1116	if 'F' in at[ct+1]:
1117	phaseVary.append(pfx+'Afrac:'+str(i))
1118	if 'X' in at[ct+1]:
1119	try: #patch for sytsym name changes
1120	xId,xCoef = G2spc.GetCSxinel(at[cs])
1121	except KeyError:
1122	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1123	at[cs] = Sytsym
1124	xId,xCoef = G2spc.GetCSxinel(at[cs])
1125	xId,xCoef = G2spc.GetCSxinel(at[cs])
1126	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1127	equivs = [[],[],[]]
1128	for j in range(3):
1129	if xId[j] > 0:
1130	phaseVary.append(names[j])
1131	equivs[xId[j]-1].append([names[j],xCoef[j]])
1132	for equiv in equivs:
1133	if len(equiv) > 1:
1134	name = equiv[0][0]
1135	coef = equiv[0][1]
1136	for eqv in equiv[1:]:
1137	eqv[1] /= coef
1138	G2mv.StoreEquivalence(name,(eqv,))
1139	if 'U' in at[ct+1]:
1140	if at[cia] == 'I':
1141	phaseVary.append(pfx+'AUiso:'+str(i))
1142	else:
1143	try: #patch for sytsym name changes
1144	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1145	except KeyError:
1146	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1147	at[cs] = Sytsym
1148	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1149	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1150	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1151	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1152	equivs = [[],[],[],[],[],[]]
1153	for j in range(6):
1154	if uId[j] > 0:
1155	phaseVary.append(names[j])
1156	equivs[uId[j]-1].append([names[j],uCoef[j]])
1157	for equiv in equivs:
1158	if len(equiv) > 1:
1159	name = equiv[0][0]
1160	coef = equiv[0][1]
1161	for eqv in equiv[1:]:
1162	eqv[1] /= coef
1163	G2mv.StoreEquivalence(name,equiv[1:])
1164	if General['Type'] in ['modulated','magnetic']:
1165	AtomSS = at[-1]['SS1']
1166	waveType = AtomSS['waveType']
1167	phaseDict[pfx+'waveType:'+str(i)] = waveType
1168	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1169	Waves = AtomSS[Stype]
1170	nx = 0
1171	for iw,wave in enumerate(Waves):
1172	if not iw:
1173	if waveType in ['ZigZag','Block']:
1174	nx = 1
1175	CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1176	else:
1177	CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1178	uId,uCoef = CSI[Stype]
1179	stiw = str(i)+':'+str(iw)
1180	if Stype == 'Spos':
1181	if waveType in ['ZigZag','Block',] and not iw:
1182	names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1183	equivs = [[],[], [],[],[]]
1184	else:
1185	names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1186	pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1187	equivs = [[],[],[], [],[],[]]
1188	elif Stype == 'Sadp':
1189	names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1190	pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1191	pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1192	pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1193	equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1194	elif Stype == 'Sfrac':
1195	equivs = [[],[]]
1196	if 'Crenel' in waveType and not iw:
1197	names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1198	else:
1199	names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1200	elif Stype == 'Smag':
1201	equivs = [[],[],[], [],[],[]]
1202	names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1203	pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1204	phaseDict.update(dict(zip(names,wave[0])))
1205	if wave[1]: #what do we do here for multiple terms in modulation constraints?
1206	for j in range(len(equivs)):
1207	if uId[j][0] > 0:
1208	phaseVary.append(names[j])
1209	equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1210	for equiv in equivs:
1211	if len(equiv) > 1:
1212	name = equiv[0][0]
1213	coef = equiv[0][1]
1214	for eqv in equiv[1:]:
1215	eqv[1] /= coef
1216	G2mv.StoreEquivalence(name,equiv[1:])
1217	maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1218	textureData = General['SH Texture']
1219	if textureData['Order'] and not seqRef:
1220	phaseDict[pfx+'SHorder'] = textureData['Order']
1221	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1222	for item in ['omega','chi','phi']:
1223	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1224	if textureData['Sample '+item][0]:
1225	phaseVary.append(pfx+'SH '+item)
1226	for item in textureData['SH Coeff'][1]:
1227	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1228	if textureData['SH Coeff'][0]:
1229	phaseVary.append(pfx+item)
1230
1231	if Print:
1232	print >>pFile,'\n Phase name: ',General['Name']
1233	print >>pFile,135*'-'
1234	PrintFFtable(FFtable)
1235	PrintBLtable(BLtable)
1236	print >>pFile,''
1237	if len(SSGtext): #if superstructure
1238	for line in SSGtext: print >>pFile,line
1239	if len(SSGtable):
1240	for item in SSGtable:
1241	line = ' %s '%(item)
1242	print >>pFile,line
1243	else:
1244	print >>pFile,' ( 1) %s'%(SSGtable[0])
1245	else:
1246	for line in SGtext: print >>pFile,line
1247	if len(SGtable):
1248	for item in SGtable:
1249	line = ' %s '%(item)
1250	print >>pFile,line
1251	else:
1252	print >>pFile,' ( 1) %s'%(SGtable[0])
1253	PrintRBObjects(resRBData,vecRBData)
1254	PrintAtoms(General,Atoms)
1255	if General['Type'] in ['modulated','magnetic']:
1256	PrintWaves(General,Atoms)
1257	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1258	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1259	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1260	if len(SSGtext): #if superstructure
1261	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1262	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1263	if not seqRef:
1264	PrintTexture(textureData)
1265	if name in RestraintDict:
1266	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1267	textureData,RestraintDict[name],pFile)
1268
1269	elif PawleyRef:
1270	if Print:
1271	print >>pFile,'\n Phase name: ',General['Name']
1272	print >>pFile,135*'-'
1273	print >>pFile,''
1274	if len(SSGtext): #if superstructure
1275	for line in SSGtext: print >>pFile,line
1276	if len(SSGtable):
1277	for item in SSGtable:
1278	line = ' %s '%(item)
1279	print >>pFile,line
1280	else:
1281	print >>pFile,' ( 1) %s'%(SSGtable[0])
1282	else:
1283	for line in SGtext: print >>pFile,line
1284	if len(SGtable):
1285	for item in SGtable:
1286	line = ' %s '%(item)
1287	print >>pFile,line
1288	else:
1289	print >>pFile,' ( 1) %s'%(SGtable[0])
1290	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1291	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1292	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1293	if len(SSGtext): #if superstructure
1294	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1295	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1296	pawleyVary = []
1297	for i,refl in enumerate(PawleyRef):
1298	phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1299	if im:
1300	pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1301	else:
1302	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1303	if refl[5+im]:
1304	pawleyVary.append(pfx+'PWLref:'+str(i))
1305	GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary) #does G2mv.StoreEquivalence
1306	phaseVary += pawleyVary
1307
1308	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave
1309
1310	def cellFill(pfx,SGData,parmDict,sigDict):
1311	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1312	from the parameter and sig dictionaries.
1313
1314	:param str pfx: parameter prefix ("n::", where n is a phase number)
1315	:param dict SGdata: a symmetry object
1316	:param dict parmDict: a dictionary of parameters
1317	:param dict sigDict: a dictionary of uncertainties on parameters
1318
1319	:returns: A,sigA where each is a list of six terms with the A terms
1320	'''
1321	if SGData['SGLaue'] in ['-1',]:
1322	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1323	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1324	elif SGData['SGLaue'] in ['2/m',]:
1325	if SGData['SGUniq'] == 'a':
1326	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1327	0,0,parmDict[pfx+'A5']]
1328	elif SGData['SGUniq'] == 'b':
1329	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1330	0,parmDict[pfx+'A4'],0]
1331	else:
1332	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1333	parmDict[pfx+'A3'],0,0]
1334	elif SGData['SGLaue'] in ['mmm',]:
1335	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1336	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1337	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1338	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1339	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1340	parmDict[pfx+'A0'],0,0]
1341	elif SGData['SGLaue'] in ['3R', '3mR']:
1342	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1343	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1344	elif SGData['SGLaue'] in ['m3m','m3']:
1345	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1346
1347	try:
1348	if SGData['SGLaue'] in ['-1',]:
1349	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1350	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1351	elif SGData['SGLaue'] in ['2/m',]:
1352	if SGData['SGUniq'] == 'a':
1353	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1354	0,0,sigDict[pfx+'A5']]
1355	elif SGData['SGUniq'] == 'b':
1356	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1357	0,sigDict[pfx+'A4'],0]
1358	else:
1359	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1360	sigDict[pfx+'A3'],0,0]
1361	elif SGData['SGLaue'] in ['mmm',]:
1362	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1363	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1364	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1365	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1366	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1367	elif SGData['SGLaue'] in ['3R', '3mR']:
1368	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1369	elif SGData['SGLaue'] in ['m3m','m3']:
1370	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1371	except KeyError:
1372	sigA = [0,0,0,0,0,0]
1373	return A,sigA
1374
1375	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1376	'needs a doc string'
1377	if phaseRest:
1378	Amat = G2lat.cell2AB(cell)[0]
1379	cx,ct,cs = AtPtrs[:3]
1380	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1381	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1382	['ChemComp','Sites'],['Texture','HKLs']]
1383	for name,rest in names:
1384	itemRest = phaseRest[name]
1385	if itemRest[rest] and itemRest['Use']:
1386	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1387	if name in ['Bond','Angle','Plane','Chiral']:
1388	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1389	for indx,ops,obs,esd in itemRest[rest]:
1390	try:
1391	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1392	AtName = ''
1393	for i,Aname in enumerate(AtNames):
1394	AtName += Aname
1395	if ops[i] == '1':
1396	AtName += '-'
1397	else:
1398	AtName += '+('+ops[i]+')-'
1399	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1400	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1401	if name == 'Bond':
1402	calc = G2mth.getRestDist(XYZ,Amat)
1403	elif name == 'Angle':
1404	calc = G2mth.getRestAngle(XYZ,Amat)
1405	elif name == 'Plane':
1406	calc = G2mth.getRestPlane(XYZ,Amat)
1407	elif name == 'Chiral':
1408	calc = G2mth.getRestChiral(XYZ,Amat)
1409	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1410	except KeyError:
1411	del itemRest[rest]
1412	elif name in ['Torsion','Rama']:
1413	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1414	coeffDict = itemRest['Coeff']
1415	for indx,ops,cofName,esd in itemRest[rest]:
1416	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1417	AtName = ''
1418	for i,Aname in enumerate(AtNames):
1419	AtName += Aname+'+('+ops[i]+')-'
1420	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1421	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1422	if name == 'Torsion':
1423	tor = G2mth.getRestTorsion(XYZ,Amat)
1424	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1425	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1426	else:
1427	phi,psi = G2mth.getRestRama(XYZ,Amat)
1428	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1429	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])
1430	elif name == 'ChemComp':
1431	print >>pFile,' atoms mul*frac factor prod'
1432	for indx,factors,obs,esd in itemRest[rest]:
1433	try:
1434	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1435	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1436	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1437	mulfrac = mul*frac
1438	calcs = mulfracfactors
1439	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1440	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1441	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1442	except KeyError:
1443	del itemRest[rest]
1444	elif name == 'Texture' and textureData['Order']:
1445	Start = False
1446	SHCoef = textureData['SH Coeff'][1]
1447	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1448	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1449	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1450	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1451	PH = np.array(hkl)
1452	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1453	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1454	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1455	Z = ma.masked_greater(Z,0.0)
1456	num = ma.count(Z)
1457	sum = 0
1458	if num:
1459	sum = np.sum(Z)
1460	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)
1461
1462	def getCellEsd(pfx,SGData,A,covData):
1463	'needs a doc string'
1464	dpr = 180./np.pi
1465	rVsq = G2lat.calc_rVsq(A)
1466	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1467	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1468	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1469	scos = cosd(cellst[3:6])
1470	ssin = sind(cellst[3:6])
1471	scot = scos/ssin
1472	rcos = cosd(cell[3:6])
1473	rsin = sind(cell[3:6])
1474	rcot = rcos/rsin
1475	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1476	varyList = covData['varyList']
1477	covMatrix = covData['covMatrix']
1478	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1479	Ax = np.array(A)
1480	Ax[3:] /= 2.
1481	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1482	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1483	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1484	Vol = 1/np.sqrt(rVsq)
1485	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1486	R123 = Ax[0]Ax[1]Ax[2]
1487	dsasdg = np.zeros((3,6))
1488	dadg = np.zeros((6,6))
1489	for i0 in range(3): #0 1 2
1490	i1 = (i0+1)%3 #1 2 0
1491	i2 = (i1+1)%3 #2 0 1
1492	i3 = 5-i2 #3 5 4
1493	i4 = 5-i1 #4 3 5
1494	i5 = 5-i0 #5 4 3
1495	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1496	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1497	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1498	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1499	denom = np.sqrt(denmsq)
1500	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1501	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1502	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1503	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1504	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1505	dadg[i5][i5] = -Ax[i0]/denom
1506	for i0 in range(3):
1507	i1 = (i0+1)%3
1508	i2 = (i1+1)%3
1509	i3 = 5-i2
1510	for ij in range(6):
1511	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1512	if ij == i0:
1513	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1514	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1515	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1516	var = np.diag(sigMat)
1517	CS = np.where(var>0.,np.sqrt(var),0.)
1518	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1519
1520	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1521	'needs a doc string'
1522
1523	def PrintAtomsAndSig(General,Atoms,atomsSig):
1524	print >>pFile,'\n Atoms:'
1525	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1526	if General['Type'] == 'magnetic':
1527	line += ' Mx My Mz'
1528	elif General['Type'] == 'macromolecular':
1529	line = ' res no residue chain '+line
1530	cx,ct,cs,cia = General['AtomPtrs']
1531	print >>pFile,line
1532	print >>pFile,135*'-'
1533	fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1534	cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1535	noFXsig = {cx:[10' ','%10s'],cx+1:[10' ','%10s'],cx+2:[10' ','%10s'],cx+3:[8' ','%8s']}
1536	for atyp in General['AtomTypes']: #zero composition
1537	General['NoAtoms'][atyp] = 0.
1538	for i,at in enumerate(Atoms):
1539	General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5]) #new composition
1540	if General['Type'] == 'macromolecular':
1541	name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1542	valstr = ' values: '
1543	sigstr = ' sig : '
1544	else:
1545	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1546	valstr = ' values:'
1547	sigstr = ' sig :'
1548	for ind in range(cx,cx+4):
1549	sigind = str(i)+':'+str(ind)
1550	valstr += fmt[ind]%(at[ind])
1551	if sigind in atomsSig:
1552	sigstr += fmt[ind]%(atomsSig[sigind])
1553	else:
1554	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1555	if at[cia] == 'I':
1556	valstr += fmt[cia+1]%(at[cia+1])
1557	if '%d:%d'%(i,cia+1) in atomsSig:
1558	sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1559	else:
1560	sigstr += 8*' '
1561	else:
1562	valstr += 8*' '
1563	sigstr += 8*' '
1564	for ind in range(cia+2,cia+8):
1565	sigind = str(i)+':'+str(ind)
1566	valstr += fmt[ind]%(at[ind])
1567	if sigind in atomsSig:
1568	sigstr += fmt[ind]%(atomsSig[sigind])
1569	else:
1570	sigstr += 8*' '
1571	print >>pFile,name
1572	print >>pFile,valstr
1573	print >>pFile,sigstr
1574
1575	def PrintWavesAndSig(General,Atoms,wavesSig):
1576	cx,ct,cs,cia = General['AtomPtrs']
1577	print >>pFile,'\n Modulation waves'
1578	names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1579	'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1580	'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1581	print >>pFile,135*'-'
1582	for i,at in enumerate(Atoms):
1583	AtomSS = at[-1]['SS1']
1584	waveType = AtomSS['waveType']
1585	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1586	Waves = AtomSS[Stype]
1587	if len(Waves):
1588	print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:' \
1589	%(at[ct-1],at[cs],Stype,waveType)
1590	line = ''
1591	for iw,wave in enumerate(Waves):
1592	stiw = ':'+str(i)+':'+str(iw)
1593	namstr = ' names :'
1594	valstr = ' values:'
1595	sigstr = ' esds :'
1596	if Stype == 'Spos':
1597	nt = 6
1598	ot = 0
1599	if waveType in ['ZigZag','Block',] and not iw:
1600	nt = 5
1601	ot = 6
1602	for j in range(nt):
1603	name = names['Spos'][j+ot]
1604	namstr += '%12s'%(name)
1605	valstr += '%12.4f'%(wave[0][j])
1606	if name+stiw in wavesSig:
1607	sigstr += '%12.4f'%(wavesSig[name+stiw])
1608	else:
1609	sigstr += 12*' '
1610	elif Stype == 'Sfrac':
1611	ot = 0
1612	if 'Crenel' in waveType and not iw:
1613	ot = 2
1614	for j in range(2):
1615	name = names['Sfrac'][j+ot]
1616	namstr += '%12s'%(names['Sfrac'][j+ot])
1617	valstr += '%12.4f'%(wave[0][j])
1618	if name+stiw in wavesSig:
1619	sigstr += '%12.4f'%(wavesSig[name+stiw])
1620	else:
1621	sigstr += 12*' '
1622	elif Stype == 'Sadp':
1623	for j in range(12):
1624	name = names['Sadp'][j]
1625	namstr += '%10s'%(names['Sadp'][j])
1626	valstr += '%10.6f'%(wave[0][j])
1627	if name+stiw in wavesSig:
1628	sigstr += '%10.6f'%(wavesSig[name+stiw])
1629	else:
1630	sigstr += 10*' '
1631	elif Stype == 'Smag':
1632	for j in range(6):
1633	name = names['Smag'][j]
1634	namstr += '%12s'%(names['Smag'][j])
1635	valstr += '%12.4f'%(wave[0][j])
1636	if name+stiw in wavesSig:
1637	sigstr += '%12.4f'%(wavesSig[name+stiw])
1638	else:
1639	sigstr += 12*' '
1640
1641	print >>pFile,namstr
1642	print >>pFile,valstr
1643	print >>pFile,sigstr
1644
1645
1646	def PrintRBObjPOAndSig(rbfx,rbsx):
1647	namstr = ' names :'
1648	valstr = ' values:'
1649	sigstr = ' esds :'
1650	for i,px in enumerate(['Px:','Py:','Pz:']):
1651	name = pfx+rbfx+px+rbsx
1652	namstr += '%12s'%('Pos '+px[1])
1653	valstr += '%12.5f'%(parmDict[name])
1654	if name in sigDict:
1655	sigstr += '%12.5f'%(sigDict[name])
1656	else:
1657	sigstr += 12*' '
1658	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1659	name = pfx+rbfx+po+rbsx
1660	namstr += '%12s'%('Ori '+po[1])
1661	valstr += '%12.5f'%(parmDict[name])
1662	if name in sigDict:
1663	sigstr += '%12.5f'%(sigDict[name])
1664	else:
1665	sigstr += 12*' '
1666	print >>pFile,namstr
1667	print >>pFile,valstr
1668	print >>pFile,sigstr
1669
1670	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1671	namstr = ' names :'
1672	valstr = ' values:'
1673	sigstr = ' esds :'
1674	if 'N' not in TLS:
1675	print >>pFile,' Thermal motion:'
1676	if 'T' in TLS:
1677	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1678	name = pfx+rbfx+pt+rbsx
1679	namstr += '%12s'%(pt[:3])
1680	valstr += '%12.5f'%(parmDict[name])
1681	if name in sigDict:
1682	sigstr += '%12.5f'%(sigDict[name])
1683	else:
1684	sigstr += 12*' '
1685	print >>pFile,namstr
1686	print >>pFile,valstr
1687	print >>pFile,sigstr
1688	if 'L' in TLS:
1689	namstr = ' names :'
1690	valstr = ' values:'
1691	sigstr = ' esds :'
1692	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1693	name = pfx+rbfx+pt+rbsx
1694	namstr += '%12s'%(pt[:3])
1695	valstr += '%12.3f'%(parmDict[name])
1696	if name in sigDict:
1697	sigstr += '%12.3f'%(sigDict[name])
1698	else:
1699	sigstr += 12*' '
1700	print >>pFile,namstr
1701	print >>pFile,valstr
1702	print >>pFile,sigstr
1703	if 'S' in TLS:
1704	namstr = ' names :'
1705	valstr = ' values:'
1706	sigstr = ' esds :'
1707	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1708	name = pfx+rbfx+pt+rbsx
1709	namstr += '%12s'%(pt[:3])
1710	valstr += '%12.4f'%(parmDict[name])
1711	if name in sigDict:
1712	sigstr += '%12.4f'%(sigDict[name])
1713	else:
1714	sigstr += 12*' '
1715	print >>pFile,namstr
1716	print >>pFile,valstr
1717	print >>pFile,sigstr
1718	if 'U' in TLS:
1719	name = pfx+rbfx+'U:'+rbsx
1720	namstr = ' names :'+'%12s'%('Uiso')
1721	valstr = ' values:'+'%12.5f'%(parmDict[name])
1722	if name in sigDict:
1723	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1724	else:
1725	sigstr = ' esds :'+12*' '
1726	print >>pFile,namstr
1727	print >>pFile,valstr
1728	print >>pFile,sigstr
1729
1730	def PrintRBObjTorAndSig(rbsx):
1731	namstr = ' names :'
1732	valstr = ' values:'
1733	sigstr = ' esds :'
1734	nTors = len(RBObj['Torsions'])
1735	if nTors:
1736	print >>pFile,' Torsions:'
1737	for it in range(nTors):
1738	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1739	namstr += '%12s'%('Tor'+str(it))
1740	valstr += '%12.4f'%(parmDict[name])
1741	if name in sigDict:
1742	sigstr += '%12.4f'%(sigDict[name])
1743	print >>pFile,namstr
1744	print >>pFile,valstr
1745	print >>pFile,sigstr
1746
1747	def PrintSHtextureAndSig(textureData,SHtextureSig):
1748	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1749	names = ['omega','chi','phi']
1750	namstr = ' names :'
1751	ptstr = ' values:'
1752	sigstr = ' esds :'
1753	for name in names:
1754	namstr += '%12s'%(name)
1755	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1756	if 'Sample '+name in SHtextureSig:
1757	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1758	else:
1759	sigstr += 12*' '
1760	print >>pFile,namstr
1761	print >>pFile,ptstr
1762	print >>pFile,sigstr
1763	print >>pFile,'\n Texture coefficients:'
1764	SHcoeff = textureData['SH Coeff'][1]
1765	SHkeys = SHcoeff.keys()
1766	nCoeff = len(SHcoeff)
1767	nBlock = nCoeff/10+1
1768	iBeg = 0
1769	iFin = min(iBeg+10,nCoeff)
1770	for block in range(nBlock):
1771	namstr = ' names :'
1772	ptstr = ' values:'
1773	sigstr = ' esds :'
1774	for name in SHkeys[iBeg:iFin]:
1775	namstr += '%12s'%(name)
1776	ptstr += '%12.3f'%(SHcoeff[name])
1777	if name in SHtextureSig:
1778	sigstr += '%12.3f'%(SHtextureSig[name])
1779	else:
1780	sigstr += 12*' '
1781	print >>pFile,namstr
1782	print >>pFile,ptstr
1783	print >>pFile,sigstr
1784	iBeg += 10
1785	iFin = min(iBeg+10,nCoeff)
1786
1787	print >>pFile,'\n Phases:'
1788	for phase in Phases:
1789	print >>pFile,' Result for phase: ',phase
1790	Phase = Phases[phase]
1791	General = Phase['General']
1792	SGData = General['SGData']
1793	Atoms = Phase['Atoms']
1794	if Atoms and not General.get('doPawley'):
1795	cx,ct,cs,cia = General['AtomPtrs']
1796	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1797	cell = General['Cell']
1798	pId = Phase['pId']
1799	pfx = str(pId)+'::'
1800	if cell[0]:
1801	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1802	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1803	print >>pFile,' Reciprocal metric tensor: '
1804	ptfmt = "%15.9f"
1805	names = ['A11','A22','A33','A12','A13','A23']
1806	namstr = ' names :'
1807	ptstr = ' values:'
1808	sigstr = ' esds :'
1809	for name,a,siga in zip(names,A,sigA):
1810	namstr += '%15s'%(name)
1811	ptstr += ptfmt%(a)
1812	if siga:
1813	sigstr += ptfmt%(siga)
1814	else:
1815	sigstr += 15*' '
1816	print >>pFile,namstr
1817	print >>pFile,ptstr
1818	print >>pFile,sigstr
1819	cell[1:7] = G2lat.A2cell(A)
1820	cell[7] = G2lat.calc_V(A)
1821	print >>pFile,' New unit cell:'
1822	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1823	names = ['a','b','c','alpha','beta','gamma','Volume']
1824	namstr = ' names :'
1825	ptstr = ' values:'
1826	sigstr = ' esds :'
1827	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1828	namstr += '%12s'%(name)
1829	ptstr += fmt%(a)
1830	if siga:
1831	sigstr += fmt%(siga)
1832	else:
1833	sigstr += 12*' '
1834	print >>pFile,namstr
1835	print >>pFile,ptstr
1836	print >>pFile,sigstr
1837	ik = 6 #for Pawley stuff below
1838	if General['Type'] in ['modulated','magnetic']:
1839	ik = 7
1840	Vec,vRef,maxH = General['SuperVec']
1841	if vRef:
1842	print >>pFile,' New modulation vector:'
1843	namstr = ' names :'
1844	ptstr = ' values:'
1845	sigstr = ' esds :'
1846	for var in ['mV0','mV1','mV2']:
1847	namstr += '%12s'%(pfx+var)
1848	ptstr += '%12.6f'%(parmDict[pfx+var])
1849	if pfx+var in sigDict:
1850	sigstr += '%12.6f'%(sigDict[pfx+var])
1851	else:
1852	sigstr += 12*' '
1853	print >>pFile,namstr
1854	print >>pFile,ptstr
1855	print >>pFile,sigstr
1856
1857	General['Mass'] = 0.
1858	if Phase['General'].get('doPawley'):
1859	pawleyRef = Phase['Pawley ref']
1860	for i,refl in enumerate(pawleyRef):
1861	key = pfx+'PWLref:'+str(i)
1862	refl[ik] = parmDict[key]
1863	if key in sigDict:
1864	refl[ik+1] = sigDict[key]
1865	else:
1866	refl[ik+1] = 0
1867	else:
1868	VRBIds = RBIds['Vector']
1869	RRBIds = RBIds['Residue']
1870	RBModels = Phase['RBModels']
1871	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1872	jrb = VRBIds.index(RBObj['RBId'])
1873	rbsx = str(irb)+':'+str(jrb)
1874	print >>pFile,' Vector rigid body parameters:'
1875	PrintRBObjPOAndSig('RBV',rbsx)
1876	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1877	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1878	jrb = RRBIds.index(RBObj['RBId'])
1879	rbsx = str(irb)+':'+str(jrb)
1880	print >>pFile,' Residue rigid body parameters:'
1881	PrintRBObjPOAndSig('RBR',rbsx)
1882	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1883	PrintRBObjTorAndSig(rbsx)
1884	atomsSig = {}
1885	wavesSig = {}
1886	cx,ct,cs,cia = General['AtomPtrs']
1887	for i,at in enumerate(Atoms):
1888	names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1889	cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1890	cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i)}
1891	for ind in range(cx,cx+4):
1892	at[ind] = parmDict[names[ind]]
1893	if ind in range(cx,cx+3):
1894	name = names[ind].replace('A','dA')
1895	else:
1896	name = names[ind]
1897	if name in sigDict:
1898	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1899	if at[cia] == 'I':
1900	at[cia+1] = parmDict[names[cia+1]]
1901	if names[cia+1] in sigDict:
1902	atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
1903	else:
1904	for ind in range(cia+2,cia+8):
1905	at[ind] = parmDict[names[ind]]
1906	if names[ind] in sigDict:
1907	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1908	ind = General['AtomTypes'].index(at[ct])
1909	General['Mass'] += General['AtomMass'][ind]at[cx+3]at[cx+5]
1910	if General['Type'] in ['modulated','magnetic']:
1911	AtomSS = at[-1]['SS1']
1912	waveType = AtomSS['waveType']
1913	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1914	Waves = AtomSS[Stype]
1915	for iw,wave in enumerate(Waves):
1916	stiw = str(i)+':'+str(iw)
1917	if Stype == 'Spos':
1918	if waveType in ['ZigZag','Block',] and not iw:
1919	names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
1920	else:
1921	names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
1922	'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
1923	elif Stype == 'Sadp':
1924	names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
1925	'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
1926	'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
1927	'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
1928	elif Stype == 'Sfrac':
1929	if 'Crenel' in waveType and not iw:
1930	names = ['Fzero:'+stiw,'Fwid:'+stiw]
1931	else:
1932	names = ['Fsin:'+stiw,'Fcos:'+stiw]
1933	elif Stype == 'Smag':
1934	names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
1935	'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
1936	for iname,name in enumerate(names):
1937	AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
1938	if pfx+name in sigDict:
1939	wavesSig[name] = sigDict[pfx+name]
1940
1941	PrintAtomsAndSig(General,Atoms,atomsSig)
1942	if General['Type'] in ['modulated','magnetic']:
1943	PrintWavesAndSig(General,Atoms,wavesSig)
1944
1945
1946	textureData = General['SH Texture']
1947	if textureData['Order']:
1948	SHtextureSig = {}
1949	for name in ['omega','chi','phi']:
1950	aname = pfx+'SH '+name
1951	textureData['Sample '+name][1] = parmDict[aname]
1952	if aname in sigDict:
1953	SHtextureSig['Sample '+name] = sigDict[aname]
1954	for name in textureData['SH Coeff'][1]:
1955	aname = pfx+name
1956	textureData['SH Coeff'][1][name] = parmDict[aname]
1957	if aname in sigDict:
1958	SHtextureSig[name] = sigDict[aname]
1959	PrintSHtextureAndSig(textureData,SHtextureSig)
1960	if phase in RestraintDict and not Phase['General'].get('doPawley'):
1961	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1962	textureData,RestraintDict[phase],pFile)
1963
1964	################################################################################
1965	##### Histogram & Phase data
1966	################################################################################
1967
1968	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1969	'''Loads the HAP histogram/phase information into dicts
1970
1971	:param dict Phases: phase information
1972	:param dict Histograms: Histogram information
1973	:param bool Print: prints information as it is read
1974	:param file pFile: file object to print to (the default, None causes printing to the console)
1975	:param bool resetRefList: Should the contents of the Reflection List be initialized
1976	on loading. The default, True, initializes the Reflection List as it is loaded.
1977
1978	:returns: (hapVary,hapDict,controlDict)
1979	* hapVary: list of refined variables
1980	* hapDict: dict with refined variables and their values
1981	* controlDict: dict with fixed parameters
1982	'''
1983
1984	def PrintSize(hapData):
1985	if hapData[0] in ['isotropic','uniaxial']:
1986	line = '\n Size model : %9s'%(hapData[0])
1987	line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1988	if hapData[0] == 'uniaxial':
1989	line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1990	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1991	print >>pFile,line
1992	else:
1993	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1994	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1995	Snames = ['S11','S22','S33','S12','S13','S23']
1996	ptlbls = ' names :'
1997	ptstr = ' values:'
1998	varstr = ' refine:'
1999	for i,name in enumerate(Snames):
2000	ptlbls += '%12s' % (name)
2001	ptstr += '%12.6f' % (hapData[4][i])
2002	varstr += '%12s' % (str(hapData[5][i]))
2003	print >>pFile,ptlbls
2004	print >>pFile,ptstr
2005	print >>pFile,varstr
2006
2007	def PrintMuStrain(hapData,SGData):
2008	if hapData[0] in ['isotropic','uniaxial']:
2009	line = '\n Mustrain model: %9s'%(hapData[0])
2010	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2011	if hapData[0] == 'uniaxial':
2012	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2013	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2014	print >>pFile,line
2015	else:
2016	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
2017	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2018	Snames = G2spc.MustrainNames(SGData)
2019	ptlbls = ' names :'
2020	ptstr = ' values:'
2021	varstr = ' refine:'
2022	for i,name in enumerate(Snames):
2023	ptlbls += '%12s' % (name)
2024	ptstr += '%12.6f' % (hapData[4][i])
2025	varstr += '%12s' % (str(hapData[5][i]))
2026	print >>pFile,ptlbls
2027	print >>pFile,ptstr
2028	print >>pFile,varstr
2029
2030	def PrintHStrain(hapData,SGData):
2031	print >>pFile,'\n Hydrostatic/elastic strain: '
2032	Hsnames = G2spc.HStrainNames(SGData)
2033	ptlbls = ' names :'
2034	ptstr = ' values:'
2035	varstr = ' refine:'
2036	for i,name in enumerate(Hsnames):
2037	ptlbls += '%12s' % (name)
2038	ptstr += '%12.4g' % (hapData[0][i])
2039	varstr += '%12s' % (str(hapData[1][i]))
2040	print >>pFile,ptlbls
2041	print >>pFile,ptstr
2042	print >>pFile,varstr
2043
2044	def PrintSHPO(hapData):
2045	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
2046	str(hapData[4])+' Refine? '+str(hapData[2])
2047	ptlbls = ' names :'
2048	ptstr = ' values:'
2049	for item in hapData[5]:
2050	ptlbls += '%12s'%(item)
2051	ptstr += '%12.3f'%(hapData[5][item])
2052	print >>pFile,ptlbls
2053	print >>pFile,ptstr
2054
2055	def PrintBabinet(hapData):
2056	print >>pFile,'\n Babinet form factor modification: '
2057	ptlbls = ' names :'
2058	ptstr = ' values:'
2059	varstr = ' refine:'
2060	for name in ['BabA','BabU']:
2061	ptlbls += '%12s' % (name)
2062	ptstr += '%12.6f' % (hapData[name][0])
2063	varstr += '%12s' % (str(hapData[name][1]))
2064	print >>pFile,ptlbls
2065	print >>pFile,ptstr
2066	print >>pFile,varstr
2067
2068	hapDict = {}
2069	hapVary = []
2070	controlDict = {}
2071	poType = {}
2072	poAxes = {}
2073	spAxes = {}
2074	spType = {}
2075
2076	for phase in Phases:
2077	HistoPhase = Phases[phase]['Histograms']
2078	SGData = Phases[phase]['General']['SGData']
2079	cell = Phases[phase]['General']['Cell'][1:7]
2080	A = G2lat.cell2A(cell)
2081	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
2082	SSGData = Phases[phase]['General']['SSGData']
2083	Vec,x,maxH = Phases[phase]['General']['SuperVec']
2084	pId = Phases[phase]['pId']
2085	histoList = HistoPhase.keys()
2086	histoList.sort()
2087	for histogram in histoList:
2088	try:
2089	Histogram = Histograms[histogram]
2090	except KeyError:
2091	#skip if histogram not included e.g. in a sequential refinement
2092	continue
2093	hapData = HistoPhase[histogram]
2094	hId = Histogram['hId']
2095	if 'PWDR' in histogram:
2096	limits = Histogram['Limits'][1]
2097	inst = Histogram['Instrument Parameters'][0]
2098	Zero = inst['Zero'][0]
2099	if 'C' in inst['Type'][1]:
2100	try:
2101	wave = inst['Lam'][1]
2102	except KeyError:
2103	wave = inst['Lam1'][1]
2104	dmin = wave/(2.0*sind(limits[1]/2.0))
2105	elif 'T' in inst['Type'][0]:
2106	dmin = limits[0]/inst['difC'][1]
2107	pfx = str(pId)+':'+str(hId)+':'
2108	for item in ['Scale','Extinction']:
2109	hapDict[pfx+item] = hapData[item][0]
2110	if hapData[item][1]:
2111	hapVary.append(pfx+item)
2112	names = G2spc.HStrainNames(SGData)
2113	HSvals = []
2114	for i,name in enumerate(names):
2115	hapDict[pfx+name] = hapData['HStrain'][0][i]
2116	HSvals.append(hapDict[pfx+name])
2117	if hapData['HStrain'][1][i]:
2118	hapVary.append(pfx+name)
2119	DIJS = G2spc.HStrainVals(HSvals,SGData)
2120	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2121	if hapData['Pref.Ori.'][0] == 'MD':
2122	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2123	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2124	if hapData['Pref.Ori.'][2]:
2125	hapVary.append(pfx+'MD')
2126	else: #'SH' spherical harmonics
2127	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2128	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2129	controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2130	controlDict[pfx+'SHhkl'] = []
2131	try: #patch for old Pref.Ori. items
2132	controlDict[pfx+'SHtoler'] = 0.1
2133	if hapData['Pref.Ori.'][6][0] != '':
2134	controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2135	controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2136	except IndexError:
2137	pass
2138	for item in hapData['Pref.Ori.'][5]:
2139	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2140	if hapData['Pref.Ori.'][2]:
2141	hapVary.append(pfx+item)
2142	for item in ['Mustrain','Size']:
2143	controlDict[pfx+item+'Type'] = hapData[item][0]
2144	hapDict[pfx+item+';mx'] = hapData[item][1][2]
2145	if hapData[item][2][2]:
2146	hapVary.append(pfx+item+';mx')
2147	if hapData[item][0] in ['isotropic','uniaxial']:
2148	hapDict[pfx+item+';i'] = hapData[item][1][0]
2149	if hapData[item][2][0]:
2150	hapVary.append(pfx+item+';i')
2151	if hapData[item][0] == 'uniaxial':
2152	controlDict[pfx+item+'Axis'] = hapData[item][3]
2153	hapDict[pfx+item+';a'] = hapData[item][1][1]
2154	if hapData[item][2][1]:
2155	hapVary.append(pfx+item+';a')
2156	else: #generalized for mustrain or ellipsoidal for size
2157	Nterms = len(hapData[item][4])
2158	if item == 'Mustrain':
2159	names = G2spc.MustrainNames(SGData)
2160	pwrs = []
2161	for name in names:
2162	h,k,l = name[1:]
2163	pwrs.append([int(h),int(k),int(l)])
2164	controlDict[pfx+'MuPwrs'] = pwrs
2165	for i in range(Nterms):
2166	sfx = ';'+str(i)
2167	hapDict[pfx+item+sfx] = hapData[item][4][i]
2168	if hapData[item][5][i]:
2169	hapVary.append(pfx+item+sfx)
2170	for bab in ['BabA','BabU']:
2171	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2172	if hapData['Babinet'][bab][1]:
2173	hapVary.append(pfx+bab)
2174
2175	if Print:
2176	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2177	print >>pFile,135*'-'
2178	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2179	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2180	if hapData['Pref.Ori.'][0] == 'MD':
2181	Ax = hapData['Pref.Ori.'][3]
2182	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2183	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2184	else: #'SH' for spherical harmonics
2185	PrintSHPO(hapData['Pref.Ori.'])
2186	print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2187	PrintSize(hapData['Size'])
2188	PrintMuStrain(hapData['Mustrain'],SGData)
2189	PrintHStrain(hapData['HStrain'],SGData)
2190	if hapData['Babinet']['BabA'][0]:
2191	PrintBabinet(hapData['Babinet'])
2192	if resetRefList:
2193	refList = []
2194	Uniq = []
2195	Phi = []
2196	if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
2197	ifSuper = True
2198	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2199	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2200	for h,k,l,m,d in HKLd:
2201	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2202	mul *= 2 # for powder overlap of Friedel pairs
2203	if m or not ext:
2204	if 'C' in inst['Type'][0]:
2205	pos = G2lat.Dsp2pos(inst,d)
2206	if limits[0] < pos < limits[1]:
2207	refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
2208	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2209	Uniq.append(uniq)
2210	Phi.append(phi)
2211	elif 'T' in inst['Type'][0]:
2212	pos = G2lat.Dsp2pos(inst,d)
2213	if limits[0] < pos < limits[1]:
2214	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2215	refList.append([h,k,l,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])
2216	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2217	Uniq.append(uniq)
2218	Phi.append(phi)
2219	else:
2220	ifSuper = False
2221	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2222	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2223	for h,k,l,d in HKLd:
2224	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2225	mul *= 2 # for powder overlap of Friedel pairs
2226	if ext:
2227	continue
2228	if 'C' in inst['Type'][0]:
2229	pos = G2lat.Dsp2pos(inst,d)
2230	if limits[0] < pos < limits[1]:
2231	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])
2232	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2233	Uniq.append(uniq)
2234	Phi.append(phi)
2235	elif 'T' in inst['Type'][0]:
2236	pos = G2lat.Dsp2pos(inst,d)
2237	if limits[0] < pos < limits[1]:
2238	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2239	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])
2240	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2241	Uniq.append(uniq)
2242	Phi.append(phi)
2243	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2244	elif 'HKLF' in histogram:
2245	inst = Histogram['Instrument Parameters'][0]
2246	hId = Histogram['hId']
2247	hfx = ':%d:'%(hId)
2248	for item in inst:
2249	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2250	hapDict[hfx+item] = inst[item][1]
2251	pfx = str(pId)+':'+str(hId)+':'
2252	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2253	if hapData['Scale'][1]:
2254	hapVary.append(pfx+'Scale')
2255
2256	extApprox,extType,extParms = hapData['Extinction']
2257	controlDict[pfx+'EType'] = extType
2258	controlDict[pfx+'EApprox'] = extApprox
2259	if 'C' in inst['Type'][0]:
2260	controlDict[pfx+'Tbar'] = extParms['Tbar']
2261	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2262	if 'Primary' in extType:
2263	Ekey = ['Ep',]
2264	elif 'I & II' in extType:
2265	Ekey = ['Eg','Es']
2266	elif 'Secondary Type II' == extType:
2267	Ekey = ['Es',]
2268	elif 'Secondary Type I' == extType:
2269	Ekey = ['Eg',]
2270	else: #'None'
2271	Ekey = []
2272	for eKey in Ekey:
2273	hapDict[pfx+eKey] = extParms[eKey][0]
2274	if extParms[eKey][1]:
2275	hapVary.append(pfx+eKey)
2276	for bab in ['BabA','BabU']:
2277	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2278	if hapData['Babinet'][bab][1]:
2279	hapVary.append(pfx+bab)
2280	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2281	if len(Twins) == 1:
2282	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2283	if hapData.get('Flack',[0,False])[1]:
2284	hapVary.append(pfx+'Flack')
2285	sumTwFr = 0.
2286	controlDict[pfx+'TwinLaw'] = []
2287	controlDict[pfx+'TwinInv'] = []
2288	NTL = 0
2289	for it,twin in enumerate(Twins):
2290	if 'bool' in str(type(twin[0])):
2291	controlDict[pfx+'TwinInv'].append(twin[0])
2292	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2293	else:
2294	NTL += 1
2295	controlDict[pfx+'TwinInv'].append(False)
2296	controlDict[pfx+'TwinLaw'].append(twin[0])
2297	if it:
2298	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2299	sumTwFr += twin[1]
2300	else:
2301	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2302	controlDict[pfx+'TwinNMN'] = twin[1][1]
2303	if Twins[0][1][1]:
2304	hapVary.append(pfx+'TwinFr:'+str(it))
2305	controlDict[pfx+'NTL'] = NTL
2306	#need to make constraint on TwinFr
2307	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2308	if len(Twins) > 1: #force sum to unity
2309	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2310	if Print:
2311	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2312	print >>pFile,135*'-'
2313	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2314	if extType != 'None':
2315	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox)
2316	text = ' Parameters :'
2317	for eKey in Ekey:
2318	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2319	print >>pFile,text
2320	if hapData['Babinet']['BabA'][0]:
2321	PrintBabinet(hapData['Babinet'])
2322	if not SGData['SGInv'] and len(Twins) == 1:
2323	print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2324	if len(Twins) > 1:
2325	for it,twin in enumerate(Twins):
2326	if 'bool' in str(type(twin[0])):
2327	print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1]
2328	else:
2329	print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1]
2330
2331	Histogram['Reflection Lists'] = phase
2332
2333	return hapVary,hapDict,controlDict
2334
2335	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2336	'needs a doc string'
2337
2338	def PrintSizeAndSig(hapData,sizeSig):
2339	line = '\n Size model: %9s'%(hapData[0])
2340	refine = False
2341	if hapData[0] in ['isotropic','uniaxial']:
2342	line += ' equatorial:%12.5f'%(hapData[1][0])
2343	if sizeSig[0][0]:
2344	line += ', sig:%8.4f'%(sizeSig[0][0])
2345	refine = True
2346	if hapData[0] == 'uniaxial':
2347	line += ' axial:%12.4f'%(hapData[1][1])
2348	if sizeSig[0][1]:
2349	refine = True
2350	line += ', sig:%8.4f'%(sizeSig[0][1])
2351	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2352	if sizeSig[0][2]:
2353	refine = True
2354	line += ', sig:%8.4f'%(sizeSig[0][2])
2355	if refine:
2356	print >>pFile,line
2357	else:
2358	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2359	if sizeSig[0][2]:
2360	refine = True
2361	line += ', sig:%8.4f'%(sizeSig[0][2])
2362	Snames = ['S11','S22','S33','S12','S13','S23']
2363	ptlbls = ' name :'
2364	ptstr = ' value :'
2365	sigstr = ' sig :'
2366	for i,name in enumerate(Snames):
2367	ptlbls += '%12s' % (name)
2368	ptstr += '%12.6f' % (hapData[4][i])
2369	if sizeSig[1][i]:
2370	refine = True
2371	sigstr += '%12.6f' % (sizeSig[1][i])
2372	else:
2373	sigstr += 12*' '
2374	if refine:
2375	print >>pFile,line
2376	print >>pFile,ptlbls
2377	print >>pFile,ptstr
2378	print >>pFile,sigstr
2379
2380	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2381	line = '\n Mustrain model: %9s'%(hapData[0])
2382	refine = False
2383	if hapData[0] in ['isotropic','uniaxial']:
2384	line += ' equatorial:%12.1f'%(hapData[1][0])
2385	if mustrainSig[0][0]:
2386	line += ', sig:%8.1f'%(mustrainSig[0][0])
2387	refine = True
2388	if hapData[0] == 'uniaxial':
2389	line += ' axial:%12.1f'%(hapData[1][1])
2390	if mustrainSig[0][1]:
2391	line += ', sig:%8.1f'%(mustrainSig[0][1])
2392	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2393	if mustrainSig[0][2]:
2394	refine = True
2395	line += ', sig:%8.3f'%(mustrainSig[0][2])
2396	if refine:
2397	print >>pFile,line
2398	else:
2399	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2400	if mustrainSig[0][2]:
2401	refine = True
2402	line += ', sig:%8.3f'%(mustrainSig[0][2])
2403	Snames = G2spc.MustrainNames(SGData)
2404	ptlbls = ' name :'
2405	ptstr = ' value :'
2406	sigstr = ' sig :'
2407	for i,name in enumerate(Snames):
2408	ptlbls += '%12s' % (name)
2409	ptstr += '%12.1f' % (hapData[4][i])
2410	if mustrainSig[1][i]:
2411	refine = True
2412	sigstr += '%12.1f' % (mustrainSig[1][i])
2413	else:
2414	sigstr += 12*' '
2415	if refine:
2416	print >>pFile,line
2417	print >>pFile,ptlbls
2418	print >>pFile,ptstr
2419	print >>pFile,sigstr
2420
2421	def PrintHStrainAndSig(hapData,strainSig,SGData):
2422	Hsnames = G2spc.HStrainNames(SGData)
2423	ptlbls = ' name :'
2424	ptstr = ' value :'
2425	sigstr = ' sig :'
2426	refine = False
2427	for i,name in enumerate(Hsnames):
2428	ptlbls += '%12s' % (name)
2429	ptstr += '%12.4g' % (hapData[0][i])
2430	if name in strainSig:
2431	refine = True
2432	sigstr += '%12.4g' % (strainSig[name])
2433	else:
2434	sigstr += 12*' '
2435	if refine:
2436	print >>pFile,'\n Hydrostatic/elastic strain: '
2437	print >>pFile,ptlbls
2438	print >>pFile,ptstr
2439	print >>pFile,sigstr
2440
2441	def PrintSHPOAndSig(pfx,hapData,POsig):
2442	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2443	ptlbls = ' names :'
2444	ptstr = ' values:'
2445	sigstr = ' sig :'
2446	for item in hapData[5]:
2447	ptlbls += '%12s'%(item)
2448	ptstr += '%12.3f'%(hapData[5][item])
2449	if pfx+item in POsig:
2450	sigstr += '%12.3f'%(POsig[pfx+item])
2451	else:
2452	sigstr += 12*' '
2453	print >>pFile,ptlbls
2454	print >>pFile,ptstr
2455	print >>pFile,sigstr
2456
2457	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2458	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2459	text = ''
2460	for item in hapData[2]:
2461	if pfx+item in ScalExtSig:
2462	text += ' %s: '%(item)
2463	text += '%12.2e'%(hapData[2][item][0])
2464	if pfx+item in ScalExtSig:
2465	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2466	print >>pFile,text
2467
2468	def PrintBabinetAndSig(pfx,hapData,BabSig):
2469	print >>pFile,'\n Babinet form factor modification: '
2470	ptlbls = ' names :'
2471	ptstr = ' values:'
2472	sigstr = ' sig :'
2473	for item in hapData:
2474	ptlbls += '%12s'%(item)
2475	ptstr += '%12.3f'%(hapData[item][0])
2476	if pfx+item in BabSig:
2477	sigstr += '%12.3f'%(BabSig[pfx+item])
2478	else:
2479	sigstr += 12*' '
2480	print >>pFile,ptlbls
2481	print >>pFile,ptstr
2482	print >>pFile,sigstr
2483
2484	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2485	print >>pFile,'\n Twin Law fractions : '
2486	ptlbls = ' names :'
2487	ptstr = ' values:'
2488	sigstr = ' sig :'
2489	for it,item in enumerate(twinData):
2490	ptlbls += ' twin #%d'%(it)
2491	if it:
2492	ptstr += '%12.3f'%(item[1])
2493	else:
2494	ptstr += '%12.3f'%(item[1][0])
2495	if pfx+'TwinFr:'+str(it) in TwinSig:
2496	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2497	else:
2498	sigstr += 12*' '
2499	print >>pFile,ptlbls
2500	print >>pFile,ptstr
2501	print >>pFile,sigstr
2502
2503
2504	PhFrExtPOSig = {}
2505	SizeMuStrSig = {}
2506	ScalExtSig = {}
2507	BabSig = {}
2508	TwinFrSig = {}
2509	wtFrSum = {}
2510	for phase in Phases:
2511	HistoPhase = Phases[phase]['Histograms']
2512	General = Phases[phase]['General']
2513	SGData = General['SGData']
2514	pId = Phases[phase]['pId']
2515	histoList = HistoPhase.keys()
2516	histoList.sort()
2517	for histogram in histoList:
2518	try:
2519	Histogram = Histograms[histogram]
2520	except KeyError:
2521	#skip if histogram not included e.g. in a sequential refinement
2522	continue
2523	hapData = HistoPhase[histogram]
2524	hId = Histogram['hId']
2525	pfx = str(pId)+':'+str(hId)+':'
2526	if hId not in wtFrSum:
2527	wtFrSum[hId] = 0.
2528	if 'PWDR' in histogram:
2529	for item in ['Scale','Extinction']:
2530	hapData[item][0] = parmDict[pfx+item]
2531	if pfx+item in sigDict:
2532	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2533	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2534	if hapData['Pref.Ori.'][0] == 'MD':
2535	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2536	if pfx+'MD' in sigDict:
2537	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2538	else: #'SH' spherical harmonics
2539	for item in hapData['Pref.Ori.'][5]:
2540	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2541	if pfx+item in sigDict:
2542	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2543	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2544	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2545	pfx+'HStrain':{}})
2546	for item in ['Mustrain','Size']:
2547	hapData[item][1][2] = parmDict[pfx+item+';mx']
2548	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2549	if pfx+item+';mx' in sigDict:
2550	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2551	if hapData[item][0] in ['isotropic','uniaxial']:
2552	hapData[item][1][0] = parmDict[pfx+item+';i']
2553	if item == 'Size':
2554	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2555	if pfx+item+';i' in sigDict:
2556	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2557	if hapData[item][0] == 'uniaxial':
2558	hapData[item][1][1] = parmDict[pfx+item+';a']
2559	if item == 'Size':
2560	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2561	if pfx+item+';a' in sigDict:
2562	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2563	else: #generalized for mustrain or ellipsoidal for size
2564	Nterms = len(hapData[item][4])
2565	for i in range(Nterms):
2566	sfx = ':'+str(i)
2567	hapData[item][4][i] = parmDict[pfx+item+sfx]
2568	if pfx+item+sfx in sigDict:
2569	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2570	names = G2spc.HStrainNames(SGData)
2571	for i,name in enumerate(names):
2572	hapData['HStrain'][0][i] = parmDict[pfx+name]
2573	if pfx+name in sigDict:
2574	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2575	for name in ['BabA','BabU']:
2576	hapData['Babinet'][name][0] = parmDict[pfx+name]
2577	if pfx+name in sigDict:
2578	BabSig[pfx+name] = sigDict[pfx+name]
2579
2580	elif 'HKLF' in histogram:
2581	for item in ['Scale','Flack']:
2582	if parmDict.get(pfx+item):
2583	hapData[item][0] = parmDict[pfx+item]
2584	if pfx+item in sigDict:
2585	ScalExtSig[pfx+item] = sigDict[pfx+item]
2586	for item in ['Ep','Eg','Es']:
2587	if parmDict.get(pfx+item):
2588	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2589	if pfx+item in sigDict:
2590	ScalExtSig[pfx+item] = sigDict[pfx+item]
2591	for item in ['BabA','BabU']:
2592	hapData['Babinet'][item][0] = parmDict[pfx+item]
2593	if pfx+item in sigDict:
2594	BabSig[pfx+item] = sigDict[pfx+item]
2595	if 'Twins' in hapData:
2596	it = 1
2597	sumTwFr = 0.
2598	while True:
2599	try:
2600	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2601	if pfx+'TwinFr:'+str(it) in sigDict:
2602	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2603	if it:
2604	sumTwFr += hapData['Twins'][it][1]
2605	it += 1
2606	except KeyError:
2607	break
2608	hapData['Twins'][0][1][0] = 1.-sumTwFr
2609
2610	if Print:
2611	for phase in Phases:
2612	HistoPhase = Phases[phase]['Histograms']
2613	General = Phases[phase]['General']
2614	SGData = General['SGData']
2615	pId = Phases[phase]['pId']
2616	histoList = HistoPhase.keys()
2617	histoList.sort()
2618	for histogram in histoList:
2619	try:
2620	Histogram = Histograms[histogram]
2621	except KeyError:
2622	#skip if histogram not included e.g. in a sequential refinement
2623	continue
2624	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2625	print >>pFile,130*'-'
2626	hapData = HistoPhase[histogram]
2627	hId = Histogram['hId']
2628	Histogram['Residuals'][str(pId)+'::Name'] = phase
2629	pfx = str(pId)+':'+str(hId)+':'
2630	hfx = ':%s:'%(hId)
2631	if 'PWDR' in histogram:
2632	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2633	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2634	print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2635
2636	if pfx+'Scale' in PhFrExtPOSig:
2637	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2638	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2639	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2640	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2641	if pfx+'Extinction' in PhFrExtPOSig:
2642	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2643	if hapData['Pref.Ori.'][0] == 'MD':
2644	if pfx+'MD' in PhFrExtPOSig:
2645	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2646	else:
2647	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2648	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2649	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2650	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2651	if len(BabSig):
2652	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2653
2654	elif 'HKLF' in histogram:
2655	Inst = Histogram['Instrument Parameters'][0]
2656	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)' \
2657	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2658	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2659	if FFtables != None and 'N' not in Inst['Type'][0]:
2660	PrintFprime(FFtables,hfx,pFile)
2661	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2662	if pfx+'Scale' in ScalExtSig:
2663	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2664	if hapData['Extinction'][0] != 'None':
2665	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2666	if len(BabSig):
2667	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2668	if pfx+'Flack' in ScalExtSig:
2669	print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2670	if len(TwinFrSig):
2671	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2672
2673	################################################################################
2674	##### Histogram data
2675	################################################################################
2676
2677	def GetHistogramData(Histograms,Print=True,pFile=None):
2678	'needs a doc string'
2679
2680	def GetBackgroundParms(hId,Background):
2681	Back = Background[0]
2682	DebyePeaks = Background[1]
2683	bakType,bakFlag = Back[:2]
2684	backVals = Back[3:]
2685	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2686	backDict = dict(zip(backNames,backVals))
2687	backVary = []
2688	if bakFlag:
2689	backVary = backNames
2690	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2691	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2692	debyeDict = {}
2693	debyeList = []
2694	for i in range(DebyePeaks['nDebye']):
2695	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2696	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2697	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2698	debyeVary = []
2699	for item in debyeList:
2700	if item[1]:
2701	debyeVary.append(item[0])
2702	backDict.update(debyeDict)
2703	backVary += debyeVary
2704	peakDict = {}
2705	peakList = []
2706	for i in range(DebyePeaks['nPeaks']):
2707	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2708	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2709	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2710	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2711	peakVary = []
2712	for item in peakList:
2713	if item[1]:
2714	peakVary.append(item[0])
2715	backDict.update(peakDict)
2716	backVary += peakVary
2717	return bakType,backDict,backVary
2718
2719	def GetInstParms(hId,Inst):
2720	dataType = Inst['Type'][0]
2721	instDict = {}
2722	insVary = []
2723	pfx = ':'+str(hId)+':'
2724	insKeys = Inst.keys()
2725	insKeys.sort()
2726	for item in insKeys:
2727	insName = pfx+item
2728	instDict[insName] = Inst[item][1]
2729	if len(Inst[item]) > 2 and Inst[item][2]:
2730	insVary.append(insName)
2731	if 'C' in dataType:
2732	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2733	return dataType,instDict,insVary
2734
2735	def GetSampleParms(hId,Sample):
2736	sampVary = []
2737	hfx = ':'+str(hId)+':'
2738	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2739	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2740	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2741	if key in Sample:
2742	sampDict[hfx+key] = Sample[key]
2743	Type = Sample['Type']
2744	if 'Bragg' in Type: #Bragg-Brentano
2745	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2746	sampDict[hfx+item] = Sample[item][0]
2747	if Sample[item][1]:
2748	sampVary.append(hfx+item)
2749	elif 'Debye' in Type: #Debye-Scherrer
2750	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2751	sampDict[hfx+item] = Sample[item][0]
2752	if Sample[item][1]:
2753	sampVary.append(hfx+item)
2754	return Type,sampDict,sampVary
2755
2756	def PrintBackground(Background):
2757	Back = Background[0]
2758	DebyePeaks = Background[1]
2759	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2760	line = ' Coefficients: '
2761	for i,back in enumerate(Back[3:]):
2762	line += '%10.3f'%(back)
2763	if i and not i%10:
2764	line += '\n'+15*' '
2765	print >>pFile,line
2766	if DebyePeaks['nDebye']:
2767	print >>pFile,'\n Debye diffuse scattering coefficients'
2768	parms = ['DebyeA','DebyeR','DebyeU']
2769	line = ' names : '
2770	for parm in parms:
2771	line += '%8s refine?'%(parm)
2772	print >>pFile,line
2773	for j,term in enumerate(DebyePeaks['debyeTerms']):
2774	line = ' term'+'%2d'%(j)+':'
2775	for i in range(3):
2776	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2777	print >>pFile,line
2778	if DebyePeaks['nPeaks']:
2779	print >>pFile,'\n Single peak coefficients'
2780	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2781	line = ' names : '
2782	for parm in parms:
2783	line += '%8s refine?'%(parm)
2784	print >>pFile,line
2785	for j,term in enumerate(DebyePeaks['peaksList']):
2786	line = ' peak'+'%2d'%(j)+':'
2787	for i in range(4):
2788	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2789	print >>pFile,line
2790
2791	def PrintInstParms(Inst):
2792	print >>pFile,'\n Instrument Parameters:'
2793	insKeys = Inst.keys()
2794	insKeys.sort()
2795	iBeg = 0
2796	Ok = True
2797	while Ok:
2798	ptlbls = ' name :'
2799	ptstr = ' value :'
2800	varstr = ' refine:'
2801	iFin = min(iBeg+9,len(insKeys))
2802	for item in insKeys[iBeg:iFin]:
2803	if item not in ['Type','Source']:
2804	ptlbls += '%12s' % (item)
2805	ptstr += '%12.6f' % (Inst[item][1])
2806	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2807	varstr += 12*' '
2808	else:
2809	varstr += '%12s' % (str(bool(Inst[item][2])))
2810	print >>pFile,ptlbls
2811	print >>pFile,ptstr
2812	print >>pFile,varstr
2813	iBeg = iFin
2814	if iBeg == len(insKeys):
2815	Ok = False
2816	else:
2817	print >>pFile,'\n'
2818
2819	def PrintSampleParms(Sample):
2820	print >>pFile,'\n Sample Parameters:'
2821	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2822	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2823	ptlbls = ' name :'
2824	ptstr = ' value :'
2825	varstr = ' refine:'
2826	if 'Bragg' in Sample['Type']:
2827	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2828	ptlbls += '%14s'%(item)
2829	ptstr += '%14.4f'%(Sample[item][0])
2830	varstr += '%14s'%(str(bool(Sample[item][1])))
2831
2832	elif 'Debye' in Type: #Debye-Scherrer
2833	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2834	ptlbls += '%14s'%(item)
2835	ptstr += '%14.4f'%(Sample[item][0])
2836	varstr += '%14s'%(str(bool(Sample[item][1])))
2837
2838	print >>pFile,ptlbls
2839	print >>pFile,ptstr
2840	print >>pFile,varstr
2841
2842	histDict = {}
2843	histVary = []
2844	controlDict = {}
2845	histoList = Histograms.keys()
2846	histoList.sort()
2847	for histogram in histoList:
2848	if 'PWDR' in histogram:
2849	Histogram = Histograms[histogram]
2850	hId = Histogram['hId']
2851	pfx = ':'+str(hId)+':'
2852	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2853	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2854	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2855	for excl in controlDict[pfx+'Exclude']:
2856	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2857	if controlDict[pfx+'Exclude']:
2858	ma.mask_rows(Histogram['Data'])
2859	Background = Histogram['Background']
2860	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2861	controlDict[pfx+'bakType'] = Type
2862	histDict.update(bakDict)
2863	histVary += bakVary
2864
2865	Inst = Histogram['Instrument Parameters'][0]
2866	Type,instDict,insVary = GetInstParms(hId,Inst)
2867	controlDict[pfx+'histType'] = Type
2868	if 'XC' in Type:
2869	if pfx+'Lam1' in instDict:
2870	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2871	else:
2872	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2873	histDict.update(instDict)
2874	histVary += insVary
2875
2876	Sample = Histogram['Sample Parameters']
2877	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2878	controlDict[pfx+'instType'] = Type
2879	histDict.update(sampDict)
2880	histVary += sampVary
2881
2882
2883	if Print:
2884	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2885	print >>pFile,135*'-'
2886	Units = {'C':' deg','T':' msec'}
2887	units = Units[controlDict[pfx+'histType'][2]]
2888	Limits = controlDict[pfx+'Limits']
2889	print >>pFile,' Instrument type: ',Sample['Type']
2890	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2891	if len(controlDict[pfx+'Exclude']):
2892	excls = controlDict[pfx+'Exclude']
2893	for excl in excls:
2894	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2895	PrintSampleParms(Sample)
2896	PrintInstParms(Inst)
2897	PrintBackground(Background)
2898	elif 'HKLF' in histogram:
2899	Histogram = Histograms[histogram]
2900	hId = Histogram['hId']
2901	pfx = ':'+str(hId)+':'
2902	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2903	Inst = Histogram['Instrument Parameters'][0]
2904	controlDict[pfx+'histType'] = Inst['Type'][0]
2905	if 'X' in Inst['Type'][0]:
2906	histDict[pfx+'Lam'] = Inst['Lam'][1]
2907	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2908	return histVary,histDict,controlDict
2909
2910	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
2911	'needs a doc string'
2912
2913	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2914	Back = Background[0]
2915	DebyePeaks = Background[1]
2916	lenBack = len(Back[3:])
2917	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2918	for i in range(lenBack):
2919	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2920	if pfx+'Back;'+str(i) in sigDict:
2921	backSig[i] = sigDict[pfx+'Back;'+str(i)]
2922	if DebyePeaks['nDebye']:
2923	for i in range(DebyePeaks['nDebye']):
2924	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
2925	for j,name in enumerate(names):
2926	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2927	if name in sigDict:
2928	backSig[lenBack+3*i+j] = sigDict[name]
2929	if DebyePeaks['nPeaks']:
2930	for i in range(DebyePeaks['nPeaks']):
2931	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2932	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2933	for j,name in enumerate(names):
2934	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2935	if name in sigDict:
2936	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2937	return backSig
2938
2939	def SetInstParms(pfx,Inst,parmDict,sigDict):
2940	instSig = {}
2941	insKeys = Inst.keys()
2942	insKeys.sort()
2943	for item in insKeys:
2944	insName = pfx+item
2945	Inst[item][1] = parmDict[insName]
2946	if insName in sigDict:
2947	instSig[item] = sigDict[insName]
2948	else:
2949	instSig[item] = 0
2950	return instSig
2951
2952	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2953	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2954	sampSig = [0 for i in range(5)]
2955	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2956	Sample[item][0] = parmDict[pfx+item]
2957	if pfx+item in sigDict:
2958	sampSig[i] = sigDict[pfx+item]
2959	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2960	sampSig = [0 for i in range(4)]
2961	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2962	Sample[item][0] = parmDict[pfx+item]
2963	if pfx+item in sigDict:
2964	sampSig[i] = sigDict[pfx+item]
2965	return sampSig
2966
2967	def PrintBackgroundSig(Background,backSig):
2968	Back = Background[0]
2969	DebyePeaks = Background[1]
2970	lenBack = len(Back[3:])
2971	valstr = ' value : '
2972	sigstr = ' sig : '
2973	refine = False
2974	for i,back in enumerate(Back[3:]):
2975	valstr += '%10.4g'%(back)
2976	if Back[1]:
2977	refine = True
2978	sigstr += '%10.4g'%(backSig[i])
2979	else:
2980	sigstr += 10*' '
2981	if refine:
2982	print >>pFile,'\n Background function: ',Back[0]
2983	print >>pFile,valstr
2984	print >>pFile,sigstr
2985	if DebyePeaks['nDebye']:
2986	ifAny = False
2987	ptfmt = "%12.3f"
2988	names = ' names :'
2989	ptstr = ' values:'
2990	sigstr = ' esds :'
2991	for item in sigDict:
2992	if 'Debye' in item:
2993	ifAny = True
2994	names += '%12s'%(item)
2995	ptstr += ptfmt%(parmDict[item])
2996	sigstr += ptfmt%(sigDict[item])
2997	if ifAny:
2998	print >>pFile,'\n Debye diffuse scattering coefficients'
2999	print >>pFile,names
3000	print >>pFile,ptstr
3001	print >>pFile,sigstr
3002	if DebyePeaks['nPeaks']:
3003	print >>pFile,'\n Single peak coefficients:'
3004	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3005	line = ' peak no. '
3006	for parm in parms:
3007	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3008	print >>pFile,line
3009	for ip in range(DebyePeaks['nPeaks']):
3010	ptstr = ' %4d '%(ip)
3011	for parm in parms:
3012	fmt = '%14.3f'
3013	efmt = '%12.3f'
3014	if parm == 'BkPkpos':
3015	fmt = '%14.4f'
3016	efmt = '%12.4f'
3017	name = pfx+parm+';%d'%(ip)
3018	ptstr += fmt%(parmDict[name])
3019	if name in sigDict:
3020	ptstr += efmt%(sigDict[name])
3021	else:
3022	ptstr += 12*' '
3023	print >>pFile,ptstr
3024	sumBk = np.array(Histogram['sumBk'])
3025	print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g' \
3026	%(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3027
3028	def PrintInstParmsSig(Inst,instSig):
3029	refine = False
3030	insKeys = instSig.keys()
3031	insKeys.sort()
3032	iBeg = 0
3033	Ok = True
3034	while Ok:
3035	ptlbls = ' names :'
3036	ptstr = ' value :'
3037	sigstr = ' sig :'
3038	iFin = min(iBeg+9,len(insKeys))
3039	for name in insKeys[iBeg:iFin]:
3040	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3041	ptlbls += '%12s' % (name)
3042	ptstr += '%12.6f' % (Inst[name][1])
3043	if instSig[name]:
3044	refine = True
3045	sigstr += '%12.6f' % (instSig[name])
3046	else:
3047	sigstr += 12*' '
3048	if refine:
3049	print >>pFile,'\n Instrument Parameters:'
3050	print >>pFile,ptlbls
3051	print >>pFile,ptstr
3052	print >>pFile,sigstr
3053	iBeg = iFin
3054	if iBeg == len(insKeys):
3055	Ok = False
3056
3057	def PrintSampleParmsSig(Sample,sampleSig):
3058	ptlbls = ' names :'
3059	ptstr = ' values:'
3060	sigstr = ' sig :'
3061	refine = False
3062	if 'Bragg' in Sample['Type']:
3063	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3064	ptlbls += '%14s'%(item)
3065	ptstr += '%14.4f'%(Sample[item][0])
3066	if sampleSig[i]:
3067	refine = True
3068	sigstr += '%14.4f'%(sampleSig[i])
3069	else:
3070	sigstr += 14*' '
3071
3072	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3073	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3074	ptlbls += '%14s'%(item)
3075	ptstr += '%14.4f'%(Sample[item][0])
3076	if sampleSig[i]:
3077	refine = True
3078	sigstr += '%14.4f'%(sampleSig[i])
3079	else:
3080	sigstr += 14*' '
3081
3082	if refine:
3083	print >>pFile,'\n Sample Parameters:'
3084	print >>pFile,ptlbls
3085	print >>pFile,ptstr
3086	print >>pFile,sigstr
3087
3088	histoList = Histograms.keys()
3089	histoList.sort()
3090	for histogram in histoList:
3091	if 'PWDR' in histogram:
3092	Histogram = Histograms[histogram]
3093	hId = Histogram['hId']
3094	pfx = ':'+str(hId)+':'
3095	Background = Histogram['Background']
3096	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3097
3098	Inst = Histogram['Instrument Parameters'][0]
3099	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3100
3101	Sample = Histogram['Sample Parameters']
3102	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3103
3104	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3105	print >>pFile,135*'-'
3106	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3107	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3108	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3109	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3110	if Print:
3111	print >>pFile,' Instrument type: ',Sample['Type']
3112	if FFtables != None and 'N' not in Inst['Type'][0]:
3113	PrintFprime(FFtables,pfx,pFile)
3114	PrintSampleParmsSig(Sample,sampSig)
3115	PrintInstParmsSig(Inst,instSig)
3116	PrintBackgroundSig(Background,backSig)
3117

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

Download in other formats: