Context Navigation

source: trunk/GSASIIpwd.py @ 4347

Last change on this file since 4347 was 4347, checked in by vondreele, 3 years ago
restrict G(R) plots to Rmax - sphere enclosed by big box model in RMCProfile
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 146.2 KB

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

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

Download in other formats: