Context Navigation

source: trunk/GSASIIpwd.py @ 3848

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

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

Download in other formats: