Context Navigation

source: trunk/GSASIIpwd.py @ 3774

Last change on this file since 3774 was 3774, checked in by vondreele, 6 years ago
fix super indexing problem in transposeHKLF fix reflection generation for incommensurate mag case in G2lattice & G2pwd clean up non Fourier modulation calcs & remove analytic derivative stuff (now numeric) fix uij derivative bug work on incommensurate magnetic sturcture factors - not working yet clean up of testDeriv - better choices for delt & force reflection regeneration before each derivative test
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 113.2 KB

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

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

Download in other formats: