Context Navigation

source: trunk/GSASIIpwd.py @ 3442

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

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

Download in other formats: