Context Navigation

source: trunk/GSASIIstruct.py @ 693

Last change on this file since 693 was 693, checked in by vondreele, 13 years ago
LS refinement results now printed to a separate file "name.lst"; it is overwritten each time refine or sequential refine is run. Console only contains error messages.
Property svn:keywords set to `Date Author Revision URL Id`
File size: 166.9 KB

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

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

Download in other formats: