Context Navigation

source: trunk/GSASIIpwd.py @ 4046

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

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

Download in other formats: