Context Navigation

source: trunk/GSASIIpwd.py @ 2849

Last change on this file since 2849 was 2849, checked in by vondreele, 6 years ago
change TOF instrument Calibration weighting
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 111.2 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-06-01 14:44:40 +0000 (Thu, 01 Jun 2017) $
10	# $Author: vondreele $
11	# $Revision: 2849 $
12	# $URL: trunk/GSASIIpwd.py $
13	# $Id: GSASIIpwd.py 2849 2017-06-01 14:44:40Z vondreele $
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: 2849 $")
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	'needs a doc string'
684	DX = xdata[1]-xdata[0]
685	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
686	x = np.linspace(pos-fmin,pos+fmin,256)
687	dx = x[1]-x[0]
688	Norm = norm.pdf(x,loc=pos,scale=widths[0])
689	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
690	arg = [pos,shl/57.2958,dx,]
691	FCJ = fcjde.pdf(x,*arg,loc=pos)
692	if len(np.nonzero(FCJ)[0])>5:
693	z = np.column_stack([Norm,Cauchy,FCJ]).T
694	Z = fft.fft(z)
695	Df = fft.ifft(Z.prod(axis=0)).real
696	else:
697	z = np.column_stack([Norm,Cauchy]).T
698	Z = fft.fft(z)
699	Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
700	Df /= np.sum(Df)
701	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
702	return intensDf(xdata)DX/dx
703
704	def getBackground(pfx,parmDict,bakType,dataType,xdata):
705	'needs a doc string'
706	if 'T' in dataType:
707	q = 2.np.piparmDict[pfx+'difC']/xdata
708	elif 'C' in dataType:
709	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
710	q = npT2q(xdata,wave)
711	yb = np.zeros_like(xdata)
712	nBak = 0
713	cw = np.diff(xdata)
714	cw = np.append(cw,cw[-1])
715	sumBk = [0.,0.,0]
716	while True:
717	key = pfx+'Back;'+str(nBak)
718	if key in parmDict:
719	nBak += 1
720	else:
721	break
722	#empirical functions
723	if bakType in ['chebyschev','cosine']:
724	dt = xdata[-1]-xdata[0]
725	for iBak in range(nBak):
726	key = pfx+'Back;'+str(iBak)
727	if bakType == 'chebyschev':
728	ybi = parmDict[key](-1.+2.(xdata-xdata[0])/dt)**iBak
729	elif bakType == 'cosine':
730	ybi = parmDict[key]npcosd(180.xdata*iBak/xdata[-1])
731	yb += ybi
732	sumBk[0] = np.sum(yb)
733	elif bakType in ['Q^2 power series','Q^-2 power series']:
734	QT = 1.
735	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
736	for iBak in range(nBak-1):
737	key = pfx+'Back;'+str(iBak+1)
738	if '-2' in bakType:
739	QT = (iBak+1)q**-2
740	else:
741	QT = q*2/(iBak+1)
742	yb += QT*parmDict[key]
743	sumBk[0] = np.sum(yb)
744	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
745	if nBak == 1:
746	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
747	elif nBak == 2:
748	dX = xdata[-1]-xdata[0]
749	T2 = (xdata-xdata[0])/dX
750	T1 = 1.0-T2
751	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
752	else:
753	if bakType == 'lin interpolate':
754	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
755	elif bakType == 'inv interpolate':
756	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
757	elif bakType == 'log interpolate':
758	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
759	bakPos[0] = xdata[0]
760	bakPos[-1] = xdata[-1]
761	bakVals = np.zeros(nBak)
762	for i in range(nBak):
763	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
764	bakInt = si.interp1d(bakPos,bakVals,'linear')
765	yb = bakInt(ma.getdata(xdata))
766	sumBk[0] = np.sum(yb)
767	#Debye function
768	if pfx+'difC' in parmDict:
769	ff = 1.
770	else:
771	try:
772	wave = parmDict[pfx+'Lam']
773	except KeyError:
774	wave = parmDict[pfx+'Lam1']
775	SQ = (q/(4.np.pi))*2
776	FF = G2elem.GetFormFactorCoeff('Si')[0]
777	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
778	iD = 0
779	while True:
780	try:
781	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
782	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
783	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
784	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
785	yb += ybi
786	sumBk[1] += np.sum(ybi)
787	iD += 1
788	except KeyError:
789	break
790	#peaks
791	iD = 0
792	while True:
793	try:
794	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
795	pkI = parmDict[pfx+'BkPkint;'+str(iD)]
796	pkS = parmDict[pfx+'BkPksig;'+str(iD)]
797	pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
798	if 'C' in dataType:
799	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
800	else: #'T'OF
801	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
802	iBeg = np.searchsorted(xdata,pkP-fmin)
803	iFin = np.searchsorted(xdata,pkP+fmax)
804	lenX = len(xdata)
805	if not iBeg:
806	iFin = np.searchsorted(xdata,pkP+fmax)
807	elif iBeg == lenX:
808	iFin = iBeg
809	else:
810	iFin = np.searchsorted(xdata,pkP+fmax)
811	if 'C' in dataType:
812	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
813	yb[iBeg:iFin] += ybi
814	else: #'T'OF
815	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
816	yb[iBeg:iFin] += ybi
817	sumBk[2] += np.sum(ybi)
818	iD += 1
819	except KeyError:
820	break
821	except ValueError:
822	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
823	break
824	return yb,sumBk
825
826	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
827	'needs a doc string'
828	if 'T' in dataType:
829	q = 2.np.piparmDict[hfx+'difC']/xdata
830	elif 'C' in dataType:
831	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
832	q = 2.np.pinpsind(xdata/2.)/wave
833	nBak = 0
834	while True:
835	key = hfx+'Back;'+str(nBak)
836	if key in parmDict:
837	nBak += 1
838	else:
839	break
840	dydb = np.zeros(shape=(nBak,len(xdata)))
841	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
842	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
843	cw = np.diff(xdata)
844	cw = np.append(cw,cw[-1])
845
846	if bakType in ['chebyschev','cosine']:
847	dt = xdata[-1]-xdata[0]
848	for iBak in range(nBak):
849	if bakType == 'chebyschev':
850	dydb[iBak] = (-1.+2.(xdata-xdata[0])/dt)*iBak
851	elif bakType == 'cosine':
852	dydb[iBak] = npcosd(180.xdataiBak/xdata[-1])
853	elif bakType in ['Q^2 power series','Q^-2 power series']:
854	QT = 1.
855	dydb[0] = np.ones_like(xdata)
856	for iBak in range(nBak-1):
857	if '-2' in bakType:
858	QT = (iBak+1)q**-2
859	else:
860	QT = q*2/(iBak+1)
861	dydb[iBak+1] = QT
862	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
863	if nBak == 1:
864	dydb[0] = np.ones_like(xdata)
865	elif nBak == 2:
866	dX = xdata[-1]-xdata[0]
867	T2 = (xdata-xdata[0])/dX
868	T1 = 1.0-T2
869	dydb = [T1,T2]
870	else:
871	if bakType == 'lin interpolate':
872	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
873	elif bakType == 'inv interpolate':
874	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
875	elif bakType == 'log interpolate':
876	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
877	bakPos[0] = xdata[0]
878	bakPos[-1] = xdata[-1]
879	for i,pos in enumerate(bakPos):
880	if i == 0:
881	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
882	elif i == len(bakPos)-1:
883	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
884	else:
885	dydb[i] = np.where(xdata>bakPos[i],
886	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
887	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
888	if hfx+'difC' in parmDict:
889	ff = 1.
890	else:
891	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
892	q = npT2q(xdata,wave)
893	SQ = (q/(4np.pi))*2
894	FF = G2elem.GetFormFactorCoeff('Si')[0]
895	ff = np.array(G2elem.ScatFac(FF,SQ)[0])np.pi*2 #needs pi^2~10. for cw data (why?)
896	iD = 0
897	while True:
898	try:
899	if hfx+'difC' in parmDict:
900	q = 2np.piparmDict[hfx+'difC']/xdata
901	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
902	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
903	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
904	sqr = np.sin(qdbR)/(qdbR)
905	cqr = np.cos(q*dbR)
906	temp = np.exp(-dbUq*2)
907	dyddb[3iD] = ffsqr*temp
908	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(dbR)
909	dyddb[3iD+2] = -ffdbAsqrtempq*2
910	iD += 1
911	except KeyError:
912	break
913	iD = 0
914	while True:
915	try:
916	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
917	pkI = parmDict[hfx+'BkPkint;'+str(iD)]
918	pkS = parmDict[hfx+'BkPksig;'+str(iD)]
919	pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
920	if 'C' in dataType:
921	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
922	else: #'T'OF
923	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
924	iBeg = np.searchsorted(xdata,pkP-fmin)
925	iFin = np.searchsorted(xdata,pkP+fmax)
926	lenX = len(xdata)
927	if not iBeg:
928	iFin = np.searchsorted(xdata,pkP+fmax)
929	elif iBeg == lenX:
930	iFin = iBeg
931	else:
932	iFin = np.searchsorted(xdata,pkP+fmax)
933	if 'C' in dataType:
934	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
935	dydpk[4iD][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdp
936	dydpk[4iD+1][iBeg:iFin] += 100.cw[iBeg:iFin]*Df
937	dydpk[4iD+2][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFds
938	dydpk[4iD+3][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdg
939	else: #'T'OF
940	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
941	dydpk[4iD][iBeg:iFin] += pkIdFdp
942	dydpk[4*iD+1][iBeg:iFin] += Df
943	dydpk[4iD+2][iBeg:iFin] += pkIdFds
944	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
945	iD += 1
946	except KeyError:
947	break
948	except ValueError:
949	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
950	break
951	return dydb,dyddb,dydpk
952
953	#use old fortran routine
954	def getFCJVoigt3(pos,sig,gam,shl,xdata):
955	'needs a doc string'
956	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
957	# Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
958	Df /= np.sum(Df)
959	return Df
960
961	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
962	'needs a doc string'
963
964	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
965	# Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
966	return Df,dFdp,dFds,dFdg,dFdsh
967
968	def getPsVoigt(pos,sig,gam,xdata):
969	'needs a doc string'
970
971	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
972	Df /= np.sum(Df)
973	return Df
974
975	def getdPsVoigt(pos,sig,gam,xdata):
976	'needs a doc string'
977
978	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
979	return Df,dFdp,dFds,dFdg
980
981	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
982	'needs a doc string'
983	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
984	Df /= np.sum(Df)
985	return Df
986
987	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
988	'needs a doc string'
989	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
990	return Df,dFdp,dFda,dFdb,dFds,dFdg
991
992	def ellipseSize(H,Sij,GB):
993	'needs a doc string'
994	HX = np.inner(H.T,GB)
995	lenHX = np.sqrt(np.sum(HX**2))
996	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
997	R = np.inner(HX/lenHX,Rsize)*Esize #want column length for hkl in crystal
998	lenR = np.sqrt(np.sum(R**2))
999	return lenR
1000
1001	def ellipseSizeDerv(H,Sij,GB):
1002	'needs a doc string'
1003	lenR = ellipseSize(H,Sij,GB)
1004	delt = 0.001
1005	dRdS = np.zeros(6)
1006	for i in range(6):
1007	Sij[i] -= delt
1008	lenM = ellipseSize(H,Sij,GB)
1009	Sij[i] += 2.*delt
1010	lenP = ellipseSize(H,Sij,GB)
1011	Sij[i] -= delt
1012	dRdS[i] = (lenP-lenM)/(2.*delt)
1013	return lenR,dRdS
1014
1015	def getHKLpeak(dmin,SGData,A,Inst=None):
1016	'needs a doc string'
1017	HKL = G2lat.GenHLaue(dmin,SGData,A)
1018	HKLs = []
1019	for h,k,l,d in HKL:
1020	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1021	if not ext:
1022	if Inst == None:
1023	HKLs.append([h,k,l,d,0,-1])
1024	else:
1025	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1026	return HKLs
1027
1028	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1029	'needs a doc string'
1030	HKLs = []
1031	vec = np.array(Vec)
1032	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1033	dvec = 1./(maxH*vstar+1./dmin)
1034	HKL = G2lat.GenHLaue(dvec,SGData,A)
1035	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1036	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1037	for h,k,l,d in HKL:
1038	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1039	if not ext and d >= dmin:
1040	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1041	for dH in SSdH:
1042	if dH:
1043	DH = SSdH[dH]
1044	H = [h+DH[0],k+DH[1],l+DH[2]]
1045	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1046	if d >= dmin:
1047	HKLM = np.array([h,k,l,dH])
1048	if G2spc.checkSSextc(HKLM,SSGData):
1049	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
1050	# GSASIIpath.IPyBreak()
1051	return G2lat.sortHKLd(HKLs,True,True,True)
1052
1053	def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1054	'needs a doc string'
1055
1056	yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1057	yc = np.zeros_like(yb)
1058	cw = np.diff(xdata)
1059	cw = np.append(cw,cw[-1])
1060	if 'C' in dataType:
1061	shl = max(parmDict['SH/L'],0.002)
1062	Ka2 = False
1063	if 'Lam1' in parmDict.keys():
1064	Ka2 = True
1065	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1066	kRatio = parmDict['I(L2)/I(L1)']
1067	iPeak = 0
1068	while True:
1069	try:
1070	pos = parmDict['pos'+str(iPeak)]
1071	tth = (pos-parmDict['Zero'])
1072	intens = parmDict['int'+str(iPeak)]
1073	sigName = 'sig'+str(iPeak)
1074	if sigName in varyList:
1075	sig = parmDict[sigName]
1076	else:
1077	sig = G2mth.getCWsig(parmDict,tth)
1078	sig = max(sig,0.001) #avoid neg sigma^2
1079	gamName = 'gam'+str(iPeak)
1080	if gamName in varyList:
1081	gam = parmDict[gamName]
1082	else:
1083	gam = G2mth.getCWgam(parmDict,tth)
1084	gam = max(gam,0.001) #avoid neg gamma
1085	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1086	iBeg = np.searchsorted(xdata,pos-fmin)
1087	iFin = np.searchsorted(xdata,pos+fmin)
1088	if not iBeg+iFin: #peak below low limit
1089	iPeak += 1
1090	continue
1091	elif not iBeg-iFin: #peak above high limit
1092	return yb+yc
1093	yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1094	if Ka2:
1095	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1096	iBeg = np.searchsorted(xdata,pos2-fmin)
1097	iFin = np.searchsorted(xdata,pos2+fmin)
1098	if iBeg-iFin:
1099	yc[iBeg:iFin] += intenskRatiogetFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1100	iPeak += 1
1101	except KeyError: #no more peaks to process
1102	return yb+yc
1103	else:
1104	Pdabc = parmDict['Pdabc']
1105	difC = parmDict['difC']
1106	iPeak = 0
1107	while True:
1108	try:
1109	pos = parmDict['pos'+str(iPeak)]
1110	tof = pos-parmDict['Zero']
1111	dsp = tof/difC
1112	intens = parmDict['int'+str(iPeak)]
1113	alpName = 'alp'+str(iPeak)
1114	if alpName in varyList:
1115	alp = parmDict[alpName]
1116	else:
1117	if len(Pdabc):
1118	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1119	else:
1120	alp = G2mth.getTOFalpha(parmDict,dsp)
1121	alp = max(0.0001,alp)
1122	betName = 'bet'+str(iPeak)
1123	if betName in varyList:
1124	bet = parmDict[betName]
1125	else:
1126	if len(Pdabc):
1127	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1128	else:
1129	bet = G2mth.getTOFbeta(parmDict,dsp)
1130	bet = max(0.0001,bet)
1131	sigName = 'sig'+str(iPeak)
1132	if sigName in varyList:
1133	sig = parmDict[sigName]
1134	else:
1135	sig = G2mth.getTOFsig(parmDict,dsp)
1136	gamName = 'gam'+str(iPeak)
1137	if gamName in varyList:
1138	gam = parmDict[gamName]
1139	else:
1140	gam = G2mth.getTOFgamma(parmDict,dsp)
1141	gam = max(gam,0.001) #avoid neg gamma
1142	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1143	iBeg = np.searchsorted(xdata,pos-fmin)
1144	iFin = np.searchsorted(xdata,pos+fmax)
1145	lenX = len(xdata)
1146	if not iBeg:
1147	iFin = np.searchsorted(xdata,pos+fmax)
1148	elif iBeg == lenX:
1149	iFin = iBeg
1150	else:
1151	iFin = np.searchsorted(xdata,pos+fmax)
1152	if not iBeg+iFin: #peak below low limit
1153	iPeak += 1
1154	continue
1155	elif not iBeg-iFin: #peak above high limit
1156	return yb+yc
1157	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1158	iPeak += 1
1159	except KeyError: #no more peaks to process
1160	return yb+yc
1161
1162	def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1163	'needs a doc string'
1164	# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1165	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1166	dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1167	if 'Back;0' in varyList: #background derivs are in front if present
1168	dMdv[0:len(dMdb)] = dMdb
1169	names = ['DebyeA','DebyeR','DebyeU']
1170	for name in varyList:
1171	if 'Debye' in name:
1172	parm,id = name.split(';')
1173	ip = names.index(parm)
1174	dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1175	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1176	for name in varyList:
1177	if 'BkPk' in name:
1178	parm,id = name.split(';')
1179	ip = names.index(parm)
1180	dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1181	cw = np.diff(xdata)
1182	cw = np.append(cw,cw[-1])
1183	if 'C' in dataType:
1184	shl = max(parmDict['SH/L'],0.002)
1185	Ka2 = False
1186	if 'Lam1' in parmDict.keys():
1187	Ka2 = True
1188	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1189	kRatio = parmDict['I(L2)/I(L1)']
1190	iPeak = 0
1191	while True:
1192	try:
1193	pos = parmDict['pos'+str(iPeak)]
1194	tth = (pos-parmDict['Zero'])
1195	intens = parmDict['int'+str(iPeak)]
1196	sigName = 'sig'+str(iPeak)
1197	if sigName in varyList:
1198	sig = parmDict[sigName]
1199	dsdU = dsdV = dsdW = 0
1200	else:
1201	sig = G2mth.getCWsig(parmDict,tth)
1202	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1203	sig = max(sig,0.001) #avoid neg sigma
1204	gamName = 'gam'+str(iPeak)
1205	if gamName in varyList:
1206	gam = parmDict[gamName]
1207	dgdX = dgdY = 0
1208	else:
1209	gam = G2mth.getCWgam(parmDict,tth)
1210	dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1211	gam = max(gam,0.001) #avoid neg gamma
1212	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1213	iBeg = np.searchsorted(xdata,pos-fmin)
1214	iFin = np.searchsorted(xdata,pos+fmin)
1215	if not iBeg+iFin: #peak below low limit
1216	iPeak += 1
1217	continue
1218	elif not iBeg-iFin: #peak above high limit
1219	break
1220	dMdpk = np.zeros(shape=(6,len(xdata)))
1221	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1222	for i in range(1,5):
1223	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intens*dMdipk[i]
1224	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk[0]
1225	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1226	if Ka2:
1227	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1228	iBeg = np.searchsorted(xdata,pos2-fmin)
1229	iFin = np.searchsorted(xdata,pos2+fmin)
1230	if iBeg-iFin:
1231	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1232	for i in range(1,5):
1233	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intenskRatiodMdipk2[i]
1234	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]kRatio*dMdipk2[0]
1235	dMdpk[5][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk2[0]
1236	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1237	for parmName in ['pos','int','sig','gam']:
1238	try:
1239	idx = varyList.index(parmName+str(iPeak))
1240	dMdv[idx] = dervDict[parmName]
1241	except ValueError:
1242	pass
1243	if 'U' in varyList:
1244	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1245	if 'V' in varyList:
1246	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1247	if 'W' in varyList:
1248	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1249	if 'X' in varyList:
1250	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1251	if 'Y' in varyList:
1252	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1253	if 'SH/L' in varyList:
1254	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1255	if 'I(L2)/I(L1)' in varyList:
1256	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1257	iPeak += 1
1258	except KeyError: #no more peaks to process
1259	break
1260	else:
1261	Pdabc = parmDict['Pdabc']
1262	difC = parmDict['difC']
1263	iPeak = 0
1264	while True:
1265	try:
1266	pos = parmDict['pos'+str(iPeak)]
1267	tof = pos-parmDict['Zero']
1268	dsp = tof/difC
1269	intens = parmDict['int'+str(iPeak)]
1270	alpName = 'alp'+str(iPeak)
1271	if alpName in varyList:
1272	alp = parmDict[alpName]
1273	else:
1274	if len(Pdabc):
1275	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1276	dada0 = 0
1277	else:
1278	alp = G2mth.getTOFalpha(parmDict,dsp)
1279	dada0 = G2mth.getTOFalphaDeriv(dsp)
1280	betName = 'bet'+str(iPeak)
1281	if betName in varyList:
1282	bet = parmDict[betName]
1283	else:
1284	if len(Pdabc):
1285	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1286	dbdb0 = dbdb1 = dbdb2 = 0
1287	else:
1288	bet = G2mth.getTOFbeta(parmDict,dsp)
1289	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1290	sigName = 'sig'+str(iPeak)
1291	if sigName in varyList:
1292	sig = parmDict[sigName]
1293	dsds0 = dsds1 = dsds2 = dsds3 = 0
1294	else:
1295	sig = G2mth.getTOFsig(parmDict,dsp)
1296	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1297	gamName = 'gam'+str(iPeak)
1298	if gamName in varyList:
1299	gam = parmDict[gamName]
1300	dsdX = dsdY = 0
1301	else:
1302	gam = G2mth.getTOFgamma(parmDict,dsp)
1303	dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1304	gam = max(gam,0.001) #avoid neg gamma
1305	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1306	iBeg = np.searchsorted(xdata,pos-fmin)
1307	lenX = len(xdata)
1308	if not iBeg:
1309	iFin = np.searchsorted(xdata,pos+fmax)
1310	elif iBeg == lenX:
1311	iFin = iBeg
1312	else:
1313	iFin = np.searchsorted(xdata,pos+fmax)
1314	if not iBeg+iFin: #peak below low limit
1315	iPeak += 1
1316	continue
1317	elif not iBeg-iFin: #peak above high limit
1318	break
1319	dMdpk = np.zeros(shape=(7,len(xdata)))
1320	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1321	for i in range(1,6):
1322	dMdpk[i][iBeg:iFin] += intenscw[iBeg:iFin]dMdipk[i]
1323	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1324	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1325	for parmName in ['pos','int','alp','bet','sig','gam']:
1326	try:
1327	idx = varyList.index(parmName+str(iPeak))
1328	dMdv[idx] = dervDict[parmName]
1329	except ValueError:
1330	pass
1331	if 'alpha' in varyList:
1332	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1333	if 'beta-0' in varyList:
1334	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1335	if 'beta-1' in varyList:
1336	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1337	if 'beta-q' in varyList:
1338	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1339	if 'sig-0' in varyList:
1340	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1341	if 'sig-1' in varyList:
1342	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1343	if 'sig-2' in varyList:
1344	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1345	if 'sig-q' in varyList:
1346	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1347	if 'X' in varyList:
1348	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1349	if 'Y' in varyList:
1350	dMdv[varyList.index('Y')] += dsdY*dervDict['gam'] #problem here
1351	iPeak += 1
1352	except KeyError: #no more peaks to process
1353	break
1354	return dMdv
1355
1356	def Dict2Values(parmdict, varylist):
1357	'''Use before call to leastsq to setup list of values for the parameters
1358	in parmdict, as selected by key in varylist'''
1359	return [parmdict[key] for key in varylist]
1360
1361	def Values2Dict(parmdict, varylist, values):
1362	''' Use after call to leastsq to update the parameter dictionary with
1363	values corresponding to keys in varylist'''
1364	parmdict.update(zip(varylist,values))
1365
1366	def SetBackgroundParms(Background):
1367	'needs a doc string'
1368	if len(Background) == 1: # fix up old backgrounds
1369	Background.append({'nDebye':0,'debyeTerms':[]})
1370	bakType,bakFlag = Background[0][:2]
1371	backVals = Background[0][3:]
1372	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1373	Debye = Background[1] #also has background peaks stuff
1374	backDict = dict(zip(backNames,backVals))
1375	backVary = []
1376	if bakFlag:
1377	backVary = backNames
1378
1379	backDict['nDebye'] = Debye['nDebye']
1380	debyeDict = {}
1381	debyeList = []
1382	for i in range(Debye['nDebye']):
1383	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1384	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1385	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1386	debyeVary = []
1387	for item in debyeList:
1388	if item[1]:
1389	debyeVary.append(item[0])
1390	backDict.update(debyeDict)
1391	backVary += debyeVary
1392
1393	backDict['nPeaks'] = Debye['nPeaks']
1394	peaksDict = {}
1395	peaksList = []
1396	for i in range(Debye['nPeaks']):
1397	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1398	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1399	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1400	peaksVary = []
1401	for item in peaksList:
1402	if item[1]:
1403	peaksVary.append(item[0])
1404	backDict.update(peaksDict)
1405	backVary += peaksVary
1406	return bakType,backDict,backVary
1407
1408	def DoCalibInst(IndexPeaks,Inst):
1409
1410	def SetInstParms():
1411	dataType = Inst['Type'][0]
1412	insVary = []
1413	insNames = []
1414	insVals = []
1415	for parm in Inst:
1416	insNames.append(parm)
1417	insVals.append(Inst[parm][1])
1418	if parm in ['Lam','difC','difA','difB','Zero',]:
1419	if Inst[parm][2]:
1420	insVary.append(parm)
1421	instDict = dict(zip(insNames,insVals))
1422	return dataType,instDict,insVary
1423
1424	def GetInstParms(parmDict,Inst,varyList):
1425	for name in Inst:
1426	Inst[name][1] = parmDict[name]
1427
1428	def InstPrint(Inst,sigDict):
1429	print 'Instrument Parameters:'
1430	if 'C' in Inst['Type'][0]:
1431	ptfmt = "%12.6f"
1432	else:
1433	ptfmt = "%12.3f"
1434	ptlbls = 'names :'
1435	ptstr = 'values:'
1436	sigstr = 'esds :'
1437	for parm in Inst:
1438	if parm in ['Lam','difC','difA','difB','Zero',]:
1439	ptlbls += "%s" % (parm.center(12))
1440	ptstr += ptfmt % (Inst[parm][1])
1441	if parm in sigDict:
1442	sigstr += ptfmt % (sigDict[parm])
1443	else:
1444	sigstr += 12*' '
1445	print ptlbls
1446	print ptstr
1447	print sigstr
1448
1449	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1450	parmDict.update(zip(varyList,values))
1451	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1452
1453	peakPos = []
1454	peakDsp = []
1455	peakWt = []
1456	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1457	if peak[2] and peak[3] and sig > 0.:
1458	peakPos.append(peak[0])
1459	peakDsp.append(peak[-1]) #d-calc
1460	# peakWt.append(peak[-1]2/sig2) #weight by d**2
1461	peakWt.append(1./(sig*peak[-1])) #
1462	peakPos = np.array(peakPos)
1463	peakDsp = np.array(peakDsp)
1464	peakWt = np.array(peakWt)
1465	dataType,insDict,insVary = SetInstParms()
1466	parmDict = {}
1467	parmDict.update(insDict)
1468	varyList = insVary
1469	if not len(varyList):
1470	print '** ERROR - nothing to refine! **'
1471	return False
1472	while True:
1473	begin = time.time()
1474	values = np.array(Dict2Values(parmDict, varyList))
1475	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1476	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1477	ncyc = int(result[2]['nfev']/2)
1478	runtime = time.time()-begin
1479	chisq = np.sum(result[2]['fvec']**2)
1480	Values2Dict(parmDict, varyList, result[0])
1481	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
1482	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1483	print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1484	print 'chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF)
1485	try:
1486	sig = np.sqrt(np.diag(result[1])*GOF)
1487	if np.any(np.isnan(sig)):
1488	print '* Least squares aborted - some invalid esds possible *'
1489	break #refinement succeeded - finish up!
1490	except ValueError: #result[1] is None on singular matrix
1491	print '** Refinement failed - singular matrix **'
1492
1493	sigDict = dict(zip(varyList,sig))
1494	GetInstParms(parmDict,Inst,varyList)
1495	InstPrint(Inst,sigDict)
1496	return True
1497
1498	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1499	'''Called to perform a peak fit, refining the selected items in the peak
1500	table as well as selected items in the background.
1501
1502	:param str FitPgm: type of fit to perform. At present this is ignored.
1503	:param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1504	four values followed by a refine flag where the values are: position, intensity,
1505	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1506	tree entry, dict item "peaks"
1507	:param list Background: describes the background. List with two items.
1508	Item 0 specifies a background model and coefficients. Item 1 is a dict.
1509	From the Histogram/Background tree entry.
1510	:param list Limits: min and max x-value to use
1511	:param dict Inst: Instrument parameters
1512	:param dict Inst2: more Instrument parameters
1513	:param numpy.array data: a 5xn array. data[0] is the x-values,
1514	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1515	calc, background and difference intensities, respectively.
1516	:param list prevVaryList: Used in sequential refinements to override the
1517	variable list. Defaults as an empty list.
1518	:param bool oneCycle: True if only one cycle of fitting should be performed
1519	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1520	and derivType = controls['deriv type']. If None default values are used.
1521	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1522	Defaults to None, which means no updates are done.
1523	'''
1524	def GetBackgroundParms(parmList,Background):
1525	iBak = 0
1526	while True:
1527	try:
1528	bakName = 'Back;'+str(iBak)
1529	Background[0][iBak+3] = parmList[bakName]
1530	iBak += 1
1531	except KeyError:
1532	break
1533	iDb = 0
1534	while True:
1535	names = ['DebyeA;','DebyeR;','DebyeU;']
1536	try:
1537	for i,name in enumerate(names):
1538	val = parmList[name+str(iDb)]
1539	Background[1]['debyeTerms'][iDb][2*i] = val
1540	iDb += 1
1541	except KeyError:
1542	break
1543	iDb = 0
1544	while True:
1545	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1546	try:
1547	for i,name in enumerate(names):
1548	val = parmList[name+str(iDb)]
1549	Background[1]['peaksList'][iDb][2*i] = val
1550	iDb += 1
1551	except KeyError:
1552	break
1553
1554	def BackgroundPrint(Background,sigDict):
1555	print 'Background coefficients for',Background[0][0],'function'
1556	ptfmt = "%12.5f"
1557	ptstr = 'value: '
1558	sigstr = 'esd : '
1559	for i,back in enumerate(Background[0][3:]):
1560	ptstr += ptfmt % (back)
1561	if Background[0][1]:
1562	prm = 'Back;'+str(i)
1563	if prm in sigDict:
1564	sigstr += ptfmt % (sigDict[prm])
1565	else:
1566	sigstr += " "*12
1567	if len(ptstr) > 75:
1568	print ptstr
1569	if Background[0][1]: print sigstr
1570	ptstr = 'value: '
1571	sigstr = 'esd : '
1572	if len(ptstr) > 8:
1573	print ptstr
1574	if Background[0][1]: print sigstr
1575
1576	if Background[1]['nDebye']:
1577	parms = ['DebyeA;','DebyeR;','DebyeU;']
1578	print 'Debye diffuse scattering coefficients'
1579	ptfmt = "%12.5f"
1580	print ' term DebyeA esd DebyeR esd DebyeU esd'
1581	for term in range(Background[1]['nDebye']):
1582	line = ' term %d'%(term)
1583	for ip,name in enumerate(parms):
1584	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1585	if name+str(term) in sigDict:
1586	line += ptfmt%(sigDict[name+str(term)])
1587	else:
1588	line += " "*12
1589	print line
1590	if Background[1]['nPeaks']:
1591	print 'Coefficients for Background Peaks'
1592	ptfmt = "%15.3f"
1593	for j,pl in enumerate(Background[1]['peaksList']):
1594	names = 'peak %3d:'%(j+1)
1595	ptstr = 'values :'
1596	sigstr = 'esds :'
1597	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1598	val = pl[2*i]
1599	prm = lbl+";"+str(j)
1600	names += '%15s'%(prm)
1601	ptstr += ptfmt%(val)
1602	if prm in sigDict:
1603	sigstr += ptfmt%(sigDict[prm])
1604	else:
1605	sigstr += " "*15
1606	print names
1607	print ptstr
1608	print sigstr
1609
1610	def SetInstParms(Inst):
1611	dataType = Inst['Type'][0]
1612	insVary = []
1613	insNames = []
1614	insVals = []
1615	for parm in Inst:
1616	insNames.append(parm)
1617	insVals.append(Inst[parm][1])
1618	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1619	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1620	insVary.append(parm)
1621	instDict = dict(zip(insNames,insVals))
1622	instDict['X'] = max(instDict['X'],0.01)
1623	instDict['Y'] = max(instDict['Y'],0.01)
1624	if 'SH/L' in instDict:
1625	instDict['SH/L'] = max(instDict['SH/L'],0.002)
1626	return dataType,instDict,insVary
1627
1628	def GetInstParms(parmDict,Inst,varyList,Peaks):
1629	for name in Inst:
1630	Inst[name][1] = parmDict[name]
1631	iPeak = 0
1632	while True:
1633	try:
1634	sigName = 'sig'+str(iPeak)
1635	pos = parmDict['pos'+str(iPeak)]
1636	if sigName not in varyList:
1637	if 'C' in Inst['Type'][0]:
1638	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1639	else:
1640	dsp = G2lat.Pos2dsp(Inst,pos)
1641	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1642	gamName = 'gam'+str(iPeak)
1643	if gamName not in varyList:
1644	if 'C' in Inst['Type'][0]:
1645	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1646	else:
1647	dsp = G2lat.Pos2dsp(Inst,pos)
1648	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1649	iPeak += 1
1650	except KeyError:
1651	break
1652
1653	def InstPrint(Inst,sigDict):
1654	print 'Instrument Parameters:'
1655	ptfmt = "%12.6f"
1656	ptlbls = 'names :'
1657	ptstr = 'values:'
1658	sigstr = 'esds :'
1659	for parm in Inst:
1660	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1661	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1662	ptlbls += "%s" % (parm.center(12))
1663	ptstr += ptfmt % (Inst[parm][1])
1664	if parm in sigDict:
1665	sigstr += ptfmt % (sigDict[parm])
1666	else:
1667	sigstr += 12*' '
1668	print ptlbls
1669	print ptstr
1670	print sigstr
1671
1672	def SetPeaksParms(dataType,Peaks):
1673	peakNames = []
1674	peakVary = []
1675	peakVals = []
1676	if 'C' in dataType:
1677	names = ['pos','int','sig','gam']
1678	else:
1679	names = ['pos','int','alp','bet','sig','gam']
1680	for i,peak in enumerate(Peaks):
1681	for j,name in enumerate(names):
1682	peakVals.append(peak[2*j])
1683	parName = name+str(i)
1684	peakNames.append(parName)
1685	if peak[2*j+1]:
1686	peakVary.append(parName)
1687	return dict(zip(peakNames,peakVals)),peakVary
1688
1689	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1690	if 'C' in Inst['Type'][0]:
1691	names = ['pos','int','sig','gam']
1692	else: #'T'
1693	names = ['pos','int','alp','bet','sig','gam']
1694	for i,peak in enumerate(Peaks):
1695	pos = parmDict['pos'+str(i)]
1696	if 'difC' in Inst:
1697	dsp = pos/Inst['difC'][1]
1698	for j in range(len(names)):
1699	parName = names[j]+str(i)
1700	if parName in varyList:
1701	peak[2*j] = parmDict[parName]
1702	elif 'alpha' in parName:
1703	peak[2*j] = parmDict['alpha']/dsp
1704	elif 'beta' in parName:
1705	peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1706	elif 'sig' in parName:
1707	if 'C' in Inst['Type'][0]:
1708	peak[2*j] = G2mth.getCWsig(parmDict,pos)
1709	else:
1710	peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1711	elif 'gam' in parName:
1712	if 'C' in Inst['Type'][0]:
1713	peak[2*j] = G2mth.getCWgam(parmDict,pos)
1714	else:
1715	peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1716
1717	def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1718	print 'Peak coefficients:'
1719	if 'C' in dataType:
1720	names = ['pos','int','sig','gam']
1721	else: #'T'
1722	names = ['pos','int','alp','bet','sig','gam']
1723	head = 13*' '
1724	for name in names:
1725	if name in ['alp','bet']:
1726	head += name.center(8)+'esd'.center(8)
1727	else:
1728	head += name.center(10)+'esd'.center(10)
1729	head += 'bins'.center(8)
1730	print head
1731	if 'C' in dataType:
1732	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1733	else:
1734	ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1735	for i,peak in enumerate(Peaks):
1736	ptstr = ':'
1737	for j in range(len(names)):
1738	name = names[j]
1739	parName = name+str(i)
1740	ptstr += ptfmt[name] % (parmDict[parName])
1741	if parName in varyList:
1742	ptstr += ptfmt[name] % (sigDict[parName])
1743	else:
1744	if name in ['alp','bet']:
1745	ptstr += 8*' '
1746	else:
1747	ptstr += 10*' '
1748	ptstr += '%9.2f'%(ptsperFW[i])
1749	print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1750
1751	def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1752	parmdict.update(zip(varylist,values))
1753	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1754
1755	def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):
1756	parmdict.update(zip(varylist,values))
1757	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1758	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
1759	if dlg:
1760	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1761	if not GoOn:
1762	return -M #abort!!
1763	return M
1764
1765	if controls:
1766	Ftol = controls['min dM/M']
1767	else:
1768	Ftol = 0.0001
1769	if oneCycle:
1770	Ftol = 1.0
1771	x,y,w,yc,yb,yd = data #these are numpy arrays!
1772	yc *= 0. #set calcd ones to zero
1773	yb *= 0.
1774	yd *= 0.
1775	cw = x[1:]-x[:-1]
1776	xBeg = np.searchsorted(x,Limits[0])
1777	xFin = np.searchsorted(x,Limits[1])+1
1778	bakType,bakDict,bakVary = SetBackgroundParms(Background)
1779	dataType,insDict,insVary = SetInstParms(Inst)
1780	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1781	parmDict = {}
1782	parmDict.update(bakDict)
1783	parmDict.update(insDict)
1784	parmDict.update(peakDict)
1785	parmDict['Pdabc'] = [] #dummy Pdabc
1786	parmDict.update(Inst2) #put in real one if there
1787	if prevVaryList:
1788	varyList = prevVaryList[:]
1789	else:
1790	varyList = bakVary+insVary+peakVary
1791	fullvaryList = varyList[:]
1792	while True:
1793	begin = time.time()
1794	values = np.array(Dict2Values(parmDict, varyList))
1795	Rvals = {}
1796	badVary = []
1797	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1798	args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1799	ncyc = int(result[2]['nfev']/2)
1800	runtime = time.time()-begin
1801	chisq = np.sum(result[2]['fvec']**2)
1802	Values2Dict(parmDict, varyList, result[0])
1803	Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]y[xBeg:xFin]2))100. #to %
1804	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
1805	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1806	if ncyc:
1807	print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1808	print 'Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1809	sig = [0]*len(varyList)
1810	if len(varyList) == 0: break # if nothing was refined
1811	try:
1812	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1813	if np.any(np.isnan(sig)):
1814	print '* Least squares aborted - some invalid esds possible *'
1815	break #refinement succeeded - finish up!
1816	except ValueError: #result[1] is None on singular matrix
1817	print '** Refinement failed - singular matrix **'
1818	Ipvt = result[2]['ipvt']
1819	for i,ipvt in enumerate(Ipvt):
1820	if not np.sum(result[2]['fjac'],axis=1)[i]:
1821	print 'Removing parameter: ',varyList[ipvt-1]
1822	badVary.append(varyList[ipvt-1])
1823	del(varyList[ipvt-1])
1824	break
1825	else: # nothing removed
1826	break
1827	if dlg: dlg.Destroy()
1828	sigDict = dict(zip(varyList,sig))
1829	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1830	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1831	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1832	GetBackgroundParms(parmDict,Background)
1833	if bakVary: BackgroundPrint(Background,sigDict)
1834	GetInstParms(parmDict,Inst,varyList,Peaks)
1835	if insVary: InstPrint(Inst,sigDict)
1836	GetPeaksParms(Inst,parmDict,Peaks,varyList)
1837	binsperFWHM = []
1838	for peak in Peaks:
1839	FWHM = getFWHM(peak[0],Inst)
1840	try:
1841	binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
1842	except IndexError:
1843	binsperFWHM.append(0.)
1844	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
1845	if len(binsperFWHM):
1846	if min(binsperFWHM) < 1.:
1847	print '* Warning: calculated peak widths are too narrow to refine profile coefficients *'
1848	if 'T' in Inst['Type'][0]:
1849	print ' Manually increase sig-0, 1, or 2 in Instrument Parameters'
1850	else:
1851	print ' Manually increase W in Instrument Parameters'
1852	elif min(binsperFWHM) < 4.:
1853	print '* Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement *'
1854	print '* recommended is 6-10; you have %.2f *'%(min(binsperFWHM))
1855	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1856
1857	def calcIncident(Iparm,xdata):
1858	'needs a doc string'
1859
1860	def IfunAdv(Iparm,xdata):
1861	Itype = Iparm['Itype']
1862	Icoef = Iparm['Icoeff']
1863	DYI = np.ones((12,xdata.shape[0]))
1864	YI = np.ones_like(xdata)*Icoef[0]
1865
1866	x = xdata/1000. #expressions are in ms
1867	if Itype == 'Exponential':
1868	for i in [1,3,5,7,9]:
1869	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1870	YI += Icoef[i]*Eterm
1871	DYI[i] *= Eterm
1872	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1873	elif 'Maxwell'in Itype:
1874	Eterm = np.exp(-Icoef[2]/x**2)
1875	DYI[1] = Eterm/x**5
1876	DYI[2] = -Icoef[1]DYI[1]/x*2
1877	YI += (Icoef[1]Eterm/x*5)
1878	if 'Exponential' in Itype:
1879	for i in range(3,11,2):
1880	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1881	YI += Icoef[i]*Eterm
1882	DYI[i] *= Eterm
1883	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1884	else: #Chebyschev
1885	T = (2./x)-1.
1886	Ccof = np.ones((12,xdata.shape[0]))
1887	Ccof[1] = T
1888	for i in range(2,12):
1889	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
1890	for i in range(1,10):
1891	YI += Ccof[i]*Icoef[i+2]
1892	DYI[i+2] =Ccof[i]
1893	return YI,DYI
1894
1895	Iesd = np.array(Iparm['Iesd'])
1896	Icovar = Iparm['Icovar']
1897	YI,DYI = IfunAdv(Iparm,xdata)
1898	YI = np.where(YI>0,YI,1.)
1899	WYI = np.zeros_like(xdata)
1900	vcov = np.zeros((12,12))
1901	k = 0
1902	for i in range(12):
1903	for j in range(i,12):
1904	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
1905	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
1906	k += 1
1907	M = np.inner(vcov,DYI.T)
1908	WYI = np.sum(M*DYI,axis=0)
1909	WYI = np.where(WYI>0.,WYI,0.)
1910	return YI,WYI
1911
1912	################################################################################
1913	# Reflectometry calculations
1914	################################################################################
1915
1916	def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
1917	print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
1918
1919	class RandomDisplacementBounds(object):
1920	"""random displacement with bounds"""
1921	def __init__(self, xmin, xmax, stepsize=0.5):
1922	self.xmin = xmin
1923	self.xmax = xmax
1924	self.stepsize = stepsize
1925
1926	def __call__(self, x):
1927	"""take a random step but ensure the new position is within the bounds"""
1928	while True:
1929	# this could be done in a much more clever way, but it will work for example purposes
1930	steps = self.xmax-self.xmin
1931	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
1932	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
1933	break
1934	return xnew
1935
1936	def GetModelParms():
1937	parmDict = {}
1938	varyList = []
1939	values = []
1940	bounds = []
1941	parmDict['dQ type'] = data['dQ type']
1942	parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig
1943	for parm in ['Scale','FltBack']:
1944	parmDict[parm] = data[parm][0]
1945	if data[parm][1]:
1946	varyList.append(parm)
1947	values.append(data[parm][0])
1948	bounds.append(Bounds[parm])
1949	parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
1950	parmDict['nLayers'] = len(parmDict['Layer Seq'])
1951	for ilay,layer in enumerate(data['Layers']):
1952	name = layer['Name']
1953	cid = str(ilay)+';'
1954	parmDict[cid+'Name'] = name
1955	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1956	parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
1957	if layer.get(parm,[0.,False])[1]:
1958	varyList.append(cid+parm)
1959	value = layer[parm][0]
1960	values.append(value)
1961	if value:
1962	bound = [value*Bfac,value/Bfac]
1963	else:
1964	bound = [0.,10.]
1965	bounds.append(bound)
1966	if name not in ['vacuum','unit scatter']:
1967	parmDict[cid+'rho'] = Substances[name]['Scatt density']
1968	parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
1969	return parmDict,varyList,values,bounds
1970
1971	def SetModelParms():
1972	line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1973	if 'Scale' in varyList:
1974	data['Scale'][0] = parmDict['Scale']
1975	line += ' esd: %.4g'%(sigDict['Scale'])
1976	print line
1977	line = ' Flat background: %15.4g'%(parmDict['FltBack'])
1978	if 'FltBack' in varyList:
1979	data['FltBack'][0] = parmDict['FltBack']
1980	line += ' esd: %15.3g'%(sigDict['FltBack'])
1981	print line
1982	for ilay,layer in enumerate(data['Layers']):
1983	name = layer['Name']
1984	print ' Parameters for layer: %d %s'%(ilay,name)
1985	cid = str(ilay)+';'
1986	line = ' '
1987	line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
1988	line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
1989	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
1990	if parm in layer:
1991	if parm == 'Rough':
1992	layer[parm][0] = abs(parmDict[cid+parm]) #make positive
1993	else:
1994	layer[parm][0] = parmDict[cid+parm]
1995	line += ' %s: %.3f'%(parm,layer[parm][0])
1996	if cid+parm in varyList:
1997	line += ' esd: %.3g'%(sigDict[cid+parm])
1998	print line
1999	print line2
2000
2001	def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2002	parmDict.update(zip(varyList,values))
2003	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2004	return M
2005
2006	def sumREFD(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 np.sum(M**2)
2010
2011	def getREFD(Q,Qsig,parmDict):
2012	Ic = np.ones_like(Q)*parmDict['FltBack']
2013	Scale = parmDict['Scale']
2014	Nlayers = parmDict['nLayers']
2015	Res = parmDict['Res']
2016	depth = np.zeros(Nlayers)
2017	rho = np.zeros(Nlayers)
2018	irho = np.zeros(Nlayers)
2019	sigma = np.zeros(Nlayers)
2020	for ilay,lay in enumerate(parmDict['Layer Seq']):
2021	cid = str(lay)+';'
2022	depth[ilay] = parmDict[cid+'Thick']
2023	sigma[ilay] = parmDict[cid+'Rough']
2024	if parmDict[cid+'Name'] == u'unit scatter':
2025	rho[ilay] = parmDict[cid+'DenMul']
2026	irho[ilay] = parmDict[cid+'iDenMul']
2027	elif 'vacuum' != parmDict[cid+'Name']:
2028	rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
2029	irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
2030	if cid+'Mag SLD' in parmDict:
2031	rho[ilay] += parmDict[cid+'Mag SLD']
2032	if parmDict['dQ type'] == 'None':
2033	AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
2034	elif 'const' in parmDict['dQ type']:
2035	AB = SmearAbeles(0.5Q,QRes,depth,rho,irho,sigma[1:])
2036	else: #dQ/Q in data
2037	AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
2038	Ic += AB*Scale
2039	return Ic
2040
2041	def estimateT0(takestep):
2042	Mmax = -1.e-10
2043	Mmin = 1.e10
2044	for i in range(100):
2045	x0 = takestep(values)
2046	M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2047	Mmin = min(M,Mmin)
2048	MMax = max(M,Mmax)
2049	return 1.5*(MMax-Mmin)
2050
2051	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2052	if data.get('2% weight'):
2053	wt = 1./(0.02Io)*2
2054	Qmin = Limits[1][0]
2055	Qmax = Limits[1][1]
2056	wtFactor = ProfDict['wtFactor']
2057	Bfac = data['Toler']
2058	Ibeg = np.searchsorted(Q,Qmin)
2059	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
2060	Ic[:] = 0
2061	Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
2062	'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
2063	parmDict,varyList,values,bounds = GetModelParms()
2064	Msg = 'Failed to converge'
2065	if varyList:
2066	if data['Minimizer'] == 'LMLS':
2067	result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
2068	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2069	parmDict.update(zip(varyList,result[0]))
2070	chisq = np.sum(result[2]['fvec']**2)
2071	ncalc = result[2]['nfev']
2072	covM = result[1]
2073	newVals = result[0]
2074	elif data['Minimizer'] == 'Basin Hopping':
2075	xyrng = np.array(bounds).T
2076	take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
2077	T0 = estimateT0(take_step)
2078	print ' Estimated temperature: %.3g'%(T0)
2079	result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
2080	interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
2081	'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
2082	chisq = result.fun
2083	ncalc = result.nfev
2084	newVals = result.x
2085	covM = []
2086	elif data['Minimizer'] == 'MC/SA Anneal':
2087	xyrng = np.array(bounds).T
2088	result = G2mth.anneal(sumREFD, values,
2089	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
2090	schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
2091	boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
2092	ranRange=0.20,autoRan=False,dlg=None)
2093	newVals = result[0]
2094	parmDict.update(zip(varyList,newVals))
2095	chisq = result[1]
2096	ncalc = result[3]
2097	covM = []
2098	print ' MC/SA final temperature: %.4g'%(result[2])
2099	elif data['Minimizer'] == 'L-BFGS-B':
2100	result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol,
2101	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
2102	parmDict.update(zip(varyList,result['x']))
2103	chisq = result.fun
2104	ncalc = result.nfev
2105	newVals = result.x
2106	covM = []
2107	else: #nothing varied
2108	M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
2109	chisq = np.sum(M**2)
2110	ncalc = 0
2111	covM = []
2112	sig = []
2113	sigDict = {}
2114	result = []
2115	Rvals = {}
2116	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]Io[Ibeg:Ifin]2))100. #to %
2117	Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2
2118	Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
2119	Ib[Ibeg:Ifin] = parmDict['FltBack']
2120	try:
2121	if not len(varyList):
2122	Msg += ' - nothing refined'
2123	raise ValueError
2124	Nans = np.isnan(newVals)
2125	if np.any(Nans):
2126	name = varyList[Nans.nonzero(True)[0]]
2127	Msg += ' Nan result for '+name+'!'
2128	raise ValueError
2129	Negs = np.less_equal(newVals,0.)
2130	if np.any(Negs):
2131	indx = Negs.nonzero()
2132	name = varyList[indx[0][0]]
2133	if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
2134	Msg += ' negative coefficient for '+name+'!'
2135	raise ValueError
2136	if len(covM):
2137	sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
2138	covMatrix = covM*Rvals['GOF']
2139	else:
2140	sig = np.zeros(len(varyList))
2141	covMatrix = []
2142	sigDict = dict(zip(varyList,sig))
2143	print ' Results of reflectometry data modelling fit:'
2144	print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2145	print 'Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2146	SetModelParms()
2147	return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
2148	except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans
2149	print Msg
2150	return False,0,0,0,0,0,0,Msg
2151
2152	def makeSLDprofile(data,Substances):
2153
2154	sq2 = np.sqrt(2.)
2155	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2156	Nlayers = len(laySeq)
2157	laySeq = np.array(laySeq,dtype=int)
2158	interfaces = np.zeros(Nlayers)
2159	rho = np.zeros(Nlayers)
2160	sigma = np.zeros(Nlayers)
2161	name = data['Layers'][0]['Name']
2162	thick = 0.
2163	for ilay,lay in enumerate(laySeq):
2164	layer = data['Layers'][lay]
2165	name = layer['Name']
2166	if 'Thick' in layer:
2167	thick += layer['Thick'][0]
2168	interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
2169	if 'Rough' in layer:
2170	sigma[ilay] = max(0.001,layer['Rough'][0])
2171	if name != 'vacuum':
2172	if name == 'unit scatter':
2173	rho[ilay] = np.sqrt(layer['DenMul'][0]2+layer['iDenMul'][0]2)
2174	else:
2175	rrho = Substances[name]['Scatt density']
2176	irho = Substances[name]['XImag density']
2177	rho[ilay] = np.sqrt(rrho2+irho2)*layer['DenMul'][0]
2178	if 'Mag SLD' in layer:
2179	rho[ilay] += layer['Mag SLD'][0]
2180	name = data['Layers'][-1]['Name']
2181	x = np.linspace(-0.15thick,1.15thick,1000,endpoint=True)
2182	xr = np.flipud(x)
2183	interfaces[-1] = x[-1]
2184	y = np.ones_like(x)*rho[0]
2185	iBeg = 0
2186	for ilayer in range(Nlayers-1):
2187	delt = rho[ilayer+1]-rho[ilayer]
2188	iPos = np.searchsorted(x,interfaces[ilayer])
2189	y[iBeg:] += (delt/2.)sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2sigma[ilayer+1]))
2190	iBeg = iPos
2191	return x,xr,y
2192
2193	def REFDModelFxn(Profile,Inst,Limits,Substances,data):
2194
2195	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
2196	Qmin = Limits[1][0]
2197	Qmax = Limits[1][1]
2198	iBeg = np.searchsorted(Q,Qmin)
2199	iFin = np.searchsorted(Q,Qmax)+1 #include last point
2200	Ib[:] = data['FltBack'][0]
2201	Ic[:] = 0
2202	Scale = data['Scale'][0]
2203	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
2204	Nlayers = len(laySeq)
2205	depth = np.zeros(Nlayers)
2206	rho = np.zeros(Nlayers)
2207	irho = np.zeros(Nlayers)
2208	sigma = np.zeros(Nlayers)
2209	for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
2210	layer = data['Layers'][lay]
2211	name = layer['Name']
2212	if 'Thick' in layer: #skips first & last layers
2213	depth[ilay] = layer['Thick'][0]
2214	if 'Rough' in layer: #skips first layer
2215	sigma[ilay] = layer['Rough'][0]
2216	if 'unit scatter' == name:
2217	rho[ilay] = layer['DenMul'][0]
2218	irho[ilay] = layer['iDenMul'][0]
2219	else:
2220	rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
2221	irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
2222	if 'Mag SLD' in layer:
2223	rho[ilay] += layer['Mag SLD'][0]
2224	if data['dQ type'] == 'None':
2225	AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
2226	elif 'const' in data['dQ type']:
2227	res = data['Resolution'][0]/(100.*sateln2)
2228	AB = SmearAbeles(0.5Q[iBeg:iFin],resQ[iBeg:iFin],depth,rho,irho,sigma[1:])
2229	else: #dQ/Q in data
2230	AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
2231	Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
2232
2233	def abeles(kz, depth, rho, irho=0, sigma=0):
2234	"""
2235	Optical matrix form of the reflectivity calculation.
2236	O.S. Heavens, Optical Properties of Thin Solid Films
2237
2238	Reflectometry as a function of kz for a set of slabs.
2239
2240	:param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
2241	This is :math:`\\tfrac12 Q_z`.
2242	:param depth: float[m] (Ang).
2243	thickness of each layer. The thickness of the incident medium
2244	and substrate are ignored.
2245	:param rho: float[n,k] (1e-6/Ang^2)
2246	Real scattering length density for each layer for each kz
2247	:param irho: float[n,k] (1e-6/Ang^2)
2248	Imaginary scattering length density for each layer for each kz
2249	Note: absorption cross section mu = 2 irho/lambda for neutrons
2250	:param sigma: float[m-1] (Ang)
2251	interfacial roughness. This is the roughness between a layer
2252	and the previous layer. The sigma array should have m-1 entries.
2253
2254	Slabs are ordered with the surface SLD at index 0 and substrate at
2255	index -1, or reversed if kz < 0.
2256	"""
2257	def calc(kz, depth, rho, irho, sigma):
2258	if len(kz) == 0: return kz
2259
2260	# Complex index of refraction is relative to the incident medium.
2261	# We can get the same effect using kz_rel^2 = kz^2 + 4pirho_o
2262	# in place of kz^2, and ignoring rho_o
2263	kz_sq = kz*2 + 4e-6np.pi*rho[:,0]
2264	k = kz
2265
2266	# According to Heavens, the initial matrix should be [ 1 F; F 1],
2267	# which we do by setting B=I and M0 to [1 F; F 1]. An extra matrix
2268	# multiply versus some coding convenience.
2269	B11 = 1
2270	B22 = 1
2271	B21 = 0
2272	B12 = 0
2273	for i in range(0, len(depth)-1):
2274	k_next = np.sqrt(kz_sq - 4e-6np.pi(rho[:,i+1] + 1j*irho[:,i+1]))
2275	F = (k - k_next) / (k + k_next)
2276	F = np.exp(-2kk_nextsigma[i]**2)
2277	#print "==== layer",i
2278	#print "kz:", kz
2279	#print "k:", k
2280	#print "k_next:",k_next
2281	#print "F:",F
2282	#print "rho:",rho[:,i+1]
2283	#print "irho:",irho[:,i+1]
2284	#print "d:",depth[i],"sigma:",sigma[i]
2285	M11 = np.exp(1jkdepth[i]) if i>0 else 1
2286	M22 = np.exp(-1jkdepth[i]) if i>0 else 1
2287	M21 = F*M11
2288	M12 = F*M22
2289	C1 = B11M11 + B21M12
2290	C2 = B11M21 + B21M22
2291	B11 = C1
2292	B21 = C2
2293	C1 = B12M11 + B22M12
2294	C2 = B12M21 + B22M22
2295	B12 = C1
2296	B22 = C2
2297	k = k_next
2298
2299	r = B12/B11
2300	return np.absolute(r)**2
2301
2302	if np.isscalar(kz): kz = np.asarray([kz], 'd')
2303
2304	m = len(depth)
2305
2306	# Make everything into arrays
2307	depth = np.asarray(depth,'d')
2308	rho = np.asarray(rho,'d')
2309	irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
2310	sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
2311
2312	# Repeat rho,irho columns as needed
2313	if len(rho.shape) == 1:
2314	rho = rho[None,:]
2315	irho = irho[None,:]
2316
2317	return calc(kz, depth, rho, irho, sigma)
2318
2319	def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
2320	y = abeles(kz, depth, rho, irho, sigma)
2321	s = dq/2.
2322	y += 0.1354(abeles(kz+2s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
2323	y += 0.24935(abeles(kz-5s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma))
2324	y += 0.4111(abeles(kz-4s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma))
2325	y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
2326	y += 0.80074(abeles(kz-2s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
2327	y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
2328	y *= 0.137023
2329	return y
2330
2331	def makeRefdFFT(Limits,Profile):
2332	print 'make fft'
2333	Q,Io = Profile[:2]
2334	Qmin = Limits[1][0]
2335	Qmax = Limits[1][1]
2336	iBeg = np.searchsorted(Q,Qmin)
2337	iFin = np.searchsorted(Q,Qmax)+1 #include last point
2338	Qf = np.linspace(0.,Q[iFin-1],5000)
2339	QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
2340	If = QI(Qf)Qf*4
2341	R = np.linspace(0.,5000.,5000)
2342	F = fft.rfft(If)
2343	return R,F
2344
2345
2346	################################################################################
2347	# Stacking fault simulation codes
2348	################################################################################
2349
2350	def GetStackParms(Layers):
2351
2352	Parms = []
2353	#cell parms
2354	if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
2355	Parms.append('cellA')
2356	Parms.append('cellC')
2357	else:
2358	Parms.append('cellA')
2359	Parms.append('cellB')
2360	Parms.append('cellC')
2361	if Layers['Laue'] != 'mmm':
2362	Parms.append('cellG')
2363	#Transition parms
2364	for iY in range(len(Layers['Layers'])):
2365	for iX in range(len(Layers['Layers'])):
2366	Parms.append('TransP;%d;%d'%(iY,iX))
2367	Parms.append('TransX;%d;%d'%(iY,iX))
2368	Parms.append('TransY;%d;%d'%(iY,iX))
2369	Parms.append('TransZ;%d;%d'%(iY,iX))
2370	return Parms
2371
2372	def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
2373	'''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
2374
2375	:param dict Layers: dict with following items
2376
2377	::
2378
2379	{'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
2380	'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
2381	'Layers':[],'Stacking':[],'Transitions':[]}
2382
2383	:param str ctrls: controls string to be written on DIFFaX controls.dif file
2384	:param float scale: scale factor
2385	:param dict background: background parameters
2386	:param list limits: min/max 2-theta to be calculated
2387	:param dict inst: instrument parameters dictionary
2388	:param list profile: powder pattern data
2389
2390	Note that parameters all updated in place
2391	'''
2392	import atmdata
2393	path = sys.path
2394	for name in path:
2395	if 'bin' in name:
2396	DIFFaX = name+'/DIFFaX.exe'
2397	print ' Execute ',DIFFaX
2398	break
2399	# make form factor file that DIFFaX wants - atom types are GSASII style
2400	sf = open('data.sfc','w')
2401	sf.write('GSASII special form factor file for DIFFaX\n\n')
2402	atTypes = Layers['AtInfo'].keys()
2403	if 'H' not in atTypes:
2404	atTypes.insert(0,'H')
2405	for atType in atTypes:
2406	if atType == 'H':
2407	blen = -.3741
2408	else:
2409	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2410	Adat = atmdata.XrayFF[atType]
2411	text = '%4s'%(atType.ljust(4))
2412	for i in range(4):
2413	text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
2414	text += '%11.6f%11.6f'%(Adat['fc'],blen)
2415	text += '%3d\n'%(Adat['Z'])
2416	sf.write(text)
2417	sf.close()
2418	#make DIFFaX control.dif file - future use GUI to set some of these flags
2419	cf = open('control.dif','w')
2420	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
2421	x0 = profile[0]
2422	iBeg = np.searchsorted(x0,limits[0])
2423	iFin = np.searchsorted(x0,limits[1])+1
2424	if iFin-iBeg > 20000:
2425	iFin = iBeg+20000
2426	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2427	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2428	cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
2429	else:
2430	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
2431	inst = {'Type':['XSC','XSC',]}
2432	cf.close()
2433	#make DIFFaX data file
2434	df = open('GSASII-DIFFaX.dat','w')
2435	df.write('INSTRUMENTAL\n')
2436	if 'X' in inst['Type'][0]:
2437	df.write('X-RAY\n')
2438	elif 'N' in inst['Type'][0]:
2439	df.write('NEUTRON\n')
2440	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
2441	df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
2442	U = ateln2*inst['U'][1]/10000.
2443	V = ateln2*inst['V'][1]/10000.
2444	W = ateln2*inst['W'][1]/10000.
2445	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
2446	HW = np.sqrt(np.mean(HWHM))
2447	# df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
2448	if 'Mean' in Layers['selInst']:
2449	df.write('GAUSSIAN %.6f TRIM\n'%(HW)) #fast option - might not really matter
2450	elif 'Gaussian' in Layers['selInst']:
2451	df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W)) #slow - make a GUI option?
2452	else:
2453	df.write('None\n')
2454	else:
2455	df.write('0.10\nNone\n')
2456	df.write('STRUCTURAL\n')
2457	a,b,c = Layers['Cell'][1:4]
2458	gam = Layers['Cell'][6]
2459	df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
2460	laue = Layers['Laue']
2461	if laue == '2/m(ab)':
2462	laue = '2/m(1)'
2463	elif laue == '2/m(c)':
2464	laue = '2/m(2)'
2465	if 'unknown' in Layers['Laue']:
2466	df.write('%s %.3f\n'%(laue,Layers['Toler']))
2467	else:
2468	df.write('%s\n'%(laue))
2469	df.write('%d\n'%(len(Layers['Layers'])))
2470	if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
2471	df.write('%.1f %.1f\n'%(Layers['Width'][0][0]10000.,Layers['Width'][0][0]10000.)) #mum to A
2472	layerNames = []
2473	for layer in Layers['Layers']:
2474	layerNames.append(layer['Name'])
2475	for il,layer in enumerate(Layers['Layers']):
2476	if layer['SameAs']:
2477	df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
2478	continue
2479	df.write('LAYER %d\n'%(il+1))
2480	if '-1' in layer['Symm']:
2481	df.write('CENTROSYMMETRIC\n')
2482	else:
2483	df.write('NONE\n')
2484	for ia,atom in enumerate(layer['Atoms']):
2485	[name,atype,x,y,z,frac,Uiso] = atom
2486	if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
2487	frac /= 2.
2488	df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
2489	df.write('STACKING\n')
2490	df.write('%s\n'%(Layers['Stacking'][0]))
2491	if 'recursive' in Layers['Stacking'][0]:
2492	df.write('%s\n'%Layers['Stacking'][1])
2493	else:
2494	if 'list' in Layers['Stacking'][1]:
2495	Slen = len(Layers['Stacking'][2])
2496	iB = 0
2497	iF = 0
2498	while True:
2499	iF += 68
2500	if iF >= Slen:
2501	break
2502	iF = min(iF,Slen)
2503	df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
2504	iB = iF
2505	else:
2506	df.write('%s\n'%Layers['Stacking'][1])
2507	df.write('TRANSITIONS\n')
2508	for iY in range(len(Layers['Layers'])):
2509	sumPx = 0.
2510	for iX in range(len(Layers['Layers'])):
2511	p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
2512	p = round(p,3)
2513	df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
2514	sumPx += p
2515	if sumPx != 1.0: #this has to be picky since DIFFaX is.
2516	print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
2517	df.close()
2518	os.remove('data.sfc')
2519	os.remove('control.dif')
2520	os.remove('GSASII-DIFFaX.dat')
2521	return
2522	df.close()
2523	time0 = time.time()
2524	try:
2525	subp.call(DIFFaX)
2526	except OSError:
2527	print ' DIFFax.exe is not available for this platform - under development'
2528	print ' DIFFaX time = %.2fs'%(time.time()-time0)
2529	if os.path.exists('GSASII-DIFFaX.spc'):
2530	Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
2531	iFin = iBeg+Xpat.shape[1]
2532	bakType,backDict,backVary = SetBackgroundParms(background)
2533	backDict['Lam1'] = G2mth.getWave(inst)
2534	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
2535	profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
2536	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
2537	rv = st.poisson(profile[3][iBeg:iFin])
2538	profile[1][iBeg:iFin] = rv.rvs()
2539	Z = np.ones_like(profile[3][iBeg:iFin])
2540	Z[1::2] *= -1
2541	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2542	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2543	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2544	#cleanup files..
2545	os.remove('GSASII-DIFFaX.spc')
2546	elif os.path.exists('GSASII-DIFFaX.sadp'):
2547	Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
2548	Sadp = np.reshape(Sadp,(256,-1))
2549	Layers['Sadp']['Img'] = Sadp
2550	os.remove('GSASII-DIFFaX.sadp')
2551	os.remove('data.sfc')
2552	os.remove('control.dif')
2553	os.remove('GSASII-DIFFaX.dat')
2554
2555	def SetPWDRscan(inst,limits,profile):
2556
2557	wave = G2mth.getMeanWave(inst)
2558	x0 = profile[0]
2559	iBeg = np.searchsorted(x0,limits[0])
2560	iFin = np.searchsorted(x0,limits[1])
2561	if iFin-iBeg > 20000:
2562	iFin = iBeg+20000
2563	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2564	pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2565	return iFin-iBeg
2566
2567	def SetStackingSF(Layers,debug):
2568	# Load scattering factors into DIFFaX arrays
2569	import atmdata
2570	atTypes = Layers['AtInfo'].keys()
2571	aTypes = []
2572	for atype in atTypes:
2573	aTypes.append('%4s'%(atype.ljust(4)))
2574	SFdat = []
2575	for atType in atTypes:
2576	Adat = atmdata.XrayFF[atType]
2577	SF = np.zeros(9)
2578	SF[:8:2] = Adat['fa']
2579	SF[1:8:2] = Adat['fb']
2580	SF[8] = Adat['fc']
2581	SFdat.append(SF)
2582	SFdat = np.array(SFdat)
2583	pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2584
2585	def SetStackingClay(Layers,Type):
2586	# Controls
2587	rand.seed()
2588	ranSeed = rand.randint(1,2**16-1)
2589	try:
2590	laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2591	'6/m','6/mmm'].index(Layers['Laue'])+1
2592	except ValueError: #for 'unknown'
2593	laueId = -1
2594	if 'SADP' in Type:
2595	planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2596	lmax = int(Layers['Sadp']['Lmax'])
2597	else:
2598	planeId = 0
2599	lmax = 0
2600	# Sequences
2601	StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2602	try:
2603	StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2604	except ValueError:
2605	StkParm = -1
2606	if StkParm == 2: #list
2607	StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2608	Nstk = len(StkSeq)
2609	else:
2610	Nstk = 1
2611	StkSeq = [0,]
2612	if StkParm == -1:
2613	StkParm = int(Layers['Stacking'][1])
2614	Wdth = Layers['Width'][0]
2615	mult = 1
2616	controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2617	LaueSym = Layers['Laue'].ljust(12)
2618	pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2619	return laueId,controls
2620
2621	def SetCellAtoms(Layers):
2622	Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2623	# atoms in layers
2624	atTypes = Layers['AtInfo'].keys()
2625	AtomXOU = []
2626	AtomTp = []
2627	LayerSymm = []
2628	LayerNum = []
2629	layerNames = []
2630	Natm = 0
2631	Nuniq = 0
2632	for layer in Layers['Layers']:
2633	layerNames.append(layer['Name'])
2634	for il,layer in enumerate(Layers['Layers']):
2635	if layer['SameAs']:
2636	LayerNum.append(layerNames.index(layer['SameAs'])+1)
2637	continue
2638	else:
2639	LayerNum.append(il+1)
2640	Nuniq += 1
2641	if '-1' in layer['Symm']:
2642	LayerSymm.append(1)
2643	else:
2644	LayerSymm.append(0)
2645	for ia,atom in enumerate(layer['Atoms']):
2646	[name,atype,x,y,z,frac,Uiso] = atom
2647	Natm += 1
2648	AtomTp.append('%4s'%(atype.ljust(4)))
2649	Ta = atTypes.index(atype)+1
2650	AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2651	AtomXOU = np.array(AtomXOU)
2652	Nlayers = len(layerNames)
2653	pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2654	return Nlayers
2655
2656	def SetStackingTrans(Layers,Nlayers):
2657	# Transitions
2658	TransX = []
2659	TransP = []
2660	for Ytrans in Layers['Transitions']:
2661	TransP.append([trans[0] for trans in Ytrans]) #get just the numbers
2662	TransX.append([trans[1:4] for trans in Ytrans]) #get just the numbers
2663	TransP = np.array(TransP,dtype='float').T
2664	TransX = np.array(TransX,dtype='float')
2665	# GSASIIpath.IPyBreak()
2666	pyx.pygettrans(Nlayers,TransP,TransX)
2667
2668	def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2669	# Scattering factors
2670	SetStackingSF(Layers,debug)
2671	# Controls & sequences
2672	laueId,controls = SetStackingClay(Layers,'PWDR')
2673	# cell & atoms
2674	Nlayers = SetCellAtoms(Layers)
2675	Volume = Layers['Cell'][7]
2676	# Transitions
2677	SetStackingTrans(Layers,Nlayers)
2678	# PWDR scan
2679	Nsteps = SetPWDRscan(inst,limits,profile)
2680	# result as Spec
2681	x0 = profile[0]
2682	profile[3] = np.zeros(len(profile[0]))
2683	profile[4] = np.zeros(len(profile[0]))
2684	profile[5] = np.zeros(len(profile[0]))
2685	iBeg = np.searchsorted(x0,limits[0])
2686	iFin = np.searchsorted(x0,limits[1])+1
2687	if iFin-iBeg > 20000:
2688	iFin = iBeg+20000
2689	Nspec = 20001
2690	spec = np.zeros(Nspec,dtype='double')
2691	time0 = time.time()
2692	pyx.pygetspc(controls,Nspec,spec)
2693	print ' GETSPC time = %.2fs'%(time.time()-time0)
2694	time0 = time.time()
2695	U = ateln2*inst['U'][1]/10000.
2696	V = ateln2*inst['V'][1]/10000.
2697	W = ateln2*inst['W'][1]/10000.
2698	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
2699	HW = np.sqrt(np.mean(HWHM))
2700	BrdSpec = np.zeros(Nsteps)
2701	if 'Mean' in Layers['selInst']:
2702	pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2703	elif 'Gaussian' in Layers['selInst']:
2704	pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2705	else:
2706	BrdSpec = spec[:Nsteps]
2707	BrdSpec /= Volume
2708	iFin = iBeg+Nsteps
2709	bakType,backDict,backVary = SetBackgroundParms(background)
2710	backDict['Lam1'] = G2mth.getWave(inst)
2711	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
2712	profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2713	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
2714	try:
2715	rv = st.poisson(profile[3][iBeg:iFin])
2716	profile[1][iBeg:iFin] = rv.rvs()
2717	except ValueError:
2718	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
2719	Z = np.ones_like(profile[3][iBeg:iFin])
2720	Z[1::2] *= -1
2721	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2722	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2723	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2724	print ' Broadening time = %.2fs'%(time.time()-time0)
2725
2726	def CalcStackingSADP(Layers,debug):
2727
2728	# Scattering factors
2729	SetStackingSF(Layers,debug)
2730	# Controls & sequences
2731	laueId,controls = SetStackingClay(Layers,'SADP')
2732	# cell & atoms
2733	Nlayers = SetCellAtoms(Layers)
2734	# Transitions
2735	SetStackingTrans(Layers,Nlayers)
2736	# result as Sadp
2737	Nspec = 20001
2738	spec = np.zeros(Nspec,dtype='double')
2739	time0 = time.time()
2740	hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2741	Sapd = np.zeros((256,256))
2742	iB = 0
2743	for i in range(hkLim):
2744	iF = iB+Nblk
2745	p1 = 127+int(i*Incr)
2746	p2 = 128-int(i*Incr)
2747	if Nblk == 128:
2748	if i:
2749	Sapd[128:,p1] = spec[iB:iF]
2750	Sapd[:128,p1] = spec[iF:iB:-1]
2751	Sapd[128:,p2] = spec[iB:iF]
2752	Sapd[:128,p2] = spec[iF:iB:-1]
2753	else:
2754	if i:
2755	Sapd[:,p1] = spec[iB:iF]
2756	Sapd[:,p2] = spec[iB:iF]
2757	iB += Nblk
2758	Layers['Sadp']['Img'] = Sapd
2759	print ' GETSAD time = %.2fs'%(time.time()-time0)
2760	# GSASIIpath.IPyBreak()
2761
2762	#testing data
2763	NeedTestData = True
2764	def TestData():
2765	'needs a doc string'
2766	# global NeedTestData
2767	global bakType
2768	bakType = 'chebyschev'
2769	global xdata
2770	xdata = np.linspace(4.0,40.0,36000)
2771	global parmDict0
2772	parmDict0 = {
2773	'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2774	'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2775	'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2776	'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2777	'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2778	'Back0':5.384,'Back1':-0.015,'Back2':.004,
2779	}
2780	global parmDict1
2781	parmDict1 = {
2782	'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2783	'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2784	'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2785	'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2786	'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2787	'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2788	'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2789	'Back0':36.897,'Back1':-0.508,'Back2':.006,
2790	'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2791	}
2792	global parmDict2
2793	parmDict2 = {
2794	'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2795	'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2796	'Back0':5.,'Back1':-0.02,'Back2':.004,
2797	# 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2798	}
2799	global varyList
2800	varyList = []
2801
2802	def test0():
2803	if NeedTestData: TestData()
2804	gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2805	gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])
2806	gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2807	fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2808	fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])
2809	fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2810
2811	def test1():
2812	if NeedTestData: TestData()
2813	time0 = time.time()
2814	for i in range(100):
2815	getPeakProfile(parmDict1,xdata,varyList,bakType)
2816	print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2817
2818	def test2(name,delt):
2819	if NeedTestData: TestData()
2820	varyList = [name,]
2821	xdata = np.linspace(5.6,5.8,400)
2822	hplot = plotter.add('derivatives test for '+name).gca()
2823	hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2824	y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2825	parmDict2[name] += delt
2826	y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2827	hplot.plot(xdata,(y1-y0)/delt,'r+')
2828
2829	def test3(name,delt):
2830	if NeedTestData: TestData()
2831	names = ['pos','sig','gam','shl']
2832	idx = names.index(name)
2833	myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2834	xdata = np.linspace(5.6,5.8,800)
2835	dx = xdata[1]-xdata[0]
2836	hplot = plotter.add('derivatives test for '+name).gca()
2837	hplot.plot(xdata,100.dxgetdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2838	y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2839	myDict[name] += delt
2840	y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2841	hplot.plot(xdata,(y1-y0)/delt,'r+')
2842
2843	if __name__ == '__main__':
2844	import GSASIItestplot as plot
2845	global plotter
2846	plotter = plot.PlotNotebook()
2847	# test0()
2848	# for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2849	for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2850	['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2851	test2(name,shft)
2852	for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2853	test3(name,shft)
2854	print "OK"
2855	plotter.StartEventLoop()

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

Download in other formats: