Context Navigation

source: trunk/GSASIIstruct.py @ 425

Last change on this file since 425 was 425, checked in by vondreele, 12 years ago
get As from Dijs in seqdata plotting
Property svn:keywords set to `Date Author Revision URL Id`
File size: 119.9 KB

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

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

Download in other formats: