Context Navigation

source: trunk/GSASIIpwd.py @ 2357

Last change on this file since 2357 was 2357, checked in by vondreele, 9 years ago
remove import .Pnn file cleanup CalcRDF
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 87.7 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: 2016-07-01 19:07:22 +0000 (Fri, 01 Jul 2016) $
10	# $Author: vondreele $
11	# $Revision: 2357 $
12	# $URL: trunk/GSASIIpwd.py $
13	# $Id: GSASIIpwd.py 2357 2016-07-01 19:07:22Z vondreele $
14	########### SVN repository information ###################
15	import sys
16	import math
17	import time
18	import os
19	import subprocess as subp
20
21	import numpy as np
22	import scipy as sp
23	import numpy.linalg as nl
24	import numpy.ma as ma
25	import random as rand
26	from numpy.fft import ifft, fft, fftshift
27	import scipy.interpolate as si
28	import scipy.stats as st
29	import scipy.optimize as so
30
31	import GSASIIpath
32	GSASIIpath.SetVersionNumber("$Revision: 2357 $")
33	import GSASIIlattice as G2lat
34	import GSASIIspc as G2spc
35	import GSASIIElem as G2elem
36	import GSASIIgrid as G2gd
37	import GSASIIIO as G2IO
38	import GSASIImath as G2mth
39	import pypowder as pyd
40	try:
41	import pydiffax as pyx
42	except ImportError:
43	print 'pydiffax is not available for this platform - under develpment'
44
45
46	# trig functions in degrees
47	sind = lambda x: math.sin(x*math.pi/180.)
48	asind = lambda x: 180.*math.asin(x)/math.pi
49	tand = lambda x: math.tan(x*math.pi/180.)
50	atand = lambda x: 180.*math.atan(x)/math.pi
51	atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
52	cosd = lambda x: math.cos(x*math.pi/180.)
53	acosd = lambda x: 180.*math.acos(x)/math.pi
54	rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
55	#numpy versions
56	npsind = lambda x: np.sin(x*np.pi/180.)
57	npasind = lambda x: 180.*np.arcsin(x)/math.pi
58	npcosd = lambda x: np.cos(x*math.pi/180.)
59	npacosd = lambda x: 180.*np.arccos(x)/math.pi
60	nptand = lambda x: np.tan(x*math.pi/180.)
61	npatand = lambda x: 180.*np.arctan(x)/np.pi
62	npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
63	npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
64	npT2q = lambda tth,wave: 2.0np.pinpT2stl(tth,wave)
65	ateln2 = 8.0*math.log(2.0)
66
67	#GSASII pdf calculation routines
68
69	def Transmission(Geometry,Abs,Diam):
70	'''
71	Calculate sample transmission
72
73	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
74	:param float Abs: absorption coeff in cm-1
75	:param float Diam: sample thickness/diameter in mm
76	'''
77	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
78	MuR = Abs*Diam/20.0
79	if MuR <= 3.0:
80	T0 = 16/(3.*math.pi)
81	T1 = -0.045780
82	T2 = -0.02489
83	T3 = 0.003045
84	T = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
85	if T < -20.:
86	return 2.06e-9
87	else:
88	return math.exp(T)
89	else:
90	T1 = 1.433902
91	T2 = 0.013869+0.337894
92	T3 = 1.933433+1.163198
93	T4 = 0.044365-0.04259
94	T = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
95	return T/100.
96	elif 'plate' in Geometry:
97	MuR = Abs*Diam/10.
98	return math.exp(-MuR)
99	elif 'Bragg' in Geometry:
100	return 0.0
101
102	def SurfaceRough(SRA,SRB,Tth):
103	''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
104	:param float SRA: Suortti surface roughness parameter
105	:param float SRB: Suortti surface roughness parameter
106	:param float Tth: 2-theta(deg) - can be numpy array
107
108	'''
109	sth = npsind(Tth/2.)
110	T1 = np.exp(-SRB/sth)
111	T2 = SRA+(1.-SRA)*np.exp(-SRB)
112	return (SRA+(1.-SRA)*T1)/T2
113
114	def SurfaceRoughDerv(SRA,SRB,Tth):
115	''' Suortti surface roughness correction derivatives
116	:param float SRA: Suortti surface roughness parameter (dimensionless)
117	:param float SRB: Suortti surface roughness parameter (dimensionless)
118	:param float Tth: 2-theta(deg) - can be numpy array
119	:return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
120	'''
121	sth = npsind(Tth/2.)
122	T1 = np.exp(-SRB/sth)
123	T2 = SRA+(1.-SRA)*np.exp(-SRB)
124	Trans = (SRA+(1.-SRA)*T1)/T2
125	dydSRA = ((1.-T1)T2-(1.-np.exp(-SRB))Trans)/T2**2
126	dydSRB = ((SRA-1.)T1T2/sth-Trans(SRA-T2))/T2*2
127	return [dydSRA,dydSRB]
128
129	def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
130	'''Calculate sample absorption
131	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
132	:param float MuR: absorption coeff * sample thickness/2 or radius
133	:param Tth: 2-theta scattering angle - can be numpy array
134	:param float Phi: flat plate tilt angle - future
135	:param float Psi: flat plate tilt axis - future
136	'''
137
138	def muRunder3(MuR,Sth2):
139	T0 = 16.0/(3.*np.pi)
140	T1 = (25.99978-0.01911Sth20.25)np.exp(-0.024551*Sth2)+ \
141	0.109561*np.sqrt(Sth2)-26.04556
142	T2 = -0.02489-0.39499Sth2+1.219077Sth2**1.5- \
143	1.31268Sth22+0.871081Sth2*2.5-0.2327Sth2**3
144	T3 = 0.003045+0.018167Sth2-0.03305Sth2**2
145	Trns = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
146	return np.exp(Trns)
147
148	def muRover3(MuR,Sth2):
149	T1 = 1.433902+11.07504Sth2-8.77629Sth2*Sth2+ \
150	10.02088Sth23-3.36778Sth2**4
151	T2 = (0.013869-0.01249Sth2)np.exp(3.27094*Sth2)+ \
152	(0.337894+13.77317Sth2)/(1.0+11.53544Sth2)**1.555039
153	T3 = 1.933433/(1.0+23.12967Sth2)*1.686715- \
154	0.13576*np.sqrt(Sth2)+1.163198
155	T4 = 0.044365-0.04259/(1.0+0.41051Sth2)*148.4202
156	Trns = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
157	return Trns/100.
158
159	Sth2 = npsind(Tth/2.0)**2
160	Cth2 = 1.-Sth2
161	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
162	if 'array' in str(type(MuR)):
163	MuRSTh2 = np.concatenate((MuR,Sth2))
164	AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
165	return AbsCr
166	else:
167	if MuR <= 3.0:
168	return muRunder3(MuR,Sth2)
169	else:
170	return muRover3(MuR,Sth2)
171	elif 'Bragg' in Geometry:
172	return 1.0
173	elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
174	# and only defined for 2theta < 90
175	MuT = 2.*MuR
176	T1 = np.exp(-MuT)
177	T2 = np.exp(-MuT/npcosd(Tth))
178	Tb = MuT-MuT/npcosd(Tth)
179	return (T2-T1)/Tb
180	elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
181	MuT = 2.*MuR
182	cth = npcosd(Tth/2.0)
183	return np.exp(-MuT/cth)/cth
184
185	def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
186	'needs a doc string'
187	dA = 0.001
188	AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
189	if MuR:
190	AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
191	return (AbsP-AbsM)/(2.0*dA)
192	else:
193	return (AbsP-1.)/dA
194
195	def Polarization(Pola,Tth,Azm=0.0):
196	""" Calculate angle dependent x-ray polarization correction (not scaled correctly!)
197
198	:param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
199	:param Azm: azimuthal angle e.g. 0.0 in plane of polarization
200	:param Tth: 2-theta scattering angle - can be numpy array
201	which (if either) of these is "right"?
202	:return: (pola, dpdPola)
203	* pola = ((1-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
204	(1-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
205	* dpdPola: derivative needed for least squares
206
207	"""
208	pola = ((1.0-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
209	(1.0-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
210	dpdPola = -npsind(Tth)*2(npsind(Azm)2-npcosd(Azm)2)
211	return pola,dpdPola
212
213	def Oblique(ObCoeff,Tth):
214	'currently assumes detector is normal to beam'
215	if ObCoeff:
216	return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
217	else:
218	return 1.0
219
220	def Ruland(RulCoff,wave,Q,Compton):
221	'needs a doc string'
222	C = 2.9978e8
223	D = 1.5e-3
224	hmc = 0.024262734687
225	sinth2 = (Qwave/(4.0np.pi))**2
226	dlam = (wave*2)Compton*Q/C
227	dlam_c = 2.0hmcsinth2-Dwave*2
228	return 1.0/((1.0+dlam/RulCoff)(1.0+(np.pidlam_c/(dlam+RulCoff))**2))
229
230	def LorchWeight(Q):
231	'needs a doc string'
232	return np.sin(np.pi(Q[-1]-Q)/(2.0Q[-1]))
233
234	def GetAsfMean(ElList,Sthl2):
235	'''Calculate various scattering factor terms for PDF calcs
236
237	:param dict ElList: element dictionary contains scattering factor coefficients, etc.
238	:param np.array Sthl2: numpy array of sin theta/lambda squared values
239	:returns: mean(f^2), mean(f)^2, mean(compton)
240	'''
241	sumNoAtoms = 0.0
242	FF = np.zeros_like(Sthl2)
243	FF2 = np.zeros_like(Sthl2)
244	CF = np.zeros_like(Sthl2)
245	for El in ElList:
246	sumNoAtoms += ElList[El]['FormulaNo']
247	for El in ElList:
248	el = ElList[El]
249	ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])2+el['fpp']2
250	cf = G2elem.ComptonFac(el,Sthl2)
251	FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
252	FF2 += ff2*el['FormulaNo']/sumNoAtoms
253	CF += cf*el['FormulaNo']/sumNoAtoms
254	return FF2,FF**2,CF
255
256	def GetNumDensity(ElList,Vol):
257	'needs a doc string'
258	sumNoAtoms = 0.0
259	for El in ElList:
260	sumNoAtoms += ElList[El]['FormulaNo']
261	return sumNoAtoms/Vol
262
263	def CalcPDF(data,inst,xydata):
264	'needs a doc string'
265	auxPlot = []
266	import copy
267	import scipy.fftpack as ft
268	#subtract backgrounds - if any
269	xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
270	if data['Sample Bkg.']['Name']:
271	xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
272	data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
273	if data['Container']['Name']:
274	xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
275	if data['Container Bkg.']['Name']:
276	xycontainer += (xydata['Container Bkg.'][1][1]+
277	data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
278	xydata['IofQ'][1][1] += xycontainer
279	#get element data & absorption coeff.
280	ElList = data['ElList']
281	Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
282	#Apply angle dependent corrections
283	Tth = xydata['Sample'][1][0]
284	dt = (Tth[1]-Tth[0])
285	MuR = Abs*data['Diam']/20.0
286	xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
287	if 'X' in inst['Type'][0]:
288	xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
289	if data['DetType'] == 'Image plate':
290	xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
291	XY = xydata['IofQ'][1]
292	#convert to Q
293	if 'C' in inst['Type'][0]:
294	hc = 12.397639
295	wave = G2mth.getWave(inst)
296	keV = hc/wave
297	minQ = npT2q(Tth[0],wave)
298	maxQ = npT2q(Tth[-1],wave)
299	Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
300	dq = Qpoints[1]-Qpoints[0]
301	XY[0] = npT2q(XY[0],wave)
302	elif 'T' in inst['Type'][0]:
303	difC = inst['difC'][1]
304	minQ = 2.np.pidifC/Tth[-1]
305	maxQ = 2.np.pidifC/Tth[0]
306	Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
307	dq = Qpoints[1]-Qpoints[0]
308	XY[0] = 2.np.pidifC/XY[0]
309	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])
310	Qdata -= np.min(Qdata)*data['BackRatio']
311
312	qLimits = data['QScaleLim']
313	minQ = np.searchsorted(Qpoints,qLimits[0])
314	maxQ = np.searchsorted(Qpoints,qLimits[1])
315	newdata = []
316	xydata['IofQ'][1][0] = Qpoints
317	xydata['IofQ'][1][1] = Qdata
318	for item in xydata['IofQ'][1]:
319	newdata.append(item[:maxQ])
320	xydata['IofQ'][1] = newdata
321
322	xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
323	FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0np.pi))*2) #these are <f^2>,<f>^2,Cf
324	Q = xydata['SofQ'][1][0]
325	auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
326	ruland = Ruland(data['Ruland'],wave,Q,CF)
327	auxPlot.append([Q,ruland,'Ruland'])
328	CF *= ruland
329	auxPlot.append([Q,CF,'CF-Corr'])
330	scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
331	xydata['SofQ'][1][1] *= scale
332	xydata['SofQ'][1][1] -= CF
333	xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
334	scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
335	xydata['SofQ'][1][1] *= scale
336	xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
337	xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
338	if data['Lorch']:
339	xydata['FofQ'][1][1] *= LorchWeight(Q)
340	xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
341	nR = len(xydata['GofR'][1][1])
342	xydata['GofR'][1][1] = -dqnp.imag(ft.fft(xydata['FofQ'][1][1],4nR)[:nR])
343	xydata['GofR'][1][0] = 0.5np.pinp.linspace(0,nR,nR)/qLimits[1]
344	return auxPlot
345
346	def MakeRDF(RDFcontrols,background,inst,pwddata):
347	import scipy.fftpack as ft
348	auxPlot = []
349	if 'C' in inst['Type'][0]:
350	Tth = pwddata[0]
351	wave = G2mth.getWave(inst)
352	minQ = npT2q(Tth[0],wave)
353	maxQ = npT2q(Tth[-1],wave)
354	powQ = npT2q(Tth,wave)
355	elif 'T' in inst['Type'][0]:
356	TOF = pwddata[0]
357	difC = inst['difC'][1]
358	minQ = 2.np.pidifC/TOF[-1]
359	maxQ = 2.np.pidifC/TOF[0]
360	powQ = 2.np.pidifC/TOF
361	piDQ = np.pi/(maxQ-minQ)
362	Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
363	if RDFcontrols['UseObsCalc']:
364	Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
365	else:
366	Qdata = si.griddata(powQ,pwddata[1],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
367	Qdata = np.sin((Qpoints-minQ)piDQ)/piDQ
368	Qdata = 0.5np.sqrt(Qpoints) #Qbin normalization
369	auxPlot.append([Qpoints,Qdata,'interp1d:'+RDFcontrols['Smooth']])
370	# GSASIIpath.IPyBreak()
371	dq = Qpoints[1]-Qpoints[0]
372	nR = len(Qdata)
373	DofR = -dqnp.imag(ft.fft(Qdata,16nR)[:nR])
374	R = 0.5np.pinp.linspace(0,nR,nR)/(4.*maxQ)
375	iFin = np.searchsorted(R,RDFcontrols['maxR'])
376	auxPlot.append([R[:iFin],DofR[:iFin],'D(R)'])
377	return auxPlot
378
379	################################################################################
380	#GSASII peak fitting routines: Finger, Cox & Jephcoat model
381	################################################################################
382
383	def factorize(num):
384	''' Provide prime number factors for integer num
385	:returns: dictionary of prime factors (keys) & power for each (data)
386	'''
387	factors = {}
388	orig = num
389
390	# we take advantage of the fact that (i +1)2 = i2 + 2*i +1
391	i, sqi = 2, 4
392	while sqi <= num:
393	while not num%i:
394	num /= i
395	factors[i] = factors.get(i, 0) + 1
396
397	sqi += 2*i + 1
398	i += 1
399
400	if num != 1 and num != orig:
401	factors[num] = factors.get(num, 0) + 1
402
403	if factors:
404	return factors
405	else:
406	return {num:1} #a prime number!
407
408	def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
409	''' Provide list of optimal data sizes for FFT calculations
410
411	:param int nmin: minimum data size >= 1
412	:param int nmax: maximum data size > nmin
413	:param int thresh: maximum prime factor allowed
414	:Returns: list of data sizes where the maximum prime factor is < thresh
415	'''
416	plist = []
417	nmin = max(1,nmin)
418	nmax = max(nmin+1,nmax)
419	for p in range(nmin,nmax):
420	if max(factorize(p).keys()) < thresh:
421	plist.append(p)
422	return plist
423
424	np.seterr(divide='ignore')
425
426	# Normal distribution
427
428	# loc = mu, scale = std
429	_norm_pdf_C = 1./math.sqrt(2*math.pi)
430	class norm_gen(st.rv_continuous):
431	'needs a doc string'
432
433	def pdf(self,x,args,*kwds):
434	loc,scale=kwds['loc'],kwds['scale']
435	x = (x-loc)/scale
436	return np.exp(-x*2/2.0) _norm_pdf_C / scale
437
438	norm = norm_gen(name='norm',longname='A normal',extradoc="""
439
440	Normal distribution
441
442	The location (loc) keyword specifies the mean.
443	The scale (scale) keyword specifies the standard deviation.
444
445	normal.pdf(x) = exp(-x*2/2)/sqrt(2pi)
446	""")
447
448	## Cauchy
449
450	# median = loc
451
452	class cauchy_gen(st.rv_continuous):
453	'needs a doc string'
454
455	def pdf(self,x,args,*kwds):
456	loc,scale=kwds['loc'],kwds['scale']
457	x = (x-loc)/scale
458	return 1.0/np.pi/(1.0+x*x) / scale
459
460	cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
461
462	Cauchy distribution
463
464	cauchy.pdf(x) = 1/(pi(1+x*2))
465
466	This is the t distribution with one degree of freedom.
467	""")
468
469
470	#GSASII peak fitting routine: Finger, Cox & Jephcoat model
471
472
473	class fcjde_gen(st.rv_continuous):
474	"""
475	Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
476	Ref: J. Appl. Cryst. (1994) 27, 892-900.
477
478	:param x: array -1 to 1
479	:param t: 2-theta position of peak
480	:param s: sum(S/L,H/L); S: sample height, H: detector opening,
481	L: sample to detector opening distance
482	:param dx: 2-theta step size in deg
483
484	:returns: for fcj.pdf
485
486	* T = x*dx+t
487	* s = S/L+H/L
488	* if x < 0::
489
490	fcj.pdf = [1/sqrt({cos(T)2/cos(t)2}-1) - 1/s]/\|cos(T)\|
491
492	* if x >= 0: fcj.pdf = 0
493	"""
494	def _pdf(self,x,t,s,dx):
495	T = dx*x+t
496	ax2 = abs(npcosd(T))
497	ax = ax2**2
498	bx = npcosd(t)**2
499	bx = np.where(ax>bx,bx,ax)
500	fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
501	fx = np.where(fx > 0.,fx,0.0)
502	return fx
503
504	def pdf(self,x,args,*kwds):
505	loc=kwds['loc']
506	return self._pdf(x-loc,*args)
507
508	fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
509
510	def getWidthsCW(pos,sig,gam,shl):
511	'''Compute the peak widths used for computing the range of a peak
512	for constant wavelength data. On low-angle side, 10 FWHM are used,
513	on high-angle side 15 are used, low angle side extended by axial divergence
514	(for peaks above 90 deg, these are reversed.)
515	'''
516	widths = [np.sqrt(sig)/100.,gam/100.]
517	fwhm = 2.355*widths[0]+widths[1]
518	fmin = 10.(fwhm+shlabs(npcosd(pos)))
519	fmax = 15.0*fwhm
520	if pos > 90:
521	fmin,fmax = [fmax,fmin]
522	return widths,fmin,fmax
523
524	def getWidthsTOF(pos,alp,bet,sig,gam):
525	'''Compute the peak widths used for computing the range of a peak
526	for constant wavelength data. 10 FWHM are used on both sides each
527	extended by exponential coeff.
528	'''
529	lnf = 3.3 # =log(0.001/2)
530	widths = [np.sqrt(sig),gam]
531	fwhm = 2.355widths[0]+2.widths[1]
532	fmin = 10.fwhm(1.+1./alp)
533	fmax = 10.fwhm(1.+1./bet)
534	return widths,fmin,fmax
535
536	def getFWHM(pos,Inst):
537	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
538	via getgamFW(g,s).
539
540	:param pos: float peak position in deg 2-theta or tof in musec
541	:param Inst: dict instrument parameters
542
543	:returns float: total FWHM of pseudoVoigt in deg or musec
544	'''
545
546	sig = lambda Th,U,V,W: np.sqrt(max(0.001,Utand(Th)2+Vtand(Th)+W))
547	sigTOF = lambda dsp,S0,S1,S2,Sq: S0+S1dsp2+S2dsp4+Sq/dsp2
548	gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))
549	gamTOF = lambda dsp,X,Y: Xdsp+Ydsp**2
550	if 'C' in Inst['Type'][0]:
551	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
552	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])
553	return getgamFW(g,s)/100. #returns FWHM in deg
554	else:
555	dsp = pos/Inst['difC'][0]
556	s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
557	g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
558	return getgamFW(g,s)
559
560	def getgamFW(g,s):
561	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
562	lambda fxn needs FWHM for both Gaussian & Lorentzian components
563
564	:param g: float Lorentzian gamma = FWHM(L)
565	:param s: float Gaussian sig
566
567	:returns float: total FWHM of pseudoVoigt
568	'''
569	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.)
570	return gamFW(2.35482s,g) #sqrt(8ln2)sig = FWHM(G)
571
572	def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
573	'needs a doc string'
574	DX = xdata[1]-xdata[0]
575	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
576	x = np.linspace(pos-fmin,pos+fmin,256)
577	dx = x[1]-x[0]
578	Norm = norm.pdf(x,loc=pos,scale=widths[0])
579	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
580	arg = [pos,shl/57.2958,dx,]
581	FCJ = fcjde.pdf(x,*arg,loc=pos)
582	if len(np.nonzero(FCJ)[0])>5:
583	z = np.column_stack([Norm,Cauchy,FCJ]).T
584	Z = fft(z)
585	Df = ifft(Z.prod(axis=0)).real
586	else:
587	z = np.column_stack([Norm,Cauchy]).T
588	Z = fft(z)
589	Df = fftshift(ifft(Z.prod(axis=0))).real
590	Df /= np.sum(Df)
591	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
592	return intensDf(xdata)DX/dx
593
594	def getBackground(pfx,parmDict,bakType,dataType,xdata):
595	'needs a doc string'
596	if 'T' in dataType:
597	q = 2.np.piparmDict[pfx+'difC']/xdata
598	elif 'C' in dataType:
599	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
600	q = 2.np.pinpsind(xdata/2.)/wave
601	yb = np.zeros_like(xdata)
602	nBak = 0
603	cw = np.diff(xdata)
604	cw = np.append(cw,cw[-1])
605	sumBk = [0.,0.,0]
606	while True:
607	key = pfx+'Back;'+str(nBak)
608	if key in parmDict:
609	nBak += 1
610	else:
611	break
612	#empirical functions
613	if bakType in ['chebyschev','cosine']:
614	dt = xdata[-1]-xdata[0]
615	for iBak in range(nBak):
616	key = pfx+'Back;'+str(iBak)
617	if bakType == 'chebyschev':
618	ybi = parmDict[key](2.(xdata-xdata[0])/dt-1.)**iBak
619	elif bakType == 'cosine':
620	ybi = parmDict[key]npcosd(xdataiBak)
621	yb += ybi
622	sumBk[0] = np.sum(yb)
623	elif bakType in ['Q^2 power series','Q^-2 power series']:
624	QT = 1.
625	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
626	for iBak in range(nBak-1):
627	key = pfx+'Back;'+str(iBak+1)
628	if '-2' in bakType:
629	QT = (iBak+1)q**-2
630	else:
631	QT = q*2/(iBak+1)
632	yb += QT*parmDict[key]
633	sumBk[0] = np.sum(yb)
634	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
635	if nBak == 1:
636	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
637	elif nBak == 2:
638	dX = xdata[-1]-xdata[0]
639	T2 = (xdata-xdata[0])/dX
640	T1 = 1.0-T2
641	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
642	else:
643	if bakType == 'lin interpolate':
644	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
645	elif bakType == 'inv interpolate':
646	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
647	elif bakType == 'log interpolate':
648	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
649	bakPos[0] = xdata[0]
650	bakPos[-1] = xdata[-1]
651	bakVals = np.zeros(nBak)
652	for i in range(nBak):
653	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
654	bakInt = si.interp1d(bakPos,bakVals,'linear')
655	yb = bakInt(ma.getdata(xdata))
656	sumBk[0] = np.sum(yb)
657	#Debye function
658	if pfx+'difC' in parmDict:
659	ff = 1.
660	else:
661	try:
662	wave = parmDict[pfx+'Lam']
663	except KeyError:
664	wave = parmDict[pfx+'Lam1']
665	SQ = (q/(4.np.pi))*2
666	FF = G2elem.GetFormFactorCoeff('Si')[0]
667	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
668	iD = 0
669	while True:
670	try:
671	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
672	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
673	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
674	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
675	yb += ybi
676	sumBk[1] += np.sum(ybi)
677	iD += 1
678	except KeyError:
679	break
680	#peaks
681	iD = 0
682	while True:
683	try:
684	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
685	pkI = parmDict[pfx+'BkPkint;'+str(iD)]
686	pkS = parmDict[pfx+'BkPksig;'+str(iD)]
687	pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
688	if 'C' in dataType:
689	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
690	else: #'T'OF
691	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
692	iBeg = np.searchsorted(xdata,pkP-fmin)
693	iFin = np.searchsorted(xdata,pkP+fmax)
694	lenX = len(xdata)
695	if not iBeg:
696	iFin = np.searchsorted(xdata,pkP+fmax)
697	elif iBeg == lenX:
698	iFin = iBeg
699	else:
700	iFin = np.searchsorted(xdata,pkP+fmax)
701	if 'C' in dataType:
702	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
703	yb[iBeg:iFin] += ybi
704	else: #'T'OF
705	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
706	yb[iBeg:iFin] += ybi
707	sumBk[2] += np.sum(ybi)
708	iD += 1
709	except KeyError:
710	break
711	except ValueError:
712	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
713	break
714	return yb,sumBk
715
716	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
717	'needs a doc string'
718	if 'T' in dataType:
719	q = 2.np.piparmDict[hfx+'difC']/xdata
720	elif 'C' in dataType:
721	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
722	q = 2.np.pinpsind(xdata/2.)/wave
723	nBak = 0
724	while True:
725	key = hfx+'Back;'+str(nBak)
726	if key in parmDict:
727	nBak += 1
728	else:
729	break
730	dydb = np.zeros(shape=(nBak,len(xdata)))
731	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
732	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
733	cw = np.diff(xdata)
734	cw = np.append(cw,cw[-1])
735
736	if bakType in ['chebyschev','cosine']:
737	dt = xdata[-1]-xdata[0]
738	for iBak in range(nBak):
739	if bakType == 'chebyschev':
740	dydb[iBak] = (2.(xdata-xdata[0])/dt-1.)*iBak
741	elif bakType == 'cosine':
742	dydb[iBak] = npcosd(xdata*iBak)
743	elif bakType in ['Q^2 power series','Q^-2 power series']:
744	QT = 1.
745	dydb[0] = np.ones_like(xdata)
746	for iBak in range(nBak-1):
747	if '-2' in bakType:
748	QT = (iBak+1)q**-2
749	else:
750	QT = q*2/(iBak+1)
751	dydb[iBak+1] = QT
752	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
753	if nBak == 1:
754	dydb[0] = np.ones_like(xdata)
755	elif nBak == 2:
756	dX = xdata[-1]-xdata[0]
757	T2 = (xdata-xdata[0])/dX
758	T1 = 1.0-T2
759	dydb = [T1,T2]
760	else:
761	if bakType == 'lin interpolate':
762	bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
763	elif bakType == 'inv interpolate':
764	bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
765	elif bakType == 'log interpolate':
766	bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
767	bakPos[0] = xdata[0]
768	bakPos[-1] = xdata[-1]
769	for i,pos in enumerate(bakPos):
770	if i == 0:
771	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
772	elif i == len(bakPos)-1:
773	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
774	else:
775	dydb[i] = np.where(xdata>bakPos[i],
776	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
777	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
778	if hfx+'difC' in parmDict:
779	ff = 1.
780	else:
781	try:
782	wave = parmDict[hfx+'Lam']
783	except KeyError:
784	wave = parmDict[hfx+'Lam1']
785	q = 4.0np.pinpsind(xdata/2.0)/wave
786	SQ = (q/(4np.pi))*2
787	FF = G2elem.GetFormFactorCoeff('Si')[0]
788	ff = np.array(G2elem.ScatFac(FF,SQ)[0])
789	iD = 0
790	while True:
791	try:
792	if hfx+'difC' in parmDict:
793	q = 2np.piparmDict[hfx+'difC']/xdata
794	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
795	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
796	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
797	sqr = np.sin(qdbR)/(qdbR)
798	cqr = np.cos(q*dbR)
799	temp = np.exp(-dbUq*2)
800	dyddb[3iD] = ffsqrtemp/(np.picw)
801	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(np.pidbRcw)
802	dyddb[3iD+2] = -ffdbAsqrtempq2/(np.picw)
803	iD += 1 #ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
804	except KeyError:
805	break
806	iD = 0
807	while True:
808	try:
809	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
810	pkI = parmDict[hfx+'BkPkint;'+str(iD)]
811	pkS = parmDict[hfx+'BkPksig;'+str(iD)]
812	pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
813	if 'C' in dataType:
814	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
815	else: #'T'OF
816	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
817	iBeg = np.searchsorted(xdata,pkP-fmin)
818	iFin = np.searchsorted(xdata,pkP+fmax)
819	lenX = len(xdata)
820	if not iBeg:
821	iFin = np.searchsorted(xdata,pkP+fmax)
822	elif iBeg == lenX:
823	iFin = iBeg
824	else:
825	iFin = np.searchsorted(xdata,pkP+fmax)
826	if 'C' in dataType:
827	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
828	dydpk[4iD][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdp
829	dydpk[4iD+1][iBeg:iFin] += 100.cw[iBeg:iFin]*Df
830	dydpk[4iD+2][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFds
831	dydpk[4iD+3][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdg
832	else: #'T'OF
833	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
834	dydpk[4iD][iBeg:iFin] += pkIdFdp
835	dydpk[4*iD+1][iBeg:iFin] += Df
836	dydpk[4iD+2][iBeg:iFin] += pkIdFds
837	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
838	iD += 1
839	except KeyError:
840	break
841	except ValueError:
842	print '** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **'
843	break
844	return dydb,dyddb,dydpk
845
846	#use old fortran routine
847	def getFCJVoigt3(pos,sig,gam,shl,xdata):
848	'needs a doc string'
849
850	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
851	# Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
852	Df /= np.sum(Df)
853	return Df
854
855	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
856	'needs a doc string'
857
858	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
859	# Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
860	sumDf = np.sum(Df)
861	return Df,dFdp,dFds,dFdg,dFdsh
862
863	def getPsVoigt(pos,sig,gam,xdata):
864	'needs a doc string'
865
866	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
867	Df /= np.sum(Df)
868	return Df
869
870	def getdPsVoigt(pos,sig,gam,xdata):
871	'needs a doc string'
872
873	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
874	sumDf = np.sum(Df)
875	return Df,dFdp,dFds,dFdg
876
877	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
878	'needs a doc string'
879	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
880	Df /= np.sum(Df)
881	return Df
882
883	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
884	'needs a doc string'
885	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
886	sumDf = np.sum(Df)
887	return Df,dFdp,dFda,dFdb,dFds,dFdg
888
889	def ellipseSize(H,Sij,GB):
890	'needs a doc string'
891	HX = np.inner(H.T,GB)
892	lenHX = np.sqrt(np.sum(HX**2))
893	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
894	R = np.inner(HX/lenHX,Rsize)*Esize #want column length for hkl in crystal
895	lenR = np.sqrt(np.sum(R**2))
896	return lenR
897
898	def ellipseSizeDerv(H,Sij,GB):
899	'needs a doc string'
900	lenR = ellipseSize(H,Sij,GB)
901	delt = 0.001
902	dRdS = np.zeros(6)
903	for i in range(6):
904	Sij[i] -= delt
905	lenM = ellipseSize(H,Sij,GB)
906	Sij[i] += 2.*delt
907	lenP = ellipseSize(H,Sij,GB)
908	Sij[i] -= delt
909	dRdS[i] = (lenP-lenM)/(2.*delt)
910	return lenR,dRdS
911
912	def getHKLpeak(dmin,SGData,A,Inst=None):
913	'needs a doc string'
914	HKL = G2lat.GenHLaue(dmin,SGData,A)
915	HKLs = []
916	for h,k,l,d in HKL:
917	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
918	if not ext:
919	if Inst == None:
920	HKLs.append([h,k,l,d,0,-1])
921	else:
922	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
923	return HKLs
924
925	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
926	'needs a doc string'
927	HKLs = []
928	vec = np.array(Vec)
929	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
930	dvec = 1./(maxH*vstar+1./dmin)
931	HKL = G2lat.GenHLaue(dvec,SGData,A)
932	SSdH = [vec*h for h in range(-maxH,maxH+1)]
933	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
934	for h,k,l,d in HKL:
935	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
936	if not ext and d >= dmin:
937	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
938	for dH in SSdH:
939	if dH:
940	DH = SSdH[dH]
941	H = [h+DH[0],k+DH[1],l+DH[2]]
942	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
943	if d >= dmin:
944	HKLM = np.array([h,k,l,dH])
945	if G2spc.checkSSextc(HKLM,SSGData):
946	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
947	# GSASIIpath.IPyBreak()
948	return G2lat.sortHKLd(HKLs,True,True,True)
949
950	def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
951	'needs a doc string'
952
953	yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
954	yc = np.zeros_like(yb)
955	cw = np.diff(xdata)
956	cw = np.append(cw,cw[-1])
957	if 'C' in dataType:
958	shl = max(parmDict['SH/L'],0.002)
959	Ka2 = False
960	if 'Lam1' in parmDict.keys():
961	Ka2 = True
962	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
963	kRatio = parmDict['I(L2)/I(L1)']
964	iPeak = 0
965	while True:
966	try:
967	pos = parmDict['pos'+str(iPeak)]
968	tth = (pos-parmDict['Zero'])
969	intens = parmDict['int'+str(iPeak)]
970	sigName = 'sig'+str(iPeak)
971	if sigName in varyList:
972	sig = parmDict[sigName]
973	else:
974	sig = G2mth.getCWsig(parmDict,tth)
975	sig = max(sig,0.001) #avoid neg sigma^2
976	gamName = 'gam'+str(iPeak)
977	if gamName in varyList:
978	gam = parmDict[gamName]
979	else:
980	gam = G2mth.getCWgam(parmDict,tth)
981	gam = max(gam,0.001) #avoid neg gamma
982	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
983	iBeg = np.searchsorted(xdata,pos-fmin)
984	iFin = np.searchsorted(xdata,pos+fmin)
985	if not iBeg+iFin: #peak below low limit
986	iPeak += 1
987	continue
988	elif not iBeg-iFin: #peak above high limit
989	return yb+yc
990	yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
991	if Ka2:
992	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
993	iBeg = np.searchsorted(xdata,pos2-fmin)
994	iFin = np.searchsorted(xdata,pos2+fmin)
995	if iBeg-iFin:
996	yc[iBeg:iFin] += intenskRatiogetFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
997	iPeak += 1
998	except KeyError: #no more peaks to process
999	return yb+yc
1000	else:
1001	Pdabc = parmDict['Pdabc']
1002	difC = parmDict['difC']
1003	iPeak = 0
1004	while True:
1005	try:
1006	pos = parmDict['pos'+str(iPeak)]
1007	tof = pos-parmDict['Zero']
1008	dsp = tof/difC
1009	intens = parmDict['int'+str(iPeak)]
1010	alpName = 'alp'+str(iPeak)
1011	if alpName in varyList:
1012	alp = parmDict[alpName]
1013	else:
1014	if len(Pdabc):
1015	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1016	else:
1017	alp = G2mth.getTOFalpha(parmDict,dsp)
1018	alp = max(0.0001,alp)
1019	betName = 'bet'+str(iPeak)
1020	if betName in varyList:
1021	bet = parmDict[betName]
1022	else:
1023	if len(Pdabc):
1024	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1025	else:
1026	bet = G2mth.getTOFbeta(parmDict,dsp)
1027	bet = max(0.0001,bet)
1028	sigName = 'sig'+str(iPeak)
1029	if sigName in varyList:
1030	sig = parmDict[sigName]
1031	else:
1032	sig = G2mth.getTOFsig(parmDict,dsp)
1033	gamName = 'gam'+str(iPeak)
1034	if gamName in varyList:
1035	gam = parmDict[gamName]
1036	else:
1037	gam = G2mth.getTOFgamma(parmDict,dsp)
1038	gam = max(gam,0.001) #avoid neg gamma
1039	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1040	iBeg = np.searchsorted(xdata,pos-fmin)
1041	iFin = np.searchsorted(xdata,pos+fmax)
1042	lenX = len(xdata)
1043	if not iBeg:
1044	iFin = np.searchsorted(xdata,pos+fmax)
1045	elif iBeg == lenX:
1046	iFin = iBeg
1047	else:
1048	iFin = np.searchsorted(xdata,pos+fmax)
1049	if not iBeg+iFin: #peak below low limit
1050	iPeak += 1
1051	continue
1052	elif not iBeg-iFin: #peak above high limit
1053	return yb+yc
1054	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1055	iPeak += 1
1056	except KeyError: #no more peaks to process
1057	return yb+yc
1058
1059	def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1060	'needs a doc string'
1061	# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1062	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1063	dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1064	if 'Back;0' in varyList: #background derivs are in front if present
1065	dMdv[0:len(dMdb)] = dMdb
1066	names = ['DebyeA','DebyeR','DebyeU']
1067	for name in varyList:
1068	if 'Debye' in name:
1069	parm,id = name.split(';')
1070	ip = names.index(parm)
1071	dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
1072	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1073	for name in varyList:
1074	if 'BkPk' in name:
1075	parm,id = name.split(';')
1076	ip = names.index(parm)
1077	dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1078	cw = np.diff(xdata)
1079	cw = np.append(cw,cw[-1])
1080	if 'C' in dataType:
1081	shl = max(parmDict['SH/L'],0.002)
1082	Ka2 = False
1083	if 'Lam1' in parmDict.keys():
1084	Ka2 = True
1085	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1086	kRatio = parmDict['I(L2)/I(L1)']
1087	iPeak = 0
1088	while True:
1089	try:
1090	pos = parmDict['pos'+str(iPeak)]
1091	tth = (pos-parmDict['Zero'])
1092	intens = parmDict['int'+str(iPeak)]
1093	sigName = 'sig'+str(iPeak)
1094	if sigName in varyList:
1095	sig = parmDict[sigName]
1096	dsdU = dsdV = dsdW = 0
1097	else:
1098	sig = G2mth.getCWsig(parmDict,tth)
1099	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1100	sig = max(sig,0.001) #avoid neg sigma
1101	gamName = 'gam'+str(iPeak)
1102	if gamName in varyList:
1103	gam = parmDict[gamName]
1104	dgdX = dgdY = 0
1105	else:
1106	gam = G2mth.getCWgam(parmDict,tth)
1107	dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1108	gam = max(gam,0.001) #avoid neg gamma
1109	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1110	iBeg = np.searchsorted(xdata,pos-fmin)
1111	iFin = np.searchsorted(xdata,pos+fmin)
1112	if not iBeg+iFin: #peak below low limit
1113	iPeak += 1
1114	continue
1115	elif not iBeg-iFin: #peak above high limit
1116	break
1117	dMdpk = np.zeros(shape=(6,len(xdata)))
1118	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1119	for i in range(1,5):
1120	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intens*dMdipk[i]
1121	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk[0]
1122	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1123	if Ka2:
1124	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1125	iBeg = np.searchsorted(xdata,pos2-fmin)
1126	iFin = np.searchsorted(xdata,pos2+fmin)
1127	if iBeg-iFin:
1128	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1129	for i in range(1,5):
1130	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intenskRatiodMdipk2[i]
1131	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]kRatio*dMdipk2[0]
1132	dMdpk[5][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk2[0]
1133	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1134	for parmName in ['pos','int','sig','gam']:
1135	try:
1136	idx = varyList.index(parmName+str(iPeak))
1137	dMdv[idx] = dervDict[parmName]
1138	except ValueError:
1139	pass
1140	if 'U' in varyList:
1141	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1142	if 'V' in varyList:
1143	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1144	if 'W' in varyList:
1145	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1146	if 'X' in varyList:
1147	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1148	if 'Y' in varyList:
1149	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1150	if 'SH/L' in varyList:
1151	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1152	if 'I(L2)/I(L1)' in varyList:
1153	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1154	iPeak += 1
1155	except KeyError: #no more peaks to process
1156	break
1157	else:
1158	Pdabc = parmDict['Pdabc']
1159	difC = parmDict['difC']
1160	iPeak = 0
1161	while True:
1162	try:
1163	pos = parmDict['pos'+str(iPeak)]
1164	tof = pos-parmDict['Zero']
1165	dsp = tof/difC
1166	intens = parmDict['int'+str(iPeak)]
1167	alpName = 'alp'+str(iPeak)
1168	if alpName in varyList:
1169	alp = parmDict[alpName]
1170	else:
1171	if len(Pdabc):
1172	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1173	dad0 = 0
1174	else:
1175	alp = G2mth.getTOFalpha(parmDict,dsp)
1176	dada0 = G2mth.getTOFalphaDeriv(dsp)
1177	betName = 'bet'+str(iPeak)
1178	if betName in varyList:
1179	bet = parmDict[betName]
1180	else:
1181	if len(Pdabc):
1182	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1183	dbdb0 = dbdb1 = dbdb2 = 0
1184	else:
1185	bet = G2mth.getTOFbeta(parmDict,dsp)
1186	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1187	sigName = 'sig'+str(iPeak)
1188	if sigName in varyList:
1189	sig = parmDict[sigName]
1190	dsds0 = dsds1 = dsds2 = dsds3 = 0
1191	else:
1192	sig = G2mth.getTOFsig(parmDict,dsp)
1193	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1194	gamName = 'gam'+str(iPeak)
1195	if gamName in varyList:
1196	gam = parmDict[gamName]
1197	dsdX = dsdY = 0
1198	else:
1199	gam = G2mth.getTOFgamma(parmDict,dsp)
1200	dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1201	gam = max(gam,0.001) #avoid neg gamma
1202	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1203	iBeg = np.searchsorted(xdata,pos-fmin)
1204	lenX = len(xdata)
1205	if not iBeg:
1206	iFin = np.searchsorted(xdata,pos+fmax)
1207	elif iBeg == lenX:
1208	iFin = iBeg
1209	else:
1210	iFin = np.searchsorted(xdata,pos+fmax)
1211	if not iBeg+iFin: #peak below low limit
1212	iPeak += 1
1213	continue
1214	elif not iBeg-iFin: #peak above high limit
1215	break
1216	dMdpk = np.zeros(shape=(7,len(xdata)))
1217	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1218	for i in range(1,6):
1219	dMdpk[i][iBeg:iFin] += intenscw[iBeg:iFin]dMdipk[i]
1220	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1221	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1222	for parmName in ['pos','int','alp','bet','sig','gam']:
1223	try:
1224	idx = varyList.index(parmName+str(iPeak))
1225	dMdv[idx] = dervDict[parmName]
1226	except ValueError:
1227	pass
1228	if 'alpha' in varyList:
1229	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1230	if 'beta-0' in varyList:
1231	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1232	if 'beta-1' in varyList:
1233	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1234	if 'beta-q' in varyList:
1235	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1236	if 'sig-0' in varyList:
1237	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1238	if 'sig-1' in varyList:
1239	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1240	if 'sig-2' in varyList:
1241	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1242	if 'sig-q' in varyList:
1243	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1244	if 'X' in varyList:
1245	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1246	if 'Y' in varyList:
1247	dMdv[varyList.index('Y')] += dsdY*dervDict['gam'] #problem here
1248	iPeak += 1
1249	except KeyError: #no more peaks to process
1250	break
1251	return dMdv
1252
1253	def Dict2Values(parmdict, varylist):
1254	'''Use before call to leastsq to setup list of values for the parameters
1255	in parmdict, as selected by key in varylist'''
1256	return [parmdict[key] for key in varylist]
1257
1258	def Values2Dict(parmdict, varylist, values):
1259	''' Use after call to leastsq to update the parameter dictionary with
1260	values corresponding to keys in varylist'''
1261	parmdict.update(zip(varylist,values))
1262
1263	def SetBackgroundParms(Background):
1264	'needs a doc string'
1265	if len(Background) == 1: # fix up old backgrounds
1266	BackGround.append({'nDebye':0,'debyeTerms':[]})
1267	bakType,bakFlag = Background[0][:2]
1268	backVals = Background[0][3:]
1269	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1270	Debye = Background[1] #also has background peaks stuff
1271	backDict = dict(zip(backNames,backVals))
1272	backVary = []
1273	if bakFlag:
1274	backVary = backNames
1275
1276	backDict['nDebye'] = Debye['nDebye']
1277	debyeDict = {}
1278	debyeList = []
1279	for i in range(Debye['nDebye']):
1280	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1281	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1282	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1283	debyeVary = []
1284	for item in debyeList:
1285	if item[1]:
1286	debyeVary.append(item[0])
1287	backDict.update(debyeDict)
1288	backVary += debyeVary
1289
1290	backDict['nPeaks'] = Debye['nPeaks']
1291	peaksDict = {}
1292	peaksList = []
1293	for i in range(Debye['nPeaks']):
1294	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1295	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1296	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1297	peaksVary = []
1298	for item in peaksList:
1299	if item[1]:
1300	peaksVary.append(item[0])
1301	backDict.update(peaksDict)
1302	backVary += peaksVary
1303	return bakType,backDict,backVary
1304
1305	def DoCalibInst(IndexPeaks,Inst):
1306
1307	def SetInstParms():
1308	dataType = Inst['Type'][0]
1309	insVary = []
1310	insNames = []
1311	insVals = []
1312	for parm in Inst:
1313	insNames.append(parm)
1314	insVals.append(Inst[parm][1])
1315	if parm in ['Lam','difC','difA','difB','Zero',]:
1316	if Inst[parm][2]:
1317	insVary.append(parm)
1318	instDict = dict(zip(insNames,insVals))
1319	return dataType,instDict,insVary
1320
1321	def GetInstParms(parmDict,Inst,varyList):
1322	for name in Inst:
1323	Inst[name][1] = parmDict[name]
1324
1325	def InstPrint(Inst,sigDict):
1326	print 'Instrument Parameters:'
1327	if 'C' in Inst['Type'][0]:
1328	ptfmt = "%12.6f"
1329	else:
1330	ptfmt = "%12.3f"
1331	ptlbls = 'names :'
1332	ptstr = 'values:'
1333	sigstr = 'esds :'
1334	for parm in Inst:
1335	if parm in ['Lam','difC','difA','difB','Zero',]:
1336	ptlbls += "%s" % (parm.center(12))
1337	ptstr += ptfmt % (Inst[parm][1])
1338	if parm in sigDict:
1339	sigstr += ptfmt % (sigDict[parm])
1340	else:
1341	sigstr += 12*' '
1342	print ptlbls
1343	print ptstr
1344	print sigstr
1345
1346	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1347	parmDict.update(zip(varyList,values))
1348	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1349
1350	peakPos = []
1351	peakDsp = []
1352	peakWt = []
1353	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1354	if peak[2] and peak[3] and sig > 0.:
1355	peakPos.append(peak[0])
1356	peakDsp.append(peak[-1]) #d-calc
1357	peakWt.append(peak[-1]2/sig2) #weight by d**2
1358	peakPos = np.array(peakPos)
1359	peakDsp = np.array(peakDsp)
1360	peakWt = np.array(peakWt)
1361	dataType,insDict,insVary = SetInstParms()
1362	parmDict = {}
1363	parmDict.update(insDict)
1364	varyList = insVary
1365	if not len(varyList):
1366	print '** ERROR - nothing to refine! **'
1367	return False
1368	while True:
1369	begin = time.time()
1370	values = np.array(Dict2Values(parmDict, varyList))
1371	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1372	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1373	ncyc = int(result[2]['nfev']/2)
1374	runtime = time.time()-begin
1375	chisq = np.sum(result[2]['fvec']**2)
1376	Values2Dict(parmDict, varyList, result[0])
1377	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
1378	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1379	print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1380	print 'chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF)
1381	try:
1382	sig = np.sqrt(np.diag(result[1])*GOF)
1383	if np.any(np.isnan(sig)):
1384	print '* Least squares aborted - some invalid esds possible *'
1385	break #refinement succeeded - finish up!
1386	except ValueError: #result[1] is None on singular matrix
1387	print '** Refinement failed - singular matrix **'
1388
1389	sigDict = dict(zip(varyList,sig))
1390	GetInstParms(parmDict,Inst,varyList)
1391	InstPrint(Inst,sigDict)
1392	return True
1393
1394	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1395	'''Called to perform a peak fit, refining the selected items in the peak
1396	table as well as selected items in the background.
1397
1398	:param str FitPgm: type of fit to perform. At present this is ignored.
1399	:param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1400	four values followed by a refine flag where the values are: position, intensity,
1401	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1402	tree entry, dict item "peaks"
1403	:param list Background: describes the background. List with two items.
1404	Item 0 specifies a background model and coefficients. Item 1 is a dict.
1405	From the Histogram/Background tree entry.
1406	:param list Limits: min and max x-value to use
1407	:param dict Inst: Instrument parameters
1408	:param dict Inst2: more Instrument parameters
1409	:param numpy.array data: a 5xn array. data[0] is the x-values,
1410	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1411	calc, background and difference intensities, respectively.
1412	:param list prevVaryList: Used in sequential refinements to override the
1413	variable list. Defaults as an empty list.
1414	:param bool oneCycle: True if only one cycle of fitting should be performed
1415	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1416	and derivType = controls['deriv type']. If None default values are used.
1417	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1418	Defaults to None, which means no updates are done.
1419	'''
1420	def GetBackgroundParms(parmList,Background):
1421	iBak = 0
1422	while True:
1423	try:
1424	bakName = 'Back;'+str(iBak)
1425	Background[0][iBak+3] = parmList[bakName]
1426	iBak += 1
1427	except KeyError:
1428	break
1429	iDb = 0
1430	while True:
1431	names = ['DebyeA;','DebyeR;','DebyeU;']
1432	try:
1433	for i,name in enumerate(names):
1434	val = parmList[name+str(iDb)]
1435	Background[1]['debyeTerms'][iDb][2*i] = val
1436	iDb += 1
1437	except KeyError:
1438	break
1439	iDb = 0
1440	while True:
1441	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1442	try:
1443	for i,name in enumerate(names):
1444	val = parmList[name+str(iDb)]
1445	Background[1]['peaksList'][iDb][2*i] = val
1446	iDb += 1
1447	except KeyError:
1448	break
1449
1450	def BackgroundPrint(Background,sigDict):
1451	print 'Background coefficients for',Background[0][0],'function'
1452	ptfmt = "%12.5f"
1453	ptstr = 'value: '
1454	sigstr = 'esd : '
1455	for i,back in enumerate(Background[0][3:]):
1456	ptstr += ptfmt % (back)
1457	if Background[0][1]:
1458	prm = 'Back;'+str(i)
1459	if prm in sigDict:
1460	sigstr += ptfmt % (sigDict[prm])
1461	else:
1462	sigstr += " "*12
1463	if len(ptstr) > 75:
1464	print ptstr
1465	if Background[0][1]: print sigstr
1466	ptstr = 'value: '
1467	sigstr = 'esd : '
1468	if len(ptstr) > 8:
1469	print ptstr
1470	if Background[0][1]: print sigstr
1471
1472	if Background[1]['nDebye']:
1473	parms = ['DebyeA;','DebyeR;','DebyeU;']
1474	print 'Debye diffuse scattering coefficients'
1475	ptfmt = "%12.5f"
1476	print ' term DebyeA esd DebyeR esd DebyeU esd'
1477	for term in range(Background[1]['nDebye']):
1478	line = ' term %d'%(term)
1479	for ip,name in enumerate(parms):
1480	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1481	if name+str(term) in sigDict:
1482	line += ptfmt%(sigDict[name+str(term)])
1483	else:
1484	line += " "*12
1485	print line
1486	if Background[1]['nPeaks']:
1487	print 'Coefficients for Background Peaks'
1488	ptfmt = "%15.3f"
1489	for j,pl in enumerate(Background[1]['peaksList']):
1490	names = 'peak %3d:'%(j+1)
1491	ptstr = 'values :'
1492	sigstr = 'esds :'
1493	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1494	val = pl[2*i]
1495	prm = lbl+";"+str(j)
1496	names += '%15s'%(prm)
1497	ptstr += ptfmt%(val)
1498	if prm in sigDict:
1499	sigstr += ptfmt%(sigDict[prm])
1500	else:
1501	sigstr += " "*15
1502	print names
1503	print ptstr
1504	print sigstr
1505
1506	def SetInstParms(Inst):
1507	dataType = Inst['Type'][0]
1508	insVary = []
1509	insNames = []
1510	insVals = []
1511	for parm in Inst:
1512	insNames.append(parm)
1513	insVals.append(Inst[parm][1])
1514	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1515	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1516	insVary.append(parm)
1517	instDict = dict(zip(insNames,insVals))
1518	instDict['X'] = max(instDict['X'],0.01)
1519	instDict['Y'] = max(instDict['Y'],0.01)
1520	if 'SH/L' in instDict:
1521	instDict['SH/L'] = max(instDict['SH/L'],0.002)
1522	return dataType,instDict,insVary
1523
1524	def GetInstParms(parmDict,Inst,varyList,Peaks):
1525	for name in Inst:
1526	Inst[name][1] = parmDict[name]
1527	iPeak = 0
1528	while True:
1529	try:
1530	sigName = 'sig'+str(iPeak)
1531	pos = parmDict['pos'+str(iPeak)]
1532	if sigName not in varyList:
1533	if 'C' in Inst['Type'][0]:
1534	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1535	else:
1536	dsp = G2lat.Pos2dsp(Inst,pos)
1537	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1538	gamName = 'gam'+str(iPeak)
1539	if gamName not in varyList:
1540	if 'C' in Inst['Type'][0]:
1541	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1542	else:
1543	dsp = G2lat.Pos2dsp(Inst,pos)
1544	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1545	iPeak += 1
1546	except KeyError:
1547	break
1548
1549	def InstPrint(Inst,sigDict):
1550	print 'Instrument Parameters:'
1551	ptfmt = "%12.6f"
1552	ptlbls = 'names :'
1553	ptstr = 'values:'
1554	sigstr = 'esds :'
1555	for parm in Inst:
1556	if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1557	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1558	ptlbls += "%s" % (parm.center(12))
1559	ptstr += ptfmt % (Inst[parm][1])
1560	if parm in sigDict:
1561	sigstr += ptfmt % (sigDict[parm])
1562	else:
1563	sigstr += 12*' '
1564	print ptlbls
1565	print ptstr
1566	print sigstr
1567
1568	def SetPeaksParms(dataType,Peaks):
1569	peakNames = []
1570	peakVary = []
1571	peakVals = []
1572	if 'C' in dataType:
1573	names = ['pos','int','sig','gam']
1574	else:
1575	names = ['pos','int','alp','bet','sig','gam']
1576	for i,peak in enumerate(Peaks):
1577	for j,name in enumerate(names):
1578	peakVals.append(peak[2*j])
1579	parName = name+str(i)
1580	peakNames.append(parName)
1581	if peak[2*j+1]:
1582	peakVary.append(parName)
1583	return dict(zip(peakNames,peakVals)),peakVary
1584
1585	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1586	if 'C' in Inst['Type'][0]:
1587	names = ['pos','int','sig','gam']
1588	else: #'T'
1589	names = ['pos','int','alp','bet','sig','gam']
1590	for i,peak in enumerate(Peaks):
1591	pos = parmDict['pos'+str(i)]
1592	if 'difC' in Inst:
1593	dsp = pos/Inst['difC'][1]
1594	for j in range(len(names)):
1595	parName = names[j]+str(i)
1596	if parName in varyList:
1597	peak[2*j] = parmDict[parName]
1598	elif 'alpha' in parName:
1599	peak[2*j] = parmDict['alpha']/dsp
1600	elif 'beta' in parName:
1601	peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1602	elif 'sig' in parName:
1603	if 'C' in Inst['Type'][0]:
1604	peak[2*j] = G2mth.getCWsig(parmDict,pos)
1605	else:
1606	peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1607	elif 'gam' in parName:
1608	if 'C' in Inst['Type'][0]:
1609	peak[2*j] = G2mth.getCWgam(parmDict,pos)
1610	else:
1611	peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1612
1613	def PeaksPrint(dataType,parmDict,sigDict,varyList):
1614	print 'Peak coefficients:'
1615	if 'C' in dataType:
1616	names = ['pos','int','sig','gam']
1617	else: #'T'
1618	names = ['pos','int','alp','bet','sig','gam']
1619	head = 13*' '
1620	for name in names:
1621	if name in ['alp','bet']:
1622	head += name.center(8)+'esd'.center(8)
1623	else:
1624	head += name.center(10)+'esd'.center(10)
1625	print head
1626	if 'C' in dataType:
1627	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1628	else:
1629	ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1630	for i,peak in enumerate(Peaks):
1631	ptstr = ':'
1632	for j in range(len(names)):
1633	name = names[j]
1634	parName = name+str(i)
1635	ptstr += ptfmt[name] % (parmDict[parName])
1636	if parName in varyList:
1637	ptstr += ptfmt[name] % (sigDict[parName])
1638	else:
1639	if name in ['alp','bet']:
1640	ptstr += 8*' '
1641	else:
1642	ptstr += 10*' '
1643	print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1644
1645	def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1646	parmdict.update(zip(varylist,values))
1647	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1648
1649	def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):
1650	parmdict.update(zip(varylist,values))
1651	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1652	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
1653	if dlg:
1654	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1655	if not GoOn:
1656	return -M #abort!!
1657	return M
1658
1659	if controls:
1660	Ftol = controls['min dM/M']
1661	derivType = controls['deriv type']
1662	else:
1663	Ftol = 0.0001
1664	derivType = 'analytic'
1665	if oneCycle:
1666	Ftol = 1.0
1667	x,y,w,yc,yb,yd = data #these are numpy arrays!
1668	yc *= 0. #set calcd ones to zero
1669	yb *= 0.
1670	yd *= 0.
1671	xBeg = np.searchsorted(x,Limits[0])
1672	xFin = np.searchsorted(x,Limits[1])
1673	bakType,bakDict,bakVary = SetBackgroundParms(Background)
1674	dataType,insDict,insVary = SetInstParms(Inst)
1675	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1676	parmDict = {}
1677	parmDict.update(bakDict)
1678	parmDict.update(insDict)
1679	parmDict.update(peakDict)
1680	parmDict['Pdabc'] = [] #dummy Pdabc
1681	parmDict.update(Inst2) #put in real one if there
1682	if prevVaryList:
1683	varyList = prevVaryList[:]
1684	else:
1685	varyList = bakVary+insVary+peakVary
1686	fullvaryList = varyList[:]
1687	while True:
1688	begin = time.time()
1689	values = np.array(Dict2Values(parmDict, varyList))
1690	Rvals = {}
1691	badVary = []
1692	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1693	args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1694	ncyc = int(result[2]['nfev']/2)
1695	runtime = time.time()-begin
1696	chisq = np.sum(result[2]['fvec']**2)
1697	Values2Dict(parmDict, varyList, result[0])
1698	Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]y[xBeg:xFin]2))100. #to %
1699	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
1700	print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1701	if ncyc:
1702	print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1703	print 'Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1704	sig = [0]*len(varyList)
1705	if len(varyList) == 0: break # if nothing was refined
1706	try:
1707	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1708	if np.any(np.isnan(sig)):
1709	print '* Least squares aborted - some invalid esds possible *'
1710	break #refinement succeeded - finish up!
1711	except ValueError: #result[1] is None on singular matrix
1712	print '** Refinement failed - singular matrix **'
1713	Ipvt = result[2]['ipvt']
1714	for i,ipvt in enumerate(Ipvt):
1715	if not np.sum(result[2]['fjac'],axis=1)[i]:
1716	print 'Removing parameter: ',varyList[ipvt-1]
1717	badVary.append(varyList[ipvt-1])
1718	del(varyList[ipvt-1])
1719	break
1720	else: # nothing removed
1721	break
1722	if dlg: dlg.Destroy()
1723	sigDict = dict(zip(varyList,sig))
1724	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1725	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1726	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1727	GetBackgroundParms(parmDict,Background)
1728	if bakVary: BackgroundPrint(Background,sigDict)
1729	GetInstParms(parmDict,Inst,varyList,Peaks)
1730	if insVary: InstPrint(Inst,sigDict)
1731	GetPeaksParms(Inst,parmDict,Peaks,varyList)
1732	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList)
1733	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1734
1735	def calcIncident(Iparm,xdata):
1736	'needs a doc string'
1737
1738	def IfunAdv(Iparm,xdata):
1739	Itype = Iparm['Itype']
1740	Icoef = Iparm['Icoeff']
1741	DYI = np.ones((12,xdata.shape[0]))
1742	YI = np.ones_like(xdata)*Icoef[0]
1743
1744	x = xdata/1000. #expressions are in ms
1745	if Itype == 'Exponential':
1746	for i in [1,3,5,7,9]:
1747	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1748	YI += Icoef[i]*Eterm
1749	DYI[i] *= Eterm
1750	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1751	elif 'Maxwell'in Itype:
1752	Eterm = np.exp(-Icoef[2]/x**2)
1753	DYI[1] = Eterm/x**5
1754	DYI[2] = -Icoef[1]DYI[1]/x*2
1755	YI += (Icoef[1]Eterm/x*5)
1756	if 'Exponential' in Itype:
1757	for i in range(3,12,2):
1758	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
1759	YI += Icoef[i]*Eterm
1760	DYI[i] *= Eterm
1761	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
1762	else: #Chebyschev
1763	T = (2./x)-1.
1764	Ccof = np.ones((12,xdata.shape[0]))
1765	Ccof[1] = T
1766	for i in range(2,12):
1767	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
1768	for i in range(1,10):
1769	YI += Ccof[i]*Icoef[i+2]
1770	DYI[i+2] =Ccof[i]
1771	return YI,DYI
1772
1773	Iesd = np.array(Iparm['Iesd'])
1774	Icovar = Iparm['Icovar']
1775	YI,DYI = IfunAdv(Iparm,xdata)
1776	YI = np.where(YI>0,YI,1.)
1777	WYI = np.zeros_like(xdata)
1778	vcov = np.zeros((12,12))
1779	k = 0
1780	for i in range(12):
1781	for j in range(i,12):
1782	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
1783	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
1784	k += 1
1785	M = np.inner(vcov,DYI.T)
1786	WYI = np.sum(M*DYI,axis=0)
1787	WYI = np.where(WYI>0.,WYI,0.)
1788	return YI,WYI
1789
1790	################################################################################
1791	# Stacking fault simulation codes
1792	################################################################################
1793
1794	def GetStackParms(Layers):
1795
1796	Parms = []
1797	#cell parms
1798	cell = Layers['Cell']
1799	if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
1800	Parms.append('cellA')
1801	Parms.append('cellC')
1802	else:
1803	Parms.append('cellA')
1804	Parms.append('cellB')
1805	Parms.append('cellC')
1806	if Layers['Laue'] != 'mmm':
1807	Parms.append('cellG')
1808	#Transition parms
1809	Trans = Layers['Transitions']
1810	for iY in range(len(Layers['Layers'])):
1811	for iX in range(len(Layers['Layers'])):
1812	Parms.append('TransP;%d;%d'%(iY,iX))
1813	Parms.append('TransX;%d;%d'%(iY,iX))
1814	Parms.append('TransY;%d;%d'%(iY,iX))
1815	Parms.append('TransZ;%d;%d'%(iY,iX))
1816	return Parms
1817
1818	def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
1819	'''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
1820
1821	param: Layers dict: 'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
1822	'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
1823	'Layers':[],'Stacking':[],'Transitions':[]}
1824	param: ctrls string: controls string to be written on DIFFaX controls.dif file
1825	param: scale float: scale factor
1826	param: background dict: background parameters
1827	param: limits list: min/max 2-theta to be calculated
1828	param: inst dict: instrument parameters dictionary
1829	param: profile list: powder pattern data
1830
1831	all updated in place
1832	'''
1833	import atmdata
1834	path = sys.path
1835	for name in path:
1836	if 'bin' in name:
1837	DIFFaX = name+'/DIFFaX.exe'
1838	print ' Execute ',DIFFaX
1839	break
1840	# make form factor file that DIFFaX wants - atom types are GSASII style
1841	sf = open('data.sfc','w')
1842	sf.write('GSASII special form factor file for DIFFaX\n\n')
1843	atTypes = Layers['AtInfo'].keys()
1844	if 'H' not in atTypes:
1845	atTypes.insert(0,'H')
1846	for atType in atTypes:
1847	if atType == 'H':
1848	blen = -.3741
1849	else:
1850	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
1851	Adat = atmdata.XrayFF[atType]
1852	text = '%4s'%(atType.ljust(4))
1853	for i in range(4):
1854	text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
1855	text += '%11.6f%11.6f'%(Adat['fc'],blen)
1856	text += '%3d\n'%(Adat['Z'])
1857	sf.write(text)
1858	sf.close()
1859	#make DIFFaX control.dif file - future use GUI to set some of these flags
1860	cf = open('control.dif','w')
1861	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
1862	x0 = profile[0]
1863	iBeg = np.searchsorted(x0,limits[0])
1864	iFin = np.searchsorted(x0,limits[1])
1865	if iFin-iBeg > 20000:
1866	iFin = iBeg+20000
1867	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
1868	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1869	cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
1870	else:
1871	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
1872	inst = {'Type':['XSC','XSC',]}
1873	cf.close()
1874	#make DIFFaX data file
1875	df = open('GSASII-DIFFaX.dat','w')
1876	df.write('INSTRUMENTAL\n')
1877	if 'X' in inst['Type'][0]:
1878	df.write('X-RAY\n')
1879	elif 'N' in inst['Type'][0]:
1880	df.write('NEUTRON\n')
1881	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
1882	df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
1883	U = ateln2*inst['U'][1]/10000.
1884	V = ateln2*inst['V'][1]/10000.
1885	W = ateln2*inst['W'][1]/10000.
1886	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
1887	HW = np.sqrt(np.mean(HWHM))
1888	# df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
1889	if 'Mean' in Layers['selInst']:
1890	df.write('GAUSSIAN %.6f TRIM\n'%(HW)) #fast option - might not really matter
1891	elif 'Gaussian' in Layers['selInst']:
1892	df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W)) #slow - make a GUI option?
1893	else:
1894	df.write('None\n')
1895	else:
1896	df.write('0.10\nNone\n')
1897	df.write('STRUCTURAL\n')
1898	a,b,c = Layers['Cell'][1:4]
1899	gam = Layers['Cell'][6]
1900	df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
1901	laue = Layers['Laue']
1902	if laue == '2/m(ab)':
1903	laue = '2/m(1)'
1904	elif laue == '2/m(c)':
1905	laue = '2/m(2)'
1906	if 'unknown' in Layers['Laue']:
1907	df.write('%s %.3f\n'%(laue,Layers['Toler']))
1908	else:
1909	df.write('%s\n'%(laue))
1910	df.write('%d\n'%(len(Layers['Layers'])))
1911	if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
1912	df.write('%.1f %.1f\n'%(Layers['Width'][0][0]10000.,Layers['Width'][0][0]10000.)) #mum to A
1913	layerNames = []
1914	for layer in Layers['Layers']:
1915	layerNames.append(layer['Name'])
1916	for il,layer in enumerate(Layers['Layers']):
1917	if layer['SameAs']:
1918	df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
1919	continue
1920	df.write('LAYER %d\n'%(il+1))
1921	if '-1' in layer['Symm']:
1922	df.write('CENTROSYMMETRIC\n')
1923	else:
1924	df.write('NONE\n')
1925	for ia,atom in enumerate(layer['Atoms']):
1926	[name,atype,x,y,z,frac,Uiso] = atom
1927	if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
1928	frac /= 2.
1929	df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
1930	df.write('STACKING\n')
1931	df.write('%s\n'%(Layers['Stacking'][0]))
1932	if 'recursive' in Layers['Stacking'][0]:
1933	df.write('%s\n'%Layers['Stacking'][1])
1934	else:
1935	if 'list' in Layers['Stacking'][1]:
1936	Slen = len(Layers['Stacking'][2])
1937	iB = 0
1938	iF = 0
1939	while True:
1940	iF += 68
1941	if iF >= Slen:
1942	break
1943	iF = min(iF,Slen)
1944	df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
1945	iB = iF
1946	else:
1947	df.write('%s\n'%Layers['Stacking'][1])
1948	df.write('TRANSITIONS\n')
1949	for iY in range(len(Layers['Layers'])):
1950	sumPx = 0.
1951	for iX in range(len(Layers['Layers'])):
1952	p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
1953	p = round(p,3)
1954	df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
1955	sumPx += p
1956	if sumPx != 1.0: #this has to be picky since DIFFaX is.
1957	print 'ERROR - Layer probabilities sum to ',sumPx,' DIFFaX will insist it = 1.0'
1958	df.close()
1959	os.remove('data.sfc')
1960	os.remove('control.dif')
1961	os.remove('GSASII-DIFFaX.dat')
1962	return
1963	df.close()
1964	time0 = time.time()
1965	try:
1966	subp.call(DIFFaX)
1967	except OSError:
1968	print ' DIFFax.exe is not available for this platform - under development'
1969	print ' DIFFaX time = %.2fs'%(time.time()-time0)
1970	if os.path.exists('GSASII-DIFFaX.spc'):
1971	Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
1972	iFin = iBeg+Xpat.shape[1]
1973	bakType,backDict,backVary = SetBackgroundParms(background)
1974	backDict['Lam1'] = G2mth.getWave(inst)
1975	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
1976	profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
1977	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
1978	rv = st.poisson(profile[3][iBeg:iFin])
1979	profile[1][iBeg:iFin] = rv.rvs()
1980	Z = np.ones_like(profile[3][iBeg:iFin])
1981	Z[1::2] *= -1
1982	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
1983	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
1984	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
1985	#cleanup files..
1986	os.remove('GSASII-DIFFaX.spc')
1987	elif os.path.exists('GSASII-DIFFaX.sadp'):
1988	Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
1989	Sadp = np.reshape(Sadp,(256,-1))
1990	Layers['Sadp']['Img'] = Sadp
1991	os.remove('GSASII-DIFFaX.sadp')
1992	os.remove('data.sfc')
1993	os.remove('control.dif')
1994	os.remove('GSASII-DIFFaX.dat')
1995
1996	def SetPWDRscan(inst,limits,profile):
1997
1998	wave = G2mth.getMeanWave(inst)
1999	x0 = profile[0]
2000	iBeg = np.searchsorted(x0,limits[0])
2001	iFin = np.searchsorted(x0,limits[1])
2002	if iFin-iBeg > 20000:
2003	iFin = iBeg+20000
2004	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
2005	pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
2006	return iFin-iBeg
2007
2008	def SetStackingSF(Layers,debug):
2009	# Load scattering factors into DIFFaX arrays
2010	import atmdata
2011	atTypes = Layers['AtInfo'].keys()
2012	aTypes = []
2013	for atype in atTypes:
2014	aTypes.append('%4s'%(atype.ljust(4)))
2015	SFdat = []
2016	for atType in atTypes:
2017	if atType == 'H':
2018	blen = -.3741
2019	else:
2020	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
2021	Adat = atmdata.XrayFF[atType]
2022	SF = np.zeros(9)
2023	SF[:8:2] = Adat['fa']
2024	SF[1:8:2] = Adat['fb']
2025	SF[8] = Adat['fc']
2026	SFdat.append(SF)
2027	SFdat = np.array(SFdat)
2028	pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
2029
2030	def SetStackingClay(Layers,Type):
2031	# Controls
2032	rand.seed()
2033	ranSeed = rand.randint(1,2**16-1)
2034	try:
2035	laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
2036	'6/m','6/mmm'].index(Layers['Laue'])+1
2037	except ValueError: #for 'unknown'
2038	laueId = -1
2039	if 'SADP' in Type:
2040	planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
2041	lmax = int(Layers['Sadp']['Lmax'])
2042	else:
2043	planeId = 0
2044	lmax = 0
2045	# Sequences
2046	StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
2047	try:
2048	StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
2049	except ValueError:
2050	StkParm = -1
2051	if StkParm == 2: #list
2052	StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
2053	Nstk = len(StkSeq)
2054	else:
2055	Nstk = 1
2056	StkSeq = [0,]
2057	if StkParm == -1:
2058	StkParm = int(Layers['Stacking'][1])
2059	Wdth = Layers['Width'][0]
2060	mult = 1
2061	controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
2062	LaueSym = Layers['Laue'].ljust(12)
2063	pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
2064	return laueId,controls
2065
2066	def SetCellAtoms(Layers):
2067	Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
2068	# atoms in layers
2069	atTypes = Layers['AtInfo'].keys()
2070	AtomXOU = []
2071	AtomTp = []
2072	LayerSymm = []
2073	LayerNum = []
2074	layerNames = []
2075	Natm = 0
2076	Nuniq = 0
2077	for layer in Layers['Layers']:
2078	layerNames.append(layer['Name'])
2079	for il,layer in enumerate(Layers['Layers']):
2080	if layer['SameAs']:
2081	LayerNum.append(layerNames.index(layer['SameAs'])+1)
2082	continue
2083	else:
2084	LayerNum.append(il+1)
2085	Nuniq += 1
2086	if '-1' in layer['Symm']:
2087	LayerSymm.append(1)
2088	else:
2089	LayerSymm.append(0)
2090	for ia,atom in enumerate(layer['Atoms']):
2091	[name,atype,x,y,z,frac,Uiso] = atom
2092	Natm += 1
2093	AtomTp.append('%4s'%(atype.ljust(4)))
2094	Ta = atTypes.index(atype)+1
2095	AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
2096	AtomXOU = np.array(AtomXOU)
2097	Nlayers = len(layerNames)
2098	pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
2099	return Nlayers
2100
2101	def SetStackingTrans(Layers,Nlayers):
2102	# Transitions
2103	TransX = []
2104	TransP = []
2105	for Ytrans in Layers['Transitions']:
2106	TransP.append([trans[0] for trans in Ytrans]) #get just the numbers
2107	TransX.append([trans[1:4] for trans in Ytrans]) #get just the numbers
2108	TransP = np.array(TransP,dtype='float').T
2109	TransX = np.array(TransX,dtype='float')
2110	# GSASIIpath.IPyBreak()
2111	pyx.pygettrans(Nlayers,TransP,TransX)
2112
2113	def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
2114	# Scattering factors
2115	SetStackingSF(Layers,debug)
2116	# Controls & sequences
2117	laueId,controls = SetStackingClay(Layers,'PWDR')
2118	# cell & atoms
2119	Nlayers = SetCellAtoms(Layers)
2120	Volume = Layers['Cell'][7]
2121	# Transitions
2122	SetStackingTrans(Layers,Nlayers)
2123	# PWDR scan
2124	Nsteps = SetPWDRscan(inst,limits,profile)
2125	# result as Spec
2126	x0 = profile[0]
2127	profile[3] = np.zeros(len(profile[0]))
2128	profile[4] = np.zeros(len(profile[0]))
2129	profile[5] = np.zeros(len(profile[0]))
2130	iBeg = np.searchsorted(x0,limits[0])
2131	iFin = np.searchsorted(x0,limits[1])
2132	if iFin-iBeg > 20000:
2133	iFin = iBeg+20000
2134	Nspec = 20001
2135	spec = np.zeros(Nspec,dtype='double')
2136	time0 = time.time()
2137	pyx.pygetspc(controls,Nspec,spec)
2138	print ' GETSPC time = %.2fs'%(time.time()-time0)
2139	time0 = time.time()
2140	U = ateln2*inst['U'][1]/10000.
2141	V = ateln2*inst['V'][1]/10000.
2142	W = ateln2*inst['W'][1]/10000.
2143	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
2144	HW = np.sqrt(np.mean(HWHM))
2145	BrdSpec = np.zeros(Nsteps)
2146	if 'Mean' in Layers['selInst']:
2147	pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
2148	elif 'Gaussian' in Layers['selInst']:
2149	pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
2150	else:
2151	BrdSpec = spec[:Nsteps]
2152	BrdSpec /= Volume
2153	iFin = iBeg+Nsteps
2154	bakType,backDict,backVary = SetBackgroundParms(background)
2155	backDict['Lam1'] = G2mth.getWave(inst)
2156	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
2157	profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
2158	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
2159	rv = st.poisson(profile[3][iBeg:iFin])
2160	profile[1][iBeg:iFin] = rv.rvs()
2161	Z = np.ones_like(profile[3][iBeg:iFin])
2162	Z[1::2] *= -1
2163	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
2164	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
2165	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
2166	print ' Broadening time = %.2fs'%(time.time()-time0)
2167
2168	def CalcStackingSADP(Layers,debug):
2169
2170	# Scattering factors
2171	SetStackingSF(Layers,debug)
2172	# Controls & sequences
2173	laueId,controls = SetStackingClay(Layers,'SADP')
2174	# cell & atoms
2175	Nlayers = SetCellAtoms(Layers)
2176	# Transitions
2177	SetStackingTrans(Layers,Nlayers)
2178	# result as Sadp
2179	mirror = laueId in [-1,2,3,7,8,9,10]
2180	Nspec = 20001
2181	spec = np.zeros(Nspec,dtype='double')
2182	time0 = time.time()
2183	hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
2184	Sapd = np.zeros((256,256))
2185	iB = 0
2186	for i in range(hkLim):
2187	iF = iB+Nblk
2188	p1 = 127+int(i*Incr)
2189	p2 = 128-int(i*Incr)
2190	if Nblk == 128:
2191	if i:
2192	Sapd[128:,p1] = spec[iB:iF]
2193	Sapd[:128,p1] = spec[iF:iB:-1]
2194	Sapd[128:,p2] = spec[iB:iF]
2195	Sapd[:128,p2] = spec[iF:iB:-1]
2196	else:
2197	if i:
2198	Sapd[:,p1] = spec[iB:iF]
2199	Sapd[:,p2] = spec[iB:iF]
2200	iB += Nblk
2201	Layers['Sadp']['Img'] = Sapd
2202	print ' GETSAD time = %.2fs'%(time.time()-time0)
2203	# GSASIIpath.IPyBreak()
2204
2205	#testing data
2206	NeedTestData = True
2207	def TestData():
2208	'needs a doc string'
2209	# global NeedTestData
2210	NeedTestData = False
2211	global bakType
2212	bakType = 'chebyschev'
2213	global xdata
2214	xdata = np.linspace(4.0,40.0,36000)
2215	global parmDict0
2216	parmDict0 = {
2217	'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
2218	'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
2219	'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
2220	'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
2221	'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
2222	'Back0':5.384,'Back1':-0.015,'Back2':.004,
2223	}
2224	global parmDict1
2225	parmDict1 = {
2226	'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
2227	'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
2228	'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
2229	'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
2230	'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
2231	'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
2232	'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
2233	'Back0':36.897,'Back1':-0.508,'Back2':.006,
2234	'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2235	}
2236	global parmDict2
2237	parmDict2 = {
2238	'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
2239	'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
2240	'Back0':5.,'Back1':-0.02,'Back2':.004,
2241	# 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
2242	}
2243	global varyList
2244	varyList = []
2245
2246	def test0():
2247	if NeedTestData: TestData()
2248	msg = 'test '
2249	gplot = plotter.add('FCJ-Voigt, 11BM').gca()
2250	gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])
2251	gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
2252	fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
2253	fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])
2254	fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
2255
2256	def test1():
2257	if NeedTestData: TestData()
2258	time0 = time.time()
2259	for i in range(100):
2260	y = getPeakProfile(parmDict1,xdata,varyList,bakType)
2261	print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
2262
2263	def test2(name,delt):
2264	if NeedTestData: TestData()
2265	varyList = [name,]
2266	xdata = np.linspace(5.6,5.8,400)
2267	hplot = plotter.add('derivatives test for '+name).gca()
2268	hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
2269	y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2270	parmDict2[name] += delt
2271	y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
2272	hplot.plot(xdata,(y1-y0)/delt,'r+')
2273
2274	def test3(name,delt):
2275	if NeedTestData: TestData()
2276	names = ['pos','sig','gam','shl']
2277	idx = names.index(name)
2278	myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
2279	xdata = np.linspace(5.6,5.8,800)
2280	dx = xdata[1]-xdata[0]
2281	hplot = plotter.add('derivatives test for '+name).gca()
2282	hplot.plot(xdata,100.dxgetdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
2283	y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2284	myDict[name] += delt
2285	y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
2286	hplot.plot(xdata,(y1-y0)/delt,'r+')
2287
2288	if __name__ == '__main__':
2289	import GSASIItestplot as plot
2290	global plotter
2291	plotter = plot.PlotNotebook()
2292	# test0()
2293	# for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
2294	for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
2295	['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
2296	test2(name,shft)
2297	for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
2298	test3(name,shft)
2299	print "OK"
2300	plotter.StartEventLoop()

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

Download in other formats: