Context Navigation

source: trunk/GSASIIstruct.py @ 578

Last change on this file since 578 was 578, checked in by vondreele, 12 years ago
fix "hard" Hessian singular matrix problems - now deletes parameters fix Pawley overlap reflection constraint problem
Property svn:keywords set to `Date Author Revision URL Id`
File size: 147.8 KB

Line
1	#GSASIIstructure - structure computation routines
2	########### SVN repository information ###################
3	# $Date: 2012-04-27 17:14:21 +0000 (Fri, 27 Apr 2012) $
4	# $Author: vondreele $
5	# $Revision: 578 $
6	# $URL: trunk/GSASIIstruct.py $
7	# $Id: GSASIIstruct.py 578 2012-04-27 17:14:21Z vondreele $
8	########### SVN repository information ###################
9	import sys
10	import os
11	import os.path as ospath
12	import time
13	import math
14	import cPickle
15	import numpy as np
16	import numpy.linalg as nl
17	import scipy.optimize as so
18	import GSASIIpath
19	import GSASIIElem as G2el
20	import GSASIIlattice as G2lat
21	import GSASIIspc as G2spc
22	import GSASIIpwd as G2pwd
23	import GSASIImapvars as G2mv
24	import GSASIImath as G2mth
25
26	sind = lambda x: np.sin(x*np.pi/180.)
27	cosd = lambda x: np.cos(x*np.pi/180.)
28	tand = lambda x: np.tan(x*np.pi/180.)
29	asind = lambda x: 180.*np.arcsin(x)/np.pi
30	acosd = lambda x: 180.*np.arccos(x)/np.pi
31	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
32
33	ateln2 = 8.0*math.log(2.0)
34
35
36	def GetControls(GPXfile):
37	''' Returns dictionary of control items found in GSASII gpx file
38	input:
39	GPXfile = .gpx full file name
40	return:
41	Controls = dictionary of control items
42	'''
43	Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
44	'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
45	fl = open(GPXfile,'rb')
46	while True:
47	try:
48	data = cPickle.load(fl)
49	except EOFError:
50	break
51	datum = data[0]
52	if datum[0] == 'Controls':
53	Controls.update(datum[1])
54	fl.close()
55	return Controls
56
57	def GetConstraints(GPXfile):
58	'''Read the constraints from the GPX file and interpret them
59	'''
60	constList = []
61	fl = open(GPXfile,'rb')
62	while True:
63	try:
64	data = cPickle.load(fl)
65	except EOFError:
66	break
67	datum = data[0]
68	if datum[0] == 'Constraints':
69	constDict = datum[1]
70	for item in constDict:
71	constList += constDict[item]
72	fl.close()
73	constDict,fixedList,ignored = ProcessConstraints(constList)
74	if ignored:
75	print ignored,'old-style Constraints were rejected'
76	return constDict,fixedList
77
78	def ProcessConstraints(constList):
79	"interpret constraints"
80	constDict = []
81	fixedList = []
82	ignored = 0
83	for item in constList:
84	if item[-1] == 'h':
85	# process a hold
86	fixedList.append('0')
87	constDict.append({item[0][1]:0.0})
88	elif item[-1] == 'f':
89	# process a new variable
90	fixedList.append(None)
91	constDict.append({})
92	for term in item[:-3]:
93	constDict[-1][term[1]] = term[0]
94	#constFlag[-1] = ['Vary']
95	elif item[-1] == 'c':
96	# process a contraint relationship
97	fixedList.append(str(item[-3]))
98	constDict.append({})
99	for term in item[:-3]:
100	constDict[-1][term[1]] = term[0]
101	#constFlag[-1] = ['VaryFree']
102	elif item[-1] == 'e':
103	# process an equivalence
104	firstmult = None
105	eqlist = []
106	for term in item[:-3]:
107	if term[0] == 0: term[0] = 1.0
108	if firstmult is None:
109	firstmult,firstvar = term
110	else:
111	eqlist.append([term[1],firstmult/term[0]])
112	G2mv.StoreEquivalence(firstvar,eqlist)
113	else:
114	ignored += 1
115	return constDict,fixedList,ignored
116
117	def CheckConstraints(GPXfile):
118	'''Load constraints and related info and return any error or warning messages'''
119	# init constraints
120	G2mv.InitVars()
121	# get variables
122	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
123	if not Phases:
124	return 'Error: No Phases!',''
125	if not Histograms:
126	return 'Error: no diffraction data',''
127	Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
128	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
129	histVary,histDict,controlDict = GetHistogramData(Histograms)
130	varyList = phaseVary+hapVary+histVary
131	constrDict,fixedList = GetConstraints(GPXfile)
132	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
133
134	def GetPhaseNames(GPXfile):
135	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
136	input:
137	GPXfile = gpx full file name
138	return:
139	PhaseNames = list of phase names
140	'''
141	fl = open(GPXfile,'rb')
142	PhaseNames = []
143	while True:
144	try:
145	data = cPickle.load(fl)
146	except EOFError:
147	break
148	datum = data[0]
149	if 'Phases' == datum[0]:
150	for datus in data[1:]:
151	PhaseNames.append(datus[0])
152	fl.close()
153	return PhaseNames
154
155	def GetAllPhaseData(GPXfile,PhaseName):
156	''' Returns the entire dictionary for PhaseName from GSASII gpx file
157	input:
158	GPXfile = gpx full file name
159	PhaseName = phase name
160	return:
161	phase dictionary
162	'''
163	fl = open(GPXfile,'rb')
164	General = {}
165	Atoms = []
166	while True:
167	try:
168	data = cPickle.load(fl)
169	except EOFError:
170	break
171	datum = data[0]
172	if 'Phases' == datum[0]:
173	for datus in data[1:]:
174	if datus[0] == PhaseName:
175	break
176	fl.close()
177	return datus[1]
178
179	def GetHistograms(GPXfile,hNames):
180	""" Returns a dictionary of histograms found in GSASII gpx file
181	input:
182	GPXfile = .gpx full file name
183	hNames = list of histogram names
184	return:
185	Histograms = dictionary of histograms (types = PWDR & HKLF)
186	"""
187	fl = open(GPXfile,'rb')
188	Histograms = {}
189	while True:
190	try:
191	data = cPickle.load(fl)
192	except EOFError:
193	break
194	datum = data[0]
195	hist = datum[0]
196	if hist in hNames:
197	if 'PWDR' in hist[:4]:
198	PWDRdata = {}
199	PWDRdata['Data'] = datum[1][1] #powder data arrays
200	PWDRdata[data[2][0]] = data[2][1] #Limits
201	PWDRdata[data[3][0]] = data[3][1] #Background
202	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
203	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
204	try:
205	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
206	except IndexError:
207	PWDRdata['Reflection lists'] = {}
208
209	Histograms[hist] = PWDRdata
210	elif 'HKLF' in hist[:4]:
211	HKLFdata = []
212	datum = data[0]
213	HKLFdata = datum[1:][0]
214	Histograms[hist] = HKLFdata
215	fl.close()
216	return Histograms
217
218	def GetHistogramNames(GPXfile,hType):
219	""" Returns a list of histogram names found in GSASII gpx file
220	input:
221	GPXfile = .gpx full file name
222	hType = list ['PWDR','HKLF']
223	return:
224	HistogramNames = list of histogram names (types = PWDR & HKLF)
225	"""
226	fl = open(GPXfile,'rb')
227	HistogramNames = []
228	while True:
229	try:
230	data = cPickle.load(fl)
231	except EOFError:
232	break
233	datum = data[0]
234	if datum[0][:4] in hType:
235	HistogramNames.append(datum[0])
236	fl.close()
237	return HistogramNames
238
239	def GetUsedHistogramsAndPhases(GPXfile):
240	''' Returns all histograms that are found in any phase
241	and any phase that uses a histogram
242	input:
243	GPXfile = .gpx full file name
244	return:
245	Histograms = dictionary of histograms as {name:data,...}
246	Phases = dictionary of phases that use histograms
247	'''
248	phaseNames = GetPhaseNames(GPXfile)
249	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
250	allHistograms = GetHistograms(GPXfile,histoList)
251	phaseData = {}
252	for name in phaseNames:
253	phaseData[name] = GetAllPhaseData(GPXfile,name)
254	Histograms = {}
255	Phases = {}
256	for phase in phaseData:
257	Phase = phaseData[phase]
258	if Phase['Histograms']:
259	if phase not in Phases:
260	pId = phaseNames.index(phase)
261	Phase['pId'] = pId
262	Phases[phase] = Phase
263	for hist in Phase['Histograms']:
264	if hist not in Histograms:
265	Histograms[hist] = allHistograms[hist]
266	#future restraint, etc. histograms here
267	hId = histoList.index(hist)
268	Histograms[hist]['hId'] = hId
269	return Histograms,Phases
270
271	def getBackupName2(GPXfile,makeBack=True): #not work correctly
272	GPXpath,GPXname = ospath.split(GPXfile)
273	if GPXpath == '': GPXpath = '.'
274	Name = ospath.splitext(GPXname)[0]
275	files = os.listdir(GPXpath)
276	last = 0
277	for name in files:
278	name = name.split('.')
279	if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
280	if makeBack:
281	last = max(last,int(name[-2].strip('bak'))+1)
282	else:
283	last = max(last,int(name[-2].strip('bak')))
284	GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
285	return GPXback
286
287	def getBackupName(GPXfile,makeBack): #recovered old one
288	GPXpath,GPXname = ospath.split(GPXfile)
289	if GPXpath == '': GPXpath = '.'
290	Name = ospath.splitext(GPXname)[0]
291	files = os.listdir(GPXpath)
292	last = 0
293	for name in files:
294	name = name.split('.')
295	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
296	if makeBack:
297	last = max(last,int(name[1].strip('bak'))+1)
298	else:
299	last = max(last,int(name[1].strip('bak')))
300	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
301	return GPXback
302
303	def GPXBackup(GPXfile,makeBack=True):
304	import distutils.file_util as dfu
305	GPXback = getBackupName(GPXfile,makeBack)
306	dfu.copy_file(GPXfile,GPXback)
307	return GPXback
308
309	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
310	''' Updates gpxfile from all histograms that are found in any phase
311	and any phase that used a histogram
312	input:
313	GPXfile = .gpx full file name
314	Histograms = dictionary of histograms as {name:data,...}
315	Phases = dictionary of phases that use histograms
316	CovData = dictionary of refined variables, varyList, & covariance matrix
317	makeBack = True if new backup of .gpx file is to be made; else use the last one made
318	'''
319
320	GPXback = GPXBackup(GPXfile,makeBack)
321	print '\n',135*'-'
322	print 'Read from file:',GPXback
323	print 'Save to file :',GPXfile
324	infile = open(GPXback,'rb')
325	outfile = open(GPXfile,'wb')
326	while True:
327	try:
328	data = cPickle.load(infile)
329	except EOFError:
330	break
331	datum = data[0]
332	# print 'read: ',datum[0]
333	if datum[0] == 'Phases':
334	for iphase in range(len(data)):
335	if data[iphase][0] in Phases:
336	phaseName = data[iphase][0]
337	data[iphase][1].update(Phases[phaseName])
338	elif datum[0] == 'Covariance':
339	data[0][1] = CovData
340	try:
341	histogram = Histograms[datum[0]]
342	# print 'found ',datum[0]
343	data[0][1][1] = histogram['Data']
344	for datus in data[1:]:
345	# print ' read: ',datus[0]
346	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
347	datus[1] = histogram[datus[0]]
348	except KeyError:
349	pass
350
351	cPickle.dump(data,outfile,1)
352	infile.close()
353	outfile.close()
354	print 'GPX file save successful'
355
356	def SetSeqResult(GPXfile,Histograms,SeqResult):
357	GPXback = GPXBackup(GPXfile)
358	print '\n',135*'-'
359	print 'Read from file:',GPXback
360	print 'Save to file :',GPXfile
361	infile = open(GPXback,'rb')
362	outfile = open(GPXfile,'wb')
363	while True:
364	try:
365	data = cPickle.load(infile)
366	except EOFError:
367	break
368	datum = data[0]
369	if datum[0] == 'Sequental results':
370	data[0][1] = SeqResult
371	try:
372	histogram = Histograms[datum[0]]
373	data[0][1][1] = histogram['Data']
374	for datus in data[1:]:
375	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
376	datus[1] = histogram[datus[0]]
377	except KeyError:
378	pass
379
380	cPickle.dump(data,outfile,1)
381	infile.close()
382	outfile.close()
383	print 'GPX file save successful'
384
385	def GetPWDRdata(GPXfile,PWDRname):
386	''' Returns powder data from GSASII gpx file
387	input:
388	GPXfile = .gpx full file name
389	PWDRname = powder histogram name as obtained from GetHistogramNames
390	return:
391	PWDRdata = powder data dictionary with:
392	Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
393
394	'''
395	fl = open(GPXfile,'rb')
396	PWDRdata = {}
397	while True:
398	try:
399	data = cPickle.load(fl)
400	except EOFError:
401	break
402	datum = data[0]
403	if datum[0] == PWDRname:
404	PWDRdata['Data'] = datum[1][1] #powder data arrays
405	PWDRdata[data[2][0]] = data[2][1] #Limits
406	PWDRdata[data[3][0]] = data[3][1] #Background
407	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
408	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
409	try:
410	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
411	except IndexError:
412	PWDRdata['Reflection lists'] = {}
413	fl.close()
414	return PWDRdata
415
416	def GetHKLFdata(GPXfile,HKLFname):
417	''' Returns single crystal data from GSASII gpx file
418	input:
419	GPXfile = .gpx full file name
420	HKLFname = single crystal histogram name as obtained from GetHistogramNames
421	return:
422	HKLFdata = single crystal data list of reflections: for each reflection:
423	HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
424	'''
425	fl = open(GPXfile,'rb')
426	HKLFdata = []
427	while True:
428	try:
429	data = cPickle.load(fl)
430	except EOFError:
431	break
432	datum = data[0]
433	if datum[0] == HKLFname:
434	HKLFdata = datum[1:][0]
435	fl.close()
436	return HKLFdata
437
438	def ShowBanner():
439	print 80''
440	print ' General Structure Analysis System-II Crystal Structure Refinement'
441	print ' by Robert B. Von Dreele & Brian H. Toby'
442	print ' Argonne National Laboratory(C), 2010'
443	print ' This product includes software developed by the UChicago Argonne, LLC,'
444	print ' as Operator of Argonne National Laboratory.'
445	print 80'','\n'
446
447	def ShowControls(Controls):
448	print ' Least squares controls:'
449	print ' Refinement type: ',Controls['deriv type']
450	if 'Hessian' in Controls['deriv type']:
451	print ' Maximum number of cycles:',Controls['max cyc']
452	else:
453	print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
454	print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
455
456	def GetFFtable(General):
457	''' returns a dictionary of form factor data for atom types found in General
458	input:
459	General = dictionary of phase info.; includes AtomTypes
460	return:
461	FFtable = dictionary of form factor data; key is atom type
462	'''
463	atomTypes = General['AtomTypes']
464	FFtable = {}
465	for El in atomTypes:
466	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
467	for item in FFs:
468	if item['Symbol'] == El.upper():
469	FFtable[El] = item
470	return FFtable
471
472	def GetBLtable(General):
473	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
474	input:
475	General = dictionary of phase info.; includes AtomTypes & Isotopes
476	return:
477	BLtable = dictionary of scattering length data; key is atom type
478	'''
479	atomTypes = General['AtomTypes']
480	BLtable = {}
481	isotopes = General['Isotopes']
482	isotope = General['Isotope']
483	for El in atomTypes:
484	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
485	return BLtable
486
487	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
488	if SGLaue in ['-1','2/m','mmm']:
489	return #no Pawley symmetry required constraints
490	eqvDict = {}
491	for i,varyI in enumerate(pawleyVary):
492	eqvDict[varyI] = []
493	refI = int(varyI.split(':')[-1])
494	ih,ik,il = PawleyRef[refI][:3]
495	for varyJ in pawleyVary[i+1:]:
496	refJ = int(varyJ.split(':')[-1])
497	jh,jk,jl = PawleyRef[refJ][:3]
498	if SGLaue in ['4/m','4/mmm']:
499	isum = ih2+ik2
500	jsum = jh2+jk2
501	if abs(il) == abs(jl) and isum == jsum:
502	eqvDict[varyI].append(varyJ)
503	elif SGLaue in ['3R','3mR']:
504	isum = ih2+ik2+il**2
505	jsum = jh2+jk2jl*2
506	isum2 = ihik+ihil+ik*il
507	jsum2 = jhjk+jhjl+jk*jl
508	if isum == jsum and isum2 == jsum2:
509	eqvDict[varyI].append(varyJ)
510	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
511	isum = ih2+ik2+ih*ik
512	jsum = jh2+jk2+jh*jk
513	if abs(il) == abs(jl) and isum == jsum:
514	eqvDict[varyI].append(varyJ)
515	elif SGLaue in ['m3','m3m']:
516	isum = ih2+ik2+il**2
517	jsum = jh2+jk2+jl**2
518	if isum == jsum:
519	eqvDict[varyI].append(varyJ)
520	for item in pawleyVary:
521	if eqvDict[item]:
522	for item2 in pawleyVary:
523	if item2 in eqvDict[item]:
524	eqvDict[item2] = []
525	G2mv.StoreEquivalence(item,eqvDict[item])
526
527	def cellVary(pfx,SGData):
528	if SGData['SGLaue'] in ['-1',]:
529	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
530	elif SGData['SGLaue'] in ['2/m',]:
531	if SGData['SGUniq'] == 'a':
532	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
533	elif SGData['SGUniq'] == 'b':
534	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
535	else:
536	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
537	elif SGData['SGLaue'] in ['mmm',]:
538	return [pfx+'A0',pfx+'A1',pfx+'A2']
539	elif SGData['SGLaue'] in ['4/m','4/mmm']:
540	return [pfx+'A0',pfx+'A2']
541	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
542	return [pfx+'A0',pfx+'A2']
543	elif SGData['SGLaue'] in ['3R', '3mR']:
544	return [pfx+'A0',pfx+'A3']
545	elif SGData['SGLaue'] in ['m3m','m3']:
546	return [pfx+'A0',]
547
548	################################################################################
549	##### Phase data
550	################################################################################
551
552	def GetPhaseData(PhaseData,Print=True):
553
554	def PrintFFtable(FFtable):
555	print '\n X-ray scattering factors:'
556	print ' Symbol fa fb fc'
557	print 99*'-'
558	for Ename in FFtable:
559	ffdata = FFtable[Ename]
560	fa = ffdata['fa']
561	fb = ffdata['fb']
562	print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
563	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
564
565	def PrintBLtable(BLtable):
566	print '\n Neutron scattering factors:'
567	print ' Symbol isotope mass b resonant terms'
568	print 99*'-'
569	for Ename in BLtable:
570	bldata = BLtable[Ename]
571	isotope = bldata[0]
572	mass = bldata[1][0]
573	blen = bldata[1][1]
574	bres = []
575	if len(bldata[1]) > 2:
576	bres = bldata[1][2:]
577	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
578	for item in bres:
579	line += '%10.5g'%(item)
580	print line
581
582	def PrintAtoms(General,Atoms):
583	print '\n Atoms:'
584	line = ' name type refine? x y z '+ \
585	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
586	if General['Type'] == 'magnetic':
587	line += ' Mx My Mz'
588	elif General['Type'] == 'macromolecular':
589	line = ' res no residue chain '+line
590	print line
591	if General['Type'] == 'nuclear':
592	print 135*'-'
593	for i,at in enumerate(Atoms):
594	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
595	'%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
596	if at[9] == 'I':
597	line += '%8.4f'%(at[10])+48*' '
598	else:
599	line += 8*' '
600	for j in range(6):
601	line += '%8.4f'%(at[11+j])
602	print line
603
604	def PrintTexture(textureData):
605	topstr = '\n Spherical harmonics texture: Order:' + \
606	str(textureData['Order'])
607	if textureData['Order']:
608	print topstr+' Refine? '+str(textureData['SH Coeff'][0])
609	else:
610	print topstr
611	return
612	names = ['omega','chi','phi']
613	line = '\n'
614	for name in names:
615	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
616	print line
617	print '\n Texture coefficients:'
618	ptlbls = ' names :'
619	ptstr = ' values:'
620	SHcoeff = textureData['SH Coeff'][1]
621	for item in SHcoeff:
622	ptlbls += '%12s'%(item)
623	ptstr += '%12.4f'%(SHcoeff[item])
624	print ptlbls
625	print ptstr
626
627	if Print: print ' Phases:'
628	phaseVary = []
629	phaseDict = {}
630	phaseConstr = {}
631	pawleyLookup = {}
632	FFtables = {} #scattering factors - xrays
633	BLtables = {} # neutrons
634	Natoms = {}
635	AtMults = {}
636	AtIA = {}
637	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
638	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
639	for name in PhaseData:
640	General = PhaseData[name]['General']
641	pId = PhaseData[name]['pId']
642	pfx = str(pId)+'::'
643	FFtable = GetFFtable(General)
644	BLtable = GetBLtable(General)
645	FFtables.update(FFtable)
646	BLtables.update(BLtable)
647	Atoms = PhaseData[name]['Atoms']
648	try:
649	PawleyRef = PhaseData[name]['Pawley ref']
650	except KeyError:
651	PawleyRef = []
652	SGData = General['SGData']
653	SGtext = G2spc.SGPrint(SGData)
654	cell = General['Cell']
655	A = G2lat.cell2A(cell[1:7])
656	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
657	if cell[0]:
658	phaseVary += cellVary(pfx,SGData)
659	Natoms[pfx] = 0
660	if Atoms and not General.get('doPawley'):
661	if General['Type'] == 'nuclear':
662	Natoms[pfx] = len(Atoms)
663	for i,at in enumerate(Atoms):
664	phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
665	pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
666	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
667	pfx+'AI/A:'+str(i):at[9],})
668	if at[9] == 'I':
669	phaseDict[pfx+'AUiso:'+str(i)] = at[10]
670	else:
671	phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
672	pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
673	if 'F' in at[2]:
674	phaseVary.append(pfx+'Afrac:'+str(i))
675	if 'X' in at[2]:
676	xId,xCoef = G2spc.GetCSxinel(at[7])
677	delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
678	for j in range(3):
679	if xId[j] > 0:
680	phaseVary.append(delnames[j])
681	for k in range(j):
682	if xId[j] == xId[k]:
683	G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),))
684	if 'U' in at[2]:
685	if at[9] == 'I':
686	phaseVary.append(pfx+'AUiso:'+str(i))
687	else:
688	uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
689	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
690	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
691	for j in range(6):
692	if uId[j] > 0:
693	phaseVary.append(names[j])
694	for k in range(j):
695	if uId[j] == uId[k]:
696	G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
697	# elif General['Type'] == 'magnetic':
698	# elif General['Type'] == 'macromolecular':
699
700	textureData = General['SH Texture']
701	if textureData['Order']:
702	phaseDict[pfx+'SHorder'] = textureData['Order']
703	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
704	for name in ['omega','chi','phi']:
705	phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
706	if textureData['Sample '+name][0]:
707	phaseVary.append(pfx+'SH '+name)
708	for name in textureData['SH Coeff'][1]:
709	phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
710	if textureData['SH Coeff'][0]:
711	phaseVary.append(pfx+name)
712
713	if Print:
714	print '\n Phase name: ',General['Name']
715	print 135*'-'
716	PrintFFtable(FFtable)
717	PrintBLtable(BLtable)
718	print ''
719	for line in SGtext: print line
720	PrintAtoms(General,Atoms)
721	print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
722	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
723	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
724	PrintTexture(textureData)
725
726	elif PawleyRef:
727	pawleyVary = []
728	for i,refl in enumerate(PawleyRef):
729	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
730	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
731	if refl[5]:
732	pawleyVary.append(pfx+'PWLref:'+str(i))
733	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
734	phaseVary += pawleyVary
735
736	return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
737
738	def cellFill(pfx,SGData,parmDict,sigDict):
739	if SGData['SGLaue'] in ['-1',]:
740	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
741	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
742	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
743	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
744	elif SGData['SGLaue'] in ['2/m',]:
745	if SGData['SGUniq'] == 'a':
746	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
747	parmDict[pfx+'A3'],0,0]
748	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
749	sigDict[pfx+'A3'],0,0]
750	elif SGData['SGUniq'] == 'b':
751	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
752	0,parmDict[pfx+'A4'],0]
753	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
754	0,sigDict[pfx+'A4'],0]
755	else:
756	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
757	0,0,parmDict[pfx+'A5']]
758	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
759	0,0,sigDict[pfx+'A5']]
760	elif SGData['SGLaue'] in ['mmm',]:
761	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
762	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
763	elif SGData['SGLaue'] in ['4/m','4/mmm']:
764	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
765	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
766	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
767	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
768	parmDict[pfx+'A0'],0,0]
769	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
770	elif SGData['SGLaue'] in ['3R', '3mR']:
771	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
772	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
773	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
774	elif SGData['SGLaue'] in ['m3m','m3']:
775	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
776	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
777	return A,sigA
778
779	def getCellEsd(pfx,SGData,A,covData):
780	dpr = 180./np.pi
781	rVsq = G2lat.calc_rVsq(A)
782	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
783	cell = np.array(G2lat.Gmat2cell(g)) #real cell
784	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
785	scos = cosd(cellst[3:6])
786	ssin = sind(cellst[3:6])
787	scot = scos/ssin
788	rcos = cosd(cell[3:6])
789	rsin = sind(cell[3:6])
790	rcot = rcos/rsin
791	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
792	varyList = covData['varyList']
793	covMatrix = covData['covMatrix']
794	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
795	Ax = np.array(A)
796	Ax[3:] /= 2.
797	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
798	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
799	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
800	Vol = 1/np.sqrt(rVsq)
801	sigVol = Vol*3np.sqrt(srcvlsq)/2.
802	R123 = Ax[0]Ax[1]Ax[2]
803	dsasdg = np.zeros((3,6))
804	dadg = np.zeros((6,6))
805	for i0 in range(3): #0 1 2
806	i1 = (i0+1)%3 #1 2 0
807	i2 = (i1+1)%3 #2 0 1
808	i3 = 5-i2 #3 5 4
809	i4 = 5-i1 #4 3 5
810	i5 = 5-i0 #5 4 3
811	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
812	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
813	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
814	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
815	denom = np.sqrt(denmsq)
816	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
817	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
818	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
819	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
820	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
821	dadg[i5][i5] = -Ax[i0]/denom
822	for i0 in range(3):
823	i1 = (i0+1)%3
824	i2 = (i1+1)%3
825	i3 = 5-i2
826	for ij in range(6):
827	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
828	if ij == i0:
829	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
830	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
831	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
832	var = np.diag(sigMat)
833	CS = np.where(var>0.,np.sqrt(var),0.)
834	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
835	return cellSig
836
837	def SetPhaseData(parmDict,sigDict,Phases,covData):
838
839	def PrintAtomsAndSig(General,Atoms,atomsSig):
840	print '\n Atoms:'
841	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
842	if General['Type'] == 'magnetic':
843	line += ' Mx My Mz'
844	elif General['Type'] == 'macromolecular':
845	line = ' res no residue chain '+line
846	print line
847	if General['Type'] == 'nuclear':
848	print 135*'-'
849	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
850	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
851	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
852	for i,at in enumerate(Atoms):
853	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
854	valstr = ' values:'
855	sigstr = ' sig :'
856	for ind in [3,4,5,6]:
857	sigind = str(i)+':'+str(ind)
858	valstr += fmt[ind]%(at[ind])
859	if sigind in atomsSig:
860	sigstr += fmt[ind]%(atomsSig[sigind])
861	else:
862	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
863	if at[9] == 'I':
864	valstr += fmt[10]%(at[10])
865	if str(i)+':10' in atomsSig:
866	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
867	else:
868	sigstr += 8*' '
869	else:
870	valstr += 8*' '
871	sigstr += 8*' '
872	for ind in [11,12,13,14,15,16]:
873	sigind = str(i)+':'+str(ind)
874	valstr += fmt[ind]%(at[ind])
875	if sigind in atomsSig:
876	sigstr += fmt[ind]%(atomsSig[sigind])
877	else:
878	sigstr += 8*' '
879	print name
880	print valstr
881	print sigstr
882
883	def PrintSHtextureAndSig(textureData,SHtextureSig):
884	print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
885	names = ['omega','chi','phi']
886	namstr = ' names :'
887	ptstr = ' values:'
888	sigstr = ' esds :'
889	for name in names:
890	namstr += '%12s'%(name)
891	ptstr += '%12.3f'%(textureData['Sample '+name][1])
892	if 'Sample '+name in SHtextureSig:
893	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
894	else:
895	sigstr += 12*' '
896	print namstr
897	print ptstr
898	print sigstr
899	print '\n Texture coefficients:'
900	namstr = ' names :'
901	ptstr = ' values:'
902	sigstr = ' esds :'
903	SHcoeff = textureData['SH Coeff'][1]
904	for name in SHcoeff:
905	namstr += '%12s'%(name)
906	ptstr += '%12.3f'%(SHcoeff[name])
907	if name in SHtextureSig:
908	sigstr += '%12.3f'%(SHtextureSig[name])
909	else:
910	sigstr += 12*' '
911	print namstr
912	print ptstr
913	print sigstr
914
915
916	print '\n Phases:'
917	for phase in Phases:
918	print ' Result for phase: ',phase
919	Phase = Phases[phase]
920	General = Phase['General']
921	SGData = General['SGData']
922	Atoms = Phase['Atoms']
923	cell = General['Cell']
924	pId = Phase['pId']
925	pfx = str(pId)+'::'
926	if cell[0]:
927	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
928	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
929	print ' Reciprocal metric tensor: '
930	ptfmt = "%15.9f"
931	names = ['A11','A22','A33','A12','A13','A23']
932	namstr = ' names :'
933	ptstr = ' values:'
934	sigstr = ' esds :'
935	for name,a,siga in zip(names,A,sigA):
936	namstr += '%15s'%(name)
937	ptstr += ptfmt%(a)
938	if siga:
939	sigstr += ptfmt%(siga)
940	else:
941	sigstr += 15*' '
942	print namstr
943	print ptstr
944	print sigstr
945	cell[1:7] = G2lat.A2cell(A)
946	cell[7] = G2lat.calc_V(A)
947	print ' New unit cell:'
948	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
949	names = ['a','b','c','alpha','beta','gamma','Volume']
950	namstr = ' names :'
951	ptstr = ' values:'
952	sigstr = ' esds :'
953	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
954	namstr += '%12s'%(name)
955	ptstr += fmt%(a)
956	if siga:
957	sigstr += fmt%(siga)
958	else:
959	sigstr += 12*' '
960	print namstr
961	print ptstr
962	print sigstr
963
964	if Phase['General'].get('doPawley'):
965	pawleyRef = Phase['Pawley ref']
966	for i,refl in enumerate(pawleyRef):
967	key = pfx+'PWLref:'+str(i)
968	refl[6] = abs(parmDict[key]) #suppress negative Fsq
969	if key in sigDict:
970	refl[7] = sigDict[key]
971	else:
972	refl[7] = 0
973	else:
974	atomsSig = {}
975	if General['Type'] == 'nuclear':
976	for i,at in enumerate(Atoms):
977	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
978	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
979	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
980	for ind in [3,4,5,6]:
981	at[ind] = parmDict[names[ind]]
982	if ind in [3,4,5]:
983	name = names[ind].replace('A','dA')
984	else:
985	name = names[ind]
986	if name in sigDict:
987	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
988	if at[9] == 'I':
989	at[10] = parmDict[names[10]]
990	if names[10] in sigDict:
991	atomsSig[str(i)+':10'] = sigDict[names[10]]
992	else:
993	for ind in [11,12,13,14,15,16]:
994	at[ind] = parmDict[names[ind]]
995	if names[ind] in sigDict:
996	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
997	PrintAtomsAndSig(General,Atoms,atomsSig)
998
999	textureData = General['SH Texture']
1000	if textureData['Order']:
1001	SHtextureSig = {}
1002	for name in ['omega','chi','phi']:
1003	aname = pfx+'SH '+name
1004	textureData['Sample '+name][1] = parmDict[aname]
1005	if aname in sigDict:
1006	SHtextureSig['Sample '+name] = sigDict[aname]
1007	for name in textureData['SH Coeff'][1]:
1008	aname = pfx+name
1009	textureData['SH Coeff'][1][name] = parmDict[aname]
1010	if aname in sigDict:
1011	SHtextureSig[name] = sigDict[aname]
1012	PrintSHtextureAndSig(textureData,SHtextureSig)
1013
1014	################################################################################
1015	##### Histogram & Phase data
1016	################################################################################
1017
1018	def GetHistogramPhaseData(Phases,Histograms,Print=True):
1019
1020	def PrintSize(hapData):
1021	if hapData[0] in ['isotropic','uniaxial']:
1022	line = '\n Size model : %9s'%(hapData[0])
1023	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1024	if hapData[0] == 'uniaxial':
1025	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1026	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1027	print line
1028	else:
1029	print '\n Size model : %s'%(hapData[0])+ \
1030	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1031	Snames = ['S11','S22','S33','S12','S13','S23']
1032	ptlbls = ' names :'
1033	ptstr = ' values:'
1034	varstr = ' refine:'
1035	for i,name in enumerate(Snames):
1036	ptlbls += '%12s' % (name)
1037	ptstr += '%12.6f' % (hapData[4][i])
1038	varstr += '%12s' % (str(hapData[5][i]))
1039	print ptlbls
1040	print ptstr
1041	print varstr
1042
1043	def PrintMuStrain(hapData,SGData):
1044	if hapData[0] in ['isotropic','uniaxial']:
1045	line = '\n Mustrain model: %9s'%(hapData[0])
1046	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1047	if hapData[0] == 'uniaxial':
1048	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1049	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1050	print line
1051	else:
1052	print '\n Mustrain model: %s'%(hapData[0])+ \
1053	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1054	Snames = G2spc.MustrainNames(SGData)
1055	ptlbls = ' names :'
1056	ptstr = ' values:'
1057	varstr = ' refine:'
1058	for i,name in enumerate(Snames):
1059	ptlbls += '%12s' % (name)
1060	ptstr += '%12.6f' % (hapData[4][i])
1061	varstr += '%12s' % (str(hapData[5][i]))
1062	print ptlbls
1063	print ptstr
1064	print varstr
1065
1066	def PrintHStrain(hapData,SGData):
1067	print '\n Hydrostatic/elastic strain: '
1068	Hsnames = G2spc.HStrainNames(SGData)
1069	ptlbls = ' names :'
1070	ptstr = ' values:'
1071	varstr = ' refine:'
1072	for i,name in enumerate(Hsnames):
1073	ptlbls += '%12s' % (name)
1074	ptstr += '%12.6f' % (hapData[0][i])
1075	varstr += '%12s' % (str(hapData[1][i]))
1076	print ptlbls
1077	print ptstr
1078	print varstr
1079
1080	def PrintSHPO(hapData):
1081	print '\n Spherical harmonics preferred orientation: Order:' + \
1082	str(hapData[4])+' Refine? '+str(hapData[2])
1083	ptlbls = ' names :'
1084	ptstr = ' values:'
1085	for item in hapData[5]:
1086	ptlbls += '%12s'%(item)
1087	ptstr += '%12.3f'%(hapData[5][item])
1088	print ptlbls
1089	print ptstr
1090
1091	hapDict = {}
1092	hapVary = []
1093	controlDict = {}
1094	poType = {}
1095	poAxes = {}
1096	spAxes = {}
1097	spType = {}
1098
1099	for phase in Phases:
1100	HistoPhase = Phases[phase]['Histograms']
1101	SGData = Phases[phase]['General']['SGData']
1102	cell = Phases[phase]['General']['Cell'][1:7]
1103	A = G2lat.cell2A(cell)
1104	pId = Phases[phase]['pId']
1105	histoList = HistoPhase.keys()
1106	histoList.sort()
1107	for histogram in histoList:
1108	try:
1109	Histogram = Histograms[histogram]
1110	except KeyError:
1111	#skip if histogram not included e.g. in a sequential refinement
1112	continue
1113	hapData = HistoPhase[histogram]
1114	hId = Histogram['hId']
1115	limits = Histogram['Limits'][1]
1116	inst = Histogram['Instrument Parameters']
1117	inst = dict(zip(inst[3],inst[1]))
1118	Zero = inst['Zero']
1119	if 'C' in inst['Type']:
1120	try:
1121	wave = inst['Lam']
1122	except KeyError:
1123	wave = inst['Lam1']
1124	dmin = wave/(2.0*sind(limits[1]/2.0))
1125	pfx = str(pId)+':'+str(hId)+':'
1126	for item in ['Scale','Extinction']:
1127	hapDict[pfx+item] = hapData[item][0]
1128	if hapData[item][1]:
1129	hapVary.append(pfx+item)
1130	names = G2spc.HStrainNames(SGData)
1131	for i,name in enumerate(names):
1132	hapDict[pfx+name] = hapData['HStrain'][0][i]
1133	if hapData['HStrain'][1][i]:
1134	hapVary.append(pfx+name)
1135	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1136	if hapData['Pref.Ori.'][0] == 'MD':
1137	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1138	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1139	if hapData['Pref.Ori.'][2]:
1140	hapVary.append(pfx+'MD')
1141	else: #'SH' spherical harmonics
1142	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1143	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1144	for item in hapData['Pref.Ori.'][5]:
1145	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1146	if hapData['Pref.Ori.'][2]:
1147	hapVary.append(pfx+item)
1148	for item in ['Mustrain','Size']:
1149	controlDict[pfx+item+'Type'] = hapData[item][0]
1150	hapDict[pfx+item+':mx'] = hapData[item][1][2]
1151	if hapData[item][2][2]:
1152	hapVary.append(pfx+item+':mx')
1153	if hapData[item][0] in ['isotropic','uniaxial']:
1154	hapDict[pfx+item+':i'] = hapData[item][1][0]
1155	if hapData[item][2][0]:
1156	hapVary.append(pfx+item+':i')
1157	if hapData[item][0] == 'uniaxial':
1158	controlDict[pfx+item+'Axis'] = hapData[item][3]
1159	hapDict[pfx+item+':a'] = hapData[item][1][1]
1160	if hapData[item][2][1]:
1161	hapVary.append(pfx+item+':a')
1162	else: #generalized for mustrain or ellipsoidal for size
1163	Nterms = len(hapData[item][4])
1164	if item == 'Mustrain':
1165	names = G2spc.MustrainNames(SGData)
1166	pwrs = []
1167	for name in names:
1168	h,k,l = name[1:]
1169	pwrs.append([int(h),int(k),int(l)])
1170	controlDict[pfx+'MuPwrs'] = pwrs
1171	for i in range(Nterms):
1172	sfx = ':'+str(i)
1173	hapDict[pfx+item+sfx] = hapData[item][4][i]
1174	if hapData[item][5][i]:
1175	hapVary.append(pfx+item+sfx)
1176
1177	if Print:
1178	print '\n Phase: ',phase,' in histogram: ',histogram
1179	print 135*'-'
1180	print ' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1181	print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1182	if hapData['Pref.Ori.'][0] == 'MD':
1183	Ax = hapData['Pref.Ori.'][3]
1184	print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1185	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1186	else: #'SH' for spherical harmonics
1187	PrintSHPO(hapData['Pref.Ori.'])
1188	PrintSize(hapData['Size'])
1189	PrintMuStrain(hapData['Mustrain'],SGData)
1190	PrintHStrain(hapData['HStrain'],SGData)
1191	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1192	refList = []
1193	for h,k,l,d in HKLd:
1194	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1195	if ext:
1196	continue
1197	if 'C' in inst['Type']:
1198	pos = 2.0asind(wave/(2.0d))+Zero
1199	if limits[0] < pos < limits[1]:
1200	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1201	#last item should contain form factor values by atom type
1202	else:
1203	raise ValueError
1204	Histogram['Reflection Lists'][phase] = refList
1205	return hapVary,hapDict,controlDict
1206
1207	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1208
1209	def PrintSizeAndSig(hapData,sizeSig):
1210	line = '\n Size model: %9s'%(hapData[0])
1211	refine = False
1212	if hapData[0] in ['isotropic','uniaxial']:
1213	line += ' equatorial:%12.3f'%(hapData[1][0])
1214	if sizeSig[0][0]:
1215	line += ', sig:%8.3f'%(sizeSig[0][0])
1216	refine = True
1217	if hapData[0] == 'uniaxial':
1218	line += ' axial:%12.3f'%(hapData[1][1])
1219	if sizeSig[0][1]:
1220	refine = True
1221	line += ', sig:%8.3f'%(sizeSig[0][1])
1222	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1223	if sizeSig[0][2]:
1224	refine = True
1225	line += ', sig:%8.3f'%(sizeSig[0][2])
1226	if refine:
1227	print line
1228	else:
1229	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1230	if sizeSig[0][2]:
1231	refine = True
1232	line += ', sig:%8.3f'%(sizeSig[0][2])
1233	Snames = ['S11','S22','S33','S12','S13','S23']
1234	ptlbls = ' name :'
1235	ptstr = ' value :'
1236	sigstr = ' sig :'
1237	for i,name in enumerate(Snames):
1238	ptlbls += '%12s' % (name)
1239	ptstr += '%12.6f' % (hapData[4][i])
1240	if sizeSig[1][i]:
1241	refine = True
1242	sigstr += '%12.6f' % (sizeSig[1][i])
1243	else:
1244	sigstr += 12*' '
1245	if refine:
1246	print line
1247	print ptlbls
1248	print ptstr
1249	print sigstr
1250
1251	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1252	line = '\n Mustrain model: %9s'%(hapData[0])
1253	refine = False
1254	if hapData[0] in ['isotropic','uniaxial']:
1255	line += ' equatorial:%12.1f'%(hapData[1][0])
1256	if mustrainSig[0][0]:
1257	line += ', sig: %8.1f'%(mustrainSig[0][0])
1258	refine = True
1259	if hapData[0] == 'uniaxial':
1260	line += ' axial:%12.1f'%(hapData[1][1])
1261	if mustrainSig[0][1]:
1262	line += ', sig:%8.1f'%(mustrainSig[0][1])
1263	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1264	if mustrainSig[0][2]:
1265	refine = True
1266	line += ', sig:%8.3f'%(mustrainSig[0][2])
1267	if refine:
1268	print line
1269	else:
1270	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1271	if mustrainSig[0][2]:
1272	refine = True
1273	line += ', sig:%8.3f'%(mustrainSig[0][2])
1274	Snames = G2spc.MustrainNames(SGData)
1275	ptlbls = ' name :'
1276	ptstr = ' value :'
1277	sigstr = ' sig :'
1278	for i,name in enumerate(Snames):
1279	ptlbls += '%12s' % (name)
1280	ptstr += '%12.6f' % (hapData[4][i])
1281	if mustrainSig[1][i]:
1282	refine = True
1283	sigstr += '%12.6f' % (mustrainSig[1][i])
1284	else:
1285	sigstr += 12*' '
1286	if refine:
1287	print line
1288	print ptlbls
1289	print ptstr
1290	print sigstr
1291
1292	def PrintHStrainAndSig(hapData,strainSig,SGData):
1293	Hsnames = G2spc.HStrainNames(SGData)
1294	ptlbls = ' name :'
1295	ptstr = ' value :'
1296	sigstr = ' sig :'
1297	refine = False
1298	for i,name in enumerate(Hsnames):
1299	ptlbls += '%12s' % (name)
1300	ptstr += '%12.6g' % (hapData[0][i])
1301	if name in strainSig:
1302	refine = True
1303	sigstr += '%12.6g' % (strainSig[name])
1304	else:
1305	sigstr += 12*' '
1306	if refine:
1307	print '\n Hydrostatic/elastic strain: '
1308	print ptlbls
1309	print ptstr
1310	print sigstr
1311
1312	def PrintSHPOAndSig(hapData,POsig):
1313	print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1314	ptlbls = ' names :'
1315	ptstr = ' values:'
1316	sigstr = ' sig :'
1317	for item in hapData[5]:
1318	ptlbls += '%12s'%(item)
1319	ptstr += '%12.3f'%(hapData[5][item])
1320	if item in POsig:
1321	sigstr += '%12.3f'%(POsig[item])
1322	else:
1323	sigstr += 12*' '
1324	print ptlbls
1325	print ptstr
1326	print sigstr
1327
1328	for phase in Phases:
1329	HistoPhase = Phases[phase]['Histograms']
1330	SGData = Phases[phase]['General']['SGData']
1331	pId = Phases[phase]['pId']
1332	histoList = HistoPhase.keys()
1333	histoList.sort()
1334	for histogram in histoList:
1335	try:
1336	Histogram = Histograms[histogram]
1337	except KeyError:
1338	#skip if histogram not included e.g. in a sequential refinement
1339	continue
1340	print '\n Phase: ',phase,' in histogram: ',histogram
1341	print 130*'-'
1342	hapData = HistoPhase[histogram]
1343	hId = Histogram['hId']
1344	pfx = str(pId)+':'+str(hId)+':'
1345	print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1346	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1347
1348	PhFrExtPOSig = {}
1349	for item in ['Scale','Extinction']:
1350	hapData[item][0] = parmDict[pfx+item]
1351	if pfx+item in sigDict:
1352	PhFrExtPOSig[item] = sigDict[pfx+item]
1353	if hapData['Pref.Ori.'][0] == 'MD':
1354	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1355	if pfx+'MD' in sigDict:
1356	PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1357	else: #'SH' spherical harmonics
1358	for item in hapData['Pref.Ori.'][5]:
1359	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1360	if pfx+item in sigDict:
1361	PhFrExtPOSig[item] = sigDict[pfx+item]
1362	if Print:
1363	if 'Scale' in PhFrExtPOSig:
1364	print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1365	if 'Extinction' in PhFrExtPOSig:
1366	print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1367	if hapData['Pref.Ori.'][0] == 'MD':
1368	if 'MD' in PhFrExtPOSig:
1369	print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1370	else:
1371	PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1372	SizeMuStrSig = {'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1373	'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1374	'HStrain':{}}
1375	for item in ['Mustrain','Size']:
1376	hapData[item][1][2] = parmDict[pfx+item+':mx']
1377	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1378	if pfx+item+':mx' in sigDict:
1379	SizeMuStrSig[item][0][2] = sigDict[pfx+item+':mx']
1380	if hapData[item][0] in ['isotropic','uniaxial']:
1381	hapData[item][1][0] = parmDict[pfx+item+':i']
1382	if item == 'Size':
1383	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1384	if pfx+item+':i' in sigDict:
1385	SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
1386	if hapData[item][0] == 'uniaxial':
1387	hapData[item][1][1] = parmDict[pfx+item+':a']
1388	if item == 'Size':
1389	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
1390	if pfx+item+':a' in sigDict:
1391	SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
1392	else: #generalized for mustrain or ellipsoidal for size
1393	Nterms = len(hapData[item][4])
1394	for i in range(Nterms):
1395	sfx = ':'+str(i)
1396	hapData[item][4][i] = parmDict[pfx+item+sfx]
1397	if pfx+item+sfx in sigDict:
1398	SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1399	names = G2spc.HStrainNames(SGData)
1400	for i,name in enumerate(names):
1401	hapData['HStrain'][0][i] = parmDict[pfx+name]
1402	if pfx+name in sigDict:
1403	SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1404	if Print:
1405	PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1406	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1407	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1408
1409	################################################################################
1410	##### Histogram data
1411	################################################################################
1412
1413	def GetHistogramData(Histograms,Print=True):
1414
1415	def GetBackgroundParms(hId,Background):
1416	Back = Background[0]
1417	Debye = Background[1]
1418	bakType,bakFlag = Back[:2]
1419	backVals = Back[3:]
1420	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1421	backDict = dict(zip(backNames,backVals))
1422	backVary = []
1423	if bakFlag:
1424	backVary = backNames
1425	backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
1426	debyeDict = {}
1427	debyeList = []
1428	for i in range(Debye['nDebye']):
1429	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1430	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1431	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1432	debyeVary = []
1433	for item in debyeList:
1434	if item[1]:
1435	debyeVary.append(item[0])
1436	backDict.update(debyeDict)
1437	backVary += debyeVary
1438	return bakType,backDict,backVary
1439
1440	def GetInstParms(hId,Inst):
1441	insVals,insFlags,insNames = Inst[1:4]
1442	dataType = insVals[0]
1443	instDict = {}
1444	insVary = []
1445	pfx = ':'+str(hId)+':'
1446	for i,flag in enumerate(insFlags):
1447	insName = pfx+insNames[i]
1448	instDict[insName] = insVals[i]
1449	if flag:
1450	insVary.append(insName)
1451	instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1452	instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1453	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1454	return dataType,instDict,insVary
1455
1456	def GetSampleParms(hId,Sample):
1457	sampVary = []
1458	hfx = ':'+str(hId)+':'
1459	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1460	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1461	Type = Sample['Type']
1462	if 'Bragg' in Type: #Bragg-Brentano
1463	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
1464	sampDict[hfx+item] = Sample[item][0]
1465	if Sample[item][1]:
1466	sampVary.append(hfx+item)
1467	elif 'Debye' in Type: #Debye-Scherrer
1468	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1469	sampDict[hfx+item] = Sample[item][0]
1470	if Sample[item][1]:
1471	sampVary.append(hfx+item)
1472	return Type,sampDict,sampVary
1473
1474	def PrintBackground(Background):
1475	Back = Background[0]
1476	Debye = Background[1]
1477	print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
1478	line = ' Coefficients: '
1479	for i,back in enumerate(Back[3:]):
1480	line += '%10.3f'%(back)
1481	if i and not i%10:
1482	line += '\n'+15*' '
1483	print line
1484	if Debye['nDebye']:
1485	print '\n Debye diffuse scattering coefficients'
1486	parms = ['DebyeA','DebyeR','DebyeU']
1487	line = ' names :'
1488	for parm in parms:
1489	line += '%16s'%(parm)
1490	print line
1491	for j,term in enumerate(Debye['debyeTerms']):
1492	line = ' term'+'%2d'%(j)+':'
1493	for i in range(3):
1494	line += '%10.4g %5s'%(term[2i],bool(term[2i+1]))
1495	print line
1496
1497	def PrintInstParms(Inst):
1498	print '\n Instrument Parameters:'
1499	ptlbls = ' name :'
1500	ptstr = ' value :'
1501	varstr = ' refine:'
1502	instNames = Inst[3][1:]
1503	for i,name in enumerate(instNames):
1504	ptlbls += '%12s' % (name)
1505	ptstr += '%12.6f' % (Inst[1][i+1])
1506	if name in ['Lam1','Lam2','Azimuth']:
1507	varstr += 12*' '
1508	else:
1509	varstr += '%12s' % (str(bool(Inst[2][i+1])))
1510	print ptlbls
1511	print ptstr
1512	print varstr
1513
1514	def PrintSampleParms(Sample):
1515	print '\n Sample Parameters:'
1516	print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1517	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
1518	ptlbls = ' name :'
1519	ptstr = ' value :'
1520	varstr = ' refine:'
1521	if 'Bragg' in Sample['Type']:
1522	for item in ['Scale','Shift','Transparency']:
1523	ptlbls += '%14s'%(item)
1524	ptstr += '%14.4f'%(Sample[item][0])
1525	varstr += '%14s'%(str(bool(Sample[item][1])))
1526
1527	elif 'Debye' in Type: #Debye-Scherrer
1528	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1529	ptlbls += '%14s'%(item)
1530	ptstr += '%14.4f'%(Sample[item][0])
1531	varstr += '%14s'%(str(bool(Sample[item][1])))
1532
1533	print ptlbls
1534	print ptstr
1535	print varstr
1536
1537
1538	histDict = {}
1539	histVary = []
1540	controlDict = {}
1541	histoList = Histograms.keys()
1542	histoList.sort()
1543	for histogram in histoList:
1544	Histogram = Histograms[histogram]
1545	hId = Histogram['hId']
1546	pfx = ':'+str(hId)+':'
1547	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1548
1549	Background = Histogram['Background']
1550	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1551	controlDict[pfx+'bakType'] = Type
1552	histDict.update(bakDict)
1553	histVary += bakVary
1554
1555	Inst = Histogram['Instrument Parameters']
1556	Type,instDict,insVary = GetInstParms(hId,Inst)
1557	controlDict[pfx+'histType'] = Type
1558	if pfx+'Lam1' in instDict:
1559	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1560	else:
1561	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
1562	histDict.update(instDict)
1563	histVary += insVary
1564
1565	Sample = Histogram['Sample Parameters']
1566	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1567	controlDict[pfx+'instType'] = Type
1568	histDict.update(sampDict)
1569	histVary += sampVary
1570
1571	if Print:
1572	print '\n Histogram: ',histogram,' histogram Id: ',hId
1573	print 135*'-'
1574	Units = {'C':' deg','T':' msec'}
1575	units = Units[controlDict[pfx+'histType'][2]]
1576	Limits = controlDict[pfx+'Limits']
1577	print ' Instrument type: ',Sample['Type']
1578	print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
1579	PrintSampleParms(Sample)
1580	PrintInstParms(Inst)
1581	PrintBackground(Background)
1582
1583	return histVary,histDict,controlDict
1584
1585	def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1586
1587	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1588	Back = Background[0]
1589	Debye = Background[1]
1590	lenBack = len(Back[3:])
1591	backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
1592	for i in range(lenBack):
1593	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1594	if pfx+'Back:'+str(i) in sigDict:
1595	backSig[i] = sigDict[pfx+'Back:'+str(i)]
1596	if Debye['nDebye']:
1597	for i in range(Debye['nDebye']):
1598	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1599	for j,name in enumerate(names):
1600	Debye['debyeTerms'][i][2*j] = parmDict[name]
1601	if name in sigDict:
1602	backSig[lenBack+3*i+j] = sigDict[name]
1603	return backSig
1604
1605	def SetInstParms(pfx,Inst,parmDict,sigDict):
1606	insVals,insFlags,insNames = Inst[1:4]
1607	instSig = [0 for i in range(len(insVals))]
1608	for i,flag in enumerate(insFlags):
1609	insName = pfx+insNames[i]
1610	insVals[i] = parmDict[insName]
1611	if insName in sigDict:
1612	instSig[i] = sigDict[insName]
1613	return instSig
1614
1615	def SetSampleParms(pfx,Sample,parmDict,sigDict):
1616	if 'Bragg' in Sample['Type']: #Bragg-Brentano
1617	sampSig = [0 for i in range(3)]
1618	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
1619	Sample[item][0] = parmDict[pfx+item]
1620	if pfx+item in sigDict:
1621	sampSig[i] = sigDict[pfx+item]
1622	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1623	sampSig = [0 for i in range(4)]
1624	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1625	Sample[item][0] = parmDict[pfx+item]
1626	if pfx+item in sigDict:
1627	sampSig[i] = sigDict[pfx+item]
1628	return sampSig
1629
1630	def PrintBackgroundSig(Background,backSig):
1631	Back = Background[0]
1632	Debye = Background[1]
1633	lenBack = len(Back[3:])
1634	valstr = ' value : '
1635	sigstr = ' sig : '
1636	refine = False
1637	for i,back in enumerate(Back[3:]):
1638	valstr += '%10.4g'%(back)
1639	if Back[1]:
1640	refine = True
1641	sigstr += '%10.4g'%(backSig[i])
1642	else:
1643	sigstr += 10*' '
1644	if refine:
1645	print '\n Background function: ',Back[0]
1646	print valstr
1647	print sigstr
1648	if Debye['nDebye']:
1649	ifAny = False
1650	ptfmt = "%12.5f"
1651	names = ' names :'
1652	ptstr = ' values:'
1653	sigstr = ' esds :'
1654	for item in sigDict:
1655	if 'Debye' in item:
1656	ifAny = True
1657	names += '%12s'%(item)
1658	ptstr += ptfmt%(parmDict[item])
1659	sigstr += ptfmt%(sigDict[item])
1660	if ifAny:
1661	print '\n Debye diffuse scattering coefficients'
1662	print names
1663	print ptstr
1664	print sigstr
1665
1666	def PrintInstParmsSig(Inst,instSig):
1667	ptlbls = ' names :'
1668	ptstr = ' value :'
1669	sigstr = ' sig :'
1670	instNames = Inst[3][1:]
1671	refine = False
1672	for i,name in enumerate(instNames):
1673	ptlbls += '%12s' % (name)
1674	ptstr += '%12.6f' % (Inst[1][i+1])
1675	if instSig[i+1]:
1676	refine = True
1677	sigstr += '%12.6f' % (instSig[i+1])
1678	else:
1679	sigstr += 12*' '
1680	if refine:
1681	print '\n Instrument Parameters:'
1682	print ptlbls
1683	print ptstr
1684	print sigstr
1685
1686	def PrintSampleParmsSig(Sample,sampleSig):
1687	ptlbls = ' names :'
1688	ptstr = ' values:'
1689	sigstr = ' sig :'
1690	refine = False
1691	if 'Bragg' in Sample['Type']:
1692	for i,item in enumerate(['Scale','Shift','Transparency']):
1693	ptlbls += '%14s'%(item)
1694	ptstr += '%14.4f'%(Sample[item][0])
1695	if sampleSig[i]:
1696	refine = True
1697	sigstr += '%14.4f'%(sampleSig[i])
1698	else:
1699	sigstr += 14*' '
1700
1701	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1702	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1703	ptlbls += '%14s'%(item)
1704	ptstr += '%14.4f'%(Sample[item][0])
1705	if sampleSig[i]:
1706	refine = True
1707	sigstr += '%14.4f'%(sampleSig[i])
1708	else:
1709	sigstr += 14*' '
1710
1711	if refine:
1712	print '\n Sample Parameters:'
1713	print ptlbls
1714	print ptstr
1715	print sigstr
1716
1717	histoList = Histograms.keys()
1718	histoList.sort()
1719	for histogram in histoList:
1720	if 'PWDR' in histogram:
1721	Histogram = Histograms[histogram]
1722	hId = Histogram['hId']
1723	pfx = ':'+str(hId)+':'
1724	Background = Histogram['Background']
1725	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1726
1727	Inst = Histogram['Instrument Parameters']
1728	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1729
1730	Sample = Histogram['Sample Parameters']
1731	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1732
1733	print '\n Histogram: ',histogram,' histogram Id: ',hId
1734	print 135*'-'
1735	print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1736	if Print:
1737	print ' Instrument type: ',Sample['Type']
1738	PrintSampleParmsSig(Sample,sampSig)
1739	PrintInstParmsSig(Inst,instSig)
1740	PrintBackgroundSig(Background,backSig)
1741
1742	################################################################################
1743	##### Function & derivative calculations
1744	################################################################################
1745
1746	def GetAtomFXU(pfx,calcControls,parmDict):
1747	Natoms = calcControls['Natoms'][pfx]
1748	Tdata = Natoms*[' ',]
1749	Mdata = np.zeros(Natoms)
1750	IAdata = Natoms*[' ',]
1751	Fdata = np.zeros(Natoms)
1752	FFdata = []
1753	BLdata = []
1754	Xdata = np.zeros((3,Natoms))
1755	dXdata = np.zeros((3,Natoms))
1756	Uisodata = np.zeros(Natoms)
1757	Uijdata = np.zeros((6,Natoms))
1758	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1759	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1760	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1761	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1762	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1763	for iatm in range(Natoms):
1764	for key in keys:
1765	parm = pfx+key+str(iatm)
1766	if parm in parmDict:
1767	keys[key][iatm] = parmDict[parm]
1768	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1769
1770	def getFFvalues(FFtables,SQ):
1771	FFvals = {}
1772	for El in FFtables:
1773	FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
1774	return FFvals
1775
1776	def getBLvalues(BLtables):
1777	BLvals = {}
1778	for El in BLtables:
1779	BLvals[El] = BLtables[El][1][1]
1780	return BLvals
1781
1782	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1783	''' Compute structure factors for all h,k,l for phase
1784	input:
1785	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1786	G: reciprocal metric tensor
1787	pfx: phase id string
1788	SGData: space group info. dictionary output from SpcGroup
1789	calcControls:
1790	ParmDict:
1791	puts result F^2 in each ref[8] in refList
1792	'''
1793	twopi = 2.0*np.pi
1794	twopisq = 2.0np.pi*2
1795	ast = np.sqrt(np.diag(G))
1796	Mast = twopisq*np.multiply.outer(ast,ast)
1797	FFtables = calcControls['FFtables']
1798	BLtables = calcControls['BLtables']
1799	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1800	FF = np.zeros(len(Tdata))
1801	if 'N' in parmDict[hfx+'Type']:
1802	FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
1803	else:
1804	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1805	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1806	maxPos = len(SGData['SGOps'])
1807	Uij = np.array(G2lat.U6toUij(Uijdata))
1808	bij = Mast*Uij.T
1809	for refl in refList:
1810	fbs = np.array([0,0])
1811	H = refl[:3]
1812	SQ = 1./(2.refl[4])*2
1813	if not len(refl[-1]): #no form factors
1814	if 'N' in parmDict[hfx+'Type']:
1815	refl[-1] = getBLvalues(BLtables)
1816	else: #'X'
1817	refl[-1] = getFFvalues(FFtables,SQ)
1818	for i,El in enumerate(Tdata):
1819	FF[i] = refl[-1][El]
1820	SQfactor = 4.0SQtwopisq
1821	Uniq = refl[11]
1822	phi = refl[12]
1823	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1824	sinp = np.sin(phase)
1825	cosp = np.cos(phase)
1826	occ = Mdata*Fdata/len(Uniq)
1827	biso = -SQfactor*Uisodata
1828	Tiso = np.where(biso<1.,np.exp(biso),1.0)
1829	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1830	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1831	Tcorr = Tiso*Tuij
1832	fa = np.array([(FF+FP)occcospTcorr,-FPPoccsinpTcorr])
1833	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
1834	if not SGData['SGInv']:
1835	fb = np.array([(FF+FP)occsinpTcorr,FPPocccospTcorr])
1836	fbs = np.sum(np.sum(fb,axis=1),axis=1)
1837	fasq = fas**2
1838	fbsq = fbs**2 #imaginary
1839	refl[9] = np.sum(fasq)+np.sum(fbsq)
1840	refl[10] = atan2d(fbs[0],fas[0])
1841	return refList
1842
1843	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1844	twopi = 2.0*np.pi
1845	twopisq = 2.0np.pi*2
1846	ast = np.sqrt(np.diag(G))
1847	Mast = twopisq*np.multiply.outer(ast,ast)
1848	FFtables = calcControls['FFtables']
1849	BLtables = calcControls['BLtables']
1850	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1851	FF = np.zeros(len(Tdata))
1852	if 'N' in parmDict[hfx+'Type']:
1853	FP = 0.
1854	FPP = 0.
1855	else:
1856	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1857	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1858	maxPos = len(SGData['SGOps'])
1859	Uij = np.array(G2lat.U6toUij(Uijdata))
1860	bij = Mast*Uij.T
1861	dFdvDict = {}
1862	dFdfr = np.zeros((len(refList),len(Mdata)))
1863	dFdx = np.zeros((len(refList),len(Mdata),3))
1864	dFdui = np.zeros((len(refList),len(Mdata)))
1865	dFdua = np.zeros((len(refList),len(Mdata),6))
1866	for iref,refl in enumerate(refList):
1867	H = np.array(refl[:3])
1868	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
1869	for i,El in enumerate(Tdata):
1870	FF[i] = refl[-1][El]
1871	SQfactor = 8.0SQnp.pi**2
1872	Uniq = refl[11]
1873	phi = refl[12]
1874	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1875	sinp = np.sin(phase)
1876	cosp = np.cos(phase)
1877	occ = Mdata*Fdata/len(Uniq)
1878	biso = -SQfactor*Uisodata
1879	Tiso = np.where(biso<1.,np.exp(biso),1.0)
1880	# HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1881	HbH = -np.inner(H,np.inner(bij,H))
1882	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1883	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1884	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1885	Tcorr = Tiso*Tuij
1886	fot = (FF+FP)occTcorr
1887	fotp = FPPoccTcorr
1888	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
1889	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
1890
1891	fas = np.sum(np.sum(fa,axis=1),axis=1)
1892	fbs = np.sum(np.sum(fb,axis=1),axis=1)
1893	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
1894	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
1895	#sum below is over Uniq
1896	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1897	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
1898	dfadui = np.sum(-SQfactor*fa,axis=2)
1899	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1900	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
1901	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
1902	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
1903	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
1904	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
1905	if not SGData['SGInv']:
1906	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
1907	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
1908	dfbdui = np.sum(-SQfactor*fb,axis=2)
1909	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1910	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
1911	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
1912	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
1913	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
1914	#loop over atoms - each dict entry is list of derivatives for all the reflections
1915	for i in range(len(Mdata)):
1916	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1917	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1918	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1919	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1920	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1921	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1922	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1923	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1924	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1925	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1926	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1927	return dFdvDict
1928
1929	def Dict2Values(parmdict, varylist):
1930	'''Use before call to leastsq to setup list of values for the parameters
1931	in parmdict, as selected by key in varylist'''
1932	return [parmdict[key] for key in varylist]
1933
1934	def Values2Dict(parmdict, varylist, values):
1935	''' Use after call to leastsq to update the parameter dictionary with
1936	values corresponding to keys in varylist'''
1937	parmdict.update(zip(varylist,values))
1938
1939	def GetNewCellParms(parmDict,varyList):
1940	newCellDict = {}
1941	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
1942	for item in varyList:
1943	keys = item.split(':')
1944	if keys[2] in Ddict:
1945	key = keys[0]+'::'+Ddict[keys[2]]
1946	parm = keys[0]+'::'+keys[2]
1947	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1948	return newCellDict
1949
1950	def ApplyXYZshifts(parmDict,varyList):
1951	''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1952	input:
1953	parmDict - parameter dictionary
1954	varyList - list of variables
1955	returns:
1956	newAtomDict - dictitemionary of new atomic coordinate names & values;
1957	key is parameter shift name
1958	'''
1959	newAtomDict = {}
1960	for item in parmDict:
1961	if 'dA' in item:
1962	parm = ''.join(item.split('d'))
1963	parmDict[parm] += parmDict[item]
1964	newAtomDict[item] = [parm,parmDict[parm]]
1965	return newAtomDict
1966
1967	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1968	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1969	odfCor = 1.0
1970	H = refl[:3]
1971	cell = G2lat.Gmat2cell(g)
1972	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1973	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1974	phi,beta = G2lat.CrsAng(H,cell,SGData)
1975	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1976	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1977	for item in SHnames:
1978	L,M,N = eval(item.strip('C'))
1979	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1980	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1981	Lnorm = G2lat.Lnorm(L)
1982	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1983	return odfCor
1984
1985	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1986	FORPI = 12.5663706143592
1987	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1988	odfCor = 1.0
1989	dFdODF = {}
1990	dFdSA = [0,0,0]
1991	H = refl[:3]
1992	cell = G2lat.Gmat2cell(g)
1993	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1994	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1995	phi,beta = G2lat.CrsAng(H,cell,SGData)
1996	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1997	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1998	for item in SHnames:
1999	L,M,N = eval(item.strip('C'))
2000	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2001	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2002	Lnorm = G2lat.Lnorm(L)
2003	odfCor += parmDict[pfx+item]LnormKcl*Ksl
2004	dFdODF[pfx+item] = LnormKclKsl
2005	for i in range(3):
2006	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
2007	return odfCor,dFdODF,dFdSA
2008
2009	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2010	odfCor = 1.0
2011	H = refl[:3]
2012	cell = G2lat.Gmat2cell(g)
2013	Sangl = [0.,0.,0.]
2014	if 'Bragg' in calcControls[hfx+'instType']:
2015	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2016	IFCoup = True
2017	else:
2018	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2019	IFCoup = False
2020	phi,beta = G2lat.CrsAng(H,cell,SGData)
2021	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2022	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2023	for item in SHnames:
2024	L,N = eval(item.strip('C'))
2025	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2026	odfCor += parmDict[phfx+item]LnormKcsl
2027	return odfCor
2028
2029	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2030	FORPI = 12.5663706143592
2031	odfCor = 1.0
2032	dFdODF = {}
2033	H = refl[:3]
2034	cell = G2lat.Gmat2cell(g)
2035	Sangl = [0.,0.,0.]
2036	if 'Bragg' in calcControls[hfx+'instType']:
2037	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2038	IFCoup = True
2039	else:
2040	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2041	IFCoup = False
2042	phi,beta = G2lat.CrsAng(H,cell,SGData)
2043	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2044	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2045	for item in SHnames:
2046	L,N = eval(item.strip('C'))
2047	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2048	odfCor += parmDict[phfx+item]LnormKcsl
2049	dFdODF[phfx+item] = Kcsl*Lnorm
2050	return odfCor,dFdODF
2051
2052	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2053	POcorr = 1.0
2054	if calcControls[phfx+'poType'] == 'MD':
2055	MD = parmDict[phfx+'MD']
2056	if MD != 1.0:
2057	MDAxis = calcControls[phfx+'MDAxis']
2058	sumMD = 0
2059	for H in refl[11]:
2060	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2061	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2062	sumMD += A**3
2063	POcorr = sumMD/len(refl[11])
2064	else: #spherical harmonics
2065	if calcControls[pfx+'SHord']:
2066	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2067	return POcorr
2068
2069	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2070	POcorr = 1.0
2071	POderv = {}
2072	if calcControls[phfx+'poType'] == 'MD':
2073	MD = parmDict[phfx+'MD']
2074	MDAxis = calcControls[phfx+'MDAxis']
2075	sumMD = 0
2076	sumdMD = 0
2077	for H in refl[11]:
2078	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2079	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2080	sumMD += A**3
2081	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
2082	POcorr = sumMD/len(refl[11])
2083	POderv[phfx+'MD'] = sumdMD/len(refl[11])
2084	else: #spherical harmonics
2085	if calcControls[pfx+'SHord']:
2086	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2087	return POcorr,POderv
2088
2089	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2090	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
2091	if 'X' in parmDict[hfx+'Type']:
2092	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2093	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2094	if pfx+'SHorder' in parmDict:
2095	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2096	refl[13] = Icorr
2097
2098	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2099	dIdsh = 1./parmDict[hfx+'Scale']
2100	dIdsp = 1./parmDict[phfx+'Scale']
2101	if 'X' in parmDict[hfx+'Type']:
2102	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2103	dIdPola /= pola
2104	else: #'N'
2105	dIdPola = 0.0
2106	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2107	for iPO in dIdPO:
2108	dIdPO[iPO] /= POcorr
2109	dFdODF = {}
2110	dFdSA = [0,0,0]
2111	if pfx+'SHorder' in parmDict:
2112	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2113	for iSH in dFdODF:
2114	dFdODF[iSH] /= odfCor
2115	for i in range(3):
2116	dFdSA[i] /= odfCor
2117	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
2118
2119	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2120	costh = cosd(refl[5]/2.)
2121	#crystallite size
2122	if calcControls[phfx+'SizeType'] == 'isotropic':
2123	Sgam = 1.8wave/(np.piparmDict[phfx+'Size:i']*costh)
2124	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2125	H = np.array(refl[:3])
2126	P = np.array(calcControls[phfx+'SizeAxis'])
2127	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2128	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size:i']parmDict[phfx+'Size:a']*costh)
2129	Sgam = np.sqrt((sinPparmDict[phfx+'Size:a'])*2+(cosPparmDict[phfx+'Size:i'])**2)
2130	else: #ellipsoidal crystallites
2131	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2132	H = np.array(refl[:3])
2133	lenR = G2pwd.ellipseSize(H,Sij,GB)
2134	Sgam = 1.8wave/(np.picosth*lenR)
2135	#microstrain
2136	if calcControls[phfx+'MustrainType'] == 'isotropic':
2137	Mgam = 0.018parmDict[phfx+'Mustrain:i']tand(refl[5]/2.)/np.pi
2138	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2139	H = np.array(refl[:3])
2140	P = np.array(calcControls[phfx+'MustrainAxis'])
2141	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2142	Si = parmDict[phfx+'Mustrain:i']
2143	Sa = parmDict[phfx+'Mustrain:a']
2144	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
2145	else: #generalized - P.W. Stephens model
2146	pwrs = calcControls[phfx+'MuPwrs']
2147	sum = 0
2148	for i,pwr in enumerate(pwrs):
2149	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2150	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2151	gam = SgamparmDict[phfx+'Size:mx']+MgamparmDict[phfx+'Mustrain:mx']
2152	sig = (Sgam(1.-parmDict[phfx+'Size:mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain:mx']))**2
2153	sig /= ateln2
2154	return sig,gam
2155
2156	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2157	gamDict = {}
2158	sigDict = {}
2159	costh = cosd(refl[5]/2.)
2160	tanth = tand(refl[5]/2.)
2161	#crystallite size derivatives
2162	if calcControls[phfx+'SizeType'] == 'isotropic':
2163	Sgam = 1.8wave/(np.piparmDict[phfx+'Size:i']*costh)
2164	gamDict[phfx+'Size:i'] = -1.80waveparmDict[phfx+'Size:mx']/(np.pi*costh)
2165	sigDict[phfx+'Size:i'] = -3.60Sgamwave(1.-parmDict[phfx+'Size:mx'])2/(np.picosth*ateln2)
2166	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2167	H = np.array(refl[:3])
2168	P = np.array(calcControls[phfx+'SizeAxis'])
2169	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2170	Si = parmDict[phfx+'Size:i']
2171	Sa = parmDict[phfx+'Size:a']
2172	gami = (1.8wave/np.pi)/(SiSa)
2173	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
2174	Sgam = gami*sqtrm
2175	gam = Sgam/costh
2176	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
2177	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
2178	gamDict[phfx+'Size:i'] = dsi*parmDict[phfx+'Size:mx']
2179	gamDict[phfx+'Size:a'] = dsa*parmDict[phfx+'Size:mx']
2180	sigDict[phfx+'Size:i'] = 2.dsiSgam(1.-parmDict[phfx+'Size:mx'])*2/ateln2
2181	sigDict[phfx+'Size:a'] = 2.dsaSgam(1.-parmDict[phfx+'Size:mx'])*2/ateln2
2182	else: #ellipsoidal crystallites
2183	const = 1.8wave/(np.picosth)
2184	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2185	H = np.array(refl[:3])
2186	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2187	Sgam = 1.8wave/(np.picosth*lenR)
2188	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2189	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size:mx']
2190	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size:mx'])*2/ateln2
2191	gamDict[phfx+'Size:mx'] = Sgam
2192	sigDict[phfx+'Size:mx'] = -2.Sgam2(1.-parmDict[phfx+'Size:mx'])/ateln2
2193
2194	#microstrain derivatives
2195	if calcControls[phfx+'MustrainType'] == 'isotropic':
2196	Mgam = 0.018parmDict[phfx+'Mustrain:i']tand(refl[5]/2.)/np.pi
2197	gamDict[phfx+'Mustrain:i'] = 0.018tanthparmDict[phfx+'Mustrain:mx']/np.pi
2198	sigDict[phfx+'Mustrain:i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain:mx'])2/(np.piateln2)
2199	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2200	H = np.array(refl[:3])
2201	P = np.array(calcControls[phfx+'MustrainAxis'])
2202	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2203	Si = parmDict[phfx+'Mustrain:i']
2204	Sa = parmDict[phfx+'Mustrain:a']
2205	gami = 0.018SiSa*tanth/np.pi
2206	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
2207	Mgam = gami/sqtrm
2208	dsi = -gamiSicosP2/sqtrm3
2209	dsa = -gamiSasinP2/sqtrm3
2210	gamDict[phfx+'Mustrain:i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain:mx']
2211	gamDict[phfx+'Mustrain:a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain:mx']
2212	sigDict[phfx+'Mustrain:i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2213	sigDict[phfx+'Mustrain:a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2214	else: #generalized - P.W. Stephens model
2215	pwrs = calcControls[phfx+'MuPwrs']
2216	const = 0.018refl[4]2tanth
2217	sum = 0
2218	for i,pwr in enumerate(pwrs):
2219	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2220	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2221	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain:mx']
2222	sigDict[phfx+'Mustrain:'+str(i)] = \
2223	2.constterm(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2224	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2225	for i in range(len(pwrs)):
2226	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2227	gamDict[phfx+'Mustrain:mx'] = Mgam
2228	sigDict[phfx+'Mustrain:mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain:mx'])/ateln2
2229	return sigDict,gamDict
2230
2231	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2232	h,k,l = refl[:3]
2233	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2234	d = np.sqrt(dsq)
2235
2236	refl[4] = d
2237	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
2238	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2239	if 'Bragg' in calcControls[hfx+'instType']:
2240	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2241	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
2242	else: #Debye-Scherrer - simple but maybe not right
2243	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2244	return pos
2245
2246	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2247	dpr = 180./np.pi
2248	h,k,l = refl[:3]
2249	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2250	dst = np.sqrt(dstsq)
2251	pos = refl[5]-parmDict[hfx+'Zero']
2252	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
2253	dpdw = const*dst
2254	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
2255	dpdA = constwave/(2.0*dst)
2256	dpdZ = 1.0
2257	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2258	if 'Bragg' in calcControls[hfx+'instType']:
2259	dpdSh = -4.constcosd(pos/2.0)
2260	dpdTr = -constsind(pos)100.0
2261	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2262	else: #Debye-Scherrer - simple but maybe not right
2263	dpdXd = -const*cosd(pos)
2264	dpdYd = -const*sind(pos)
2265	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2266
2267	def GetHStrainShift(refl,SGData,phfx,parmDict):
2268	laue = SGData['SGLaue']
2269	uniq = SGData['SGUniq']
2270	h,k,l = refl[:3]
2271	if laue in ['m3','m3m']:
2272	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2273	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2274	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2275	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2276	elif laue in ['3R','3mR']:
2277	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2278	elif laue in ['4/m','4/mmm']:
2279	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2280	elif laue in ['mmm']:
2281	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2282	elif laue in ['2/m']:
2283	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2284	if uniq == 'a':
2285	Dij += parmDict[phfx+'D23']kl
2286	elif uniq == 'b':
2287	Dij += parmDict[phfx+'D13']hl
2288	elif uniq == 'c':
2289	Dij += parmDict[phfx+'D12']hk
2290	else:
2291	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2292	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2293	return Dijrefl[4]2tand(refl[5]/2.0)
2294
2295	def GetHStrainShiftDerv(refl,SGData,phfx):
2296	laue = SGData['SGLaue']
2297	uniq = SGData['SGUniq']
2298	h,k,l = refl[:3]
2299	if laue in ['m3','m3m']:
2300	dDijDict = {phfx+'D11':h2+k2+l**2,
2301	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2302	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2303	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2304	elif laue in ['3R','3mR']:
2305	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2306	elif laue in ['4/m','4/mmm']:
2307	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2308	elif laue in ['mmm']:
2309	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2310	elif laue in ['2/m']:
2311	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2312	if uniq == 'a':
2313	dDijDict[phfx+'D23'] = k*l
2314	elif uniq == 'b':
2315	dDijDict[phfx+'D13'] = h*l
2316	elif uniq == 'c':
2317	dDijDict[phfx+'D12'] = h*k
2318	names.append()
2319	else:
2320	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2321	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2322	for item in dDijDict:
2323	dDijDict[item] = refl[4]2tand(refl[5]/2.0)
2324	return dDijDict
2325
2326	def GetFprime(controlDict,Histograms):
2327	FFtables = controlDict['FFtables']
2328	if not FFtables:
2329	return
2330	histoList = Histograms.keys()
2331	histoList.sort()
2332	for histogram in histoList:
2333	if 'PWDR' in histogram[:4]:
2334	Histogram = Histograms[histogram]
2335	hId = Histogram['hId']
2336	hfx = ':%d:'%(hId)
2337	keV = controlDict[hfx+'keV']
2338	for El in FFtables:
2339	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2340	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2341	FFtables[El][hfx+'FP'] = FP
2342	FFtables[El][hfx+'FPP'] = FPP
2343
2344	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2345	histoList = Histograms.keys()
2346	histoList.sort()
2347	for histogram in histoList:
2348	if 'PWDR' in histogram[:4]:
2349	Histogram = Histograms[histogram]
2350	hId = Histogram['hId']
2351	hfx = ':%d:'%(hId)
2352	Limits = calcControls[hfx+'Limits']
2353	shl = max(parmDict[hfx+'SH/L'],0.002)
2354	Ka2 = False
2355	kRatio = 0.0
2356	if hfx+'Lam1' in parmDict.keys():
2357	Ka2 = True
2358	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2359	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2360	x,y,w,yc,yb,yd = Histogram['Data']
2361	xB = np.searchsorted(x,Limits[0])
2362	xF = np.searchsorted(x,Limits[1])
2363	ymb = np.array(y-yb)
2364	ycmb = np.array(yc-yb)
2365	ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)
2366	refLists = Histogram['Reflection Lists']
2367	for phase in refLists:
2368	Phase = Phases[phase]
2369	pId = Phase['pId']
2370	phfx = '%d:%d:'%(pId,hId)
2371	refList = refLists[phase]
2372	sumFo = 0.0
2373	sumdF = 0.0
2374	sumFosq = 0.0
2375	sumdFsq = 0.0
2376	for refl in refList:
2377	if 'C' in calcControls[hfx+'histType']:
2378	yp = np.zeros_like(yb)
2379	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2380	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
2381	iFin = min(np.searchsorted(x,refl[5]+fmax),xF)
2382	iFin2 = iFin
2383	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2384	if Ka2:
2385	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2386	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2387	iBeg2 = max(xB,np.searchsorted(x[xB:xF],pos2-fmin))
2388	iFin2 = min(np.searchsorted(x[xB:xF],pos2+fmax),xF)
2389	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
2390	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
2391	elif 'T' in calcControls[hfx+'histType']:
2392	print 'TOF Undefined at present'
2393	raise Exception #no TOF yet
2394	Fo = np.sqrt(np.abs(refl[8]))
2395	Fc = np.sqrt(np.abs(refl[9]))
2396	sumFo += Fo
2397	sumFosq += refl[8]**2
2398	sumdF += np.abs(Fo-Fc)
2399	sumdFsq += (refl[8]-refl[9])**2
2400	Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2401	Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2402	Histogram[phfx+'Nref'] = len(refList)
2403
2404	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2405
2406	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
2407	U = parmDict[hfx+'U']
2408	V = parmDict[hfx+'V']
2409	W = parmDict[hfx+'W']
2410	X = parmDict[hfx+'X']
2411	Y = parmDict[hfx+'Y']
2412	tanPos = tand(refl[5]/2.0)
2413	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
2414	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
2415	sig = max(0.001,sig)
2416	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
2417	gam = max(0.001,gam)
2418	return sig,gam
2419
2420	hId = Histogram['hId']
2421	hfx = ':%d:'%(hId)
2422	bakType = calcControls[hfx+'bakType']
2423	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2424	yc = np.zeros_like(yb)
2425
2426	if 'C' in calcControls[hfx+'histType']:
2427	shl = max(parmDict[hfx+'SH/L'],0.002)
2428	Ka2 = False
2429	if hfx+'Lam1' in parmDict.keys():
2430	wave = parmDict[hfx+'Lam1']
2431	Ka2 = True
2432	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2433	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2434	else:
2435	wave = parmDict[hfx+'Lam']
2436	else:
2437	print 'TOF Undefined at present'
2438	raise ValueError
2439	for phase in Histogram['Reflection Lists']:
2440	refList = Histogram['Reflection Lists'][phase]
2441	Phase = Phases[phase]
2442	pId = Phase['pId']
2443	pfx = '%d::'%(pId)
2444	phfx = '%d:%d:'%(pId,hId)
2445	hfx = ':%d:'%(hId)
2446	SGData = Phase['General']['SGData']
2447	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2448	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2449	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2450	Vst = np.sqrt(nl.det(G)) #V*
2451	if not Phase['General'].get('doPawley'):
2452	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2453	for refl in refList:
2454	if 'C' in calcControls[hfx+'histType']:
2455	h,k,l = refl[:3]
2456	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
2457	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
2458	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
2459	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
2460	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
2461	refl[13] = VstLorenz
2462	if Phase['General'].get('doPawley'):
2463	try:
2464	refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2465	except KeyError:
2466	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2467	continue
2468	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2469	iBeg = np.searchsorted(x,refl[5]-fmin)
2470	iFin = np.searchsorted(x,refl[5]+fmax)
2471	if not iBeg+iFin: #peak below low limit - skip peak
2472	continue
2473	elif not iBeg-iFin: #peak above high limit - done
2474	break
2475	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2476	if Ka2:
2477	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2478	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2479	iBeg = np.searchsorted(x,pos2-fmin)
2480	iFin = np.searchsorted(x,pos2+fmax)
2481	if not iBeg+iFin: #peak below low limit - skip peak
2482	continue
2483	elif not iBeg-iFin: #peak above high limit - done
2484	return yc,yb
2485	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here
2486	elif 'T' in calcControls[hfx+'histType']:
2487	print 'TOF Undefined at present'
2488	raise Exception #no TOF yet
2489	return yc,yb
2490
2491	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2492
2493	def cellVaryDerv(pfx,SGData,dpdA):
2494	if SGData['SGLaue'] in ['-1',]:
2495	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2496	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2497	elif SGData['SGLaue'] in ['2/m',]:
2498	if SGData['SGUniq'] == 'a':
2499	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2500	elif SGData['SGUniq'] == 'b':
2501	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2502	else:
2503	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2504	elif SGData['SGLaue'] in ['mmm',]:
2505	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2506	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2507	# return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2508	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2509	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2510	# return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2511	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2512	elif SGData['SGLaue'] in ['3R', '3mR']:
2513	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2514	elif SGData['SGLaue'] in ['m3m','m3']:
2515	# return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2516	return [[pfx+'A0',dpdA[0]]]
2517	# create a list of dependent variables and set up a dictionary to hold their derivatives
2518	dependentVars = G2mv.GetDependentVars()
2519	depDerivDict = {}
2520	for j in dependentVars:
2521	depDerivDict[j] = np.zeros(shape=(len(x)))
2522	#print 'dependent vars',dependentVars
2523	lenX = len(x)
2524	hId = Histogram['hId']
2525	hfx = ':%d:'%(hId)
2526	bakType = calcControls[hfx+'bakType']
2527	dMdv = np.zeros(shape=(len(varylist),len(x)))
2528	dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2529	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2530	bBpos =varylist.index(hfx+'Back:0')
2531	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2532	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2533	for name in varylist:
2534	if 'Debye' in name:
2535	id = int(name.split(':')[-1])
2536	parm = name[:int(name.rindex(':'))]
2537	ip = names.index(parm)
2538	dMdv[varylist.index(name)] = dMddb[3*id+ip]
2539	if 'C' in calcControls[hfx+'histType']:
2540	dx = x[1]-x[0]
2541	shl = max(parmDict[hfx+'SH/L'],0.002)
2542	Ka2 = False
2543	if hfx+'Lam1' in parmDict.keys():
2544	wave = parmDict[hfx+'Lam1']
2545	Ka2 = True
2546	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2547	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2548	else:
2549	wave = parmDict[hfx+'Lam']
2550	else:
2551	print 'TOF Undefined at present'
2552	raise ValueError
2553	for phase in Histogram['Reflection Lists']:
2554	refList = Histogram['Reflection Lists'][phase]
2555	Phase = Phases[phase]
2556	SGData = Phase['General']['SGData']
2557	pId = Phase['pId']
2558	pfx = '%d::'%(pId)
2559	phfx = '%d:%d:'%(pId,hId)
2560	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2561	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2562	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2563	if not Phase['General'].get('doPawley'):
2564	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2565	for iref,refl in enumerate(refList):
2566	if 'C' in calcControls[hfx+'histType']: #CW powder
2567	h,k,l = refl[:3]
2568	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2569	if Phase['General'].get('doPawley'):
2570	try:
2571	refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2572	except KeyError:
2573	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2574	continue
2575	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2576	iBeg = np.searchsorted(x,refl[5]-fmin)
2577	iFin = np.searchsorted(x,refl[5]+fmax)
2578	if not iBeg+iFin: #peak below low limit - skip peak
2579	continue
2580	elif not iBeg-iFin: #peak above high limit - done
2581	break
2582	pos = refl[5]
2583	tanth = tand(pos/2.0)
2584	costh = cosd(pos/2.0)
2585	lenBF = iFin-iBeg
2586	dMdpk = np.zeros(shape=(6,lenBF))
2587	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2588	for i in range(1,5):
2589	dMdpk[i] += 100.dxrefl[13]refl[9]dMdipk[i]
2590	dMdpk[0] += 100.dxrefl[13]refl[9]dMdipk[0]
2591	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2592	if Ka2:
2593	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2594	kdelt = int((pos2-refl[5])/dx)
2595	iBeg2 = min(lenX,iBeg+kdelt)
2596	iFin2 = min(lenX,iFin+kdelt)
2597	if iBeg2-iFin2:
2598	lenBF2 = iFin2-iBeg2
2599	dMdpk2 = np.zeros(shape=(6,lenBF2))
2600	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2601	for i in range(1,5):
2602	dMdpk2[i] = 100.dxrefl[13]refl[9]kRatio*dMdipk2[i]
2603	dMdpk2[0] = 100.dxrefl[13]refl[9]kRatio*dMdipk2[0]
2604	dMdpk2[5] = 100.dxrefl[13]*dMdipk2[0]
2605	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2606	if Phase['General'].get('doPawley'):
2607	try:
2608	idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2609	dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
2610	if Ka2:
2611	dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
2612	# Assuming Pawley variables not in constraints
2613	except ValueError:
2614	pass
2615	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2616	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2617	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2618	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2619	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2620	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2621	hfx+'DisplaceY':[dpdY,'pos'],}
2622	for name in names:
2623	item = names[name]
2624	if name in varylist:
2625	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2626	if Ka2:
2627	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2628	elif name in dependentVars:
2629	if Ka2:
2630	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2631	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2632	for iPO in dIdPO:
2633	if iPO in varylist:
2634	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2635	if Ka2:
2636	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2637	elif iPO in dependentVars:
2638	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2639	if Ka2:
2640	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2641	for i,name in enumerate(['omega','chi','phi']):
2642	aname = pfx+'SH '+name
2643	if aname in varylist:
2644	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2645	if Ka2:
2646	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2647	elif aname in dependentVars:
2648	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2649	if Ka2:
2650	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2651	for iSH in dFdODF:
2652	if iSH in varylist:
2653	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2654	if Ka2:
2655	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2656	elif iSH in dependentVars:
2657	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2658	if Ka2:
2659	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2660	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2661	for name,dpdA in cellDervNames:
2662	if name in varylist:
2663	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2664	if Ka2:
2665	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2666	elif name in dependentVars:
2667	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2668	if Ka2:
2669	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2670	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2671	for name in dDijDict:
2672	if name in varylist:
2673	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2674	if Ka2:
2675	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2676	elif name in dependentVars:
2677	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2678	if Ka2:
2679	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2680	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2681	for name in gamDict:
2682	if name in varylist:
2683	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2684	if Ka2:
2685	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2686	elif name in dependentVars:
2687	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2688	if Ka2:
2689	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2690	for name in sigDict:
2691	if name in varylist:
2692	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2693	if Ka2:
2694	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2695	elif name in dependentVars:
2696	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2697	if Ka2:
2698	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2699
2700	elif 'T' in calcControls[hfx+'histType']:
2701	print 'TOF Undefined at present'
2702	raise Exception #no TOF yet
2703	#do atom derivatives - for F,X & U so far
2704	corr = dervDict['int']/refl[9]
2705	if Ka2:
2706	corr2 = dervDict2['int']/refl[9]
2707	for name in varylist+dependentVars:
2708	try:
2709	aname = name.split(pfx)[1][:2]
2710	if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2711	except IndexError:
2712	continue
2713	if name in varylist:
2714	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2715	if Ka2:
2716	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2717	elif name in dependentVars:
2718	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2719	if Ka2:
2720	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2721	# now process derivatives in constraints
2722	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2723	return dMdv
2724
2725	def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2726	parmdict.update(zip(varylist,values))
2727	G2mv.Dict2Map(parmdict,varylist)
2728	Histograms,Phases = HistoPhases
2729	nvar = len(varylist)
2730	dMdv = np.empty(0)
2731	histoList = Histograms.keys()
2732	histoList.sort()
2733	for histogram in histoList:
2734	if 'PWDR' in histogram[:4]:
2735	Histogram = Histograms[histogram]
2736	hId = Histogram['hId']
2737	hfx = ':%d:'%(hId)
2738	Limits = calcControls[hfx+'Limits']
2739	x,y,w,yc,yb,yd = Histogram['Data']
2740	xB = np.searchsorted(x,Limits[0])
2741	xF = np.searchsorted(x,Limits[1])
2742	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2743	varylist,Histogram,Phases,calcControls,pawleyLookup)
2744	if len(dMdv):
2745	dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2746	else:
2747	dMdv = dMdvh
2748	return dMdv
2749
2750	def ComputePowderHessian(args):
2751	Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2752	hId = Histogram['hId']
2753	hfx = ':%d:'%(hId)
2754	Limits = calcControls[hfx+'Limits']
2755	x,y,w,yc,yb,yd = Histogram['Data']
2756	dy = y-yc
2757	xB = np.searchsorted(x,Limits[0])
2758	xF = np.searchsorted(x,Limits[1])
2759	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
2760	parmdict,x[xB:xF],
2761	varylist,Histogram,Phases,calcControls,pawleyLookup)
2762	Vec = np.sum(dMdvhnp.sqrt(w[xB:xF])dy[xB:xF],axis=1)
2763	Hess = np.inner(dMdvh,dMdvh)
2764	return Vec,Hess
2765
2766	def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2767	parmdict.update(zip(varylist,values))
2768	G2mv.Dict2Map(parmdict,varylist)
2769	Histograms,Phases = HistoPhases
2770	nvar = len(varylist)
2771	Hess = np.empty(0)
2772	histoList = Histograms.keys()
2773	histoList.sort()
2774	for histogram in histoList:
2775	if 'PWDR' in histogram[:4]:
2776	Histogram = Histograms[histogram]
2777	hId = Histogram['hId']
2778	hfx = ':%d:'%(hId)
2779	Limits = calcControls[hfx+'Limits']
2780	x,y,w,yc,yb,yd = Histogram['Data']
2781	dy = y-yc
2782	xB = np.searchsorted(x,Limits[0])
2783	xF = np.searchsorted(x,Limits[1])
2784	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2785	varylist,Histogram,Phases,calcControls,pawleyLookup)
2786	if dlg:
2787	dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2788	if len(Hess):
2789	Vec += np.sum(dMdvhnp.sqrt(w[xB:xF])dy[xB:xF],axis=1)
2790	Hess += np.inner(dMdvh,dMdvh)
2791	else:
2792	Vec = np.sum(dMdvhnp.sqrt(w[xB:xF])dy[xB:xF],axis=1)
2793	Hess = np.inner(dMdvh,dMdvh)
2794	return Vec,Hess
2795
2796	def ComputePowderProfile(args):
2797	Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2798	hId = Histogram['hId']
2799	hfx = ':%d:'%(hId)
2800	x,y,w,yc,yb,yd = Histogram['Data']
2801	Limits = calcControls[hfx+'Limits']
2802	xB = np.searchsorted(x,Limits[0])
2803	xF = np.searchsorted(x,Limits[1])
2804	yc,yb = getPowderProfile(parmdict,x[xB:xF],varylist,Histogram,Phases,
2805	calcControls,pawleyLookup)
2806	return xB,xF,yc,yb,Histogram['Reflection Lists']
2807
2808	def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2809	parmdict.update(zip(varylist,values))
2810	Values2Dict(parmdict, varylist, values)
2811	G2mv.Dict2Map(parmdict,varylist)
2812	Histograms,Phases = HistoPhases
2813	M = np.empty(0)
2814	sumwYo = 0
2815	Nobs = 0
2816	histoList = Histograms.keys()
2817	histoList.sort()
2818	for histogram in histoList:
2819	if 'PWDR' in histogram[:4]:
2820	Histogram = Histograms[histogram]
2821	hId = Histogram['hId']
2822	hfx = ':%d:'%(hId)
2823	Limits = calcControls[hfx+'Limits']
2824	x,y,w,yc,yb,yd = Histogram['Data']
2825	yc *= 0.0 #zero full calcd profiles
2826	yb *= 0.0
2827	yd *= 0.0
2828	xB = np.searchsorted(x,Limits[0])
2829	xF = np.searchsorted(x,Limits[1])
2830	Histogram['Nobs'] = xF-xB
2831	Nobs += Histogram['Nobs']
2832	Histogram['sumwYo'] = np.sum(w[xB:xF]y[xB:xF]*2)
2833	sumwYo += Histogram['sumwYo']
2834	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2835	varylist,Histogram,Phases,calcControls,pawleyLookup)
2836	yc[xB:xF] += yb[xB:xF]
2837	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2838	Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2839	wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2840	Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy*2)/Histogram['sumwYo'])100.)
2841	if dlg:
2842	dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2843	M = np.concatenate((M,wdy))
2844	Histograms['sumwYo'] = sumwYo
2845	Histograms['Nobs'] = Nobs
2846	Rwp = min(100.,np.sqrt(np.sum(M*2)/sumwYo)100.)
2847	if dlg:
2848	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2849	if not GoOn:
2850	parmDict['saved values'] = values
2851	raise Exception #Abort!!
2852	return M
2853
2854	def Refine(GPXfile,dlg):
2855	import pytexture as ptx
2856	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
2857
2858	ShowBanner()
2859	varyList = []
2860	parmDict = {}
2861	G2mv.InitVars()
2862	Controls = GetControls(GPXfile)
2863	ShowControls(Controls)
2864	calcControls = {}
2865	calcControls.update(Controls)
2866	constrDict,fixedList = GetConstraints(GPXfile)
2867	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2868	if not Phases:
2869	print ' * ERROR - you have no phases! *'
2870	print ' * Refine aborted *'
2871	raise Exception
2872	if not Histograms:
2873	print ' * ERROR - you have no data to refine with! *'
2874	print ' * Refine aborted *'
2875	raise Exception
2876	Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2877	calcControls['Natoms'] = Natoms
2878	calcControls['FFtables'] = FFtables
2879	calcControls['BLtables'] = BLtables
2880	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2881	calcControls.update(controlDict)
2882	histVary,histDict,controlDict = GetHistogramData(Histograms)
2883	calcControls.update(controlDict)
2884	varyList = phaseVary+hapVary+histVary
2885	parmDict.update(phaseDict)
2886	parmDict.update(hapDict)
2887	parmDict.update(histDict)
2888	GetFprime(calcControls,Histograms)
2889	# do constraint processing
2890	try:
2891	groups,parmlist = G2mv.GroupConstraints(constrDict)
2892	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
2893	except:
2894	print ' * ERROR - your constraints are internally inconsistent *'
2895	# traceback for debug
2896	#print 'varyList',varyList
2897	#print 'constrDict',constrDict
2898	#print 'fixedList',fixedList
2899	#import traceback
2900	#print traceback.format_exc()
2901	raise Exception(' * Refine aborted *')
2902	# # check to see which generated parameters are fully varied
2903	# msg = G2mv.SetVaryFlags(varyList)
2904	# if msg:
2905	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
2906	# print msg
2907	# raise Exception(' * Refine aborted *')
2908	#print G2mv.VarRemapShow(varyList)
2909	G2mv.Map2Dict(parmDict,varyList)
2910	Rvals = {}
2911	while True:
2912	begin = time.time()
2913	values = np.array(Dict2Values(parmDict, varyList))
2914	Ftol = Controls['min dM/M']
2915	Factor = Controls['shift factor']
2916	maxCyc = Controls['max cyc']
2917	if 'Jacobian' in Controls['deriv type']:
2918	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2919	ftol=Ftol,col_deriv=True,factor=Factor,
2920	args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2921	ncyc = int(result[2]['nfev']/2)
2922	elif 'Hessian' in Controls['deriv type']:
2923	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
2924	args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2925	ncyc = result[2]['num cyc']+1
2926	Rvals['lamMax'] = result[2]['lamMax']
2927	else: #'numeric'
2928	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2929	args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2930	ncyc = int(result[2]['nfev']/len(varyList))
2931	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2932	# for item in table: print item,table[item] #useful debug - are things shifting?
2933	runtime = time.time()-begin
2934	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
2935	Values2Dict(parmDict, varyList, result[0])
2936	G2mv.Dict2Map(parmDict,varyList)
2937
2938	Rvals['Nobs'] = Histograms['Nobs']
2939	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
2940	Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
2941	print '\n Refinement results:'
2942	print 135*'-'
2943	print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2944	print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
2945	print ' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
2946	try:
2947	covMatrix = result[1]*Rvals['GOF']
2948	sig = np.sqrt(np.diag(covMatrix))
2949	if np.any(np.isnan(sig)):
2950	print '* Least squares aborted - some invalid esds possible *'
2951	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2952	# for item in table: print item,table[item] #useful debug - are things shifting?
2953	break #refinement succeeded - finish up!
2954	except TypeError: #result[1] is None on singular matrix
2955	print '** Refinement failed - singular matrix **'
2956	if 'Hessian' in Controls['deriv type']:
2957	num = len(varyList)-1
2958	for i,val in enumerate(np.flipud(result[2]['psing'])):
2959	if val:
2960	print 'Removing parameter: ',varyList[num-i]
2961	del(varyList[num-i])
2962	else:
2963	Ipvt = result[2]['ipvt']
2964	for i,ipvt in enumerate(Ipvt):
2965	if not np.sum(result[2]['fjac'],axis=1)[i]:
2966	print 'Removing parameter: ',varyList[ipvt-1]
2967	del(varyList[ipvt-1])
2968	break
2969
2970	# print 'dependentParmList: ',G2mv.dependentParmList
2971	# print 'arrayList: ',G2mv.arrayList
2972	# print 'invarrayList: ',G2mv.invarrayList
2973	# print 'indParmList: ',G2mv.indParmList
2974	# print 'fixedDict: ',G2mv.fixedDict
2975	# print 'test1'
2976	GetFobsSq(Histograms,Phases,parmDict,calcControls)
2977	# print 'test2'
2978	sigDict = dict(zip(varyList,sig))
2979	newCellDict = GetNewCellParms(parmDict,varyList)
2980	newAtomDict = ApplyXYZshifts(parmDict,varyList)
2981	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
2982	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2983	# add the uncertainties into the esd dictionary (sigDict)
2984	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
2985	SetPhaseData(parmDict,sigDict,Phases,covData)
2986	SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2987	SetHistogramData(parmDict,sigDict,Histograms)
2988	G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
2989	SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2990
2991	#for testing purposes!!!
2992	# import cPickle
2993	# fl = open('structTestdata.dat','wb')
2994	# cPickle.dump(parmDict,fl,1)
2995	# cPickle.dump(varyList,fl,1)
2996	# for histogram in Histograms:
2997	# if 'PWDR' in histogram[:4]:
2998	# Histogram = Histograms[histogram]
2999	# cPickle.dump(Histogram,fl,1)
3000	# cPickle.dump(Phases,fl,1)
3001	# cPickle.dump(calcControls,fl,1)
3002	# cPickle.dump(pawleyLookup,fl,1)
3003	# fl.close()
3004
3005	if dlg:
3006	return Rvals['Rwp']
3007
3008	def SeqRefine(GPXfile,dlg):
3009	import pytexture as ptx
3010	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
3011
3012	ShowBanner()
3013	print ' Sequential Refinement'
3014	G2mv.InitVars()
3015	Controls = GetControls(GPXfile)
3016	ShowControls(Controls)
3017	constrDict,fixedList = GetConstraints(GPXfile)
3018	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3019	if not Phases:
3020	print ' * ERROR - you have no histograms to refine! *'
3021	print ' * Refine aborted *'
3022	raise Exception
3023	if not Histograms:
3024	print ' * ERROR - you have no data to refine with! *'
3025	print ' * Refine aborted *'
3026	raise Exception
3027	Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
3028	if 'Seq Data' in Controls:
3029	histNames = Controls['Seq Data']
3030	else:
3031	histNames = GetHistogramNames(GPXfile,['PWDR',])
3032	if 'Reverse Seq' in Controls:
3033	if Controls['Reverse Seq']:
3034	histNames.reverse()
3035	SeqResult = {'histNames':histNames}
3036	makeBack = True
3037	for ihst,histogram in enumerate(histNames):
3038	ifPrint = False
3039	if dlg:
3040	dlg.SetTitle('Residual for histogram '+str(ihst))
3041	calcControls = {}
3042	calcControls['Natoms'] = Natoms
3043	calcControls['FFtables'] = FFtables
3044	calcControls['BLtables'] = BLtables
3045	varyList = []
3046	parmDict = {}
3047	Histo = {histogram:Histograms[histogram],}
3048	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3049	calcControls.update(controlDict)
3050	histVary,histDict,controlDict = GetHistogramData(Histo,False)
3051	calcControls.update(controlDict)
3052	varyList = phaseVary+hapVary+histVary
3053	if not ihst:
3054	saveVaryList = varyList[:]
3055	for i,item in enumerate(saveVaryList):
3056	items = item.split(':')
3057	if items[1]:
3058	items[1] = ''
3059	item = ':'.join(items)
3060	saveVaryList[i] = item
3061	SeqResult['varyList'] = saveVaryList
3062	else:
3063	newVaryList = varyList[:]
3064	for i,item in enumerate(newVaryList):
3065	items = item.split(':')
3066	if items[1]:
3067	items[1] = ''
3068	item = ':'.join(items)
3069	newVaryList[i] = item
3070	if newVaryList != SeqResult['varyList']:
3071	print newVaryList
3072	print SeqResult['varyList']
3073	print '**** ERROR - variable list for this histogram does not match previous'
3074	raise Exception
3075	parmDict.update(phaseDict)
3076	parmDict.update(hapDict)
3077	parmDict.update(histDict)
3078	GetFprime(calcControls,Histo)
3079	# do constraint processing
3080	try:
3081	groups,parmlist = G2mv.GroupConstraints(constrDict)
3082	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3083	except:
3084	print ' * ERROR - your constraints are internally inconsistent *'
3085	raise Exception(' * Refine aborted *')
3086	# check to see which generated parameters are fully varied
3087	# msg = G2mv.SetVaryFlags(varyList)
3088	# if msg:
3089	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
3090	# print msg
3091	# raise Exception(' * Refine aborted *')
3092	#print G2mv.VarRemapShow(varyList)
3093	G2mv.Map2Dict(parmDict,varyList)
3094	Rvals = {}
3095	while True:
3096	begin = time.time()
3097	values = np.array(Dict2Values(parmDict, varyList))
3098	Ftol = Controls['min dM/M']
3099	Factor = Controls['shift factor']
3100	maxCyc = Controls['max cyc']
3101
3102	if 'Jacobian' in Controls['deriv type']:
3103	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3104	ftol=Ftol,col_deriv=True,factor=Factor,
3105	args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3106	ncyc = int(result[2]['nfev']/2)
3107	elif 'Hessian' in Controls['deriv type']:
3108	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3109	args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3110	ncyc = result[2]['num cyc']+1
3111	else: #'numeric'
3112	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3113	args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3114	ncyc = int(result[2]['nfev']/len(varyList))
3115
3116
3117
3118	runtime = time.time()-begin
3119	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3120	Values2Dict(parmDict, varyList, result[0])
3121	G2mv.Dict2Map(parmDict,varyList)
3122
3123	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100. #to %
3124	Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3125	Rvals['Nobs'] = Histo['Nobs']
3126	print '\n Refinement results for histogram: v'+histogram
3127	print 135*'-'
3128	print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3129	print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3130	print ' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3131	try:
3132	covMatrix = result[1]*Rvals['GOF']
3133	sig = np.sqrt(np.diag(covMatrix))
3134	if np.any(np.isnan(sig)):
3135	print '* Least squares aborted - some invalid esds possible *'
3136	ifPrint = True
3137	break #refinement succeeded - finish up!
3138	except TypeError: #result[1] is None on singular matrix
3139	print '** Refinement failed - singular matrix **'
3140	if 'Hessian' in Controls['deriv type']:
3141	num = len(varyList)-1
3142	for i,val in enumerate(np.flipud(result[2]['psing'])):
3143	if val:
3144	print 'Removing parameter: ',varyList[num-i]
3145	del(varyList[num-i])
3146	else:
3147	Ipvt = result[2]['ipvt']
3148	for i,ipvt in enumerate(Ipvt):
3149	if not np.sum(result[2]['fjac'],axis=1)[i]:
3150	print 'Removing parameter: ',varyList[ipvt-1]
3151	del(varyList[ipvt-1])
3152	break
3153
3154	GetFobsSq(Histo,Phases,parmDict,calcControls)
3155	sigDict = dict(zip(varyList,sig))
3156	newCellDict = GetNewCellParms(parmDict,varyList)
3157	newAtomDict = ApplyXYZshifts(parmDict,varyList)
3158	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3159	'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3160	SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
3161	SetHistogramData(parmDict,sigDict,Histo,ifPrint)
3162	SeqResult[histogram] = covData
3163	SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3164	makeBack = False
3165	SetSeqResult(GPXfile,Histograms,SeqResult)
3166
3167	def DistAngle(DisAglCtls,DisAglData):
3168	import numpy.ma as ma
3169
3170	def ShowBanner(name):
3171	print 80''
3172	print ' Interatomic Distances and Angles for phase '+name
3173	print 80'','\n'
3174
3175	ShowBanner(DisAglCtls['Name'])
3176	SGData = DisAglData['SGData']
3177	SGtext = G2spc.SGPrint(SGData)
3178	for line in SGtext: print line
3179	Cell = DisAglData['Cell']
3180
3181	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3182	covData = {}
3183	if 'covData' in DisAglData:
3184	covData = DisAglData['covData']
3185	covMatrix = covData['covMatrix']
3186	varyList = covData['varyList']
3187	pfx = str(DisAglData['pId'])+'::'
3188	A = G2lat.cell2A(Cell[:6])
3189	cellSig = getCellEsd(pfx,SGData,A,covData)
3190	names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3191	valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3192	line = '\n Unit cell:'
3193	for name,vals in zip(names,valEsd):
3194	line += name+vals
3195	print line
3196	else:
3197	print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3198	' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3199	'%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3200	Factor = DisAglCtls['Factors']
3201	Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3202	Units = np.array([ #is there a nicer way to make this?
3203	[-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1],
3204	[0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1],
3205	[1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]])
3206	origAtoms = DisAglData['OrigAtoms']
3207	targAtoms = DisAglData['TargAtoms']
3208	for Oatom in origAtoms:
3209	OxyzNames = ''
3210	IndBlist = []
3211	Dist = []
3212	Vect = []
3213	VectA = []
3214	angles = []
3215	for Tatom in targAtoms:
3216	Xvcov = []
3217	TxyzNames = ''
3218	if 'covData' in DisAglData:
3219	OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3220	TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3221	Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3222	result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3223	BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3224	AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3225	for Txyz,Top,Tunit in result:
3226	Dx = (Txyz-np.array(Oatom[3:6]))+Units
3227	dx = np.inner(Amat,Dx)
3228	dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3229	IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3230	if np.any(IndB):
3231	for indb in IndB:
3232	for i in range(len(indb)):
3233	if str(dx.T[indb][i]) not in IndBlist:
3234	IndBlist.append(str(dx.T[indb][i]))
3235	unit = Units[indb][i]
3236	tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3237	pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3238	sig = 0.0
3239	if len(Xvcov):
3240	sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3241	Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3242	if (Dist[-1][-1]-AsumR) <= 0.:
3243	Vect.append(dx.T[indb][i]/Dist[-1][-2])
3244	VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3245	else:
3246	Vect.append([0.,0.,0.])
3247	VectA.append([])
3248	Vect = np.array(Vect)
3249	angles = np.zeros((len(Vect),len(Vect)))
3250	angsig = np.zeros((len(Vect),len(Vect)))
3251	for i,veca in enumerate(Vect):
3252	if np.any(veca):
3253	for j,vecb in enumerate(Vect):
3254	if np.any(vecb):
3255	angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3256	line = ''
3257	for i,x in enumerate(Oatom[3:6]):
3258	if len(Xvcov):
3259	line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3260	else:
3261	line += '%12.5f'%(x)
3262	print '\n Distances & angles for ',Oatom[1],' at ',line
3263	print 80''
3264	line = ''
3265	for dist in Dist[:-1]:
3266	line += '%12s'%(dist[1].center(12))
3267	print ' To cell +(sym. op.) dist. ',line
3268	for i,dist in enumerate(Dist):
3269	line = ''
3270	for j,angle in enumerate(angles[i][0:i]):
3271	sig = angsig[i][j]
3272	if angle:
3273	if sig:
3274	line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3275	else:
3276	val = '%.3f'%(angle)
3277	line += '%12s'%(val.center(12))
3278	else:
3279	line += 12*' '
3280	if dist[5]: #sig exists!
3281	val = G2mth.ValEsd(dist[4],dist[5])
3282	else:
3283	val = '%8.4f'%(dist[4])
3284	print ' %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3285
3286	def DisAglTor(DATData):
3287	SGData = DATData['SGData']
3288	Cell = DATData['Cell']
3289
3290	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3291	covData = {}
3292	pfx = ''
3293	if 'covData' in DATData:
3294	covData = DATData['covData']
3295	covMatrix = covData['covMatrix']
3296	varyList = covData['varyList']
3297	pfx = str(DATData['pId'])+'::'
3298	Datoms = []
3299	Oatoms = []
3300	for i,atom in enumerate(DATData['Datoms']):
3301	symop = atom[-1].split('+')
3302	if len(symop) == 1:
3303	symop.append('0,0,0')
3304	symop[0] = int(symop[0])
3305	symop[1] = eval(symop[1])
3306	atom.append(symop)
3307	Datoms.append(atom)
3308	oatom = DATData['Oatoms'][i]
3309	names = ['','','']
3310	if pfx:
3311	names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3312	oatom += [names,]
3313	Oatoms.append(oatom)
3314	atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3315	if DATData['Natoms'] == 4: #torsion
3316	Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3317	print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3318	print ' NB: Atom sequence determined by selection order'
3319	return # done with torsion
3320	elif DATData['Natoms'] == 3: #angle
3321	Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3322	print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3323	print ' NB: Atom sequence determined by selection order'
3324	else: #2 atoms - distance
3325	Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3326	print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3327
3328	def BestPlane(PlaneData):
3329
3330	def ShowBanner(name):
3331	print 80''
3332	print ' Best plane result for phase '+name
3333	print 80'','\n'
3334
3335	ShowBanner(PlaneData['Name'])
3336
3337	Cell = PlaneData['Cell']
3338	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3339	Atoms = PlaneData['Atoms']
3340	sumXYZ = np.zeros(3)
3341	XYZ = []
3342	Natoms = len(Atoms)
3343	for atom in Atoms:
3344	xyz = np.array(atom[3:6])
3345	XYZ.append(xyz)
3346	sumXYZ += xyz
3347	sumXYZ /= Natoms
3348	XYZ = np.array(XYZ)-sumXYZ
3349	XYZ = np.inner(Amat,XYZ).T
3350	Zmat = np.zeros((3,3))
3351	for i,xyz in enumerate(XYZ):
3352	Zmat += np.outer(xyz.T,xyz)
3353	print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3354	Evec,Emat = nl.eig(Zmat)
3355	Evec = np.sqrt(Evec)/(Natoms-3)
3356	Order = np.argsort(Evec)
3357	XYZ = np.inner(XYZ,Emat.T).T
3358	XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3359	print ' Atoms in Cartesian best plane coordinates:'
3360	print ' Name X Y Z'
3361	for i,xyz in enumerate(XYZ):
3362	print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3363	print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])
3364
3365
3366	def main():
3367	arg = sys.argv
3368	if len(arg) > 1:
3369	GPXfile = arg[1]
3370	if not ospath.exists(GPXfile):
3371	print 'ERROR - ',GPXfile," doesn't exist!"
3372	exit()
3373	GPXpath = ospath.dirname(arg[1])
3374	Refine(GPXfile,None)
3375	else:
3376	print 'ERROR - missing filename'
3377	exit()
3378	print "Done"
3379
3380	if __name__ == '__main__':
3381	main()

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

Download in other formats: