Context Navigation

source: branch/2frame/GSASIIpwd.py @ 2928

Last change on this file since 2928 was 2928, checked in by toby, 6 years ago
fix seq. ref. window; screen positioning bug; read of multiple data files after use of defaults parms
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 111.5 KB

Line
1	#/usr/bin/env python
2	# -- coding: utf-8 --
3	'''
4	GSASII powder calculation module
5	==================================
6
7	'''
8	########### SVN repository information ###################
9	# $Date: 2017-07-12 18:08:49 +0000 (Wed, 12 Jul 2017) $
10	# $Author: toby $
11	# $Revision: 2928 $
12	# $URL: branch/2frame/GSASIIpwd.py $
13	# $Id: GSASIIpwd.py 2928 2017-07-12 18:08:49Z toby $
14	########### SVN repository information ###################
15	import sys
16	import math
17	import time
18	import os
19	import subprocess as subp
20	import copy
21
22	import numpy as np
23	import numpy.linalg as nl
24	import numpy.ma as ma
25	import random as rand
26	import numpy.fft as fft
27	import scipy.interpolate as si
28	import scipy.stats as st
29	import scipy.optimize as so
30	import scipy.special as sp
31
32	import GSASIIpath
33	GSASIIpath.SetVersionNumber("$Revision: 2928 $")
34	import GSASIIlattice as G2lat
35	import GSASIIspc as G2spc
36	import GSASIIElem as G2elem
37	import GSASIImath as G2mth
38	import pypowder as pyd
39	try:
40	import pydiffax as pyx
41	except ImportError:
42	print 'pydiffax is not available for this platform - under develpment'
43
44
45	# trig functions in degrees
46	tand = lambda x: math.tan(x*math.pi/180.)
47	atand = lambda x: 180.*math.atan(x)/math.pi
48	atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
49	cosd = lambda x: math.cos(x*math.pi/180.)
50	acosd = lambda x: 180.*math.acos(x)/math.pi
51	rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
52	#numpy versions
53	npsind = lambda x: np.sin(x*np.pi/180.)
54	npasind = lambda x: 180.*np.arcsin(x)/math.pi
55	npcosd = lambda x: np.cos(x*math.pi/180.)
56	npacosd = lambda x: 180.*np.arccos(x)/math.pi
57	nptand = lambda x: np.tan(x*math.pi/180.)
58	npatand = lambda x: 180.*np.arctan(x)/np.pi
59	npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
60	npT2stl = lambda tth, wave: 2.0npsind(tth/2.0)/wave #=d
61	npT2q = lambda tth,wave: 2.0np.pinpT2stl(tth,wave) #=2pid
62	ateln2 = 8.0*math.log(2.0)
63	sateln2 = np.sqrt(ateln2)
64	nxs = np.newaxis
65
66	################################################################################
67	#### Powder utilities
68	################################################################################
69
70	def PhaseWtSum(G2frame,histo):
71	'''
72	Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
73
74	:param G2frame: GSASII main frame structure
75	:param str histo: histogram name
76	:returns: sum(scale*mass) for phases in histo
77	'''
78	Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
79	wtSum = 0.0
80	for phase in Phases:
81	if Phases[phase]['General']['Type'] != 'magnetic':
82	if histo in Phases[phase]['Histograms']:
83	mass = Phases[phase]['General']['Mass']
84	phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
85	wtSum += mass*phFr
86	return wtSum
87
88	################################################################################
89	#GSASII pdf calculation routines
90	################################################################################
91
92	def Transmission(Geometry,Abs,Diam):
93	'''
94	Calculate sample transmission
95
96	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
97	:param float Abs: absorption coeff in cm-1
98	:param float Diam: sample thickness/diameter in mm
99	'''
100	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
101	MuR = Abs*Diam/20.0
102	if MuR <= 3.0:
103	T0 = 16/(3.*math.pi)
104	T1 = -0.045780
105	T2 = -0.02489
106	T3 = 0.003045
107	T = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
108	if T < -20.:
109	return 2.06e-9
110	else:
111	return math.exp(T)
112	else:
113	T1 = 1.433902
114	T2 = 0.013869+0.337894
115	T3 = 1.933433+1.163198
116	T4 = 0.044365-0.04259
117	T = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
118	return T/100.
119	elif 'plate' in Geometry:
120	MuR = Abs*Diam/10.
121	return math.exp(-MuR)
122	elif 'Bragg' in Geometry:
123	return 0.0
124
125	def SurfaceRough(SRA,SRB,Tth):
126	''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
127	:param float SRA: Suortti surface roughness parameter
128	:param float SRB: Suortti surface roughness parameter
129	:param float Tth: 2-theta(deg) - can be numpy array
130
131	'''
132	sth = npsind(Tth/2.)
133	T1 = np.exp(-SRB/sth)
134	T2 = SRA+(1.-SRA)*np.exp(-SRB)
135	return (SRA+(1.-SRA)*T1)/T2
136
137	def SurfaceRoughDerv(SRA,SRB,Tth):
138	''' Suortti surface roughness correction derivatives
139	:param float SRA: Suortti surface roughness parameter (dimensionless)
140	:param float SRB: Suortti surface roughness parameter (dimensionless)
141	:param float Tth: 2-theta(deg) - can be numpy array
142	:return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
143	'''
144	sth = npsind(Tth/2.)
145	T1 = np.exp(-SRB/sth)
146	T2 = SRA+(1.-SRA)*np.exp(-SRB)
147	Trans = (SRA+(1.-SRA)*T1)/T2
148	dydSRA = ((1.-T1)T2-(1.-np.exp(-SRB))Trans)/T2**2
149	dydSRB = ((SRA-1.)T1T2/sth-Trans(SRA-T2))/T2*2
150	return [dydSRA,dydSRB]
151
152	def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
153	'''Calculate sample absorption
154	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
155	:param float MuR: absorption coeff * sample thickness/2 or radius
156	:param Tth: 2-theta scattering angle - can be numpy array
157	:param float Phi: flat plate tilt angle - future
158	:param float Psi: flat plate tilt axis - future
159	'''
160
161	def muRunder3(MuR,Sth2):
162	T0 = 16.0/(3.*np.pi)
163	T1 = (25.99978-0.01911Sth20.25)np.exp(-0.024551*Sth2)+ \
164	0.109561*np.sqrt(Sth2)-26.04556
165	T2 = -0.02489-0.39499Sth2+1.219077Sth2**1.5- \
166	1.31268Sth22+0.871081Sth2*2.5-0.2327Sth2**3
167	T3 = 0.003045+0.018167Sth2-0.03305Sth2**2
168	Trns = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
169	return np.exp(Trns)
170
171	def muRover3(MuR,Sth2):
172	T1 = 1.433902+11.07504Sth2-8.77629Sth2*Sth2+ \
173	10.02088Sth23-3.36778Sth2**4
174	T2 = (0.013869-0.01249Sth2)np.exp(3.27094*Sth2)+ \
175	(0.337894+13.77317Sth2)/(1.0+11.53544Sth2)**1.555039
176	T3 = 1.933433/(1.0+23.12967Sth2)*1.686715- \
177	0.13576*np.sqrt(Sth2)+1.163198
178	T4 = 0.044365-0.04259/(1.0+0.41051Sth2)*148.4202
179	Trns = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
180	return Trns/100.
181
182	Sth2 = npsind(Tth/2.0)**2
183	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
184	if 'array' in str(type(MuR)):
185	MuRSTh2 = np.concatenate((MuR,Sth2))
186	AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
187	return AbsCr
188	else:
189	if MuR <= 3.0:
190	return muRunder3(MuR,Sth2)
191	else:
192	return muRover3(MuR,Sth2)
193	elif 'Bragg' in Geometry:
194	return 1.0
195	elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
196	# and only defined for 2theta < 90
197	MuT = 2.*MuR
198	T1 = np.exp(-MuT)
199	T2 = np.exp(-MuT/npcosd(Tth))
200	Tb = MuT-MuT/npcosd(Tth)
201	return (T2-T1)/Tb
202	elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
203	MuT = 2.*MuR
204	cth = npcosd(Tth/2.0)
205	return np.exp(-MuT/cth)/cth
206
207	def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
208	'needs a doc string'
209	dA = 0.001
210	AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
211	if MuR:
212	AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
213	return (AbsP-AbsM)/(2.0*dA)
214	else:
215	return (AbsP-1.)/dA
216
217	def Polarization(Pola,Tth,Azm=0.0):
218	""" Calculate angle dependent x-ray polarization correction (not scaled correctly!)
219
220	:param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
221	:param Azm: azimuthal angle e.g. 0.0 in plane of polarization
222	:param Tth: 2-theta scattering angle - can be numpy array
223	which (if either) of these is "right"?
224	:return: (pola, dpdPola)
225	* pola = ((1-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
226	(1-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
227	* dpdPola: derivative needed for least squares
228
229	"""
230	pola = ((1.0-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
231	(1.0-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
232	dpdPola = -npsind(Tth)*2(npsind(Azm)2-npcosd(Azm)2)
233	return pola,dpdPola
234
235	def Oblique(ObCoeff,Tth):
236	'currently assumes detector is normal to beam'
237	if ObCoeff:
238	return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
239	else:
240	return 1.0
241
242	def Ruland(RulCoff,wave,Q,Compton):
243	'needs a doc string'
244	C = 2.9978e8
245	D = 1.5e-3
246	hmc = 0.024262734687
247	sinth2 = (Qwave/(4.0np.pi))**2
248	dlam = (wave*2)Compton*Q/C
249	dlam_c = 2.0hmcsinth2-Dwave*2
250	return 1.0/((1.0+dlam/RulCoff)(1.0+(np.pidlam_c/(dlam+RulCoff))**2))
251
252	def LorchWeight(Q):
253	'needs a doc string'
254	return np.sin(np.pi(Q[-1]-Q)/(2.0Q[-1]))
255
256	def GetAsfMean(ElList,Sthl2):
257	'''Calculate various scattering factor terms for PDF calcs
258
259	:param dict ElList: element dictionary contains scattering factor coefficients, etc.
260	:param np.array Sthl2: numpy array of sin theta/lambda squared values
261	:returns: mean(f^2), mean(f)^2, mean(compton)
262	'''
263	sumNoAtoms = 0.0
264	FF = np.zeros_like(Sthl2)
265	FF2 = np.zeros_like(Sthl2)
266	CF = np.zeros_like(Sthl2)
267	for El in ElList:
268	sumNoAtoms += ElList[El]['FormulaNo']
269	for El in ElList:
270	el = ElList[El]
271	ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])2+el['fpp']2
272	cf = G2elem.ComptonFac(el,Sthl2)
273	FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
274	FF2 += ff2*el['FormulaNo']/sumNoAtoms
275	CF += cf*el['FormulaNo']/sumNoAtoms
276	return FF2,FF**2,CF
277
278	def GetNumDensity(ElList,Vol):
279	'needs a doc string'
280	sumNoAtoms = 0.0
281	for El in ElList:
282	sumNoAtoms += ElList[El]['FormulaNo']
283	return sumNoAtoms/Vol
284
285	def CalcPDF(data,inst,limits,xydata):
286	'''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
287	dict xydata; results are placed in xydata.
288	Calculation parameters are found in dicts data and inst and list limits.
289	The return value is at present an empty list.
290	'''
291	auxPlot = []
292	Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
293	Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
294	#subtract backgrounds - if any & use PWDR limits
295	# GSASIIpath.IPyBreak()
296	IofQ = copy.deepcopy(xydata['Sample'])
297	IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
298	if data['Sample Bkg.']['Name']:
299	IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
300	if data['Container']['Name']:
301	xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
302	if data['Container Bkg.']['Name']:
303	xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
304	IofQ[1][1] += xycontainer[Ibeg:Ifin]
305	data['IofQmin'] = IofQ[1][1][-1]
306	IofQ[1][1] -= data.get('Flat Bkg',0.)
307	#get element data & absorption coeff.
308	ElList = data['ElList']
309	Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
310	#Apply angle dependent corrections
311	Tth = IofQ[1][0]
312	MuR = Abs*data['Diam']/20.0
313	IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
314	if 'X' in inst['Type'][0]:
315	IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
316	if data['DetType'] == 'Image plate':
317	IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
318	XY = IofQ[1]
319	#convert to Q
320	# nQpoints = len(XY[0]) #points for Q interpolation
321	nQpoints = 5000
322	if 'C' in inst['Type'][0]:
323	wave = G2mth.getWave(inst)
324	minQ = npT2q(Tth[0],wave)
325	maxQ = npT2q(Tth[-1],wave)
326	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
327	dq = Qpoints[1]-Qpoints[0]
328	XY[0] = npT2q(XY[0],wave)
329	elif 'T' in inst['Type'][0]:
330	difC = inst['difC'][1]
331	minQ = 2.np.pidifC/Tth[-1]
332	maxQ = 2.np.pidifC/Tth[0]
333	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
334	dq = Qpoints[1]-Qpoints[0]
335	XY[0] = 2.np.pidifC/XY[0]
336	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0]) #interpolate I(Q)
337	Qdata -= np.min(Qdata)*data['BackRatio']
338
339	qLimits = data['QScaleLim']
340	minQ = np.searchsorted(Qpoints,qLimits[0])
341	maxQ = np.searchsorted(Qpoints,qLimits[1])+1
342	newdata = []
343	if len(IofQ) < 3:
344	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
345	else:
346	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
347	for item in xydata['IofQ'][1]:
348	newdata.append(item[:maxQ])
349	xydata['IofQ'][1] = newdata
350
351	xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
352	FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0np.pi))*2) #these are <f^2>,<f>^2,Cf
353	Q = xydata['SofQ'][1][0]
354	# auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
355	ruland = Ruland(data['Ruland'],wave,Q,CF)
356	# auxPlot.append([Q,ruland,'Ruland'])
357	CF *= ruland
358	# auxPlot.append([Q,CF,'CF-Corr'])
359	scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
360	xydata['SofQ'][1][1] *= scale
361	xydata['SofQ'][1][1] -= CF
362	xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
363	scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
364	xydata['SofQ'][1][1] *= scale
365	xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
366	xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
367	if data['Lorch']:
368	xydata['FofQ'][1][1] *= LorchWeight(Q)
369	xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
370	nR = len(xydata['GofR'][1][1])
371	mul = int(round(2.np.pinR/(data.get('Rmax',100.)*qLimits[1])))
372	xydata['GofR'][1][0] = 2.np.pinp.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
373	xydata['GofR'][1][1] = -dqnp.imag(fft.fft(xydata['FofQ'][1][1],mulnR)[:nR])
374	if data.get('noRing',True):
375	xydata['GofR'][1][1] = np.where(xydata['GofR'][1][0]<0.5,0.,xydata['GofR'][1][1])
376	return auxPlot
377
378	def PDFPeakFit(peaks,data):
379	rs2pi = 1./np.sqrt(2*np.pi)
380
381	def MakeParms(peaks):
382	varyList = []
383	parmDict = {'slope':peaks['Background'][1][1]}
384	if peaks['Background'][2]:
385	varyList.append('slope')
386	for i,peak in enumerate(peaks['Peaks']):
387	parmDict['PDFpos;'+str(i)] = peak[0]
388	parmDict['PDFmag;'+str(i)] = peak[1]
389	parmDict['PDFsig;'+str(i)] = peak[2]
390	if 'P' in peak[3]:
391	varyList.append('PDFpos;'+str(i))
392	if 'M' in peak[3]:
393	varyList.append('PDFmag;'+str(i))
394	if 'S' in peak[3]:
395	varyList.append('PDFsig;'+str(i))
396	return parmDict,varyList
397
398	def SetParms(peaks,parmDict,varyList):
399	if 'slope' in varyList:
400	peaks['Background'][1][1] = parmDict['slope']
401	for i,peak in enumerate(peaks['Peaks']):
402	if 'PDFpos;'+str(i) in varyList:
403	peak[0] = parmDict['PDFpos;'+str(i)]
404	if 'PDFmag;'+str(i) in varyList:
405	peak[1] = parmDict['PDFmag;'+str(i)]
406	if 'PDFsig;'+str(i) in varyList:
407	peak[2] = parmDict['PDFsig;'+str(i)]
408
409
410	def CalcPDFpeaks(parmdict,Xdata):
411	Z = parmDict['slope']*Xdata
412	ipeak = 0
413	while True:
414	try:
415	pos = parmdict['PDFpos;'+str(ipeak)]
416	mag = parmdict['PDFmag;'+str(ipeak)]
417	wid = parmdict['PDFsig;'+str(ipeak)]
418	wid2 = 2.wid*2
419	Z += magrs2pinp.exp(-(Xdata-pos)**2/wid2)/wid
420	ipeak += 1
421	except KeyError: #no more peaks to process
422	return Z
423
424	def errPDFProfile(values,xdata,ydata,parmdict,varylist):
425	parmdict.update(zip(varylist,values))
426	M = CalcPDFpeaks(parmdict,xdata)-ydata
427	return M
428
429	newpeaks = copy.copy(peaks)
430	iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
431	iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
432	X = data[1][0][iBeg:iFin]
433	Y = data[1][1][iBeg:iFin]
434	parmDict,varyList = MakeParms(peaks)
435	if not len(varyList):
436	print ' Nothing varied'
437	return newpeaks,None,None,None,None,None
438
439	Rvals = {}
440	values = np.array(Dict2Values(parmDict, varyList))
441	result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
442	args=(X,Y,parmDict,varyList))
443	chisq = np.sum(result[2]['fvec']**2)
444	Values2Dict(parmDict, varyList, result[0])
445	SetParms(peaks,parmDict,varyList)
446	Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y*2))100. #to %
447	chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values)) #reduced chi^2 = M/(Nobs-Nvar)
448	sigList = list(np.sqrt(chisq*np.diag(result[1])))
449	Z = CalcPDFpeaks(parmDict,X)
450	newpeaks['calc'] = [X,Z]
451	return newpeaks,result[0],varyList,sigList,parmDict,Rvals
452
453	def MakeRDF(RDFcontrols,background,inst,pwddata):
454	import scipy.signal as signal
455	auxPlot = []
456	if 'C' in inst['Type'][0]:
457	Tth = pwddata[0]
458	wave = G2mth.getWave(inst)
459	minQ = npT2q(Tth[0],wave)
460	maxQ = npT2q(Tth[-1],wave)
461	powQ = npT2q(Tth,wave)
462	elif 'T' in inst['Type'][0]:
463	TOF = pwddata[0]
464	difC = inst['difC'][1]
465	minQ = 2.np.pidifC/TOF[-1]
466	maxQ = 2.np.pidifC/TOF[0]
467	powQ = 2.np.pidifC/TOF
468	piDQ = np.pi/(maxQ-minQ)
469	Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
470	if RDFcontrols['UseObsCalc']:
471	Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
472	else:
473	Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
474	Qdata = np.sin((Qpoints-minQ)piDQ)/piDQ
475	Qdata = 0.5np.sqrt(Qpoints) #Qbin normalization
476	# GSASIIpath.IPyBreak()
477	dq = Qpoints[1]-Qpoints[0]
478	nR = len(Qdata)
479	R = 0.5np.pinp.linspace(0,nR,nR)/(4.*maxQ)
480	iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
481	bBut,aBut = signal.butter(4,0.01)
482	Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
483	# auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
484	# auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
485	DofR = dqnp.imag(fft.fft(Qsmooth,16nR)[:nR])
486	# DofR = dqnp.imag(ft.fft(Qsmooth,16nR)[:nR])
487	auxPlot.append([R[:iFin],DofR[:iFin],'D(R)'])
488	return auxPlot
489
490	################################################################################
491	#GSASII peak fitting routines: Finger, Cox & Jephcoat model
492	################################################################################
493
494	def factorize(num):
495	''' Provide prime number factors for integer num
496	:returns: dictionary of prime factors (keys) & power for each (data)
497	'''
498	factors = {}
499	orig = num
500
501	# we take advantage of the fact that (i +1)2 = i2 + 2*i +1
502	i, sqi = 2, 4
503	while sqi <= num:
504	while not num%i:
505	num /= i
506	factors[i] = factors.get(i, 0) + 1
507
508	sqi += 2*i + 1
509	i += 1
510
511	if num != 1 and num != orig:
512	factors[num] = factors.get(num, 0) + 1
513
514	if factors:
515	return factors
516	else:
517	return {num:1} #a prime number!
518
519	def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
520	''' Provide list of optimal data sizes for FFT calculations
521
522	:param int nmin: minimum data size >= 1
523	:param int nmax: maximum data size > nmin
524	:param int thresh: maximum prime factor allowed
525	:Returns: list of data sizes where the maximum prime factor is < thresh
526	'''
527	plist = []
528	nmin = max(1,nmin)
529	nmax = max(nmin+1,nmax)
530	for p in range(nmin,nmax):
531	if max(factorize(p).keys()) < thresh:
532	plist.append(p)
533	return plist
534
535	np.seterr(divide='ignore')
536
537	# Normal distribution
538
539	# loc = mu, scale = std
540	_norm_pdf_C = 1./math.sqrt(2*math.pi)
541	class norm_gen(st.rv_continuous):
542	'needs a doc string'
543
544	def pdf(self,x,args,*kwds):
545	loc,scale=kwds['loc'],kwds['scale']
546	x = (x-loc)/scale
547	return np.exp(-x*2/2.0) _norm_pdf_C / scale
548
549	norm = norm_gen(name='norm',longname='A normal',extradoc="""
550
551	Normal distribution
552
553	The location (loc) keyword specifies the mean.
554	The scale (scale) keyword specifies the standard deviation.
555
556	normal.pdf(x) = exp(-x*2/2)/sqrt(2pi)
557	""")
558
559	## Cauchy
560
561	# median = loc
562
563	class cauchy_gen(st.rv_continuous):
564	'needs a doc string'
565
566	def pdf(self,x,args,*kwds):
567	loc,scale=kwds['loc'],kwds['scale']
568	x = (x-loc)/scale
569	return 1.0/np.pi/(1.0+x*x) / scale
570
571	cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
572
573	Cauchy distribution
574
575	cauchy.pdf(x) = 1/(pi(1+x*2))
576
577	This is the t distribution with one degree of freedom.
578	""")
579
580
581	#GSASII peak fitting routine: Finger, Cox & Jephcoat model
582
583
584	class fcjde_gen(st.rv_continuous):
585	"""
586	Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
587	Ref: J. Appl. Cryst. (1994) 27, 892-900.
588
589	:param x: array -1 to 1
590	:param t: 2-theta position of peak
591	:param s: sum(S/L,H/L); S: sample height, H: detector opening,
592	L: sample to detector opening distance
593	:param dx: 2-theta step size in deg
594
595	:returns: for fcj.pdf
596
597	* T = x*dx+t
598	* s = S/L+H/L
599	* if x < 0::
600
601	fcj.pdf = [1/sqrt({cos(T)2/cos(t)2}-1) - 1/s]/\|cos(T)\|
602
603	* if x >= 0: fcj.pdf = 0
604	"""
605	def _pdf(self,x,t,s,dx):
606	T = dx*x+t
607	ax2 = abs(npcosd(T))
608	ax = ax2**2
609	bx = npcosd(t)**2
610	bx = np.where(ax>bx,bx,ax)
611	fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
612	fx = np.where(fx > 0.,fx,0.0)
613	return fx
614
615	def pdf(self,x,args,*kwds):
616	loc=kwds['loc']
617	return self._pdf(x-loc,*args)
618
619	fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
620
621	def getWidthsCW(pos,sig,gam,shl):
622	'''Compute the peak widths used for computing the range of a peak
623	for constant wavelength data. On low-angle side, 50 FWHM are used,
624	on high-angle side 75 are used, low angle side extended for axial divergence
625	(for peaks above 90 deg, these are reversed.)
626	'''
627	widths = [np.sqrt(sig)/100.,gam/100.]
628	fwhm = 2.355*widths[0]+widths[1]
629	fmin = 50.(fwhm+shlabs(npcosd(pos)))
630	fmax = 75.0*fwhm
631	if pos > 90:
632	fmin,fmax = [fmax,fmin]
633	return widths,fmin,fmax
634
635	def getWidthsTOF(pos,alp,bet,sig,gam):
636	'''Compute the peak widths used for computing the range of a peak
637	for constant wavelength data. 50 FWHM are used on both sides each
638	extended by exponential coeff.
639	'''
640	widths = [np.sqrt(sig),gam]
641	fwhm = 2.355widths[0]+2.widths[1]
642	fmin = 50.fwhm(1.+1./alp)
643	fmax = 50.fwhm(1.+1./bet)
644	return widths,fmin,fmax
645
646	def getFWHM(pos,Inst):
647	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
648	via getgamFW(g,s).
649
650	:param pos: float peak position in deg 2-theta or tof in musec
651	:param Inst: dict instrument parameters
652
653	:returns float: total FWHM of pseudoVoigt in deg or musec
654	'''
655
656	sig = lambda Th,U,V,W: np.sqrt(max(0.001,Utand(Th)2+Vtand(Th)+W))
657	sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1dsp2+S2dsp4+Sq/dsp2)
658	gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))
659	gamTOF = lambda dsp,X,Y: Xdsp+Ydsp**2
660	if 'C' in Inst['Type'][0]:
661	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
662	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])
663	return getgamFW(g,s)/100. #returns FWHM in deg
664	else:
665	dsp = pos/Inst['difC'][0]
666	s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
667	g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
668	return getgamFW(g,s)
669
670	def getgamFW(g,s):
671	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
672	lambda fxn needs FWHM for both Gaussian & Lorentzian components
673
674	:param g: float Lorentzian gamma = FWHM(L)
675	:param s: float Gaussian sig
676
677	:returns float: total FWHM of pseudoVoigt
678	'''
679	gamFW = lambda s,g: np.exp(np.log(s*5+2.69269s*4g+2.42843s3g*2+4.47163s*2g*3+0.07842sg4+g*5)/5.)
680	return gamFW(2.35482s,g) #sqrt(8ln2)sig = FWHM(G)
681
682	def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
683	'''Compute the Finger-Cox-Jepcoat modified Voigt function for a
684	CW powder peak by direct convolution. This version is not used.
685	'''
686	DX = xdata[1]-xdata[0]
687	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
688	x = np.linspace(pos-fmin,pos+fmin,256)
689	dx = x[1]-x[0]
690	Norm = norm.pdf(x,loc=pos,scale=widths[0])
691	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
692	arg = [pos,shl/57.2958,dx,]
693	FCJ = fcjde.pdf(x,*arg,loc=pos)
694	if len(np.nonzero(FCJ)[0])>5:
695	z = np.column_stack([Norm,Cauchy,FCJ]).T
696	Z = fft.fft(z)
697	Df = fft.ifft(Z.prod(axis=0)).real
698	else:
699	z = np.column_stack([Norm,Cauchy]).T
700	Z = fft.fft(z)
701	Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
702	Df /= np.sum(Df)
703	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
704	return intensDf(xdata)DX/dx
705
706	def getBackground(pfx,parmDict,bakType,dataType,xdata):
707	'needs a doc string'
708	if 'T' in dataType:
709	q = 2.np.piparmDict[pfx+'difC']/xdata
710	elif 'C' in dataType:
711	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
712	q = npT2q(xdata,wave)
713	yb = np.zeros_like(xdata)
714	nBak = 0
715	cw = np.diff(xdata)
716	cw = np.append(cw,cw[-1])
717	sumBk = [0.,0.,0]
718	while True:
719	key = pfx+'Back;'+str(nBak)
720	if key in parmDict:
721	nBak += 1
722	else:
723	break
724	#empirical functions
725	if bakType in ['chebyschev','cosine']:
726	dt = xdata[-1]-xdata[0]
727	for iBak in range(nBak):
728	key = pfx+'Back;'+str(iBak)
729	if bakType == 'chebyschev':
730	ybi = parmDict[key](-1.+2.(xdata-xdata[0])/dt)**iBak
731	elif bakType == 'cosine':
732	ybi = parmDict[key]npcosd(180.xdata*iBak/xdata[-1])
733	yb += ybi
734	sumBk[0] = np.sum(yb)
735	elif bakType in ['Q^2 power series','Q^-2 power series']:
736	QT = 1.
737	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
738	for iBak in range(nBak-1):
739	key = pfx+'Back;'+str(iBak+1)
740	if '-2' in bakType:
741	QT = (iBak+1)q**-2
742	else:
743	QT = q*2/(iBak+1)
744	yb += QT*parmDict[key]
745	sumBk[0] = np.sum(yb)
746	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
747	if nBak == 1:
748	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
749	elif nBak == 2:
750	dX = xdata[-1]-xdata[0]
751	T2 = (xdata-xdata[0])/dX
752	T1 = 1.0-T2
753	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
754	else:
755	if bakType == 'lin interpolate':
756	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
757	elif bakType == 'inv interpolate':
758	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
759	elif bakType == 'log interpolate':
760	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
761	bakPos[0] = xdata[0]
762	bakPos[-1] = xdata[-1]
763	bakVals = np.zeros(nBak)
764	for i in range(nBak):
765	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
766	bakInt = si.interp1d(bakPos,bakVals,'linear')
767	yb = bakInt(ma.getdata(xdata))
768	sumBk[0] = np.sum(yb)
769	#Debye function
770	if pfx+'difC' in parmDict:
771	ff = 1.
772	else:
773	try:
774	wave = parmDict[pfx+'Lam']
775	except KeyError:
776	wave = parmDict[pfx+'Lam1']
777	SQ = (q/(4.np.pi))*2
778	FF = G2elem.GetFormFactorCoeff('Si')[0]
779	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
780	iD = 0
781	while True:
782	try:
783	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
784	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
785	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
786	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
787	yb += ybi
788	sumBk[1] += np.sum(ybi)
789	iD += 1
790	except KeyError:
791	break
792	#peaks
793	iD = 0
794	while True:
795	try:
796	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
797	pkI = parmDict[pfx+'BkPkint;'+str(iD)]
798	pkS = parmDict[pfx+'BkPksig;'+str(iD)]
799	pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
800	if 'C' in dataType:
801	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
802	else: #'T'OF
803	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
804	iBeg = np.searchsorted(xdata,pkP-fmin)
805	iFin = np.searchsorted(xdata,pkP+fmax)
806	lenX = len(xdata)
807	if not iBeg:
808	iFin = np.searchsorted(xdata,pkP+fmax)
809	elif iBeg == lenX:
810	iFin = iBeg
811	else:
812	iFin = np.searchsorted(xdata,pkP+fmax)
813	if 'C' in dataType:
814	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
815	yb[iBeg:iFin] += ybi
816	else: #'T'OF
817	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
818	yb[iBeg:iFin] += ybi
819	sumBk[2] += np.sum(ybi)
820	iD += 1
821	except KeyError:
822	break
823	except ValueError:
824	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
825	break
826	return yb,sumBk
827
828	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
829	'needs a doc string'
830	if 'T' in dataType:
831	q = 2.np.piparmDict[hfx+'difC']/xdata
832	elif 'C' in dataType:
833	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
834	q = 2.np.pinpsind(xdata/2.)/wave
835	nBak = 0
836	while True:
837	key = hfx+'Back;'+str(nBak)
838	if key in parmDict:
839	nBak += 1
840	else:
841	break
842	dydb = np.zeros(shape=(nBak,len(xdata)))
843	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
844	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
845	cw = np.diff(xdata)
846	cw = np.append(cw,cw[-1])
847
848	if bakType in ['chebyschev','cosine']:
849	dt = xdata[-1]-xdata[0]
850	for iBak in range(nBak):
851	if bakType == 'chebyschev':
852	dydb[iBak] = (-1.+2.(xdata-xdata[0])/dt)*iBak
853	elif bakType == 'cosine':
854	dydb[iBak] = npcosd(180.xdataiBak/xdata[-1])
855	elif bakType in ['Q^2 power series','Q^-2 power series']:
856	QT = 1.
857	dydb[0] = np.ones_like(xdata)
858	for iBak in range(nBak-1):
859	if '-2' in bakType:
860	QT = (iBak+1)q**-2
861	else:
862	QT = q*2/(iBak+1)
863	dydb[iBak+1] = QT
864	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
865	if nBak == 1:
866	dydb[0] = np.ones_like(xdata)
867	elif nBak == 2:
868	dX = xdata[-1]-xdata[0]
869	T2 = (xdata-xdata[0])/dX
870	T1 = 1.0-T2
871	dydb = [T1,T2]
872	else:
873	if bakType == 'lin interpolate':
874	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
875	elif bakType == 'inv interpolate':
876	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
877	elif bakType == 'log interpolate':
878	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
879	bakPos[0] = xdata[0]
880	bakPos[-1] = xdata[-1]
881	for i,pos in enumerate(bakPos):
882	if i == 0:
883	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
884	elif i == len(bakPos)-1:
885	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
886	else:
887	dydb[i] = np.where(xdata>bakPos[i],
888	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
889	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
890	if hfx+'difC' in parmDict:
891	ff = 1.
892	else:
893	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
894	q = npT2q(xdata,wave)
895	SQ = (q/(4np.pi))*2
896	FF = G2elem.GetFormFactorCoeff('Si')[0]
897	ff = np.array(G2elem.ScatFac(FF,SQ)[0])np.pi*2 #needs pi^2~10. for cw data (why?)
898	iD = 0
899	while True:
900	try:
901	if hfx+'difC' in parmDict:
902	q = 2np.piparmDict[hfx+'difC']/xdata
903	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
904	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
905	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
906	sqr = np.sin(qdbR)/(qdbR)
907	cqr = np.cos(q*dbR)
908	temp = np.exp(-dbUq*2)
909	dyddb[3iD] = ffsqr*temp
910	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(dbR)
911	dyddb[3iD+2] = -ffdbAsqrtempq*2
912	iD += 1
913	except KeyError:
914	break
915	iD = 0
916	while True:
917	try:
918	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
919	pkI = parmDict[hfx+'BkPkint;'+str(iD)]
920	pkS = parmDict[hfx+'BkPksig;'+str(iD)]
921	pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
922	if 'C' in dataType:
923	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
924	else: #'T'OF
925	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
926	iBeg = np.searchsorted(xdata,pkP-fmin)
927	iFin = np.searchsorted(xdata,pkP+fmax)
928	lenX = len(xdata)
929	if not iBeg:
930	iFin = np.searchsorted(xdata,pkP+fmax)
931	elif iBeg == lenX:
932	iFin = iBeg
933	else:
934	iFin = np.searchsorted(xdata,pkP+fmax)
935	if 'C' in dataType:
936	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
937	dydpk[4iD][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdp
938	dydpk[4iD+1][iBeg:iFin] += 100.cw[iBeg:iFin]*Df
939	dydpk[4iD+2][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFds
940	dydpk[4iD+3][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdg
941	else: #'T'OF
942	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
943	dydpk[4iD][iBeg:iFin] += pkIdFdp
944	dydpk[4*iD+1][iBeg:iFin] += Df
945	dydpk[4iD+2][iBeg:iFin] += pkIdFds
946	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
947	iD += 1
948	except KeyError:
949	break
950	except ValueError:
951	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
952	break
953	return dydb,dyddb,dydpk
954
955	#use old fortran routine
956	def getFCJVoigt3(pos,sig,gam,shl,xdata):
957	'''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
958	CW powder peak in external Fortran routine
959	'''
960	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
961	# Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
962	Df /= np.sum(Df)
963	return Df
964
965	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
966	'''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
967	function for a CW powder peak
968	'''
969	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
970	# Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
971	return Df,dFdp,dFds,dFdg,dFdsh
972
973	def getPsVoigt(pos,sig,gam,xdata):
974	'needs a doc string'
975
976	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
977	Df /= np.sum(Df)
978	return Df
979
980	def getdPsVoigt(pos,sig,gam,xdata):
981	'needs a doc string'
982
983	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
984	return Df,dFdp,dFds,dFdg
985
986	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
987	'needs a doc string'
988	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
989	Df /= np.sum(Df)
990	return Df
991
992	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
993	'needs a doc string'
994	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
995	return Df,dFdp,dFda,dFdb,dFds,dFdg
996
997	def ellipseSize(H,Sij,GB):
998	'needs a doc string'
999	HX = np.inner(H.T,GB)
1000	lenHX = np.sqrt(np.sum(HX**2))
1001	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
1002	R = np.inner(HX/lenHX,Rsize)*Esize #want column length for hkl in crystal
1003	lenR = np.sqrt(np.sum(R**2))
1004	return lenR
1005
1006	def ellipseSizeDerv(H,Sij,GB):
1007	'needs a doc string'
1008	lenR = ellipseSize(H,Sij,GB)
1009	delt = 0.001
1010	dRdS = np.zeros(6)
1011	for i in range(6):
1012	Sij[i] -= delt
1013	lenM = ellipseSize(H,Sij,GB)
1014	Sij[i] += 2.*delt
1015	lenP = ellipseSize(H,Sij,GB)
1016	Sij[i] -= delt
1017	dRdS[i] = (lenP-lenM)/(2.*delt)
1018	return lenR,dRdS
1019
1020	def getHKLpeak(dmin,SGData,A,Inst=None):
1021	'needs a doc string'
1022	HKL = G2lat.GenHLaue(dmin,SGData,A)
1023	HKLs = []
1024	for h,k,l,d in HKL:
1025	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1026	if not ext:
1027	if Inst == None:
1028	HKLs.append([h,k,l,d,0,-1])
1029	else:
1030	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1031	return HKLs
1032
1033	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1034	'needs a doc string'
1035	HKLs = []
1036	vec = np.array(Vec)
1037	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1038	dvec = 1./(maxH*vstar+1./dmin)
1039	HKL = G2lat.GenHLaue(dvec,SGData,A)
1040	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1041	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1042	for h,k,l,d in HKL:
1043	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1044	if not ext and d >= dmin:
1045	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1046	for dH in SSdH:
1047	if dH:
1048	DH = SSdH[dH]
1049	H = [h+DH[0],k+DH[1],l+DH[2]]
1050	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1051	if d >= dmin:
1052	HKLM = np.array([h,k,l,dH])
1053	if G2spc.checkSSextc(HKLM,SSGData):
1054	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
1055	# GSASIIpath.IPyBreak()
1056	return G2lat.sortHKLd(HKLs,True,True,True)
1057
1058	def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1059	'needs a doc string'
1060
1061	yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1062	yc = np.zeros_like(yb)
1063	cw = np.diff(xdata)
1064	cw = np.append(cw,cw[-1])
1065	if 'C' in dataType:
1066	shl = max(parmDict['SH/L'],0.002)
1067	Ka2 = False
1068	if 'Lam1' in parmDict.keys():
1069	Ka2 = True
1070	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1071	kRatio = parmDict['I(L2)/I(L1)']
1072	iPeak = 0
1073	while True:
1074	try:
1075	pos = parmDict['pos'+str(iPeak)]
1076	tth = (pos-parmDict['Zero'])
1077	intens = parmDict['int'+str(iPeak)]
1078	sigName = 'sig'+str(iPeak)
1079	if sigName in varyList:
1080	sig = parmDict[sigName]
1081	else:
1082	sig = G2mth.getCWsig(parmDict,tth)
1083	sig = max(sig,0.001) #avoid neg sigma^2
1084	gamName = 'gam'+str(iPeak)
1085	if gamName in varyList:
1086	gam = parmDict[gamName]
1087	else:
1088	gam = G2mth.getCWgam(parmDict,tth)
1089	gam = max(gam,0.001) #avoid neg gamma
1090	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1091	iBeg = np.searchsorted(xdata,pos-fmin)
1092	iFin = np.searchsorted(xdata,pos+fmin)
1093	if not iBeg+iFin: #peak below low limit
1094	iPeak += 1
1095	continue
1096	elif not iBeg-iFin: #peak above high limit
1097	return yb+yc
1098	yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1099	if Ka2:
1100	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1101	iBeg = np.searchsorted(xdata,pos2-fmin)
1102	iFin = np.searchsorted(xdata,pos2+fmin)
1103	if iBeg-iFin:
1104	yc[iBeg:iFin] += intenskRatiogetFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1105	iPeak += 1
1106	except KeyError: #no more peaks to process
1107	return yb+yc
1108	else:
1109	Pdabc = parmDict['Pdabc']
1110	difC = parmDict['difC']
1111	iPeak = 0
1112	while True:
1113	try:
1114	pos = parmDict['pos'+str(iPeak)]
1115	tof = pos-parmDict['Zero']
1116	dsp = tof/difC
1117	intens = parmDict['int'+str(iPeak)]
1118	alpName = 'alp'+str(iPeak)
1119	if alpName in varyList:
1120	alp = parmDict[alpName]
1121	else:
1122	if len(Pdabc):
1123	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1124	else:
1125	alp = G2mth.getTOFalpha(parmDict,dsp)
1126	alp = max(0.0001,alp)
1127	betName = 'bet'+str(iPeak)
1128	if betName in varyList:
1129	bet = parmDict[betName]
1130	else:
1131	if len(Pdabc):
1132	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1133	else:
1134	bet = G2mth.getTOFbeta(parmDict,dsp)
1135	bet = max(0.0001,bet)
1136	sigName = 'sig'+str(iPeak)
1137	if sigName in varyList:
1138	sig = parmDict[sigName]
1139	else:
1140	sig = G2mth.getTOFsig(parmDict,dsp)
1141	gamName = 'gam'+str(iPeak)
1142	if gamName in varyList:
1143	gam = parmDict[gamName]
1144	else:
1145	gam = G2mth.getTOFgamma(parmDict,dsp)
1146	gam = max(gam,0.001) #avoid neg gamma
1147	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1148	iBeg = np.searchsorted(xdata,pos-fmin)
1149	iFin = np.searchsorted(xdata,pos+fmax)
1150	lenX = len(xdata)
1151	if not iBeg:
1152	iFin = np.searchsorted(xdata,pos+fmax)
1153	elif iBeg == lenX:
1154	iFin = iBeg
1155	else:
1156	iFin = np.searchsorted(xdata,pos+fmax)
1157	if not iBeg+iFin: #peak below low limit
1158	iPeak += 1
1159	continue
1160	elif not iBeg-iFin: #peak above high limit
1161	return yb+yc
1162	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1163	iPeak += 1
1164	except KeyError: #no more peaks to process
1165	return yb+yc
1166
1167	def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1168	'needs a doc string'
1169	# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1170	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1171	dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1172	if 'Back;0' in varyList: #background derivs are in front if present
1173	dMdv[0:len(dMdb)] = dMdb
1174	names = ['DebyeA','DebyeR','DebyeU']
1175	for name in varyList:
1176	if 'Debye' in name:
1177	parm,id = name.split(';')
1178	ip = names.index(parm)
1179	dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1180	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1181	for name in varyList:
1182	if 'BkPk' in name:
1183	parm,id = name.split(';')
1184	ip = names.index(parm)
1185	dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1186	cw = np.diff(xdata)
1187	cw = np.append(cw,cw[-1])
1188	if 'C' in dataType:
1189	shl = max(parmDict['SH/L'],0.002)
1190	Ka2 = False
1191	if 'Lam1' in parmDict.keys():
1192	Ka2 = True
1193	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1194	kRatio = parmDict['I(L2)/I(L1)']
1195	iPeak = 0
1196	while True:
1197	try:
1198	pos = parmDict['pos'+str(iPeak)]
1199	tth = (pos-parmDict['Zero'])
1200	intens = parmDict['int'+str(iPeak)]
1201	sigName = 'sig'+str(iPeak)
1202	if sigName in varyList:
1203	sig = parmDict[sigName]
1204	dsdU = dsdV = dsdW = 0
1205	else:
1206	sig = G2mth.getCWsig(parmDict,tth)
1207	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1208	sig = max(sig,0.001) #avoid neg sigma
1209	gamName = 'gam'+str(iPeak)
1210	if gamName in varyList:
1211	gam = parmDict[gamName]
1212	dgdX = dgdY = 0
1213	else:
1214	gam = G2mth.getCWgam(parmDict,tth)
1215	dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1216	gam = max(gam,0.001) #avoid neg gamma
1217	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1218	iBeg = np.searchsorted(xdata,pos-fmin)
1219	iFin = np.searchsorted(xdata,pos+fmin)
1220	if not iBeg+iFin: #peak below low limit
1221	iPeak += 1
1222	continue
1223	elif not iBeg-iFin: #peak above high limit
1224	break
1225	dMdpk = np.zeros(shape=(6,len(xdata)))
1226	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1227	for i in range(1,5):
1228	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intens*dMdipk[i]
1229	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk[0]
1230	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1231	if Ka2:
1232	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1233	iBeg = np.searchsorted(xdata,pos2-fmin)
1234	iFin = np.searchsorted(xdata,pos2+fmin)
1235	if iBeg-iFin:
1236	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1237	for i in range(1,5):
1238	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intenskRatiodMdipk2[i]
1239	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]kRatio*dMdipk2[0]
1240	dMdpk[5][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk2[0]
1241	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1242	for parmName in ['pos','int','sig','gam']:
1243	try:
1244	idx = varyList.index(parmName+str(iPeak))
1245	dMdv[idx] = dervDict[parmName]
1246	except ValueError:
1247	pass
1248	if 'U' in varyList:
1249	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1250	if 'V' in varyList:
1251	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1252	if 'W' in varyList:
1253	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1254	if 'X' in varyList:
1255	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1256	if 'Y' in varyList:
1257	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1258	if 'SH/L' in varyList:
1259	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1260	if 'I(L2)/I(L1)' in varyList:
1261	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1262	iPeak += 1
1263	except KeyError: #no more peaks to process
1264	break
1265	else:
1266	Pdabc = parmDict['Pdabc']
1267	difC = parmDict['difC']
1268	iPeak = 0
1269	while True:
1270	try:
1271	pos = parmDict['pos'+str(iPeak)]
1272	tof = pos-parmDict['Zero']
1273	dsp = tof/difC
1274	intens = parmDict['int'+str(iPeak)]
1275	alpName = 'alp'+str(iPeak)
1276	if alpName in varyList:
1277	alp = parmDict[alpName]
1278	else:
1279	if len(Pdabc):
1280	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1281	dada0 = 0
1282	else:
1283	alp = G2mth.getTOFalpha(parmDict,dsp)
1284	dada0 = G2mth.getTOFalphaDeriv(dsp)
1285	betName = 'bet'+str(iPeak)
1286	if betName in varyList:
1287	bet = parmDict[betName]
1288	else:
1289	if len(Pdabc):
1290	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1291	dbdb0 = dbdb1 = dbdb2 = 0
1292	else:
1293	bet = G2mth.getTOFbeta(parmDict,dsp)
1294	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1295	sigName = 'sig'+str(iPeak)
1296	if sigName in varyList:
1297	sig = parmDict[sigName]
1298	dsds0 = dsds1 = dsds2 = dsds3 = 0
1299	else:
1300	sig = G2mth.getTOFsig(parmDict,dsp)
1301	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1302	gamName = 'gam'+str(iPeak)
1303	if gamName in varyList:
1304	gam = parmDict[gamName]
1305	dsdX = dsdY = 0
1306	else:
1307	gam = G2mth.getTOFgamma(parmDict,dsp)
1308	dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1309	gam = max(gam,0.001) #avoid neg gamma
1310	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1311	iBeg = np.searchsorted(xdata,pos-fmin)
1312	lenX = len(xdata)
1313	if not iBeg:
1314	iFin = np.searchsorted(xdata,pos+fmax)
1315	elif iBeg == lenX:
1316	iFin = iBeg
1317	else:
1318	iFin = np.searchsorted(xdata,pos+fmax)
1319	if not iBeg+iFin: #peak below low limit
1320	iPeak += 1
1321	continue
1322	elif not iBeg-iFin: #peak above high limit
1323	break
1324	dMdpk = np.zeros(shape=(7,len(xdata)))
1325	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1326	for i in range(1,6):
1327	dMdpk[i][iBeg:iFin] += intenscw[iBeg:iFin]dMdipk[i]
1328	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1329	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1330	for parmName in ['pos','int','alp','bet','sig','gam']:
1331	try:
1332	idx = varyList.index(parmName+str(iPeak))
1333	dMdv[idx] = dervDict[parmName]
1334	except ValueError:
1335	pass
1336	if 'alpha' in varyList:
1337	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1338	if 'beta-0' in varyList:
1339	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1340	if 'beta-1' in varyList:
1341	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1342	if 'beta-q' in varyList:
1343	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1344	if 'sig-0' in varyList:
1345	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1346	if 'sig-1' in varyList:
1347	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1348	if 'sig-2' in varyList:
1349	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1350	if 'sig-q' in varyList:
1351	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1352	if 'X' in varyList:
1353	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1354	if 'Y' in varyList:
1355	dMdv[varyList.index('Y')] += dsdY*dervDict['gam'] #problem here
1356	iPeak += 1
1357	except KeyError: #no more peaks to process
1358	break
1359	return dMdv
1360
1361	def Dict2Values(parmdict, varylist):
1362	'''Use before call to leastsq to setup list of values for the parameters
1363	in parmdict, as selected by key in varylist'''
1364	return [parmdict[key] for key in varylist]
1365
1366	def Values2Dict(parmdict, varylist, values):
1367	''' Use after call to leastsq to update the parameter dictionary with
1368	values corresponding to keys in varylist'''
1369	parmdict.update(zip(varylist,values))
1370
1371	def SetBackgroundParms(Background):
1372	'needs a doc string'
1373	if len(Background) == 1: # fix up old backgrounds
1374	Background.append({'nDebye':0,'debyeTerms':[]})
1375	bakType,bakFlag = Background[0][:2]
1376	backVals = Background[0][3:]
1377	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1378	Debye = Background[1] #also has background peaks stuff
1379	backDict = dict(zip(backNames,backVals))
1380	backVary = []
1381	if bakFlag:
1382	backVary = backNames
1383
1384	backDict['nDebye'] = Debye['nDebye']
1385	debyeDict = {}
1386	debyeList = []
1387	for i in range(Debye['nDebye']):
1388	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1389	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1390	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1391	debyeVary = []
1392	for item in debyeList:
1393	if item[1]:
1394	debyeVary.append(item[0])
1395	backDict.update(debyeDict)
1396	backVary += debyeVary
1397
1398	backDict['nPeaks'] = Debye['nPeaks']
1399	peaksDict = {}
1400	peaksList = []
1401	for i in range(Debye['nPeaks']):
1402	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1403	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1404	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1405	peaksVary = []
1406	for item in peaksList:
1407	if item[1]:
1408	peaksVary.append(item[0])
1409	backDict.update(peaksDict)
1410	backVary += peaksVary
1411	return bakType,backDict,backVary
1412
1413	def DoCalibInst(IndexPeaks,Inst):
1414
1415	def SetInstParms():
1416	dataType = Inst['Type'][0]
1417	insVary = []
1418	insNames = []
1419	insVals = []
1420	for parm in Inst:
1421	insNames.append(parm)
1422	insVals.append(Inst[parm][1])
1423	if parm in ['Lam','difC','difA','difB','Zero',]:
1424	if Inst[parm][2]:
1425	insVary.append(parm)
1426	instDict = dict(zip(insNames,insVals))
1427	return dataType,instDict,insVary
1428
1429	def GetInstParms(parmDict,Inst,varyList):
1430	for name in Inst:
1431	Inst[name][1] = parmDict[name]
1432
1433	def InstPrint(Inst,sigDict):
1434	print 'Instrument Parameters:'
1435	if 'C' in Inst['Type'][0]:
1436	ptfmt = "%12.6f"
1437	else:
1438	ptfmt = "%12.3f"
1439	ptlbls = 'names :'
1440	ptstr = 'values:'
1441	sigstr = 'esds :'
1442	for parm in Inst:
1443	if parm in ['Lam','difC','difA','difB','Zero',]:
1444	ptlbls += "%s" % (parm.center(12))
1445	ptstr += ptfmt % (Inst[parm][1])
1446	if parm in sigDict:
1447	sigstr += ptfmt % (sigDict[parm])
1448	else:
1449	sigstr += 12*' '
1450	print ptlbls
1451	print ptstr
1452	print sigstr
1453
1454	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1455	parmDict.update(zip(varyList,values))
1456	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1457
1458	peakPos = []
1459	peakDsp = []
1460	peakWt = []
1461	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1462	if peak[2] and peak[3] and sig > 0.:
1463	peakPos.append(peak[0])
1464	peakDsp.append(peak[-1]) #d-calc
1465	# peakWt.append(peak[-1]2/sig2) #weight by d**2
1466	peakWt.append(1./(sig*peak[-1])) #
1467	peakPos = np.array(peakPos)
1468	peakDsp = np.array(peakDsp)
1469	peakWt = np.array(peakWt)
1470	dataType,insDict,insVary = SetInstParms()
1471	parmDict = {}
1472	parmDict.update(insDict)
1473	varyList = insVary
1474	if not len(varyList):
1475	print '** ERROR - nothing to refine! **'
1476	return False
1477	while True:
1478	begin = time.time()
1479	values = np.array(Dict2Values(parmDict, varyList))
1480	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1481	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1482	ncyc = int(result[2]['nfev']/2)
1483	runtime = time.time()-begin
1484	chisq = np.sum(result[2]['fvec']**2)
1485	Values2Dict(parmDict, varyList, result[0])
1486	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
1487	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1488	print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1489	print 'chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF)
1490	try:
1491	sig = np.sqrt(np.diag(result[1])*GOF)
1492	if np.any(np.isnan(sig)):
1493	print '* Least squares aborted - some invalid esds possible *'
1494	break #refinement succeeded - finish up!
1495	except ValueError: #result[1] is None on singular matrix
1496	print '** Refinement failed - singular matrix **'
1497
1498	sigDict = dict(zip(varyList,sig))
1499	GetInstParms(parmDict,Inst,varyList)
1500	InstPrint(Inst,sigDict)
1501	return True
1502
1503	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1504	'''Called to perform a peak fit, refining the selected items in the peak
1505	table as well as selected items in the background.
1506
1507	:param str FitPgm: type of fit to perform. At present this is ignored.
1508	:param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1509	four values followed by a refine flag where the values are: position, intensity,
1510	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1511	tree entry, dict item "peaks"
1512	:param list Background: describes the background. List with two items.
1513	Item 0 specifies a background model and coefficients. Item 1 is a dict.
1514	From the Histogram/Background tree entry.
1515	:param list Limits: min and max x-value to use
1516	:param dict Inst: Instrument parameters
1517	:param dict Inst2: more Instrument parameters
1518	:param numpy.array data: a 5xn array. data[0] is the x-values,
1519	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1520	calc, background and difference intensities, respectively.
1521	:param list prevVaryList: Used in sequential refinements to override the
1522	variable list. Defaults as an empty list.
1523	:param bool oneCycle: True if only one cycle of fitting should be performed
1524	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1525	and derivType = controls['deriv type']. If None default values are used.
1526	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1527	Defaults to None, which means no updates are done.
1528	'''
1529	def GetBackgroundParms(parmList,Background):
1530	iBak = 0
1531	while True:
1532	try:
1533	bakName = 'Back;'+str(iBak)
1534	Background[0][iBak+3] = parmList[bakName]
1535	iBak += 1
1536	except KeyError:
1537	break
1538	iDb = 0
1539	while True:
1540	names = ['DebyeA;','DebyeR;','DebyeU;']
1541	try:
1542	for i,name in enumerate(names):
1543	val = parmList[name+str(iDb)]
1544	Background[1]['debyeTerms'][iDb][2*i] = val
1545	iDb += 1
1546	except KeyError:
1547	break
1548	iDb = 0
1549	while True:
1550	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1551	try:
1552	for i,name in enumerate(names):
1553	val = parmList[name+str(iDb)]
1554	Background[1]['peaksList'][iDb][2*i] = val
1555	iDb += 1
1556	except KeyError:
1557	break
1558
1559	def BackgroundPrint(Background,sigDict):
1560	print 'Background coefficients for',Background[0][0],'function'
1561	ptfmt = "%12.5f"
1562	ptstr = 'value: '
1563	sigstr = 'esd : '
1564	for i,back in enumerate(Background[0][3:]):
1565	ptstr += ptfmt % (back)
1566	if Background[0][1]:
1567	prm = 'Back;'+str(i)
1568	if prm in sigDict:
1569	sigstr += ptfmt % (sigDict[prm])
1570	else:
1571	sigstr += " "*12
1572	if len(ptstr) > 75:
1573	print ptstr
1574	if Background[0][1]: print sigstr
1575	ptstr = 'value: '
1576	sigstr = 'esd : '
1577	if len(ptstr) > 8:
1578	print ptstr
1579	if Background[0][1]: print sigstr
1580
1581	if Background[1]['nDebye']:
1582	parms = ['DebyeA;','DebyeR;','DebyeU;']
1583	print 'Debye diffuse scattering coefficients'
1584	ptfmt = "%12.5f"
1585	print ' term DebyeA esd DebyeR esd DebyeU esd'
1586	for term in range(Background[1]['nDebye']):
1587	line = ' term %d'%(term)
1588	for ip,name in enumerate(parms):
1589	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1590	if name+str(term) in sigDict:
1591	line += ptfmt%(sigDict[name+str(term)])
1592	else:
1593	line += " "*12
1594	print line
1595	if Background[1]['nPeaks']:
1596	print 'Coefficients for Background Peaks'
1597	ptfmt = "%15.3f"
1598	for j,pl in enumerate(Background[1]['peaksList']):
1599	names = 'peak %3d:'%(j+1)
1600	ptstr = 'values :'
1601	sigstr = 'esds :'
1602	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1603	val = pl[2*i]
1604	prm = lbl+";"+str(j)
1605	names += '%15s'%(prm)
1606	ptstr += ptfmt%(val)
1607	if prm in sigDict:
1608	sigstr += ptfmt%(sigDict[prm])
1609	else:
1610	sigstr += " "*15
1611	print names
1612	print ptstr
1613	print sigstr
1614
1615	def SetInstParms(Inst):
1616	dataType = Inst['Type'][0]
1617	insVary = []
1618	insNames = []
1619	insVals = []
1620	for parm in Inst:
1621	insNames.append(parm)
1622	insVals.append(Inst[parm][1])
1623	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1624	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1625	insVary.append(parm)
1626	instDict = dict(zip(insNames,insVals))
1627	instDict['X'] = max(instDict['X'],0.01)
1628	instDict['Y'] = max(instDict['Y'],0.01)
1629	if 'SH/L' in instDict:
1630	instDict['SH/L'] = max(instDict['SH/L'],0.002)
1631	return dataType,instDict,insVary
1632
1633	def GetInstParms(parmDict,Inst,varyList,Peaks):
1634	for name in Inst:
1635	Inst[name][1] = parmDict[name]
1636	iPeak = 0
1637	while True:
1638	try:
1639	sigName = 'sig'+str(iPeak)
1640	pos = parmDict['pos'+str(iPeak)]
1641	if sigName not in varyList:
1642	if 'C' in Inst['Type'][0]:
1643	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1644	else:
1645	dsp = G2lat.Pos2dsp(Inst,pos)
1646	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1647	gamName = 'gam'+str(iPeak)
1648	if gamName not in varyList:
1649	if 'C' in Inst['Type'][0]:
1650	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1651	else:
1652	dsp = G2lat.Pos2dsp(Inst,pos)
1653	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1654	iPeak += 1
1655	except KeyError:
1656	break
1657
1658	def InstPrint(Inst,sigDict):
1659	print 'Instrument Parameters:'
1660	ptfmt = "%12.6f"
1661	ptlbls = 'names :'
1662	ptstr = 'values:'
1663	sigstr = 'esds :'
1664	for parm in Inst:
1665	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1666	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1667	ptlbls += "%s" % (parm.center(12))
1668	ptstr += ptfmt % (Inst[parm][1])
1669	if parm in sigDict:
1670	sigstr += ptfmt % (sigDict[parm])
1671	else:
1672	sigstr += 12*' '
1673	print ptlbls
1674	print ptstr
1675	print sigstr
1676
1677	def SetPeaksParms(dataType,Peaks):
1678	peakNames = []
1679	peakVary = []
1680	peakVals = []
1681	if 'C' in dataType:
1682	names = ['pos','int','sig','gam']
1683	else:
1684	names = ['pos','int','alp','bet','sig','gam']
1685	for i,peak in enumerate(Peaks):
1686	for j,name in enumerate(names):
1687	peakVals.append(peak[2*j])
1688	parName = name+str(i)
1689	peakNames.append(parName)
1690	if peak[2*j+1]:
1691	peakVary.append(parName)
1692	return dict(zip(peakNames,peakVals)),peakVary
1693
1694	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1695	if 'C' in Inst['Type'][0]:
1696	names = ['pos','int','sig','gam']
1697	else: #'T'
1698	names = ['pos','int','alp','bet','sig','gam']
1699	for i,peak in enumerate(Peaks):
1700	pos = parmDict['pos'+str(i)]
1701	if 'difC' in Inst:
1702	dsp = pos/Inst['difC'][1]
1703	for j in range(len(names)):
1704	parName = names[j]+str(i)
1705	if parName in varyList:
1706	peak[2*j] = parmDict[parName]
1707	elif 'alpha' in parName:
1708	peak[2*j] = parmDict['alpha']/dsp
1709	elif 'beta' in parName:
1710	peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1711	elif 'sig' in parName:
1712	if 'C' in Inst['Type'][0]:
1713	peak[2*j] = G2mth.getCWsig(parmDict,pos)
1714	else:
1715	peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1716	elif 'gam' in parName:
1717	if 'C' in Inst['Type'][0]:
1718	peak[2*j] = G2mth.getCWgam(parmDict,pos)
1719	else:
1720	peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1721
1722	def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1723	print 'Peak coefficients:'
1724	if 'C' in dataType:
1725	names = ['pos','int','sig','gam']
1726	else: #'T'
1727	names = ['pos','int','alp','bet','sig','gam']
1728	head = 13*' '
1729	for name in names:
1730	if name in ['alp','bet']:
1731	head += name.center(8)+'esd'.center(8)
1732	else:
1733	head += name.center(10)+'esd'.center(10)
1734	head += 'bins'.center(8)
1735	print head
1736	if 'C' in dataType:
1737	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1738	else:
1739	ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1740	for i,peak in enumerate(Peaks):
1741	ptstr = ':'
1742	for j in range(len(names)):
1743	name = names[j]
1744	parName = name+str(i)
1745	ptstr += ptfmt[name] % (parmDict[parName])
1746	if parName in varyList:
1747	ptstr += ptfmt[name] % (sigDict[parName])
1748	else:
1749	if name in ['alp','bet']:
1750	ptstr += 8*' '
1751	else:
1752	ptstr += 10*' '
1753	ptstr += '%9.2f'%(ptsperFW[i])
1754	print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1755
1756	def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1757	parmdict.update(zip(varylist,values))
1758	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1759
1760	def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):
1761	parmdict.update(zip(varylist,values))
1762	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1763	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
1764	if dlg:
1765	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1766	if not GoOn:
1767	return -M #abort!!
1768	return M
1769
1770	if controls:
1771	Ftol = controls['min dM/M']
1772	else:
1773	Ftol = 0.0001
1774	if oneCycle:
1775	Ftol = 1.0
1776	x,y,w,yc,yb,yd = data #these are numpy arrays!
1777	yc *= 0. #set calcd ones to zero
1778	yb *= 0.
1779	yd *= 0.
1780	cw = x[1:]-x[:-1]
1781	xBeg = np.searchsorted(x,Limits[0])
1782	xFin = np.searchsorted(x,Limits[1])+1
1783	bakType,bakDict,bakVary = SetBackgroundParms(Background)
1784	dataType,insDict,insVary = SetInstParms(Inst)
1785	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1786	parmDict = {}
1787	parmDict.update(bakDict)
1788	parmDict.update(insDict)
1789	parmDict.update(peakDict)
1790	parmDict['Pdabc'] = [] #dummy Pdabc
1791	parmDict.update(Inst2) #put in real one if there
1792	if prevVaryList:
1793	varyList = prevVaryList[:]
1794	else:
1795	varyList = bakVary+insVary+peakVary
1796	fullvaryList = varyList[:]
1797	while True:
1798	begin = time.time()
1799	values = np.array(Dict2Values(parmDict, varyList))
1800	Rvals = {}
1801	badVary = []
1802	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1803	args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1804	ncyc = int(result[2]['nfev']/2)
1805	runtime = time.time()-begin
1806	chisq = np.sum(result[2]['fvec']**2)
1807	Values2Dict(parmDict, varyList, result[0])
1808	Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]y[xBeg:xFin]2))100. #to %
1809	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
1810	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1811	if ncyc:
1812	print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1813	print 'Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1814	sig = [0]*len(varyList)
1815	if len(varyList) == 0: break # if nothing was refined
1816	try:
1817	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1818	if np.any(np.isnan(sig)):
1819	print '* Least squares aborted - some invalid esds possible *'
1820	break #refinement succeeded - finish up!
1821	except ValueError: #result[1] is None on singular matrix
1822	print '** Refinement failed - singular matrix **'
1823	Ipvt = result[2]['ipvt']
1824	for i,ipvt in enumerate(Ipvt):
1825	if not np.sum(result[2]['fjac'],axis=1)[i]:
1826	print 'Removing parameter: ',varyList[ipvt-1]
1827	badVary.append(varyList[ipvt-1])
1828	del(varyList[ipvt-1])
1829	break
1830	else: # nothing removed
1831	break
1832	if dlg: dlg.Destroy()
1833	sigDict = dict(zip(varyList,sig))
1834	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1835	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1836	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1837	GetBackgroundParms(parmDict,Background)
1838	if bakVary: BackgroundPrint(Background,sigDict)
1839	GetInstParms(parmDict,Inst,varyList,Peaks)
1840	if insVary: InstPrint(Inst,sigDict)
1841	GetPeaksParms(Inst,parmDict,Peaks,varyList)
1842	binsperFWHM = []
1843	for peak in Peaks:
1844	FWHM = getFWHM(peak[0],Inst)
1845	try:
1846	binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1847	except IndexError:
1848	binsperFWHM.append(0.)
1849	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1850	if len(binsperFWHM):
1851	if min(binsperFWHM) < 1.:
1852	print '* Warning: calculated peak widths are too narrow to refine profile coefficients *'
1853	if 'T' in Inst['Type'][0]:
1854	print ' Manually increase sig-0, 1, or 2 in Instrument Parameters'
1855	else:
1856	print ' Manually increase W in Instrument Parameters'
1857	elif min(binsperFWHM) < 4.:
1858	print '* Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement *'
1859	print '* recommended is 6-10; you have %.2f *'%(min(binsperFWHM))
1860	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1861
1862	def calcIncident(Iparm,xdata):
1863	'needs a doc string'
1864
1865	def IfunAdv(Iparm,xdata):
1866	Itype = Iparm['Itype']
1867	Icoef = Iparm['Icoeff']
1868	DYI = np.ones((12,xdata.shape[0]))
1869	YI = np.ones_like(xdata)*Icoef[0]
1870
1871	x = xdata/1000. #expressions are in ms
1872	if Itype == 'Exponential':
1873	for i in [1,3,5,7,9]:
1874	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1875	YI += Icoef[i]*Eterm
1876	DYI[i] *= Eterm
1877	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1878	elif 'Maxwell'in Itype:
1879	Eterm = np.exp(-Icoef[2]/x**2)
1880	DYI[1] = Eterm/x**5
1881	DYI[2] = -Icoef[1]DYI[1]/x*2
1882	YI += (Icoef[1]Eterm/x*5)
1883	if 'Exponential' in Itype:
1884	for i in range(3,11,2):
1885	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1886	YI += Icoef[i]*Eterm
1887	DYI[i] *= Eterm
1888	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1889	else: #Chebyschev
1890	T = (2./x)-1.
1891	Ccof = np.ones((12,xdata.shape[0]))
1892	Ccof[1] = T
1893	for i in range(2,12):
1894	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
1895	for i in range(1,10):
1896	YI += Ccof[i]*Icoef[i+2]
1897	DYI[i+2] =Ccof[i]
1898	return YI,DYI
1899
1900	Iesd = np.array(Iparm['Iesd'])
1901	Icovar = Iparm['Icovar']
1902	YI,DYI = IfunAdv(Iparm,xdata)
1903	YI = np.where(YI>0,YI,1.)
1904	WYI = np.zeros_like(xdata)
1905	vcov = np.zeros((12,12))
1906	k = 0
1907	for i in range(12):
1908	for j in range(i,12):
1909	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
1910	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
1911	k += 1
1912	M = np.inner(vcov,DYI.T)
1913	WYI = np.sum(M*DYI,axis=0)
1914	WYI = np.where(WYI>0.,WYI,0.)
1915	return YI,WYI
1916
1917	################################################################################
1918	# Reflectometry calculations
1919	################################################################################
1920
1921	def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1922	print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
1923
1924	class RandomDisplacementBounds(object):
1925	"""random displacement with bounds"""
1926	def __init__(self, xmin, xmax, stepsize=0.5):
1927	self.xmin = xmin
1928	self.xmax = xmax
1929	self.stepsize = stepsize
1930
1931	def __call__(self, x):
1932	"""take a random step but ensure the new position is within the bounds"""
1933	while True:
1934	# this could be done in a much more clever way, but it will work for example purposes
1935	steps = self.xmax-self.xmin
1936	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
1937	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1938	break
1939	return xnew
1940
1941	def GetModelParms():
1942	parmDict = {}
1943	varyList = []
1944	values = []
1945	bounds = []
1946	parmDict['dQ type'] = data['dQ type']
1947	parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig
1948	for parm in ['Scale','FltBack']:
1949	parmDict[parm] = data[parm][0]
1950	if data[parm][1]:
1951	varyList.append(parm)
1952	values.append(data[parm][0])
1953	bounds.append(Bounds[parm])
1954	parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1955	parmDict['nLayers'] = len(parmDict['Layer Seq'])
1956	for ilay,layer in enumerate(data['Layers']):
1957	name = layer['Name']
1958	cid = str(ilay)+';'
1959	parmDict[cid+'Name'] = name
1960	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1961	parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1962	if layer.get(parm,[0.,False])[1]:
1963	varyList.append(cid+parm)
1964	value = layer[parm][0]
1965	values.append(value)
1966	if value:
1967	bound = [value*Bfac,value/Bfac]
1968	else:
1969	bound = [0.,10.]
1970	bounds.append(bound)
1971	if name not in ['vacuum','unit scatter']:
1972	parmDict[cid+'rho'] = Substances[name]['Scatt density']
1973	parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
1974	return parmDict,varyList,values,bounds
1975
1976	def SetModelParms():
1977	line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1978	if 'Scale' in varyList:
1979	data['Scale'][0] = parmDict['Scale']
1980	line += ' esd: %.4g'%(sigDict['Scale'])
1981	print line
1982	line = ' Flat background: %15.4g'%(parmDict['FltBack'])
1983	if 'FltBack' in varyList:
1984	data['FltBack'][0] = parmDict['FltBack']
1985	line += ' esd: %15.3g'%(sigDict['FltBack'])
1986	print line
1987	for ilay,layer in enumerate(data['Layers']):
1988	name = layer['Name']
1989	print ' Parameters for layer: %d %s'%(ilay,name)
1990	cid = str(ilay)+';'
1991	line = ' '
1992	line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
1993	line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
1994	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1995	if parm in layer:
1996	if parm == 'Rough':
1997	layer[parm][0] = abs(parmDict[cid+parm]) #make positive
1998	else:
1999	layer[parm][0] = parmDict[cid+parm]
2000	line += ' %s: %.3f'%(parm,layer[parm][0])
2001	if cid+parm in varyList:
2002	line += ' esd: %.3g'%(sigDict[cid+parm])
2003	print line
2004	print line2
2005
2006	def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2007	parmDict.update(zip(varyList,values))
2008	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2009	return M
2010
2011	def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2012	parmDict.update(zip(varyList,values))
2013	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2014	return np.sum(M**2)
2015
2016	def getREFD(Q,Qsig,parmDict):
2017	Ic = np.ones_like(Q)*parmDict['FltBack']
2018	Scale = parmDict['Scale']
2019	Nlayers = parmDict['nLayers']
2020	Res = parmDict['Res']
2021	depth = np.zeros(Nlayers)
2022	rho = np.zeros(Nlayers)
2023	irho = np.zeros(Nlayers)
2024	sigma = np.zeros(Nlayers)
2025	for ilay,lay in enumerate(parmDict['Layer Seq']):
2026	cid = str(lay)+';'
2027	depth[ilay] = parmDict[cid+'Thick']
2028	sigma[ilay] = parmDict[cid+'Rough']
2029	if parmDict[cid+'Name'] == u'unit scatter':
2030	rho[ilay] = parmDict[cid+'DenMul']
2031	irho[ilay] = parmDict[cid+'iDenMul']
2032	elif 'vacuum' != parmDict[cid+'Name']:
2033	rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2034	irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2035	if cid+'Mag SLD' in parmDict:
2036	rho[ilay] += parmDict[cid+'Mag SLD']
2037	if parmDict['dQ type'] == 'None':
2038	AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
2039	elif 'const' in parmDict['dQ type']:
2040	AB = SmearAbeles(0.5Q,QRes,depth,rho,irho,sigma[1:])
2041	else: #dQ/Q in data
2042	AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2043	Ic += AB*Scale
2044	return Ic
2045
2046	def estimateT0(takestep):
2047	Mmax = -1.e-10
2048	Mmin = 1.e10
2049	for i in range(100):
2050	x0 = takestep(values)
2051	M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2052	Mmin = min(M,Mmin)
2053	MMax = max(M,Mmax)
2054	return 1.5*(MMax-Mmin)
2055
2056	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2057	if data.get('2% weight'):
2058	wt = 1./(0.02Io)*2
2059	Qmin = Limits[1][0]
2060	Qmax = Limits[1][1]
2061	wtFactor = ProfDict['wtFactor']
2062	Bfac = data['Toler']
2063	Ibeg = np.searchsorted(Q,Qmin)
2064	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
2065	Ic[:] = 0
2066	Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2067	'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2068	parmDict,varyList,values,bounds = GetModelParms()
2069	Msg = 'Failed to converge'
2070	if varyList:
2071	if data['Minimizer'] == 'LMLS':
2072	result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2073	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2074	parmDict.update(zip(varyList,result[0]))
2075	chisq = np.sum(result[2]['fvec']**2)
2076	ncalc = result[2]['nfev']
2077	covM = result[1]
2078	newVals = result[0]
2079	elif data['Minimizer'] == 'Basin Hopping':
2080	xyrng = np.array(bounds).T
2081	take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2082	T0 = estimateT0(take_step)
2083	print ' Estimated temperature: %.3g'%(T0)
2084	result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2085	interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2086	'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2087	chisq = result.fun
2088	ncalc = result.nfev
2089	newVals = result.x
2090	covM = []
2091	elif data['Minimizer'] == 'MC/SA Anneal':
2092	xyrng = np.array(bounds).T
2093	result = G2mth.anneal(sumREFD, values,
2094	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2095	schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2096	boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2097	ranRange=0.20,autoRan=False,dlg=None)
2098	newVals = result[0]
2099	parmDict.update(zip(varyList,newVals))
2100	chisq = result[1]
2101	ncalc = result[3]
2102	covM = []
2103	print ' MC/SA final temperature: %.4g'%(result[2])
2104	elif data['Minimizer'] == 'L-BFGS-B':
2105	result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol,
2106	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2107	parmDict.update(zip(varyList,result['x']))
2108	chisq = result.fun
2109	ncalc = result.nfev
2110	newVals = result.x
2111	covM = []
2112	else: #nothing varied
2113	M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2114	chisq = np.sum(M**2)
2115	ncalc = 0
2116	covM = []
2117	sig = []
2118	sigDict = {}
2119	result = []
2120	Rvals = {}
2121	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]Io[Ibeg:Ifin]2))100. #to %
2122	Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2
2123	Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2124	Ib[Ibeg:Ifin] = parmDict['FltBack']
2125	try:
2126	if not len(varyList):
2127	Msg += ' - nothing refined'
2128	raise ValueError
2129	Nans = np.isnan(newVals)
2130	if np.any(Nans):
2131	name = varyList[Nans.nonzero(True)[0]]
2132	Msg += ' Nan result for '+name+'!'
2133	raise ValueError
2134	Negs = np.less_equal(newVals,0.)
2135	if np.any(Negs):
2136	indx = Negs.nonzero()
2137	name = varyList[indx[0][0]]
2138	if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2139	Msg += ' negative coefficient for '+name+'!'
2140	raise ValueError
2141	if len(covM):
2142	sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2143	covMatrix = covM*Rvals['GOF']
2144	else:
2145	sig = np.zeros(len(varyList))
2146	covMatrix = []
2147	sigDict = dict(zip(varyList,sig))
2148	print ' Results of reflectometry data modelling fit:'
2149	print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2150	print 'Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2151	SetModelParms()
2152	return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2153	except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans
2154	print Msg
2155	return False,0,0,0,0,0,0,Msg
2156
2157	def makeSLDprofile(data,Substances):
2158
2159	sq2 = np.sqrt(2.)
2160	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2161	Nlayers = len(laySeq)
2162	laySeq = np.array(laySeq,dtype=int)
2163	interfaces = np.zeros(Nlayers)
2164	rho = np.zeros(Nlayers)
2165	sigma = np.zeros(Nlayers)
2166	name = data['Layers'][0]['Name']
2167	thick = 0.
2168	for ilay,lay in enumerate(laySeq):
2169	layer = data['Layers'][lay]
2170	name = layer['Name']
2171	if 'Thick' in layer:
2172	thick += layer['Thick'][0]
2173	interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2174	if 'Rough' in layer:
2175	sigma[ilay] = max(0.001,layer['Rough'][0])
2176	if name != 'vacuum':
2177	if name == 'unit scatter':
2178	rho[ilay] = np.sqrt(layer['DenMul'][0]2+layer['iDenMul'][0]2)
2179	else:
2180	rrho = Substances[name]['Scatt density']
2181	irho = Substances[name]['XImag density']
2182	rho[ilay] = np.sqrt(rrho2+irho2)*layer['DenMul'][0]
2183	if 'Mag SLD' in layer:
2184	rho[ilay] += layer['Mag SLD'][0]
2185	name = data['Layers'][-1]['Name']
2186	x = np.linspace(-0.15thick,1.15thick,1000,endpoint=True)
2187	xr = np.flipud(x)
2188	interfaces[-1] = x[-1]
2189	y = np.ones_like(x)*rho[0]
2190	iBeg = 0
2191	for ilayer in range(Nlayers-1):
2192	delt = rho[ilayer+1]-rho[ilayer]
2193	iPos = np.searchsorted(x,interfaces[ilayer])
2194	y[iBeg:] += (delt/2.)sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2sigma[ilayer+1]))
2195	iBeg = iPos
2196	return x,xr,y
2197
2198	def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2199
2200	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2201	Qmin = Limits[1][0]
2202	Qmax = Limits[1][1]
2203	iBeg = np.searchsorted(Q,Qmin)
2204	iFin = np.searchsorted(Q,Qmax)+1 #include last point
2205	Ib[:] = data['FltBack'][0]
2206	Ic[:] = 0
2207	Scale = data['Scale'][0]
2208	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2209	Nlayers = len(laySeq)
2210	depth = np.zeros(Nlayers)
2211	rho = np.zeros(Nlayers)
2212	irho = np.zeros(Nlayers)
2213	sigma = np.zeros(Nlayers)
2214	for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2215	layer = data['Layers'][lay]
2216	name = layer['Name']
2217	if 'Thick' in layer: #skips first & last layers
2218	depth[ilay] = layer['Thick'][0]
2219	if 'Rough' in layer: #skips first layer
2220	sigma[ilay] = layer['Rough'][0]
2221	if 'unit scatter' == name:
2222	rho[ilay] = layer['DenMul'][0]
2223	irho[ilay] = layer['iDenMul'][0]
2224	else:
2225	rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2226	irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2227	if 'Mag SLD' in layer:
2228	rho[ilay] += layer['Mag SLD'][0]
2229	if data['dQ type'] == 'None':
2230	AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
2231	elif 'const' in data['dQ type']:
2232	res = data['Resolution'][0]/(100.*sateln2)
2233	AB = SmearAbeles(0.5Q[iBeg:iFin],resQ[iBeg:iFin],depth,rho,irho,sigma[1:])
2234	else: #dQ/Q in data
2235	AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2236	Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2237
2238	def abeles(kz, depth, rho, irho=0, sigma=0):
2239	"""
2240	Optical matrix form of the reflectivity calculation.
2241	O.S. Heavens, Optical Properties of Thin Solid Films
2242
2243	Reflectometry as a function of kz for a set of slabs.
2244
2245	:param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2246	This is :math:`\\tfrac12 Q_z`.
2247	:param depth: float[m] (Ang).
2248	thickness of each layer. The thickness of the incident medium
2249	and substrate are ignored.
2250	:param rho: float[n,k] (1e-6/Ang^2)
2251	Real scattering length density for each layer for each kz
2252	:param irho: float[n,k] (1e-6/Ang^2)
2253	Imaginary scattering length density for each layer for each kz
2254	Note: absorption cross section mu = 2 irho/lambda for neutrons
2255	:param sigma: float[m-1] (Ang)
2256	interfacial roughness. This is the roughness between a layer
2257	and the previous layer. The sigma array should have m-1 entries.
2258
2259	Slabs are ordered with the surface SLD at index 0 and substrate at
2260	index -1, or reversed if kz < 0.
2261	"""
2262	def calc(kz, depth, rho, irho, sigma):
2263	if len(kz) == 0: return kz
2264
2265	# Complex index of refraction is relative to the incident medium.
2266	# We can get the same effect using kz_rel^2 = kz^2 + 4pirho_o
2267	# in place of kz^2, and ignoring rho_o
2268	kz_sq = kz*2 + 4e-6np.pi*rho[:,0]
2269	k = kz
2270
2271	# According to Heavens, the initial matrix should be [ 1 F; F 1],
2272	# which we do by setting B=I and M0 to [1 F; F 1]. An extra matrix
2273	# multiply versus some coding convenience.
2274	B11 = 1
2275	B22 = 1
2276	B21 = 0
2277	B12 = 0
2278	for i in range(0, len(depth)-1):
2279	k_next = np.sqrt(kz_sq - 4e-6np.pi(rho[:,i+1] + 1j*irho[:,i+1]))
2280	F = (k - k_next) / (k + k_next)
2281	F = np.exp(-2kk_nextsigma[i]**2)
2282	#print "==== layer",i
2283	#print "kz:", kz
2284	#print "k:", k
2285	#print "k_next:",k_next
2286	#print "F:",F
2287	#print "rho:",rho[:,i+1]
2288	#print "irho:",irho[:,i+1]
2289	#print "d:",depth[i],"sigma:",sigma[i]
2290	M11 = np.exp(1jkdepth[i]) if i>0 else 1
2291	M22 = np.exp(-1jkdepth[i]) if i>0 else 1
2292	M21 = F*M11
2293	M12 = F*M22
2294	C1 = B11M11 + B21M12
2295	C2 = B11M21 + B21M22
2296	B11 = C1
2297	B21 = C2
2298	C1 = B12M11 + B22M12
2299	C2 = B12M21 + B22M22
2300	B12 = C1
2301	B22 = C2
2302	k = k_next
2303
2304	r = B12/B11
2305	return np.absolute(r)**2
2306
2307	if np.isscalar(kz): kz = np.asarray([kz], 'd')
2308
2309	m = len(depth)
2310
2311	# Make everything into arrays
2312	depth = np.asarray(depth,'d')
2313	rho = np.asarray(rho,'d')
2314	irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2315	sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2316
2317	# Repeat rho,irho columns as needed
2318	if len(rho.shape) == 1:
2319	rho = rho[None,:]
2320	irho = irho[None,:]
2321
2322	return calc(kz, depth, rho, irho, sigma)
2323
2324	def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2325	y = abeles(kz, depth, rho, irho, sigma)
2326	s = dq/2.
2327	y += 0.1354(abeles(kz+2s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2328	y += 0.24935(abeles(kz-5s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma))
2329	y += 0.4111(abeles(kz-4s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma))
2330	y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2331	y += 0.80074(abeles(kz-2s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2332	y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2333	y *= 0.137023
2334	return y
2335
2336	def makeRefdFFT(Limits,Profile):
2337	print 'make fft'
2338	Q,Io = Profile[:2]
2339	Qmin = Limits[1][0]
2340	Qmax = Limits[1][1]
2341	iBeg = np.searchsorted(Q,Qmin)
2342	iFin = np.searchsorted(Q,Qmax)+1 #include last point
2343	Qf = np.linspace(0.,Q[iFin-1],5000)
2344	QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2345	If = QI(Qf)Qf*4
2346	R = np.linspace(0.,5000.,5000)
2347	F = fft.rfft(If)
2348	return R,F
2349
2350
2351	################################################################################
2352	# Stacking fault simulation codes
2353	################################################################################
2354
2355	def GetStackParms(Layers):
2356
2357	Parms = []
2358	#cell parms
2359	if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2360	Parms.append('cellA')
2361	Parms.append('cellC')
2362	else:
2363	Parms.append('cellA')
2364	Parms.append('cellB')
2365	Parms.append('cellC')
2366	if Layers['Laue'] != 'mmm':
2367	Parms.append('cellG')
2368	#Transition parms
2369	for iY in range(len(Layers['Layers'])):
2370	for iX in range(len(Layers['Layers'])):
2371	Parms.append('TransP;%d;%d'%(iY,iX))
2372	Parms.append('TransX;%d;%d'%(iY,iX))
2373	Parms.append('TransY;%d;%d'%(iY,iX))
2374	Parms.append('TransZ;%d;%d'%(iY,iX))
2375	return Parms
2376
2377	def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2378	'''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2379
2380	:param dict Layers: dict with following items
2381
2382	::
2383
2384	{'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2385	'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2386	'Layers':[],'Stacking':[],'Transitions':[]}
2387
2388	:param str ctrls: controls string to be written on DIFFaX controls.dif file
2389	:param float scale: scale factor
2390	:param dict background: background parameters
2391	:param list limits: min/max 2-theta to be calculated
2392	:param dict inst: instrument parameters dictionary
2393	:param list profile: powder pattern data
2394
2395	Note that parameters all updated in place
2396	'''
2397	import atmdata
2398	path = sys.path
2399	for name in path:
2400	if 'bin' in name:
2401	DIFFaX = name+'/DIFFaX.exe'
2402	print ' Execute ',DIFFaX
2403	break
2404	# make form factor file that DIFFaX wants - atom types are GSASII style
2405	sf = open('data.sfc','w')
2406	sf.write('GSASII special form factor file for DIFFaX\n\n')
2407	atTypes = Layers['AtInfo'].keys()
2408	if 'H' not in atTypes:
2409	atTypes.insert(0,'H')
2410	for atType in atTypes:
2411	if atType == 'H':
2412	blen = -.3741
2413	else:
2414	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2415	Adat = atmdata.XrayFF[atType]
2416	text = '%4s'%(atType.ljust(4))
2417	for i in range(4):
2418	text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2419	text += '%11.6f%11.6f'%(Adat['fc'],blen)
2420	text += '%3d\n'%(Adat['Z'])
2421	sf.write(text)
2422	sf.close()
2423	#make DIFFaX control.dif file - future use GUI to set some of these flags
2424	cf = open('control.dif','w')
2425	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
2426	x0 = profile[0]
2427	iBeg = np.searchsorted(x0,limits[0])
2428	iFin = np.searchsorted(x0,limits[1])+1
2429	if iFin-iBeg > 20000:
2430	iFin = iBeg+20000
2431	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2432	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2433	cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2434	else:
2435	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2436	inst = {'Type':['XSC','XSC',]}
2437	cf.close()
2438	#make DIFFaX data file
2439	df = open('GSASII-DIFFaX.dat','w')
2440	df.write('INSTRUMENTAL\n')
2441	if 'X' in inst['Type'][0]:
2442	df.write('X-RAY\n')
2443	elif 'N' in inst['Type'][0]:
2444	df.write('NEUTRON\n')
2445	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
2446	df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2447	U = ateln2*inst['U'][1]/10000.
2448	V = ateln2*inst['V'][1]/10000.
2449	W = ateln2*inst['W'][1]/10000.
2450	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
2451	HW = np.sqrt(np.mean(HWHM))
2452	# df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2453	if 'Mean' in Layers['selInst']:
2454	df.write('GAUSSIAN %.6f TRIM\n'%(HW)) #fast option - might not really matter
2455	elif 'Gaussian' in Layers['selInst']:
2456	df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W)) #slow - make a GUI option?
2457	else:
2458	df.write('None\n')
2459	else:
2460	df.write('0.10\nNone\n')
2461	df.write('STRUCTURAL\n')
2462	a,b,c = Layers['Cell'][1:4]
2463	gam = Layers['Cell'][6]
2464	df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2465	laue = Layers['Laue']
2466	if laue == '2/m(ab)':
2467	laue = '2/m(1)'
2468	elif laue == '2/m(c)':
2469	laue = '2/m(2)'
2470	if 'unknown' in Layers['Laue']:
2471	df.write('%s %.3f\n'%(laue,Layers['Toler']))
2472	else:
2473	df.write('%s\n'%(laue))
2474	df.write('%d\n'%(len(Layers['Layers'])))
2475	if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2476	df.write('%.1f %.1f\n'%(Layers['Width'][0][0]10000.,Layers['Width'][0][0]10000.)) #mum to A
2477	layerNames = []
2478	for layer in Layers['Layers']:
2479	layerNames.append(layer['Name'])
2480	for il,layer in enumerate(Layers['Layers']):
2481	if layer['SameAs']:
2482	df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2483	continue
2484	df.write('LAYER %d\n'%(il+1))
2485	if '-1' in layer['Symm']:
2486	df.write('CENTROSYMMETRIC\n')
2487	else:
2488	df.write('NONE\n')
2489	for ia,atom in enumerate(layer['Atoms']):
2490	[name,atype,x,y,z,frac,Uiso] = atom
2491	if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2492	frac /= 2.
2493	df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2494	df.write('STACKING\n')
2495	df.write('%s\n'%(Layers['Stacking'][0]))
2496	if 'recursive' in Layers['Stacking'][0]:
2497	df.write('%s\n'%Layers['Stacking'][1])
2498	else:
2499	if 'list' in Layers['Stacking'][1]:
2500	Slen = len(Layers['Stacking'][2])
2501	iB = 0
2502	iF = 0
2503	while True:
2504	iF += 68
2505	if iF >= Slen:
2506	break
2507	iF = min(iF,Slen)
2508	df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2509	iB = iF
2510	else:
2511	df.write('%s\n'%Layers['Stacking'][1])
2512	df.write('TRANSITIONS\n')
2513	for iY in range(len(Layers['Layers'])):
2514	sumPx = 0.
2515	for iX in range(len(Layers['Layers'])):
2516	p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2517	p = round(p,3)
2518	df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2519	sumPx += p
2520	if sumPx != 1.0: #this has to be picky since DIFFaX is.
2521	print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
2522	df.close()
2523	os.remove('data.sfc')
2524	os.remove('control.dif')
2525	os.remove('GSASII-DIFFaX.dat')
2526	return
2527	df.close()
2528	time0 = time.time()
2529	try:
2530	subp.call(DIFFaX)
2531	except OSError:
2532	print ' DIFFax.exe is not available for this platform - under development'
2533	print ' DIFFaX time = %.2fs'%(time.time()-time0)
2534	if os.path.exists('GSASII-DIFFaX.spc'):
2535	Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2536	iFin = iBeg+Xpat.shape[1]
2537	bakType,backDict,backVary = SetBackgroundParms(background)
2538	backDict['Lam1'] = G2mth.getWave(inst)
2539	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
2540	profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2541	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
2542	rv = st.poisson(profile[3][iBeg:iFin])
2543	profile[1][iBeg:iFin] = rv.rvs()
2544	Z = np.ones_like(profile[3][iBeg:iFin])
2545	Z[1::2] *= -1
2546	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2547	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2548	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2549	#cleanup files..
2550	os.remove('GSASII-DIFFaX.spc')
2551	elif os.path.exists('GSASII-DIFFaX.sadp'):
2552	Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2553	Sadp = np.reshape(Sadp,(256,-1))
2554	Layers['Sadp']['Img'] = Sadp
2555	os.remove('GSASII-DIFFaX.sadp')
2556	os.remove('data.sfc')
2557	os.remove('control.dif')
2558	os.remove('GSASII-DIFFaX.dat')
2559
2560	def SetPWDRscan(inst,limits,profile):
2561
2562	wave = G2mth.getMeanWave(inst)
2563	x0 = profile[0]
2564	iBeg = np.searchsorted(x0,limits[0])
2565	iFin = np.searchsorted(x0,limits[1])
2566	if iFin-iBeg > 20000:
2567	iFin = iBeg+20000
2568	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2569	pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2570	return iFin-iBeg
2571
2572	def SetStackingSF(Layers,debug):
2573	# Load scattering factors into DIFFaX arrays
2574	import atmdata
2575	atTypes = Layers['AtInfo'].keys()
2576	aTypes = []
2577	for atype in atTypes:
2578	aTypes.append('%4s'%(atype.ljust(4)))
2579	SFdat = []
2580	for atType in atTypes:
2581	Adat = atmdata.XrayFF[atType]
2582	SF = np.zeros(9)
2583	SF[:8:2] = Adat['fa']
2584	SF[1:8:2] = Adat['fb']
2585	SF[8] = Adat['fc']
2586	SFdat.append(SF)
2587	SFdat = np.array(SFdat)
2588	pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2589
2590	def SetStackingClay(Layers,Type):
2591	# Controls
2592	rand.seed()
2593	ranSeed = rand.randint(1,2**16-1)
2594	try:
2595	laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2596	'6/m','6/mmm'].index(Layers['Laue'])+1
2597	except ValueError: #for 'unknown'
2598	laueId = -1
2599	if 'SADP' in Type:
2600	planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2601	lmax = int(Layers['Sadp']['Lmax'])
2602	else:
2603	planeId = 0
2604	lmax = 0
2605	# Sequences
2606	StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2607	try:
2608	StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2609	except ValueError:
2610	StkParm = -1
2611	if StkParm == 2: #list
2612	StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2613	Nstk = len(StkSeq)
2614	else:
2615	Nstk = 1
2616	StkSeq = [0,]
2617	if StkParm == -1:
2618	StkParm = int(Layers['Stacking'][1])
2619	Wdth = Layers['Width'][0]
2620	mult = 1
2621	controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2622	LaueSym = Layers['Laue'].ljust(12)
2623	pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2624	return laueId,controls
2625
2626	def SetCellAtoms(Layers):
2627	Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2628	# atoms in layers
2629	atTypes = Layers['AtInfo'].keys()
2630	AtomXOU = []
2631	AtomTp = []
2632	LayerSymm = []
2633	LayerNum = []
2634	layerNames = []
2635	Natm = 0
2636	Nuniq = 0
2637	for layer in Layers['Layers']:
2638	layerNames.append(layer['Name'])
2639	for il,layer in enumerate(Layers['Layers']):
2640	if layer['SameAs']:
2641	LayerNum.append(layerNames.index(layer['SameAs'])+1)
2642	continue
2643	else:
2644	LayerNum.append(il+1)
2645	Nuniq += 1
2646	if '-1' in layer['Symm']:
2647	LayerSymm.append(1)
2648	else:
2649	LayerSymm.append(0)
2650	for ia,atom in enumerate(layer['Atoms']):
2651	[name,atype,x,y,z,frac,Uiso] = atom
2652	Natm += 1
2653	AtomTp.append('%4s'%(atype.ljust(4)))
2654	Ta = atTypes.index(atype)+1
2655	AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2656	AtomXOU = np.array(AtomXOU)
2657	Nlayers = len(layerNames)
2658	pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2659	return Nlayers
2660
2661	def SetStackingTrans(Layers,Nlayers):
2662	# Transitions
2663	TransX = []
2664	TransP = []
2665	for Ytrans in Layers['Transitions']:
2666	TransP.append([trans[0] for trans in Ytrans]) #get just the numbers
2667	TransX.append([trans[1:4] for trans in Ytrans]) #get just the numbers
2668	TransP = np.array(TransP,dtype='float').T
2669	TransX = np.array(TransX,dtype='float')
2670	# GSASIIpath.IPyBreak()
2671	pyx.pygettrans(Nlayers,TransP,TransX)
2672
2673	def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2674	# Scattering factors
2675	SetStackingSF(Layers,debug)
2676	# Controls & sequences
2677	laueId,controls = SetStackingClay(Layers,'PWDR')
2678	# cell & atoms
2679	Nlayers = SetCellAtoms(Layers)
2680	Volume = Layers['Cell'][7]
2681	# Transitions
2682	SetStackingTrans(Layers,Nlayers)
2683	# PWDR scan
2684	Nsteps = SetPWDRscan(inst,limits,profile)
2685	# result as Spec
2686	x0 = profile[0]
2687	profile[3] = np.zeros(len(profile[0]))
2688	profile[4] = np.zeros(len(profile[0]))
2689	profile[5] = np.zeros(len(profile[0]))
2690	iBeg = np.searchsorted(x0,limits[0])
2691	iFin = np.searchsorted(x0,limits[1])+1
2692	if iFin-iBeg > 20000:
2693	iFin = iBeg+20000
2694	Nspec = 20001
2695	spec = np.zeros(Nspec,dtype='double')
2696	time0 = time.time()
2697	pyx.pygetspc(controls,Nspec,spec)
2698	print ' GETSPC time = %.2fs'%(time.time()-time0)
2699	time0 = time.time()
2700	U = ateln2*inst['U'][1]/10000.
2701	V = ateln2*inst['V'][1]/10000.
2702	W = ateln2*inst['W'][1]/10000.
2703	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
2704	HW = np.sqrt(np.mean(HWHM))
2705	BrdSpec = np.zeros(Nsteps)
2706	if 'Mean' in Layers['selInst']:
2707	pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2708	elif 'Gaussian' in Layers['selInst']:
2709	pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2710	else:
2711	BrdSpec = spec[:Nsteps]
2712	BrdSpec /= Volume
2713	iFin = iBeg+Nsteps
2714	bakType,backDict,backVary = SetBackgroundParms(background)
2715	backDict['Lam1'] = G2mth.getWave(inst)
2716	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
2717	profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2718	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
2719	try:
2720	rv = st.poisson(profile[3][iBeg:iFin])
2721	profile[1][iBeg:iFin] = rv.rvs()
2722	except ValueError:
2723	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2724	Z = np.ones_like(profile[3][iBeg:iFin])
2725	Z[1::2] *= -1
2726	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2727	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2728	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2729	print ' Broadening time = %.2fs'%(time.time()-time0)
2730
2731	def CalcStackingSADP(Layers,debug):
2732
2733	# Scattering factors
2734	SetStackingSF(Layers,debug)
2735	# Controls & sequences
2736	laueId,controls = SetStackingClay(Layers,'SADP')
2737	# cell & atoms
2738	Nlayers = SetCellAtoms(Layers)
2739	# Transitions
2740	SetStackingTrans(Layers,Nlayers)
2741	# result as Sadp
2742	Nspec = 20001
2743	spec = np.zeros(Nspec,dtype='double')
2744	time0 = time.time()
2745	hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2746	Sapd = np.zeros((256,256))
2747	iB = 0
2748	for i in range(hkLim):
2749	iF = iB+Nblk
2750	p1 = 127+int(i*Incr)
2751	p2 = 128-int(i*Incr)
2752	if Nblk == 128:
2753	if i:
2754	Sapd[128:,p1] = spec[iB:iF]
2755	Sapd[:128,p1] = spec[iF:iB:-1]
2756	Sapd[128:,p2] = spec[iB:iF]
2757	Sapd[:128,p2] = spec[iF:iB:-1]
2758	else:
2759	if i:
2760	Sapd[:,p1] = spec[iB:iF]
2761	Sapd[:,p2] = spec[iB:iF]
2762	iB += Nblk
2763	Layers['Sadp']['Img'] = Sapd
2764	print ' GETSAD time = %.2fs'%(time.time()-time0)
2765	# GSASIIpath.IPyBreak()
2766
2767	#testing data
2768	NeedTestData = True
2769	def TestData():
2770	'needs a doc string'
2771	# global NeedTestData
2772	global bakType
2773	bakType = 'chebyschev'
2774	global xdata
2775	xdata = np.linspace(4.0,40.0,36000)
2776	global parmDict0
2777	parmDict0 = {
2778	'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2779	'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2780	'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2781	'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2782	'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2783	'Back0':5.384,'Back1':-0.015,'Back2':.004,
2784	}
2785	global parmDict1
2786	parmDict1 = {
2787	'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2788	'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2789	'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2790	'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2791	'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2792	'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2793	'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2794	'Back0':36.897,'Back1':-0.508,'Back2':.006,
2795	'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2796	}
2797	global parmDict2
2798	parmDict2 = {
2799	'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2800	'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2801	'Back0':5.,'Back1':-0.02,'Back2':.004,
2802	# 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2803	}
2804	global varyList
2805	varyList = []
2806
2807	def test0():
2808	if NeedTestData: TestData()
2809	gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2810	gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])
2811	gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2812	fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2813	fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])
2814	fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2815
2816	def test1():
2817	if NeedTestData: TestData()
2818	time0 = time.time()
2819	for i in range(100):
2820	getPeakProfile(parmDict1,xdata,varyList,bakType)
2821	print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2822
2823	def test2(name,delt):
2824	if NeedTestData: TestData()
2825	varyList = [name,]
2826	xdata = np.linspace(5.6,5.8,400)
2827	hplot = plotter.add('derivatives test for '+name).gca()
2828	hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2829	y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2830	parmDict2[name] += delt
2831	y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2832	hplot.plot(xdata,(y1-y0)/delt,'r+')
2833
2834	def test3(name,delt):
2835	if NeedTestData: TestData()
2836	names = ['pos','sig','gam','shl']
2837	idx = names.index(name)
2838	myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2839	xdata = np.linspace(5.6,5.8,800)
2840	dx = xdata[1]-xdata[0]
2841	hplot = plotter.add('derivatives test for '+name).gca()
2842	hplot.plot(xdata,100.dxgetdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2843	y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2844	myDict[name] += delt
2845	y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2846	hplot.plot(xdata,(y1-y0)/delt,'r+')
2847
2848	if __name__ == '__main__':
2849	import GSASIItestplot as plot
2850	global plotter
2851	plotter = plot.PlotNotebook()
2852	# test0()
2853	# for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2854	for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2855	['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2856	test2(name,shft)
2857	for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2858	test3(name,shft)
2859	print "OK"
2860	plotter.StartEventLoop()

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

Download in other formats: