Context Navigation

source: trunk/GSASIIpwd.py @ 4179

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

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

Download in other formats: