Context Navigation

source: trunk/GSASIIpwd.py @ 3500

Last change on this file since 3500 was 3500, checked in by vondreele, 5 years ago
undo that mask fix - wrong thing to do.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 112.7 KB

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

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

Download in other formats: