Context Navigation

source: trunk/GSASIIstrIO.py @ 2311

Last change on this file since 2311 was 2245, checked in by vondreele, 9 years ago
add a line to the layer GUI informing lack of refinement capability G2strIO - skip faulted phases in refinements.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 139.8 KB

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

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

Download in other formats: