Context Navigation

source: trunk/GSASIIpwd.py @ 5058

Last change on this file since 5058 was 5058, checked in by vondreele, 2 years ago
provide sliders for mode displacements - not connected yet
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 186.5 KB

Line
1	#/usr/bin/env python
2	# -- coding: utf-8 --
3	'''
4	GSASIIpwd: Powder calculations module
5	==============================================
6
7	'''
8	########### SVN repository information ###################
9	# $Date: 2021-10-29 21:20:04 +0000 (Fri, 29 Oct 2021) $
10	# $Author: vondreele $
11	# $Revision: 5058 $
12	# $URL: trunk/GSASIIpwd.py $
13	# $Id: GSASIIpwd.py 5058 2021-10-29 21:20:04Z 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 os.path
21	import subprocess as subp
22	import datetime as dt
23	import copy
24
25	import numpy as np
26	import numpy.linalg as nl
27	import numpy.ma as ma
28	import random as rand
29	import numpy.fft as fft
30	import scipy.interpolate as si
31	import scipy.stats as st
32	import scipy.optimize as so
33	import scipy.special as sp
34	import scipy.signal as signal
35
36	import GSASIIpath
37	filversion = "$Revision: 5058 $"
38	GSASIIpath.SetVersionNumber("$Revision: 5058 $")
39	import GSASIIlattice as G2lat
40	import GSASIIspc as G2spc
41	import GSASIIElem as G2elem
42	import GSASIImath as G2mth
43	try:
44	import pypowder as pyd
45	except ImportError:
46	print ('pypowder is not available - profile calcs. not allowed')
47	try:
48	import pydiffax as pyx
49	except ImportError:
50	print ('pydiffax is not available for this platform')
51	import GSASIIfiles as G2fil
52
53
54	# trig functions in degrees
55	tand = lambda x: math.tan(x*math.pi/180.)
56	atand = lambda x: 180.*math.atan(x)/math.pi
57	atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
58	cosd = lambda x: math.cos(x*math.pi/180.)
59	acosd = lambda x: 180.*math.acos(x)/math.pi
60	rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
61	#numpy versions
62	npsind = lambda x: np.sin(x*np.pi/180.)
63	npasind = lambda x: 180.*np.arcsin(x)/math.pi
64	npcosd = lambda x: np.cos(x*math.pi/180.)
65	npacosd = lambda x: 180.*np.arccos(x)/math.pi
66	nptand = lambda x: np.tan(x*math.pi/180.)
67	npatand = lambda x: 180.*np.arctan(x)/np.pi
68	npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
69	npT2stl = lambda tth, wave: 2.0npsind(tth/2.0)/wave #=d
70	npT2q = lambda tth,wave: 2.0np.pinpT2stl(tth,wave) #=2pid
71	npq2T = lambda Q,wave: 2.0npasind(0.25Q*wave/np.pi)
72	ateln2 = 8.0*math.log(2.0)
73	sateln2 = np.sqrt(ateln2)
74	nxs = np.newaxis
75
76	#### Powder utilities ################################################################################
77	def PhaseWtSum(G2frame,histo):
78	'''
79	Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
80
81	:param G2frame: GSASII main frame structure
82	:param str histo: histogram name
83	:returns: sum(scale*mass) for phases in histo
84	'''
85	Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
86	wtSum = 0.0
87	for phase in Phases:
88	if Phases[phase]['General']['Type'] != 'magnetic':
89	if histo in Phases[phase]['Histograms']:
90	if not Phases[phase]['Histograms'][histo]['Use']: continue
91	mass = Phases[phase]['General']['Mass']
92	phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
93	wtSum += mass*phFr
94	return wtSum
95
96	#### GSASII pwdr & pdf calculation routines ################################################################################
97	def Transmission(Geometry,Abs,Diam):
98	'''
99	Calculate sample transmission
100
101	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
102	:param float Abs: absorption coeff in cm-1
103	:param float Diam: sample thickness/diameter in mm
104	'''
105	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
106	MuR = Abs*Diam/20.0
107	if MuR <= 3.0:
108	T0 = 16/(3.*math.pi)
109	T1 = -0.045780
110	T2 = -0.02489
111	T3 = 0.003045
112	T = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
113	if T < -20.:
114	return 2.06e-9
115	else:
116	return math.exp(T)
117	else:
118	T1 = 1.433902
119	T2 = 0.013869+0.337894
120	T3 = 1.933433+1.163198
121	T4 = 0.044365-0.04259
122	T = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
123	return T/100.
124	elif 'plate' in Geometry:
125	MuR = Abs*Diam/10.
126	return math.exp(-MuR)
127	elif 'Bragg' in Geometry:
128	return 0.0
129
130	def SurfaceRough(SRA,SRB,Tth):
131	''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
132	:param float SRA: Suortti surface roughness parameter
133	:param float SRB: Suortti surface roughness parameter
134	:param float Tth: 2-theta(deg) - can be numpy array
135
136	'''
137	sth = npsind(Tth/2.)
138	T1 = np.exp(-SRB/sth)
139	T2 = SRA+(1.-SRA)*np.exp(-SRB)
140	return (SRA+(1.-SRA)*T1)/T2
141
142	def SurfaceRoughDerv(SRA,SRB,Tth):
143	''' Suortti surface roughness correction derivatives
144	:param float SRA: Suortti surface roughness parameter (dimensionless)
145	:param float SRB: Suortti surface roughness parameter (dimensionless)
146	:param float Tth: 2-theta(deg) - can be numpy array
147	:return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
148	'''
149	sth = npsind(Tth/2.)
150	T1 = np.exp(-SRB/sth)
151	T2 = SRA+(1.-SRA)*np.exp(-SRB)
152	Trans = (SRA+(1.-SRA)*T1)/T2
153	dydSRA = ((1.-T1)T2-(1.-np.exp(-SRB))Trans)/T2**2
154	dydSRB = ((SRA-1.)T1T2/sth-Trans(SRA-T2))/T2*2
155	return [dydSRA,dydSRB]
156
157	def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
158	'''Calculate sample absorption
159	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
160	:param float MuR: absorption coeff * sample thickness/2 or radius
161	:param Tth: 2-theta scattering angle - can be numpy array
162	:param float Phi: flat plate tilt angle - future
163	:param float Psi: flat plate tilt axis - future
164	'''
165
166	def muRunder3(MuR,Sth2):
167	T0 = 16.0/(3.*np.pi)
168	T1 = (25.99978-0.01911Sth20.25)np.exp(-0.024551*Sth2)+ \
169	0.109561*np.sqrt(Sth2)-26.04556
170	T2 = -0.02489-0.39499Sth2+1.219077Sth2**1.5- \
171	1.31268Sth22+0.871081Sth2*2.5-0.2327Sth2**3
172	T3 = 0.003045+0.018167Sth2-0.03305Sth2**2
173	Trns = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
174	return np.exp(Trns)
175
176	def muRover3(MuR,Sth2):
177	T1 = 1.433902+11.07504Sth2-8.77629Sth2*Sth2+ \
178	10.02088Sth23-3.36778Sth2**4
179	T2 = (0.013869-0.01249Sth2)np.exp(3.27094*Sth2)+ \
180	(0.337894+13.77317Sth2)/(1.0+11.53544Sth2)**1.555039
181	T3 = 1.933433/(1.0+23.12967Sth2)*1.686715- \
182	0.13576*np.sqrt(Sth2)+1.163198
183	T4 = 0.044365-0.04259/(1.0+0.41051Sth2)*148.4202
184	Trns = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
185	return Trns/100.
186
187	Sth2 = npsind(Tth/2.0)**2
188	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
189	if 'array' in str(type(MuR)):
190	MuRSTh2 = np.vstack((MuR,Sth2))
191	AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
192	return AbsCr
193	else:
194	if MuR <= 3.0:
195	return muRunder3(MuR,Sth2)
196	else:
197	return muRover3(MuR,Sth2)
198	elif 'Bragg' in Geometry:
199	return 1.0
200	elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
201	# and only defined for 2theta < 90
202	MuT = 2.*MuR
203	T1 = np.exp(-MuT)
204	T2 = np.exp(-MuT/npcosd(Tth))
205	Tb = MuT-MuT/npcosd(Tth)
206	return (T2-T1)/Tb
207	elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
208	MuT = 2.*MuR
209	cth = npcosd(Tth/2.0)
210	return np.exp(-MuT/cth)/cth
211
212	def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
213	'needs a doc string'
214	dA = 0.001
215	AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
216	if MuR:
217	AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
218	return (AbsP-AbsM)/(2.0*dA)
219	else:
220	return (AbsP-1.)/dA
221
222	def Polarization(Pola,Tth,Azm=0.0):
223	""" Calculate angle dependent x-ray polarization correction (not scaled correctly!)
224
225	:param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
226	:param Azm: azimuthal angle e.g. 0.0 in plane of polarization - can be numpy array
227	:param Tth: 2-theta scattering angle - can be numpy array
228	which (if either) of these is "right"?
229	:return: (pola, dpdPola) - both 2-d arrays
230	* pola = ((1-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
231	(1-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
232	* dpdPola: derivative needed for least squares
233
234	"""
235	cazm = npcosd(Azm)**2
236	sazm = npsind(Azm)**2
237	pola = ((1.0-Pola)cazm+Polasazm)npcosd(Tth)2+(1.0-Pola)sazm+Pola*cazm
238	dpdPola = -npsind(Tth)*2(sazm-cazm)
239	return pola,dpdPola
240
241	def Oblique(ObCoeff,Tth):
242	'currently assumes detector is normal to beam'
243	if ObCoeff:
244	K = (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
245	return K
246	else:
247	return 1.0
248
249	def Ruland(RulCoff,wave,Q,Compton):
250	'needs a doc string'
251	C = 2.9978e8
252	D = 1.5e-3
253	hmc = 0.024262734687 #Compton wavelength in A
254	sinth2 = (Qwave/(4.0np.pi))**2
255	dlam = (wave*2)Compton*Q/C
256	dlam_c = 2.0hmcsinth2-Dwave*2
257	return 1.0/((1.0+dlam/RulCoff)(1.0+(np.pidlam_c/(dlam+RulCoff))**2))
258
259	def KleinNishina(wave,Q):
260	hmc = 0.024262734687 #Compton wavelength in A
261	TTh = npq2T(Q,wave)
262	P = 1./(1.+(1.-npcosd(TTh)*(hmc/wave)))
263	KN = (P*3-(Pnpsind(TTh))2+P)/(1.+npcosd(TTh)2)
264	return KN
265
266	def LorchWeight(Q):
267	'needs a doc string'
268	return np.sin(np.pi(Q[-1]-Q)/(2.0Q[-1]))
269
270	def GetAsfMean(ElList,Sthl2):
271	'''Calculate various scattering factor terms for PDF calcs
272
273	:param dict ElList: element dictionary contains scattering factor coefficients, etc.
274	:param np.array Sthl2: numpy array of sin theta/lambda squared values
275	:returns: mean(f^2), mean(f)^2, mean(compton)
276	'''
277	sumNoAtoms = 0.0
278	FF = np.zeros_like(Sthl2)
279	FF2 = np.zeros_like(Sthl2)
280	CF = np.zeros_like(Sthl2)
281	for El in ElList:
282	sumNoAtoms += ElList[El]['FormulaNo']
283	for El in ElList:
284	el = ElList[El]
285	ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])2+el['fpp']2
286	cf = G2elem.ComptonFac(el,Sthl2)
287	FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
288	FF2 += ff2*el['FormulaNo']/sumNoAtoms
289	CF += cf*el['FormulaNo']/sumNoAtoms
290	return FF2,FF**2,CF
291
292	def GetNumDensity(ElList,Vol):
293	'needs a doc string'
294	sumNoAtoms = 0.0
295	for El in ElList:
296	sumNoAtoms += ElList[El]['FormulaNo']
297	return sumNoAtoms/Vol
298
299	def CalcPDF(data,inst,limits,xydata):
300	'''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
301	dict xydata; results are placed in xydata.
302	Calculation parameters are found in dicts data and inst and list limits.
303	The return value is at present an empty list.
304	'''
305	auxPlot = []
306	if 'T' in inst['Type'][0]:
307	Ibeg = 0
308	Ifin = len(xydata['Sample'][1][0])
309	else:
310	Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
311	Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
312	#subtract backgrounds - if any & use PWDR limits
313	IofQ = copy.deepcopy(xydata['Sample'])
314	IofQ[1] = np.array([I[Ibeg:Ifin] for I in IofQ[1]])
315	if data['Sample Bkg.']['Name']:
316	IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
317	if data['Container']['Name']:
318	xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
319	if data['Container Bkg.']['Name']:
320	xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
321	IofQ[1][1] += xycontainer[Ibeg:Ifin]
322	data['IofQmin'] = IofQ[1][1][-1]
323	IofQ[1][1] -= data.get('Flat Bkg',0.)
324	#get element data & absorption coeff.
325	ElList = data['ElList']
326	Tth = IofQ[1][0] #2-theta or TOF!
327	if 'X' in inst['Type'][0]:
328	Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
329	#Apply angle dependent corrections
330	MuR = Abs*data['Diam']/20.0
331	IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
332	IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
333	if data['DetType'] == 'Area detector':
334	IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
335	elif 'T' in inst['Type'][0]: #neutron TOF normalized data - needs wavelength dependent absorption
336	wave = 2.G2lat.TOF2dsp(inst,IofQ[1][0])npsind(inst['2-theta'][1]/2.)
337	Els = ElList.keys()
338	Isotope = {El:'Nat. abund.' for El in Els}
339	GD = {'AtomTypes':ElList,'Isotope':Isotope}
340	BLtables = G2elem.GetBLtable(GD)
341	FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave)
342	Abs = np.zeros(len(wave))
343	for iel,El in enumerate(Els):
344	BL = BLtables[El][1]
345	SA = BL['SA']wave/1.798197+4.0np.piFPP[iel]*2 #+BL['SL'][1]?
346	SA *= ElList[El]['FormulaNo']/data['Form Vol']
347	Abs += SA
348	MuR = Abs*data['Diam']/2.
349	IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave)))
350	# improves look of F(Q) but no impact on G(R)
351	# bBut,aBut = signal.butter(8,.5,"lowpass")
352	# IofQ[1][1] = signal.filtfilt(bBut,aBut,IofQ[1][1])
353	XY = IofQ[1]
354	#convert to Q
355	nQpoints = 5000
356	if 'C' in inst['Type'][0]:
357	wave = G2mth.getWave(inst)
358	minQ = npT2q(Tth[0],wave)
359	maxQ = npT2q(Tth[-1],wave)
360	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
361	dq = Qpoints[1]-Qpoints[0]
362	XY[0] = npT2q(XY[0],wave)
363	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0]) #interpolate I(Q)
364	elif 'T' in inst['Type'][0]:
365	difC = inst['difC'][1]
366	minQ = 2.np.pidifC/Tth[-1]
367	maxQ = 2.np.pidifC/Tth[0]
368	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
369	dq = Qpoints[1]-Qpoints[0]
370	XY[0] = 2.np.pidifC/XY[0]
371	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1]) #interpolate I(Q)
372	Qdata -= np.min(Qdata)*data['BackRatio']
373
374	qLimits = data['QScaleLim']
375	maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1
376	minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1]))
377	qLimits = [Qpoints[minQ],Qpoints[maxQ-1]]
378	newdata = []
379	if len(IofQ) < 3:
380	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
381	else:
382	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
383	for item in xydata['IofQ'][1]:
384	newdata.append(item[:maxQ])
385	xydata['IofQ'][1] = newdata
386
387	xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
388	if 'XC' in inst['Type'][0]:
389	FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0np.pi))*2) #these are <f^2>,<f>^2,Cf
390	else: #TOF
391	CF = np.zeros(len(xydata['SofQ'][1][0]))
392	FFSq = np.ones(len(xydata['SofQ'][1][0]))
393	SqFF = np.ones(len(xydata['SofQ'][1][0]))
394	Q = xydata['SofQ'][1][0]
395	# auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
396	if 'XC' in inst['Type'][0]:
397	# CF *= KleinNishina(wave,Q)
398	ruland = Ruland(data['Ruland'],wave,Q,CF)
399	# auxPlot.append([Q,ruland,'Ruland'])
400	CF *= ruland
401	# auxPlot.append([Q,CF,'CF-Corr'])
402	scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
403	xydata['SofQ'][1][1] *= scale
404	if 'XC' in inst['Type'][0]:
405	xydata['SofQ'][1][1] -= CF
406	xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
407	scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
408	xydata['SofQ'][1][1] *= scale
409	xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
410	xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
411	if data['Lorch']:
412	xydata['FofQ'][1][1] *= LorchWeight(Q)
413	xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
414	xydata['gofr'] = copy.deepcopy(xydata['FofQ'])
415	nR = len(xydata['GofR'][1][1])
416	Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.)
417	mul = int(round(2.np.pinR/(Rmax*qLimits[1])))
418	# mul = int(round(2.np.pinR/(data.get('Rmax',100.)*qLimits[1])))
419	R = 2.np.pinp.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
420	xydata['GofR'][1][0] = R
421	xydata['gofr'][1][0] = R
422	GR = -dqnp.imag(fft.fft(xydata['FofQ'][1][1],mulnR)[:nR])*data.get('GR Scale',1.0)
423	xydata['GofR'][1][1] = GR
424	gr = GR/(np.pi*R)
425	xydata['gofr'][1][1] = gr
426	numbDen = 0.
427	if 'ElList' in data:
428	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
429	if data.get('noRing',True):
430	Rmin = data['Rmin']
431	xydata['gofr'][1][1] = np.where(R<Rmin,-4.*numbDen,xydata['gofr'][1][1])
432	xydata['GofR'][1][1] = np.where(R<Rmin,-4.Rnp.pi*numbDen,xydata['GofR'][1][1])
433	return auxPlot
434
435	def PDFPeakFit(peaks,data):
436	rs2pi = 1./np.sqrt(2*np.pi)
437
438	def MakeParms(peaks):
439	varyList = []
440	parmDict = {'slope':peaks['Background'][1][1]}
441	if peaks['Background'][2]:
442	varyList.append('slope')
443	for i,peak in enumerate(peaks['Peaks']):
444	parmDict['PDFpos;'+str(i)] = peak[0]
445	parmDict['PDFmag;'+str(i)] = peak[1]
446	parmDict['PDFsig;'+str(i)] = peak[2]
447	if 'P' in peak[3]:
448	varyList.append('PDFpos;'+str(i))
449	if 'M' in peak[3]:
450	varyList.append('PDFmag;'+str(i))
451	if 'S' in peak[3]:
452	varyList.append('PDFsig;'+str(i))
453	return parmDict,varyList
454
455	def SetParms(peaks,parmDict,varyList):
456	if 'slope' in varyList:
457	peaks['Background'][1][1] = parmDict['slope']
458	for i,peak in enumerate(peaks['Peaks']):
459	if 'PDFpos;'+str(i) in varyList:
460	peak[0] = parmDict['PDFpos;'+str(i)]
461	if 'PDFmag;'+str(i) in varyList:
462	peak[1] = parmDict['PDFmag;'+str(i)]
463	if 'PDFsig;'+str(i) in varyList:
464	peak[2] = parmDict['PDFsig;'+str(i)]
465
466
467	def CalcPDFpeaks(parmdict,Xdata):
468	Z = parmDict['slope']*Xdata
469	ipeak = 0
470	while True:
471	try:
472	pos = parmdict['PDFpos;'+str(ipeak)]
473	mag = parmdict['PDFmag;'+str(ipeak)]
474	wid = parmdict['PDFsig;'+str(ipeak)]
475	wid2 = 2.wid*2
476	Z += magrs2pinp.exp(-(Xdata-pos)**2/wid2)/wid
477	ipeak += 1
478	except KeyError: #no more peaks to process
479	return Z
480
481	def errPDFProfile(values,xdata,ydata,parmdict,varylist):
482	parmdict.update(zip(varylist,values))
483	M = CalcPDFpeaks(parmdict,xdata)-ydata
484	return M
485
486	newpeaks = copy.copy(peaks)
487	iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
488	iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
489	X = data[1][0][iBeg:iFin]
490	Y = data[1][1][iBeg:iFin]
491	parmDict,varyList = MakeParms(peaks)
492	if not len(varyList):
493	G2fil.G2Print (' Nothing varied')
494	return newpeaks,None,None,None,None,None
495
496	Rvals = {}
497	values = np.array(Dict2Values(parmDict, varyList))
498	result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
499	args=(X,Y,parmDict,varyList))
500	chisq = np.sum(result[2]['fvec']**2)
501	Values2Dict(parmDict, varyList, result[0])
502	SetParms(peaks,parmDict,varyList)
503	Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y*2))100. #to %
504	chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values)) #reduced chi^2 = M/(Nobs-Nvar)
505	sigList = list(np.sqrt(chisq*np.diag(result[1])))
506	Z = CalcPDFpeaks(parmDict,X)
507	newpeaks['calc'] = [X,Z]
508	return newpeaks,result[0],varyList,sigList,parmDict,Rvals
509
510	def MakeRDF(RDFcontrols,background,inst,pwddata):
511	auxPlot = []
512	if 'C' in inst['Type'][0] or 'B' in inst['Type'][0]:
513	Tth = pwddata[0]
514	wave = G2mth.getWave(inst)
515	minQ = npT2q(Tth[0],wave)
516	maxQ = npT2q(Tth[-1],wave)
517	powQ = npT2q(Tth,wave)
518	elif 'T' in inst['Type'][0]:
519	TOF = pwddata[0]
520	difC = inst['difC'][1]
521	minQ = 2.np.pidifC/TOF[-1]
522	maxQ = 2.np.pidifC/TOF[0]
523	powQ = 2.np.pidifC/TOF
524	piDQ = np.pi/(maxQ-minQ)
525	Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
526	if RDFcontrols['UseObsCalc'] == 'obs-calc':
527	Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
528	elif RDFcontrols['UseObsCalc'] == 'obs-back':
529	Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
530	elif RDFcontrols['UseObsCalc'] == 'calc-back':
531	Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
532	Qdata = np.sin((Qpoints-minQ)piDQ)/piDQ
533	Qdata = 0.5np.sqrt(Qpoints) #Qbin normalization
534	dq = Qpoints[1]-Qpoints[0]
535	nR = len(Qdata)
536	R = 0.5np.pinp.linspace(0,nR,nR)/(4.*maxQ)
537	iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
538	bBut,aBut = signal.butter(4,0.01)
539	Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
540	# auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
541	# auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
542	DofR = dqnp.imag(fft.fft(Qsmooth,16nR)[:nR])
543	auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']])
544	return auxPlot
545
546	# PDF optimization =============================================================
547	def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=25):
548	import scipy.optimize as opt
549	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
550	Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen)
551	xstart = Init()
552	bakMul = data['Sample Bkg.']['Mult']
553	if showFit:
554	rms = Min(xstart)
555	G2fil.G2Print(' Optimizing corrections to improve G(r) at low r')
556	if data['Sample Bkg.'].get('Refine',False):
557	# data['Flat Bkg'] = 0.
558	G2fil.G2Print(' start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format(
559	data['Ruland'],data['Sample Bkg.']['Mult'],rms))
560	else:
561	G2fil.G2Print(' start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format(
562	data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
563	if data['Sample Bkg.'].get('Refine',False):
564	res = opt.minimize(Min,xstart,bounds=([0.01,1.],[1.2bakMul,0.8bakMul]),
565	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
566	else:
567	res = opt.minimize(Min,xstart,bounds=([0.,None],[0,1],[0.01,1.]),
568	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
569	Done(res['x'])
570	if showFit:
571	if res['success']:
572	msg = 'Converged'
573	else:
574	msg = 'Not Converged'
575	if data['Sample Bkg.'].get('Refine',False):
576	G2fil.G2Print(' end: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) * {} *\n'.format(
577	data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
578	else:
579	G2fil.G2Print(' end: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) * {} *\n'.format(
580	data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
581	return res
582
583	def SetupPDFEval(data,xydata,limits,inst,numbDen):
584	Data = copy.deepcopy(data)
585	BkgMax = 1.
586	def EvalLowPDF(arg):
587	'''Objective routine -- evaluates the RMS deviations in G(r)
588	from -4(pi)#densityr for for r<Rmin
589	arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
590	the min & max values are between 0 and 1.
591	'''
592	if Data['Sample Bkg.'].get('Refine',False):
593	R,S = arg
594	Data['Sample Bkg.']['Mult'] = S
595	else:
596	F,B,R = arg
597	Data['Flat Bkg'] = BkgMax(2.F-1.)
598	Data['BackRatio'] = B
599	Data['Ruland'] = R
600	CalcPDF(Data,inst,limits,xydata)
601	# test low r computation
602	g = xydata['GofR'][1][1]
603	r = xydata['GofR'][1][0]
604	g0 = g[r < Data['Rmin']] + 4np.pir[r < Data['Rmin']]*numbDen
605	M = sum(g0**2)/len(g0)
606	return M
607	def GetCurrentVals():
608	'''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
609	'''
610	if data['Sample Bkg.'].get('Refine',False):
611	return [max(data['Ruland'],.05),data['Sample']['Mult']]
612	try:
613	F = 0.5+0.5*data['Flat Bkg']/BkgMax
614	except:
615	F = 0
616	return [F,data['BackRatio'],max(data['Ruland'],.05)]
617	def SetFinalVals(arg):
618	'''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
619	scaled, refined values and plot corrected region of G(r)
620	'''
621	if data['Sample Bkg.'].get('Refine',False):
622	R,S = arg
623	data['Sample Bkg.']['Mult'] = S
624	else:
625	F,B,R = arg
626	data['Flat Bkg'] = BkgMax(2.F-1.)
627	data['BackRatio'] = B
628	data['Ruland'] = R
629	CalcPDF(data,inst,limits,xydata)
630	EvalLowPDF(GetCurrentVals())
631	BkgMax = max(xydata['IofQ'][1][1])/50.
632	return EvalLowPDF,GetCurrentVals,SetFinalVals
633
634	#### GSASII peak fitting routines: Finger, Cox & Jephcoat model ################################################################################
635	def factorize(num):
636	''' Provide prime number factors for integer num
637	:returns: dictionary of prime factors (keys) & power for each (data)
638	'''
639	factors = {}
640	orig = num
641
642	# we take advantage of the fact that (i +1)2 = i2 + 2*i +1
643	i, sqi = 2, 4
644	while sqi <= num:
645	while not num%i:
646	num /= i
647	factors[i] = factors.get(i, 0) + 1
648
649	sqi += 2*i + 1
650	i += 1
651
652	if num != 1 and num != orig:
653	factors[num] = factors.get(num, 0) + 1
654
655	if factors:
656	return factors
657	else:
658	return {num:1} #a prime number!
659
660	def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
661	''' Provide list of optimal data sizes for FFT calculations
662
663	:param int nmin: minimum data size >= 1
664	:param int nmax: maximum data size > nmin
665	:param int thresh: maximum prime factor allowed
666	:Returns: list of data sizes where the maximum prime factor is < thresh
667	'''
668	plist = []
669	nmin = max(1,nmin)
670	nmax = max(nmin+1,nmax)
671	for p in range(nmin,nmax):
672	if max(list(factorize(p).keys())) < thresh:
673	plist.append(p)
674	return plist
675
676	np.seterr(divide='ignore')
677
678	# Normal distribution
679
680	# loc = mu, scale = std
681	_norm_pdf_C = 1./math.sqrt(2*math.pi)
682	class norm_gen(st.rv_continuous):
683	'needs a doc string'
684
685	def pdf(self,x,args,*kwds):
686	loc,scale=kwds['loc'],kwds['scale']
687	x = (x-loc)/scale
688	return np.exp(-x*2/2.0) _norm_pdf_C / scale
689
690	norm = norm_gen(name='norm',longname='A normal',extradoc="""
691
692	Normal distribution
693
694	The location (loc) keyword specifies the mean.
695	The scale (scale) keyword specifies the standard deviation.
696
697	normal.pdf(x) = exp(-x*2/2)/sqrt(2pi)
698	""")
699
700	## Cauchy
701
702	# median = loc
703
704	class cauchy_gen(st.rv_continuous):
705	'needs a doc string'
706
707	def pdf(self,x,args,*kwds):
708	loc,scale=kwds['loc'],kwds['scale']
709	x = (x-loc)/scale
710	return 1.0/np.pi/(1.0+x*x) / scale
711
712	cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
713
714	Cauchy distribution
715
716	cauchy.pdf(x) = 1/(pi(1+x*2))
717
718	This is the t distribution with one degree of freedom.
719	""")
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	"""
744	def _pdf(self,x,t,s,dx):
745	T = dx*x+t
746	ax2 = abs(npcosd(T))
747	ax = ax2**2
748	bx = npcosd(t)**2
749	bx = np.where(ax>bx,bx,ax)
750	fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
751	fx = np.where(fx > 0.,fx,0.0)
752	return fx
753
754	def pdf(self,x,args,*kwds):
755	loc=kwds['loc']
756	return self._pdf(x-loc,*args)
757
758	fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
759
760	def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
761	'''Compute the Finger-Cox-Jepcoat modified Voigt function for a
762	CW powder peak by direct convolution. This version is not used.
763	'''
764	DX = xdata[1]-xdata[0]
765	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
766	x = np.linspace(pos-fmin,pos+fmin,256)
767	dx = x[1]-x[0]
768	Norm = norm.pdf(x,loc=pos,scale=widths[0])
769	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
770	arg = [pos,shl/57.2958,dx,]
771	FCJ = fcjde.pdf(x,*arg,loc=pos)
772	if len(np.nonzero(FCJ)[0])>5:
773	z = np.column_stack([Norm,Cauchy,FCJ]).T
774	Z = fft.fft(z)
775	Df = fft.ifft(Z.prod(axis=0)).real
776	else:
777	z = np.column_stack([Norm,Cauchy]).T
778	Z = fft.fft(z)
779	Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
780	Df /= np.sum(Df)
781	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
782	return intensDf(xdata)DX/dx
783
784	#### GSASII peak fitting routine: Finger, Cox & Jephcoat model
785
786	def getWidthsCW(pos,sig,gam,shl):
787	'''Compute the peak widths used for computing the range of a peak
788	for constant wavelength data. On low-angle side, 50 FWHM are used,
789	on high-angle side 75 are used, high angle side extended for axial divergence
790	(for peaks above 90 deg, these are reversed.)
791
792	:param pos: peak position; 2-theta in degrees
793	:param sig: Gaussian peak variance in centideg^2
794	:param gam: Lorentzian peak width in centidegrees
795	:param shl: axial divergence parameter (S+H)/L
796
797	:returns: widths; [Gaussian sigma, Lorentzian gamma] in degrees, and
798	low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position
799	reversed for 2-theta > 90 deg.
800	'''
801	widths = [np.sqrt(sig)/100.,gam/100.]
802	fwhm = 2.355*widths[0]+widths[1]
803	fmin = 50.(fwhm+shlabs(npcosd(pos)))
804	fmax = 75.0*fwhm
805	if pos > 90:
806	fmin,fmax = [fmax,fmin]
807	return widths,fmin,fmax
808
809	def getWidthsTOF(pos,alp,bet,sig,gam):
810	'''Compute the peak widths used for computing the range of a peak
811	for constant wavelength data. 50 FWHM are used on both sides each
812	extended by exponential coeff.
813
814	param pos: peak position; TOF in musec
815	param alp,bet: TOF peak exponential rise & decay parameters
816	param sig: Gaussian peak variance in musec^2
817	param gam: Lorentzian peak width in musec
818
819	returns: widths; [Gaussian sigma, Lornetzian gamma] in musec
820	returns: low TOF, high TOF ends of peak; 50FWHM from position
821	'''
822	widths = [np.sqrt(sig),gam]
823	fwhm = 2.355widths[0]+2.widths[1]
824	fmin = 50.fwhm(1.+1./alp)
825	fmax = 50.fwhm(1.+1./bet)
826	return widths,fmin,fmax
827
828	def getFWHM(pos,Inst):
829	'''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83
830	via getgamFW(g,s).
831
832	:param pos: float peak position in deg 2-theta or tof in musec
833	:param Inst: dict instrument parameters
834
835	:returns float: total FWHM of pseudoVoigt in deg or musec
836	'''
837
838	sig = lambda Th,U,V,W: np.sqrt(max(0.001,Utand(Th)2+Vtand(Th)+W))
839	sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1dsp2+S2dsp*4+Sqdsp)
840	gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th)
841	gamTOF = lambda dsp,X,Y,Z: Z+Xdsp+Ydsp**2
842	alpTOF = lambda dsp,alp: alp/dsp
843	betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp4+betq/dsp2
844	alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.)
845	betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.)
846	if 'T' in Inst['Type'][0]:
847	dsp = pos/Inst['difC'][1]
848	alp = alpTOF(dsp,Inst['alpha'][1])
849	bet = betTOF(dsp,Inst['beta-0'][1],Inst['beta-1'][1],Inst['beta-q'][1])
850	s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
851	g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
852	return getgamFW(g,s)+np.log(2.0)(alp+bet)/(alpbet)
853	elif 'C' in Inst['Type'][0]:
854	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
855	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
856	return getgamFW(g,s)/100. #returns FWHM in deg
857	else: #'B'
858	alp = alpPink(pos,Inst['alpha-0'][1],Inst['alpha-1'][1])
859	bet = betPink(pos,Inst['beta-0'][1],Inst['beta-1'][1])
860	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
861	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
862	return getgamFW(g,s)/100.+np.log(2.0)(alp+bet)/(alpbet) #returns FWHM in deg
863
864	def getgamFW(g,s):
865	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
866	lambda fxn needs FWHM for both Gaussian & Lorentzian components
867
868	:param g: float Lorentzian gamma = FWHM(L)
869	:param s: float Gaussian sig
870
871	:returns float: total FWHM of pseudoVoigt
872	'''
873	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.)
874	return gamFW(2.35482s,g) #sqrt(8ln2)sig = FWHM(G)
875
876	def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None):
877	'''Computes the background from vars pulled from gpx file or tree.
878	'''
879	if 'T' in dataType:
880	q = 2.np.piparmDict[pfx+'difC']/xdata
881	else:
882	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
883	q = npT2q(xdata,wave)
884	yb = np.zeros_like(xdata)
885	nBak = 0
886	cw = np.diff(xdata)
887	cw = np.append(cw,cw[-1])
888	sumBk = [0.,0.,0]
889	while True:
890	key = pfx+'Back;'+str(nBak)
891	if key in parmDict:
892	nBak += 1
893	else:
894	break
895	#empirical functions
896	if bakType in ['chebyschev','cosine','chebyschev-1']:
897	dt = xdata[-1]-xdata[0]
898	for iBak in range(nBak):
899	key = pfx+'Back;'+str(iBak)
900	if bakType == 'chebyschev':
901	ybi = parmDict[key](-1.+2.(xdata-xdata[0])/dt)**iBak
902	elif bakType == 'chebyschev-1':
903	xpos = -1.+2.*(xdata-xdata[0])/dt
904	ybi = parmDict[key]np.cos(iBaknp.arccos(xpos))
905	elif bakType == 'cosine':
906	ybi = parmDict[key]npcosd(180.xdata*iBak/xdata[-1])
907	yb += ybi
908	sumBk[0] = np.sum(yb)
909	elif bakType in ['Q^2 power series','Q^-2 power series']:
910	QT = 1.
911	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
912	for iBak in range(nBak-1):
913	key = pfx+'Back;'+str(iBak+1)
914	if '-2' in bakType:
915	QT = (iBak+1)q**-2
916	else:
917	QT = q*2/(iBak+1)
918	yb += QT*parmDict[key]
919	sumBk[0] = np.sum(yb)
920	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
921	if nBak == 1:
922	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
923	elif nBak == 2:
924	dX = xdata[-1]-xdata[0]
925	T2 = (xdata-xdata[0])/dX
926	T1 = 1.0-T2
927	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
928	else:
929	xnomask = ma.getdata(xdata)
930	xmin,xmax = xnomask[0],xnomask[-1]
931	if bakType == 'lin interpolate':
932	bakPos = np.linspace(xmin,xmax,nBak,True)
933	elif bakType == 'inv interpolate':
934	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
935	elif bakType == 'log interpolate':
936	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
937	bakPos[0] = xmin
938	bakPos[-1] = xmax
939	bakVals = np.zeros(nBak)
940	for i in range(nBak):
941	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
942	bakInt = si.interp1d(bakPos,bakVals,'linear')
943	yb = bakInt(ma.getdata(xdata))
944	sumBk[0] = np.sum(yb)
945	#Debye function
946	if pfx+'difC' in parmDict:
947	ff = 1.
948	else:
949	try:
950	wave = parmDict[pfx+'Lam']
951	except KeyError:
952	wave = parmDict[pfx+'Lam1']
953	SQ = (q/(4.np.pi))*2
954	FF = G2elem.GetFormFactorCoeff('Si')[0]
955	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
956	iD = 0
957	while True:
958	try:
959	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
960	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
961	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
962	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
963	yb += ybi
964	sumBk[1] += np.sum(ybi)
965	iD += 1
966	except KeyError:
967	break
968	#peaks
969	iD = 0
970	while True:
971	try:
972	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
973	pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
974	pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],1.)
975	pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
976	if 'C' in dataType:
977	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
978	else: #'T'OF
979	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
980	iBeg = np.searchsorted(xdata,pkP-fmin)
981	iFin = np.searchsorted(xdata,pkP+fmax)
982	lenX = len(xdata)
983	if not iBeg:
984	iFin = np.searchsorted(xdata,pkP+fmax)
985	elif iBeg == lenX:
986	iFin = iBeg
987	else:
988	iFin = np.searchsorted(xdata,pkP+fmax)
989	if 'C' in dataType:
990	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0]
991	yb[iBeg:iFin] += ybi/cw[iBeg:iFin]
992	elif 'T' in dataType:
993	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0]
994	yb[iBeg:iFin] += ybi
995	elif 'B' in dataType:
996	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0]
997	yb[iBeg:iFin] += ybi
998	sumBk[2] += np.sum(ybi)
999	iD += 1
1000	except KeyError:
1001	break
1002	except ValueError:
1003	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
1004	break
1005	if fixback is not None:
1006	yb += parmDict[pfx+'BF mult']*fixback
1007	sumBk[0] = sum(yb)
1008	return yb,sumBk
1009
1010	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata,fixback=None):
1011	'needs a doc string'
1012	if 'T' in dataType:
1013	q = 2.np.piparmDict[hfx+'difC']/xdata
1014	else:
1015	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1016	q = 2.np.pinpsind(xdata/2.)/wave
1017	nBak = 0
1018	while True:
1019	key = hfx+'Back;'+str(nBak)
1020	if key in parmDict:
1021	nBak += 1
1022	else:
1023	break
1024	dydb = np.zeros(shape=(nBak,len(xdata)))
1025	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
1026	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
1027	dydfb = []
1028	cw = np.diff(xdata)
1029	cw = np.append(cw,cw[-1])
1030
1031	if bakType in ['chebyschev','cosine','chebyschev-1']:
1032	dt = xdata[-1]-xdata[0]
1033	for iBak in range(nBak):
1034	if bakType == 'chebyschev':
1035	dydb[iBak] = (-1.+2.(xdata-xdata[0])/dt)*iBak
1036	elif bakType == 'chebyschev-1':
1037	xpos = -1.+2.*(xdata-xdata[0])/dt
1038	dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1039	elif bakType == 'cosine':
1040	dydb[iBak] = npcosd(180.xdataiBak/xdata[-1])
1041	elif bakType in ['Q^2 power series','Q^-2 power series']:
1042	QT = 1.
1043	dydb[0] = np.ones_like(xdata)
1044	for iBak in range(nBak-1):
1045	if '-2' in bakType:
1046	QT = (iBak+1)q**-2
1047	else:
1048	QT = q*2/(iBak+1)
1049	dydb[iBak+1] = QT
1050	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1051	if nBak == 1:
1052	dydb[0] = np.ones_like(xdata)
1053	elif nBak == 2:
1054	dX = xdata[-1]-xdata[0]
1055	T2 = (xdata-xdata[0])/dX
1056	T1 = 1.0-T2
1057	dydb = [T1,T2]
1058	else:
1059	xnomask = ma.getdata(xdata)
1060	xmin,xmax = xnomask[0],xnomask[-1]
1061	if bakType == 'lin interpolate':
1062	bakPos = np.linspace(xmin,xmax,nBak,True)
1063	elif bakType == 'inv interpolate':
1064	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1065	elif bakType == 'log interpolate':
1066	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1067	bakPos[0] = xmin
1068	bakPos[-1] = xmax
1069	for i,pos in enumerate(bakPos):
1070	if i == 0:
1071	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1072	elif i == len(bakPos)-1:
1073	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1074	else:
1075	dydb[i] = np.where(xdata>bakPos[i],
1076	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1077	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1078	if hfx+'difC' in parmDict:
1079	ff = 1.
1080	else:
1081	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1082	q = npT2q(xdata,wave)
1083	SQ = (q/(4np.pi))*2
1084	FF = G2elem.GetFormFactorCoeff('Si')[0]
1085	ff = np.array(G2elem.ScatFac(FF,SQ)[0])np.pi*2 #needs pi^2~10. for cw data (why?)
1086	iD = 0
1087	while True:
1088	try:
1089	if hfx+'difC' in parmDict:
1090	q = 2np.piparmDict[hfx+'difC']/xdata
1091	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1092	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1093	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1094	sqr = np.sin(qdbR)/(qdbR)
1095	cqr = np.cos(q*dbR)
1096	temp = np.exp(-dbUq*2)
1097	dyddb[3iD] = ffsqr*temp
1098	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(dbR)
1099	dyddb[3iD+2] = -ffdbAsqrtempq*2
1100	iD += 1
1101	except KeyError:
1102	break
1103	iD = 0
1104	while True:
1105	try:
1106	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1107	pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1108	pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],1.0)
1109	pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1110	if 'C' in dataType:
1111	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1112	else: #'T' or 'B'
1113	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1114	iBeg = np.searchsorted(xdata,pkP-fmin)
1115	iFin = np.searchsorted(xdata,pkP+fmax)
1116	lenX = len(xdata)
1117	if not iBeg:
1118	iFin = np.searchsorted(xdata,pkP+fmax)
1119	elif iBeg == lenX:
1120	iFin = iBeg
1121	else:
1122	iFin = np.searchsorted(xdata,pkP+fmax)
1123	if 'C' in dataType:
1124	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1125	# dydpk[4iD][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdp
1126	# dydpk[4iD+1][iBeg:iFin] += 100.cw[iBeg:iFin]*Df
1127	# dydpk[4iD+2][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFds
1128	# dydpk[4iD+3][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdg
1129	dydpk[4iD][iBeg:iFin] += 1000.pkI*dFdp
1130	dydpk[4iD+1][iBeg:iFin] += 1000.Df
1131	dydpk[4iD+2][iBeg:iFin] += 1000.pkI*dFds
1132	dydpk[4iD+3][iBeg:iFin] += 1000.pkI*dFdg
1133	else: #'T'OF
1134	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1135	dydpk[4iD][iBeg:iFin] += pkIdFdp
1136	dydpk[4*iD+1][iBeg:iFin] += Df
1137	dydpk[4iD+2][iBeg:iFin] += pkIdFds
1138	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
1139	iD += 1
1140	except KeyError:
1141	break
1142	except ValueError:
1143	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
1144	break
1145	# fixed background from file
1146	if fixback is not None:
1147	dydfb = fixback
1148	return dydb,dyddb,dydpk,dydfb
1149
1150	#### Using old gsas fortran routines for powder peak shapes & derivatives
1151	def getFCJVoigt3(pos,sig,gam,shl,xdata):
1152	'''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1153	CW powder peak in external Fortran routine
1154
1155	param pos: peak position in degrees
1156	param sig: Gaussian variance in centideg^2
1157	param gam: Lorentzian width in centideg
1158	param shl: axial divergence parameter (S+H)/L
1159	param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1160
1161	returns: array: calculated peak function at each xdata
1162	returns: integral of peak; nominally = 1.0
1163	'''
1164	if len(xdata):
1165	cw = np.diff(xdata)
1166	cw = np.append(cw,cw[-1])
1167	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1168	return Df,np.sum(100.Dfcw)
1169	else:
1170	return 0.,1.
1171
1172	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1173	'''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1174	function for a CW powder peak
1175
1176	param pos: peak position in degrees
1177	param sig: Gaussian variance in centideg^2
1178	param gam: Lorentzian width in centideg
1179	param shl: axial divergence parameter (S+H)/L
1180	param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1181
1182	returns: arrays: function and derivatives of pos, sig, gam, & shl
1183	'''
1184	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1185	return Df,dFdp,dFds,dFdg,dFdsh
1186
1187	def getPsVoigt(pos,sig,gam,xdata):
1188	'''Compute the simple Pseudo-Voigt function for a
1189	small angle Bragg peak in external Fortran routine
1190
1191	param pos: peak position in degrees
1192	param sig: Gaussian variance in centideg^2
1193	param gam: Lorentzian width in centideg
1194
1195	returns: array: calculated peak function at each xdata
1196	returns: integral of peak; nominally = 1.0
1197	'''
1198
1199	cw = np.diff(xdata)
1200	cw = np.append(cw,cw[-1])
1201	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1202	return Df,np.sum(100.Dfcw)
1203
1204	def getdPsVoigt(pos,sig,gam,xdata):
1205	'''Compute the simple Pseudo-Voigt function derivatives for a
1206	small angle Bragg peak peak in external Fortran routine
1207
1208	param pos: peak position in degrees
1209	param sig: Gaussian variance in centideg^2
1210	param gam: Lorentzian width in centideg
1211
1212	returns: arrays: function and derivatives of pos, sig, gam, & shl
1213	'''
1214
1215	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1216	return Df,dFdp,dFds,dFdg
1217
1218	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1219	'''Compute the double exponential Pseudo-Voigt convolution function for a
1220	neutron TOF & CW pink peak in external Fortran routine
1221	'''
1222
1223	cw = np.diff(xdata)
1224	cw = np.append(cw,cw[-1])
1225	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1226	return Df,np.sum(Df*cw)
1227
1228	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1229	'''Compute the double exponential Pseudo-Voigt convolution function derivatives for a
1230	neutron TOF & CW pink peak in external Fortran routine
1231	'''
1232
1233	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1234	return Df,dFdp,dFda,dFdb,dFds,dFdg
1235
1236	def ellipseSize(H,Sij,GB):
1237	'''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady
1238	'''
1239
1240	HX = np.inner(H.T,GB)
1241	lenHX = np.sqrt(np.sum(HX**2))
1242	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
1243	R = np.inner(HX/lenHX,Rsize)*2Esize #want column length for hkl in crystal
1244	lenR = 1./np.sqrt(np.sum(R))
1245	return lenR
1246
1247	def ellipseSizeDerv(H,Sij,GB):
1248	'''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady
1249	'''
1250
1251	lenR = ellipseSize(H,Sij,GB)
1252	delt = 0.001
1253	dRdS = np.zeros(6)
1254	for i in range(6):
1255	Sij[i] -= delt
1256	lenM = ellipseSize(H,Sij,GB)
1257	Sij[i] += 2.*delt
1258	lenP = ellipseSize(H,Sij,GB)
1259	Sij[i] -= delt
1260	dRdS[i] = (lenP-lenM)/(2.*delt)
1261	return lenR,dRdS
1262
1263	def getMustrain(HKL,G,SGData,muStrData):
1264	if muStrData[0] == 'isotropic':
1265	return np.ones(HKL.shape[1])*muStrData[1][0]
1266	elif muStrData[0] == 'uniaxial':
1267	H = np.array(HKL)
1268	P = np.array(muStrData[3])
1269	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1270	Si = muStrData[1][0]
1271	Sa = muStrData[1][1]
1272	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1273	else: #generalized - P.W. Stephens model
1274	H = np.array(HKL)
1275	rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1276	Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1277	Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1278	return np.sqrt(Sum)/rdsq
1279
1280	def getCrSize(HKL,G,GB,sizeData):
1281	if sizeData[0] == 'isotropic':
1282	return np.ones(HKL.shape[1])*sizeData[1][0]
1283	elif sizeData[0] == 'uniaxial':
1284	H = np.array(HKL)
1285	P = np.array(sizeData[3])
1286	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1287	Si = sizeData[1][0]
1288	Sa = sizeData[1][1]
1289	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1290	else:
1291	Sij =[sizeData[4][i] for i in range(6)]
1292	H = np.array(HKL)
1293	return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1294
1295	def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1296	'''
1297	Generates allowed by symmetry reflections with d >= dmin
1298	NB: GenHKLf & checkMagextc return True for extinct reflections
1299
1300	:param dmin: minimum d-spacing
1301	:param SGData: space group data obtained from SpcGroup
1302	:param A: lattice parameter terms A1-A6
1303	:param Inst: instrument parameter info
1304	:returns: HKLs: np.array hkl, etc for allowed reflections
1305
1306	'''
1307	HKL = G2lat.GenHLaue(dmin,SGData,A)
1308	HKLs = []
1309	ds = []
1310	for h,k,l,d in HKL:
1311	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1312	if ext and 'MagSpGrp' in SGData:
1313	ext = G2spc.checkMagextc([h,k,l],SGData)
1314	if not ext:
1315	if nodup and int(10000*d) in ds:
1316	continue
1317	ds.append(int(10000*d))
1318	if Inst == None:
1319	HKLs.append([h,k,l,d,0,-1])
1320	else:
1321	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1322	return np.array(HKLs)
1323
1324	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1325	'needs a doc string'
1326	HKLs = []
1327	vec = np.array(Vec)
1328	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1329	dvec = 1./(maxH*vstar+1./dmin)
1330	HKL = G2lat.GenHLaue(dvec,SGData,A)
1331	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1332	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1333	ifMag = False
1334	if 'MagSpGrp' in SGData:
1335	ifMag = True
1336	for h,k,l,d in HKL:
1337	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1338	if not ext and d >= dmin:
1339	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1340	for dH in SSdH:
1341	if dH:
1342	DH = SSdH[dH]
1343	H = [h+DH[0],k+DH[1],l+DH[2]]
1344	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1345	if d >= dmin:
1346	HKLM = np.array([h,k,l,dH])
1347	if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1348	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
1349	return G2lat.sortHKLd(HKLs,True,True,True)
1350
1351	peakInstPrmMode = True
1352	'''Determines the mode used for peak fitting. When peakInstPrmMode=True peak
1353	width parameters are computed from the instrument parameters (UVW,... or
1354	alpha,... etc) unless the individual parameter is refined. This allows the
1355	instrument parameters to be refined. When peakInstPrmMode=False, the instrument
1356	parameters are not used and cannot be refined.
1357	The default is peakFitMode=True.
1358	'''
1359
1360	def setPeakInstPrmMode(normal=True):
1361	'''Determines the mode used for peak fitting. If normal=True (default)
1362	peak width parameters are computed from the instrument parameters
1363	unless the individual parameter is refined. If normal=False,
1364	peak widths are used as supplied for each peak.
1365
1366	Note that normal=True unless this routine is called. Also,
1367	instrument parameters can only be refined with normal=True.
1368
1369	:param bool normal: setting to apply to global variable
1370	:data:`peakInstPrmMode`
1371	'''
1372	global peakInstPrmMode
1373	peakInstPrmMode = normal
1374
1375	def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType):
1376	'''Computes the profile for a powder pattern for single peak fitting
1377	NB: not used for Rietveld refinement
1378	'''
1379
1380	yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
1381	yc = np.zeros_like(yb)
1382	# cw = np.diff(xdata)
1383	# cw = np.append(cw,cw[-1])
1384	if 'C' in dataType:
1385	shl = max(parmDict['SH/L'],0.002)
1386	Ka2 = False
1387	if 'Lam1' in parmDict.keys():
1388	Ka2 = True
1389	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1390	kRatio = parmDict['I(L2)/I(L1)']
1391	iPeak = 0
1392	while True:
1393	try:
1394	pos = parmDict['pos'+str(iPeak)]
1395	tth = (pos-parmDict['Zero'])
1396	intens = parmDict['int'+str(iPeak)]
1397	sigName = 'sig'+str(iPeak)
1398	if sigName in varyList or not peakInstPrmMode:
1399	sig = parmDict[sigName]
1400	else:
1401	sig = G2mth.getCWsig(parmDict,tth)
1402	sig = max(sig,0.001) #avoid neg sigma^2
1403	gamName = 'gam'+str(iPeak)
1404	if gamName in varyList or not peakInstPrmMode:
1405	gam = parmDict[gamName]
1406	else:
1407	gam = G2mth.getCWgam(parmDict,tth)
1408	gam = max(gam,0.001) #avoid neg gamma
1409	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1410	iBeg = np.searchsorted(xdata,pos-fmin)
1411	iFin = np.searchsorted(xdata,pos+fmin)
1412	if not iBeg+iFin: #peak below low limit
1413	iPeak += 1
1414	continue
1415	elif not iBeg-iFin: #peak above high limit
1416	return yb+yc
1417	fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0]
1418	yc[iBeg:iFin] += intens*fp
1419	if Ka2:
1420	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1421	iBeg = np.searchsorted(xdata,pos2-fmin)
1422	iFin = np.searchsorted(xdata,pos2+fmin)
1423	if iBeg-iFin:
1424	fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0]
1425	yc[iBeg:iFin] += intenskRatiofp2
1426	iPeak += 1
1427	except KeyError: #no more peaks to process
1428	return yb+yc
1429	elif 'B' in dataType:
1430	iPeak = 0
1431	dsp = 1.0 #for now - fix later
1432	while True:
1433	try:
1434	pos = parmDict['pos'+str(iPeak)]
1435	tth = (pos-parmDict['Zero'])
1436	intens = parmDict['int'+str(iPeak)]
1437	alpName = 'alp'+str(iPeak)
1438	if alpName in varyList or not peakInstPrmMode:
1439	alp = parmDict[alpName]
1440	else:
1441	alp = G2mth.getPinkalpha(parmDict,tth)
1442	alp = max(0.1,alp)
1443	betName = 'bet'+str(iPeak)
1444	if betName in varyList or not peakInstPrmMode:
1445	bet = parmDict[betName]
1446	else:
1447	bet = G2mth.getPinkbeta(parmDict,tth)
1448	bet = max(0.1,bet)
1449	sigName = 'sig'+str(iPeak)
1450	if sigName in varyList or not peakInstPrmMode:
1451	sig = parmDict[sigName]
1452	else:
1453	sig = G2mth.getCWsig(parmDict,tth)
1454	sig = max(sig,0.001) #avoid neg sigma^2
1455	gamName = 'gam'+str(iPeak)
1456	if gamName in varyList or not peakInstPrmMode:
1457	gam = parmDict[gamName]
1458	else:
1459	gam = G2mth.getCWgam(parmDict,tth)
1460	gam = max(gam,0.001) #avoid neg gamma
1461	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1462	iBeg = np.searchsorted(xdata,pos-fmin)
1463	iFin = np.searchsorted(xdata,pos+fmin)
1464	if not iBeg+iFin: #peak below low limit
1465	iPeak += 1
1466	continue
1467	elif not iBeg-iFin: #peak above high limit
1468	return yb+yc
1469	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0]
1470	iPeak += 1
1471	except KeyError: #no more peaks to process
1472	return yb+yc
1473	else:
1474	Pdabc = parmDict['Pdabc']
1475	difC = parmDict['difC']
1476	iPeak = 0
1477	while True:
1478	try:
1479	pos = parmDict['pos'+str(iPeak)]
1480	tof = pos-parmDict['Zero']
1481	dsp = tof/difC
1482	intens = parmDict['int'+str(iPeak)]
1483	alpName = 'alp'+str(iPeak)
1484	if alpName in varyList or not peakInstPrmMode:
1485	alp = parmDict[alpName]
1486	else:
1487	if len(Pdabc):
1488	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1489	else:
1490	alp = G2mth.getTOFalpha(parmDict,dsp)
1491	alp = max(0.1,alp)
1492	betName = 'bet'+str(iPeak)
1493	if betName in varyList or not peakInstPrmMode:
1494	bet = parmDict[betName]
1495	else:
1496	if len(Pdabc):
1497	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1498	else:
1499	bet = G2mth.getTOFbeta(parmDict,dsp)
1500	bet = max(0.0001,bet)
1501	sigName = 'sig'+str(iPeak)
1502	if sigName in varyList or not peakInstPrmMode:
1503	sig = parmDict[sigName]
1504	else:
1505	sig = G2mth.getTOFsig(parmDict,dsp)
1506	gamName = 'gam'+str(iPeak)
1507	if gamName in varyList or not peakInstPrmMode:
1508	gam = parmDict[gamName]
1509	else:
1510	gam = G2mth.getTOFgamma(parmDict,dsp)
1511	gam = max(gam,0.001) #avoid neg gamma
1512	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1513	iBeg = np.searchsorted(xdata,pos-fmin)
1514	iFin = np.searchsorted(xdata,pos+fmax)
1515	lenX = len(xdata)
1516	if not iBeg:
1517	iFin = np.searchsorted(xdata,pos+fmax)
1518	elif iBeg == lenX:
1519	iFin = iBeg
1520	else:
1521	iFin = np.searchsorted(xdata,pos+fmax)
1522	if not iBeg+iFin: #peak below low limit
1523	iPeak += 1
1524	continue
1525	elif not iBeg-iFin: #peak above high limit
1526	return yb+yc
1527	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0]
1528	iPeak += 1
1529	except KeyError: #no more peaks to process
1530	return yb+yc
1531
1532	def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType):
1533	'''Computes the profile derivatives for a powder pattern for single peak fitting
1534
1535	return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1536
1537	NB: not used for Rietveld refinement
1538	'''
1539	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1540	dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
1541	if 'Back;0' in varyList: #background derivs are in front if present
1542	dMdv[0:len(dMdb)] = dMdb
1543	names = ['DebyeA','DebyeR','DebyeU']
1544	for name in varyList:
1545	if 'Debye' in name:
1546	parm,Id = name.split(';')
1547	ip = names.index(parm)
1548	dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1549	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1550	for name in varyList:
1551	if 'BkPk' in name:
1552	parm,Id = name.split(';')
1553	ip = names.index(parm)
1554	dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1555	cw = np.diff(xdata)
1556	cw = np.append(cw,cw[-1])
1557	if 'C' in dataType:
1558	shl = max(parmDict['SH/L'],0.002)
1559	Ka2 = False
1560	if 'Lam1' in parmDict.keys():
1561	Ka2 = True
1562	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1563	kRatio = parmDict['I(L2)/I(L1)']
1564	iPeak = 0
1565	while True:
1566	try:
1567	pos = parmDict['pos'+str(iPeak)]
1568	tth = (pos-parmDict['Zero'])
1569	intens = parmDict['int'+str(iPeak)]
1570	sigName = 'sig'+str(iPeak)
1571	if sigName in varyList or not peakInstPrmMode:
1572	sig = parmDict[sigName]
1573	dsdU = dsdV = dsdW = 0
1574	else:
1575	sig = G2mth.getCWsig(parmDict,tth)
1576	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1577	sig = max(sig,0.001) #avoid neg sigma
1578	gamName = 'gam'+str(iPeak)
1579	if gamName in varyList or not peakInstPrmMode:
1580	gam = parmDict[gamName]
1581	dgdX = dgdY = dgdZ = 0
1582	else:
1583	gam = G2mth.getCWgam(parmDict,tth)
1584	dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1585	gam = max(gam,0.001) #avoid neg gamma
1586	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1587	iBeg = np.searchsorted(xdata,pos-fmin)
1588	iFin = np.searchsorted(xdata,pos+fmin)
1589	if not iBeg+iFin: #peak below low limit
1590	iPeak += 1
1591	continue
1592	elif not iBeg-iFin: #peak above high limit
1593	break
1594	dMdpk = np.zeros(shape=(6,len(xdata)))
1595	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1596	for i in range(1,5):
1597	dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1598	dMdpk[0][iBeg:iFin] += dMdipk[0]
1599	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1600	if Ka2:
1601	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1602	iBeg = np.searchsorted(xdata,pos2-fmin)
1603	iFin = np.searchsorted(xdata,pos2+fmin)
1604	if iBeg-iFin:
1605	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1606	for i in range(1,5):
1607	dMdpk[i][iBeg:iFin] += intenskRatiodMdipk2[i]
1608	dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0]
1609	dMdpk[5][iBeg:iFin] += dMdipk2[0]
1610	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1611	for parmName in ['pos','int','sig','gam']:
1612	try:
1613	idx = varyList.index(parmName+str(iPeak))
1614	dMdv[idx] = dervDict[parmName]
1615	except ValueError:
1616	pass
1617	if 'U' in varyList:
1618	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1619	if 'V' in varyList:
1620	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1621	if 'W' in varyList:
1622	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1623	if 'X' in varyList:
1624	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1625	if 'Y' in varyList:
1626	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1627	if 'Z' in varyList:
1628	dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1629	if 'SH/L' in varyList:
1630	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1631	if 'I(L2)/I(L1)' in varyList:
1632	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1633	iPeak += 1
1634	except KeyError: #no more peaks to process
1635	break
1636	elif 'B' in dataType:
1637	iPeak = 0
1638	while True:
1639	try:
1640	pos = parmDict['pos'+str(iPeak)]
1641	tth = (pos-parmDict['Zero'])
1642	intens = parmDict['int'+str(iPeak)]
1643	alpName = 'alp'+str(iPeak)
1644	if alpName in varyList or not peakInstPrmMode:
1645	alp = parmDict[alpName]
1646	dada0 = dada1 = 0.0
1647	else:
1648	alp = G2mth.getPinkalpha(parmDict,tth)
1649	dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
1650	alp = max(0.0001,alp)
1651	betName = 'bet'+str(iPeak)
1652	if betName in varyList or not peakInstPrmMode:
1653	bet = parmDict[betName]
1654	dbdb0 = dbdb1 = 0.0
1655	else:
1656	bet = G2mth.getPinkbeta(parmDict,tth)
1657	dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
1658	bet = max(0.0001,bet)
1659	sigName = 'sig'+str(iPeak)
1660	if sigName in varyList or not peakInstPrmMode:
1661	sig = parmDict[sigName]
1662	dsdU = dsdV = dsdW = 0
1663	else:
1664	sig = G2mth.getCWsig(parmDict,tth)
1665	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1666	sig = max(sig,0.001) #avoid neg sigma
1667	gamName = 'gam'+str(iPeak)
1668	if gamName in varyList or not peakInstPrmMode:
1669	gam = parmDict[gamName]
1670	dgdX = dgdY = dgdZ = 0
1671	else:
1672	gam = G2mth.getCWgam(parmDict,tth)
1673	dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1674	gam = max(gam,0.001) #avoid neg gamma
1675	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
1676	iBeg = np.searchsorted(xdata,pos-fmin)
1677	iFin = np.searchsorted(xdata,pos+fmin)
1678	if not iBeg+iFin: #peak below low limit
1679	iPeak += 1
1680	continue
1681	elif not iBeg-iFin: #peak above high limit
1682	break
1683	dMdpk = np.zeros(shape=(7,len(xdata)))
1684	dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1685	for i in range(1,6):
1686	dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]intensdMdipk[i]
1687	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1688	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
1689	for parmName in ['pos','int','alp','bet','sig','gam']:
1690	try:
1691	idx = varyList.index(parmName+str(iPeak))
1692	dMdv[idx] = dervDict[parmName]
1693	except ValueError:
1694	pass
1695	if 'U' in varyList:
1696	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1697	if 'V' in varyList:
1698	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1699	if 'W' in varyList:
1700	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1701	if 'X' in varyList:
1702	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1703	if 'Y' in varyList:
1704	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1705	if 'Z' in varyList:
1706	dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1707	if 'alpha-0' in varyList:
1708	dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
1709	if 'alpha-1' in varyList:
1710	dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
1711	if 'beta-0' in varyList:
1712	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1713	if 'beta-1' in varyList:
1714	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1715	iPeak += 1
1716	except KeyError: #no more peaks to process
1717	break
1718	else:
1719	Pdabc = parmDict['Pdabc']
1720	difC = parmDict['difC']
1721	iPeak = 0
1722	while True:
1723	try:
1724	pos = parmDict['pos'+str(iPeak)]
1725	tof = pos-parmDict['Zero']
1726	dsp = tof/difC
1727	intens = parmDict['int'+str(iPeak)]
1728	alpName = 'alp'+str(iPeak)
1729	if alpName in varyList or not peakInstPrmMode:
1730	alp = parmDict[alpName]
1731	else:
1732	if len(Pdabc):
1733	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1734	dada0 = 0
1735	else:
1736	alp = G2mth.getTOFalpha(parmDict,dsp)
1737	dada0 = G2mth.getTOFalphaDeriv(dsp)
1738	betName = 'bet'+str(iPeak)
1739	if betName in varyList or not peakInstPrmMode:
1740	bet = parmDict[betName]
1741	else:
1742	if len(Pdabc):
1743	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1744	dbdb0 = dbdb1 = dbdb2 = 0
1745	else:
1746	bet = G2mth.getTOFbeta(parmDict,dsp)
1747	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1748	sigName = 'sig'+str(iPeak)
1749	if sigName in varyList or not peakInstPrmMode:
1750	sig = parmDict[sigName]
1751	dsds0 = dsds1 = dsds2 = dsds3 = 0
1752	else:
1753	sig = G2mth.getTOFsig(parmDict,dsp)
1754	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1755	gamName = 'gam'+str(iPeak)
1756	if gamName in varyList or not peakInstPrmMode:
1757	gam = parmDict[gamName]
1758	dsdX = dsdY = dsdZ = 0
1759	else:
1760	gam = G2mth.getTOFgamma(parmDict,dsp)
1761	dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1762	gam = max(gam,0.001) #avoid neg gamma
1763	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1764	iBeg = np.searchsorted(xdata,pos-fmin)
1765	lenX = len(xdata)
1766	if not iBeg:
1767	iFin = np.searchsorted(xdata,pos+fmax)
1768	elif iBeg == lenX:
1769	iFin = iBeg
1770	else:
1771	iFin = np.searchsorted(xdata,pos+fmax)
1772	if not iBeg+iFin: #peak below low limit
1773	iPeak += 1
1774	continue
1775	elif not iBeg-iFin: #peak above high limit
1776	break
1777	dMdpk = np.zeros(shape=(7,len(xdata)))
1778	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1779	for i in range(1,6):
1780	dMdpk[i][iBeg:iFin] += intenscw[iBeg:iFin]dMdipk[i]
1781	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1782	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1783	for parmName in ['pos','int','alp','bet','sig','gam']:
1784	try:
1785	idx = varyList.index(parmName+str(iPeak))
1786	dMdv[idx] = dervDict[parmName]
1787	except ValueError:
1788	pass
1789	if 'alpha' in varyList:
1790	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1791	if 'beta-0' in varyList:
1792	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1793	if 'beta-1' in varyList:
1794	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1795	if 'beta-q' in varyList:
1796	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1797	if 'sig-0' in varyList:
1798	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1799	if 'sig-1' in varyList:
1800	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1801	if 'sig-2' in varyList:
1802	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1803	if 'sig-q' in varyList:
1804	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1805	if 'X' in varyList:
1806	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1807	if 'Y' in varyList:
1808	dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1809	if 'Z' in varyList:
1810	dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1811	iPeak += 1
1812	except KeyError: #no more peaks to process
1813	break
1814	if 'BF mult' in varyList:
1815	dMdv[varyList.index('BF mult')] = fixback
1816
1817	return dMdv
1818
1819	def Dict2Values(parmdict, varylist):
1820	'''Use before call to leastsq to setup list of values for the parameters
1821	in parmdict, as selected by key in varylist'''
1822	return [parmdict[key] for key in varylist]
1823
1824	def Values2Dict(parmdict, varylist, values):
1825	''' Use after call to leastsq to update the parameter dictionary with
1826	values corresponding to keys in varylist'''
1827	parmdict.update(zip(varylist,values))
1828
1829	def SetBackgroundParms(Background):
1830	'Loads background parameters into dicts/lists to create varylist & parmdict'
1831	if len(Background) == 1: # fix up old backgrounds
1832	Background.append({'nDebye':0,'debyeTerms':[]})
1833	bakType,bakFlag = Background[0][:2]
1834	backVals = Background[0][3:]
1835	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1836	Debye = Background[1] #also has background peaks stuff
1837	backDict = dict(zip(backNames,backVals))
1838	backVary = []
1839	if bakFlag:
1840	backVary = backNames
1841
1842	backDict['nDebye'] = Debye['nDebye']
1843	debyeDict = {}
1844	debyeList = []
1845	for i in range(Debye['nDebye']):
1846	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1847	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1848	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1849	debyeVary = []
1850	for item in debyeList:
1851	if item[1]:
1852	debyeVary.append(item[0])
1853	backDict.update(debyeDict)
1854	backVary += debyeVary
1855
1856	backDict['nPeaks'] = Debye['nPeaks']
1857	peaksDict = {}
1858	peaksList = []
1859	for i in range(Debye['nPeaks']):
1860	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1861	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1862	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1863	peaksVary = []
1864	for item in peaksList:
1865	if item[1]:
1866	peaksVary.append(item[0])
1867	backDict.update(peaksDict)
1868	backVary += peaksVary
1869	if 'background PWDR' in Background[1]:
1870	backDict['Back File'] = Background[1]['background PWDR'][0]
1871	backDict['BF mult'] = Background[1]['background PWDR'][1]
1872	if len(Background[1]['background PWDR']) > 2:
1873	if Background[1]['background PWDR'][2]:
1874	backVary += ['BF mult',]
1875	return bakType,backDict,backVary
1876
1877	def DoCalibInst(IndexPeaks,Inst):
1878
1879	def SetInstParms():
1880	dataType = Inst['Type'][0]
1881	insVary = []
1882	insNames = []
1883	insVals = []
1884	for parm in Inst:
1885	insNames.append(parm)
1886	insVals.append(Inst[parm][1])
1887	if parm in ['Lam','difC','difA','difB','Zero',]:
1888	if Inst[parm][2]:
1889	insVary.append(parm)
1890	instDict = dict(zip(insNames,insVals))
1891	return dataType,instDict,insVary
1892
1893	def GetInstParms(parmDict,Inst,varyList):
1894	for name in Inst:
1895	Inst[name][1] = parmDict[name]
1896
1897	def InstPrint(Inst,sigDict):
1898	print ('Instrument Parameters:')
1899	if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]:
1900	ptfmt = "%12.6f"
1901	else:
1902	ptfmt = "%12.3f"
1903	ptlbls = 'names :'
1904	ptstr = 'values:'
1905	sigstr = 'esds :'
1906	for parm in Inst:
1907	if parm in ['Lam','difC','difA','difB','Zero',]:
1908	ptlbls += "%s" % (parm.center(12))
1909	ptstr += ptfmt % (Inst[parm][1])
1910	if parm in sigDict:
1911	sigstr += ptfmt % (sigDict[parm])
1912	else:
1913	sigstr += 12*' '
1914	print (ptlbls)
1915	print (ptstr)
1916	print (sigstr)
1917
1918	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1919	parmDict.update(zip(varyList,values))
1920	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1921
1922	peakPos = []
1923	peakDsp = []
1924	peakWt = []
1925	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1926	if peak[2] and peak[3] and sig > 0.:
1927	peakPos.append(peak[0])
1928	peakDsp.append(peak[-1]) #d-calc
1929	# peakWt.append(peak[-1]2/sig2) #weight by d**2
1930	peakWt.append(1./(sig*peak[-1])) #
1931	peakPos = np.array(peakPos)
1932	peakDsp = np.array(peakDsp)
1933	peakWt = np.array(peakWt)
1934	dataType,insDict,insVary = SetInstParms()
1935	parmDict = {}
1936	parmDict.update(insDict)
1937	varyList = insVary
1938	if not len(varyList):
1939	G2fil.G2Print ('** ERROR - nothing to refine! **')
1940	return False
1941	while True:
1942	begin = time.time()
1943	values = np.array(Dict2Values(parmDict, varyList))
1944	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1945	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1946	ncyc = int(result[2]['nfev']/2)
1947	runtime = time.time()-begin
1948	chisq = np.sum(result[2]['fvec']**2)
1949	Values2Dict(parmDict, varyList, result[0])
1950	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
1951	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1952	G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1953	G2fil.G2Print ('chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF))
1954	try:
1955	sig = np.sqrt(np.diag(result[1])*GOF)
1956	if np.any(np.isnan(sig)):
1957	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
1958	break #refinement succeeded - finish up!
1959	except ValueError: #result[1] is None on singular matrix
1960	G2fil.G2Print ('** Refinement failed - singular matrix **')
1961
1962	sigDict = dict(zip(varyList,sig))
1963	GetInstParms(parmDict,Inst,varyList)
1964	InstPrint(Inst,sigDict)
1965	return True
1966
1967	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None):
1968	'''Called to perform a peak fit, refining the selected items in the peak
1969	table as well as selected items in the background.
1970
1971	:param str FitPgm: type of fit to perform. At present this is ignored.
1972	:param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1973	four values followed by a refine flag where the values are: position, intensity,
1974	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1975	tree entry, dict item "peaks"
1976	:param list Background: describes the background. List with two items.
1977	Item 0 specifies a background model and coefficients. Item 1 is a dict.
1978	From the Histogram/Background tree entry.
1979	:param list Limits: min and max x-value to use
1980	:param dict Inst: Instrument parameters
1981	:param dict Inst2: more Instrument parameters
1982	:param numpy.array data: a 5xn array. data[0] is the x-values,
1983	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1984	calc, background and difference intensities, respectively.
1985	:param array fixback: fixed background array; same size as data[0-5]
1986	:param list prevVaryList: Used in sequential refinements to override the
1987	variable list. Defaults as an empty list.
1988	:param bool oneCycle: True if only one cycle of fitting should be performed
1989	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1990	and derivType = controls['deriv type']. If None default values are used.
1991	:param float wtFactor: weight multiplier; = 1.0 by default
1992	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1993	Defaults to None, which means no updates are done.
1994	'''
1995	def GetBackgroundParms(parmList,Background):
1996	iBak = 0
1997	while True:
1998	try:
1999	bakName = 'Back;'+str(iBak)
2000	Background[0][iBak+3] = parmList[bakName]
2001	iBak += 1
2002	except KeyError:
2003	break
2004	iDb = 0
2005	while True:
2006	names = ['DebyeA;','DebyeR;','DebyeU;']
2007	try:
2008	for i,name in enumerate(names):
2009	val = parmList[name+str(iDb)]
2010	Background[1]['debyeTerms'][iDb][2*i] = val
2011	iDb += 1
2012	except KeyError:
2013	break
2014	iDb = 0
2015	while True:
2016	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
2017	try:
2018	for i,name in enumerate(names):
2019	val = parmList[name+str(iDb)]
2020	Background[1]['peaksList'][iDb][2*i] = val
2021	iDb += 1
2022	except KeyError:
2023	break
2024	if 'BF mult' in parmList:
2025	Background[1]['background PWDR'][1] = parmList['BF mult']
2026
2027	def BackgroundPrint(Background,sigDict):
2028	print ('Background coefficients for '+Background[0][0]+' function')
2029	ptfmt = "%12.5f"
2030	ptstr = 'value: '
2031	sigstr = 'esd : '
2032	for i,back in enumerate(Background[0][3:]):
2033	ptstr += ptfmt % (back)
2034	if Background[0][1]:
2035	prm = 'Back;'+str(i)
2036	if prm in sigDict:
2037	sigstr += ptfmt % (sigDict[prm])
2038	else:
2039	sigstr += " "*12
2040	if len(ptstr) > 75:
2041	print (ptstr)
2042	if Background[0][1]: print (sigstr)
2043	ptstr = 'value: '
2044	sigstr = 'esd : '
2045	if len(ptstr) > 8:
2046	print (ptstr)
2047	if Background[0][1]: print (sigstr)
2048
2049	if Background[1]['nDebye']:
2050	parms = ['DebyeA;','DebyeR;','DebyeU;']
2051	print ('Debye diffuse scattering coefficients')
2052	ptfmt = "%12.5f"
2053	print (' term DebyeA esd DebyeR esd DebyeU esd')
2054	for term in range(Background[1]['nDebye']):
2055	line = ' term %d'%(term)
2056	for ip,name in enumerate(parms):
2057	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
2058	if name+str(term) in sigDict:
2059	line += ptfmt%(sigDict[name+str(term)])
2060	else:
2061	line += " "*12
2062	print (line)
2063	if Background[1]['nPeaks']:
2064	print ('Coefficients for Background Peaks')
2065	ptfmt = "%15.3f"
2066	for j,pl in enumerate(Background[1]['peaksList']):
2067	names = 'peak %3d:'%(j+1)
2068	ptstr = 'values :'
2069	sigstr = 'esds :'
2070	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
2071	val = pl[2*i]
2072	prm = lbl+";"+str(j)
2073	names += '%15s'%(prm)
2074	ptstr += ptfmt%(val)
2075	if prm in sigDict:
2076	sigstr += ptfmt%(sigDict[prm])
2077	else:
2078	sigstr += " "*15
2079	print (names)
2080	print (ptstr)
2081	print (sigstr)
2082	if 'BF mult' in sigDict:
2083	print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult'])))
2084
2085	def SetInstParms(Inst):
2086	dataType = Inst['Type'][0]
2087	insVary = []
2088	insNames = []
2089	insVals = []
2090	for parm in Inst:
2091	insNames.append(parm)
2092	insVals.append(Inst[parm][1])
2093	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2094	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
2095	insVary.append(parm)
2096	instDict = dict(zip(insNames,insVals))
2097	# instDict['X'] = max(instDict['X'],0.01)
2098	# instDict['Y'] = max(instDict['Y'],0.01)
2099	if 'SH/L' in instDict:
2100	instDict['SH/L'] = max(instDict['SH/L'],0.002)
2101	return dataType,instDict,insVary
2102
2103	def GetInstParms(parmDict,Inst,varyList,Peaks):
2104	for name in Inst:
2105	Inst[name][1] = parmDict[name]
2106	iPeak = 0
2107	while True:
2108	try:
2109	sigName = 'sig'+str(iPeak)
2110	pos = parmDict['pos'+str(iPeak)]
2111	if sigName not in varyList and peakInstPrmMode:
2112	if 'T' in Inst['Type'][0]:
2113	dsp = G2lat.Pos2dsp(Inst,pos)
2114	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
2115	else:
2116	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
2117	gamName = 'gam'+str(iPeak)
2118	if gamName not in varyList and peakInstPrmMode:
2119	if 'T' in Inst['Type'][0]:
2120	dsp = G2lat.Pos2dsp(Inst,pos)
2121	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
2122	else:
2123	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
2124	iPeak += 1
2125	except KeyError:
2126	break
2127
2128	def InstPrint(Inst,sigDict):
2129	print ('Instrument Parameters:')
2130	ptfmt = "%12.6f"
2131	ptlbls = 'names :'
2132	ptstr = 'values:'
2133	sigstr = 'esds :'
2134	for parm in Inst:
2135	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2136	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
2137	ptlbls += "%s" % (parm.center(12))
2138	ptstr += ptfmt % (Inst[parm][1])
2139	if parm in sigDict:
2140	sigstr += ptfmt % (sigDict[parm])
2141	else:
2142	sigstr += 12*' '
2143	print (ptlbls)
2144	print (ptstr)
2145	print (sigstr)
2146
2147	def SetPeaksParms(dataType,Peaks):
2148	peakNames = []
2149	peakVary = []
2150	peakVals = []
2151	if 'C' in dataType:
2152	names = ['pos','int','sig','gam']
2153	else: #'T' and 'B'
2154	names = ['pos','int','alp','bet','sig','gam']
2155	for i,peak in enumerate(Peaks):
2156	for j,name in enumerate(names):
2157	peakVals.append(peak[2*j])
2158	parName = name+str(i)
2159	peakNames.append(parName)
2160	if peak[2*j+1]:
2161	peakVary.append(parName)
2162	return dict(zip(peakNames,peakVals)),peakVary
2163
2164	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
2165	if 'C' in Inst['Type'][0]:
2166	names = ['pos','int','sig','gam']
2167	else: #'T' & 'B'
2168	names = ['pos','int','alp','bet','sig','gam']
2169	for i,peak in enumerate(Peaks):
2170	pos = parmDict['pos'+str(i)]
2171	if 'difC' in Inst:
2172	dsp = pos/Inst['difC'][1]
2173	for j in range(len(names)):
2174	parName = names[j]+str(i)
2175	if parName in varyList or not peakInstPrmMode:
2176	peak[2*j] = parmDict[parName]
2177	elif 'alp' in parName:
2178	if 'T' in Inst['Type'][0]:
2179	peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
2180	else: #'B'
2181	peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
2182	elif 'bet' in parName:
2183	if 'T' in Inst['Type'][0]:
2184	peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
2185	else: #'B'
2186	peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
2187	elif 'sig' in parName:
2188	if 'T' in Inst['Type'][0]:
2189	peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
2190	else: #'C' & 'B'
2191	peak[2*j] = G2mth.getCWsig(parmDict,pos)
2192	elif 'gam' in parName:
2193	if 'T' in Inst['Type'][0]:
2194	peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
2195	else: #'C' & 'B'
2196	peak[2*j] = G2mth.getCWgam(parmDict,pos)
2197
2198	def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
2199	print ('Peak coefficients:')
2200	if 'C' in dataType:
2201	names = ['pos','int','sig','gam']
2202	else: #'T' & 'B'
2203	names = ['pos','int','alp','bet','sig','gam']
2204	head = 13*' '
2205	for name in names:
2206	if name in ['alp','bet']:
2207	head += name.center(8)+'esd'.center(8)
2208	else:
2209	head += name.center(10)+'esd'.center(10)
2210	head += 'bins'.center(8)
2211	print (head)
2212	if 'C' in dataType:
2213	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
2214	elif 'T' in dataType:
2215	ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
2216	else: #'B'
2217	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'alp':"%8.2f",'bet':"%8.4f",'sig':"%10.3f",'gam':"%10.3f"}
2218	for i,peak in enumerate(Peaks):
2219	ptstr = ':'
2220	for j in range(len(names)):
2221	name = names[j]
2222	parName = name+str(i)
2223	ptstr += ptfmt[name] % (parmDict[parName])
2224	if parName in varyList:
2225	ptstr += ptfmt[name] % (sigDict[parName])
2226	else:
2227	if name in ['alp','bet']:
2228	ptstr += 8*' '
2229	else:
2230	ptstr += 10*' '
2231	ptstr += '%9.2f'%(ptsperFW[i])
2232	print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
2233
2234	def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg):
2235	parmdict.update(zip(varylist,values))
2236	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType)
2237
2238	def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg):
2239	parmdict.update(zip(varylist,values))
2240	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata)
2241	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
2242	if dlg:
2243	dlg.Raise()
2244	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
2245	if not GoOn:
2246	return -M #abort!!
2247	return M
2248
2249	# beginning of DoPeakFit
2250	if controls:
2251	Ftol = controls['min dM/M']
2252	else:
2253	Ftol = 0.0001
2254	if oneCycle:
2255	Ftol = 1.0
2256	x,y,w,yc,yb,yd = data #these are numpy arrays - remove masks!
2257	if fixback is None:
2258	fixback = np.zeros_like(y)
2259	yc *= 0. #set calcd ones to zero
2260	yb *= 0.
2261	yd *= 0.
2262	xBeg = np.searchsorted(x,Limits[0])
2263	xFin = np.searchsorted(x,Limits[1])+1
2264	bakType,bakDict,bakVary = SetBackgroundParms(Background)
2265	dataType,insDict,insVary = SetInstParms(Inst)
2266	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2267	parmDict = {}
2268	parmDict.update(bakDict)
2269	parmDict.update(insDict)
2270	parmDict.update(peakDict)
2271	parmDict['Pdabc'] = [] #dummy Pdabc
2272	parmDict.update(Inst2) #put in real one if there
2273	if prevVaryList:
2274	varyList = prevVaryList[:]
2275	else:
2276	varyList = bakVary+insVary+peakVary
2277	fullvaryList = varyList[:]
2278	if not peakInstPrmMode:
2279	for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1',
2280	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2281	if v in varyList:
2282	raise Exception('Instrumental profile terms cannot be varied '+
2283	'after setPeakInstPrmMode(False) is used')
2284	while True:
2285	begin = time.time()
2286	values = np.array(Dict2Values(parmDict, varyList))
2287	Rvals = {}
2288	badVary = []
2289	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2290	args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2291	ncyc = int(result[2]['nfev']/2)
2292	runtime = time.time()-begin
2293	chisq = np.sum(result[2]['fvec']**2)
2294	Values2Dict(parmDict, varyList, result[0])
2295	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactorw[xBeg:xFin]y[xBeg:xFin]*2))100. #to %
2296	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
2297	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2298	if ncyc:
2299	G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2300	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2301	sig = [0]*len(varyList)
2302	if len(varyList) == 0: break # if nothing was refined
2303	try:
2304	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2305	if np.any(np.isnan(sig)):
2306	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
2307	break #refinement succeeded - finish up!
2308	except ValueError: #result[1] is None on singular matrix
2309	G2fil.G2Print ('** Refinement failed - singular matrix **')
2310	Ipvt = result[2]['ipvt']
2311	for i,ipvt in enumerate(Ipvt):
2312	if not np.sum(result[2]['fjac'],axis=1)[i]:
2313	G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2314	badVary.append(varyList[ipvt-1])
2315	del(varyList[ipvt-1])
2316	break
2317	else: # nothing removed
2318	break
2319	if dlg: dlg.Destroy()
2320	sigDict = dict(zip(varyList,sig))
2321	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0]
2322	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType)
2323	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2324	GetBackgroundParms(parmDict,Background)
2325	if bakVary: BackgroundPrint(Background,sigDict)
2326	GetInstParms(parmDict,Inst,varyList,Peaks)
2327	if insVary: InstPrint(Inst,sigDict)
2328	GetPeaksParms(Inst,parmDict,Peaks,varyList)
2329	binsperFWHM = []
2330	for peak in Peaks:
2331	FWHM = getFWHM(peak[0],Inst)
2332	try:
2333	xpk = x.searchsorted(peak[0])
2334	cw = x[xpk]-x[xpk-1]
2335	binsperFWHM.append(FWHM/cw)
2336	except IndexError:
2337	binsperFWHM.append(0.)
2338	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2339	if len(binsperFWHM):
2340	if min(binsperFWHM) < 1.:
2341	G2fil.G2Print ('* Warning: calculated peak widths are too narrow to refine profile coefficients *')
2342	if 'T' in Inst['Type'][0]:
2343	G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2344	else:
2345	G2fil.G2Print (' Manually increase W in Instrument Parameters')
2346	elif min(binsperFWHM) < 4.:
2347	G2fil.G2Print ('* Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement *')
2348	G2fil.G2Print ('* recommended is 6-10; you have %.2f *'%(min(binsperFWHM)))
2349	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2350
2351	def calcIncident(Iparm,xdata):
2352	'needs a doc string'
2353
2354	def IfunAdv(Iparm,xdata):
2355	Itype = Iparm['Itype']
2356	Icoef = Iparm['Icoeff']
2357	DYI = np.ones((12,xdata.shape[0]))
2358	YI = np.ones_like(xdata)*Icoef[0]
2359
2360	x = xdata/1000. #expressions are in ms
2361	if Itype == 'Exponential':
2362	for i in [1,3,5,7,9]:
2363	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2364	YI += Icoef[i]*Eterm
2365	DYI[i] *= Eterm
2366	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2367	elif 'Maxwell'in Itype:
2368	Eterm = np.exp(-Icoef[2]/x**2)
2369	DYI[1] = Eterm/x**5
2370	DYI[2] = -Icoef[1]DYI[1]/x*2
2371	YI += (Icoef[1]Eterm/x*5)
2372	if 'Exponential' in Itype:
2373	for i in range(3,11,2):
2374	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2375	YI += Icoef[i]*Eterm
2376	DYI[i] *= Eterm
2377	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2378	else: #Chebyschev
2379	T = (2./x)-1.
2380	Ccof = np.ones((12,xdata.shape[0]))
2381	Ccof[1] = T
2382	for i in range(2,12):
2383	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
2384	for i in range(1,10):
2385	YI += Ccof[i]*Icoef[i+2]
2386	DYI[i+2] =Ccof[i]
2387	return YI,DYI
2388
2389	Iesd = np.array(Iparm['Iesd'])
2390	Icovar = Iparm['Icovar']
2391	YI,DYI = IfunAdv(Iparm,xdata)
2392	YI = np.where(YI>0,YI,1.)
2393	WYI = np.zeros_like(xdata)
2394	vcov = np.zeros((12,12))
2395	k = 0
2396	for i in range(12):
2397	for j in range(i,12):
2398	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
2399	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
2400	k += 1
2401	M = np.inner(vcov,DYI.T)
2402	WYI = np.sum(M*DYI,axis=0)
2403	WYI = np.where(WYI>0.,WYI,0.)
2404	return YI,WYI
2405
2406	#### RMCutilities ################################################################################
2407	def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2408	inst = PWDdata['Instrument Parameters'][0]
2409	Xsb = 0.
2410	Ysb = 0.
2411	if 'T' in inst['Type'][1]:
2412	difC = inst['difC'][1]
2413	if useSamBrd[0]:
2414	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2415	Xsb = 1.e-4*difC/Size[1][0]
2416	if useSamBrd[1]:
2417	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2418	Ysb = 1.e-6difCMustrain[1][0]
2419	prms = ['Bank',
2420	'difC','difA','Zero','2-theta',
2421	'alpha','beta-0','beta-1',
2422	'sig-0','sig-1','sig-2',
2423	'Z','X','Y']
2424	fname = Name+'.inst'
2425	fl = open(fname,'w')
2426	fl.write('1\n')
2427	fl.write('%d\n'%int(inst[prms[0]][1]))
2428	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]))
2429	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2430	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))
2431	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))
2432	fl.write('\n\n\n')
2433	fl.close()
2434	else:
2435	if useSamBrd[0]:
2436	wave = G2mth.getWave(inst)
2437	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2438	Xsb = 1.8wave/(np.piSize[1][0])
2439	if useSamBrd[1]:
2440	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2441	Ysb = 0.0180*Mustrain[1][0]/np.pi
2442	prms = ['Bank',
2443	'Lam','Zero','Polariz.',
2444	'U','V','W',
2445	'X','Y']
2446	fname = Name+'.inst'
2447	fl = open(fname,'w')
2448	fl.write('1\n')
2449	fl.write('%d\n'%int(inst[prms[0]][1]))
2450	fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2451	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2452	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))
2453	fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2454	fl.close()
2455	return fname
2456
2457	def MakeBack(PWDdata,Name):
2458	Back = PWDdata['Background'][0]
2459	inst = PWDdata['Instrument Parameters'][0]
2460	if 'chebyschev-1' != Back[0]:
2461	return None
2462	Nback = Back[2]
2463	BackVals = Back[3:]
2464	fname = Name+'.back'
2465	fl = open(fname,'w')
2466	fl.write('%10d\n'%Nback)
2467	for val in BackVals:
2468	if 'T' in inst['Type'][1]:
2469	fl.write('%12.6g\n'%(float(val)))
2470	else:
2471	fl.write('%12.6g\n'%val)
2472	fl.close()
2473	return fname
2474
2475	def findDup(Atoms):
2476	Dup = []
2477	Fracs = []
2478	for iat1,at1 in enumerate(Atoms):
2479	if any([at1[0] in dup for dup in Dup]):
2480	continue
2481	else:
2482	Dup.append([at1[0],])
2483	Fracs.append([at1[6],])
2484	for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2485	if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2486	Dup[-1] += [at2[0],]
2487	Fracs[-1]+= [at2[6],]
2488	return Dup,Fracs
2489
2490	def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
2491
2492	Meta = RMCPdict['metadata']
2493	Atseq = RMCPdict['atSeq']
2494	Supercell = RMCPdict['SuperCell']
2495	generalData = Phase['General']
2496	Dups,Fracs = findDup(Phase['Atoms'])
2497	Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2498	ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2499	Sample = PWDdata['Sample Parameters']
2500	Meta['temperature'] = Sample['Temperature']
2501	Meta['pressure'] = Sample['Pressure']
2502	Cell = generalData['Cell'][1:7]
2503	Trans = np.eye(3)*np.array(Supercell)
2504	newPhase = copy.deepcopy(Phase)
2505	newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2506	newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2507	GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2508	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.
2509	newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2510	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2511	Natm = np.count_nonzero(Natm-1)
2512	Atoms = newPhase['Atoms']
2513	reset = False
2514
2515	if ifSfracs:
2516	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2517	Natm = np.count_nonzero(Natm-1)
2518	Satoms = []
2519	for i in range(len(Atoms)//Natm):
2520	ind = i*Natm
2521	Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2522	Natoms = []
2523	for satoms in Satoms:
2524	for idup,dup in enumerate(Dups):
2525	ldup = len(dup)
2526	natm = len(satoms)
2527	i = 0
2528	while i < natm:
2529	if satoms[i][0] in dup:
2530	atoms = satoms[i:i+ldup]
2531	try:
2532	atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2533	Natoms.append(atom)
2534	except IndexError: #what about vacancies?
2535	if 'Va' not in Atseq:
2536	reset = True
2537	Atseq.append('Va')
2538	RMCPdict['aTypes']['Va'] = 0.0
2539	atom = atoms[0]
2540	atom[1] = 'Va'
2541	Natoms.append(atom)
2542	i += ldup
2543	else:
2544	i += 1
2545	else:
2546	Natoms = Atoms
2547
2548	NAtype = np.zeros(len(Atseq))
2549	for atom in Natoms:
2550	NAtype[Atseq.index(atom[1])] += 1
2551	NAstr = ['%6d'%i for i in NAtype]
2552	Cell = newPhase['General']['Cell'][1:7]
2553	if os.path.exists(Name+'.his6f'):
2554	os.remove(Name+'.his6f')
2555	if os.path.exists(Name+'.neigh'):
2556	os.remove(Name+'.neigh')
2557	fname = Name+'.rmc6f'
2558	fl = open(fname,'w')
2559	fl.write('(Version 6f format configuration file)\n')
2560	for item in Meta:
2561	fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2562	fl.write('Atom types present: %s\n'%' '.join(Atseq))
2563	fl.write('Number of each atom type: %s\n'%''.join(NAstr))
2564	fl.write('Number of atoms: %d\n'%len(Natoms))
2565	fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2566	fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2567	Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2568	A,B = G2lat.cell2AB(Cell,True)
2569	fl.write('Lattice vectors (Ang):\n')
2570	for i in [0,1,2]:
2571	fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2572	fl.write('Atoms (fractional coordinates):\n')
2573	nat = 0
2574	for atm in Atseq:
2575	for iat,atom in enumerate(Natoms):
2576	if atom[1] == atm:
2577	nat += 1
2578	atcode = Atcodes[iat].split(':')
2579	cell = [0,0,0]
2580	if '+' in atcode[1]:
2581	cell = eval(atcode[1].split('+')[1])
2582	fl.write('%6d%4s [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(
2583	nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2584	fl.close()
2585	return fname,reset
2586
2587	def MakeBragg(PWDdata,Name,Phase):
2588	generalData = Phase['General']
2589	Vol = generalData['Cell'][7]
2590	Data = PWDdata['Data']
2591	Inst = PWDdata['Instrument Parameters'][0]
2592	Bank = int(Inst['Bank'][1])
2593	Sample = PWDdata['Sample Parameters']
2594	Scale = Sample['Scale'][0]
2595	if 'X' in Inst['Type'][1]:
2596	Scale *= 2.
2597	Limits = PWDdata['Limits'][1]
2598	Ibeg = np.searchsorted(Data[0],Limits[0])
2599	Ifin = np.searchsorted(Data[0],Limits[1])+1
2600	fname = Name+'.bragg'
2601	fl = open(fname,'w')
2602	fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2603	if 'T' in Inst['Type'][0]:
2604	fl.write('%12s%12s\n'%(' TOF,ms',' I(obs)'))
2605	for i in range(Ibeg,Ifin-1):
2606	fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2607	else:
2608	fl.write('%12s%12s\n'%(' 2-theta, deg',' I(obs)'))
2609	for i in range(Ibeg,Ifin-1):
2610	fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))
2611	fl.close()
2612	return fname
2613
2614	def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2615	Meta = RMCPdict['metadata']
2616	Times = RMCPdict['runTimes']
2617	Atseq = RMCPdict['atSeq']
2618	Natoms = RMCPdict['NoAtoms']
2619	sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms]))
2620	Isotope = RMCPdict['Isotope']
2621	Isotopes = RMCPdict['Isotopes']
2622	Atypes = RMCPdict['aTypes']
2623	atPairs = RMCPdict['Pairs']
2624	Files = RMCPdict['files']
2625	BraggWt = RMCPdict['histogram'][1]
2626	inst = PWDdata['Instrument Parameters'][0]
2627	try:
2628	refList = PWDdata['Reflection Lists'][Name]['RefList']
2629	except KeyError:
2630	return 'Error - missing reflection list; you must do Refine first'
2631	dMin = refList[-1][4]
2632	gsasType = 'xray2'
2633	if 'T' in inst['Type'][1]:
2634	gsasType = 'gsas3'
2635	elif 'X' in inst['Type'][1]:
2636	XFF = G2elem.GetFFtable(Atseq)
2637	Xfl = open(Name+'.xray','w')
2638	for atm in Atseq:
2639	fa = XFF[atm]['fa']
2640	fb = XFF[atm]['fb']
2641	fc = XFF[atm]['fc']
2642	Xfl.write('%2s %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2643	atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2644	Xfl.close()
2645	lenA = len(Atseq)
2646	Pairs = []
2647	Ncoeff = []
2648	Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes]
2649	for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2650	Pairs += pair
2651	for pair in Pairs:
2652	pair = pair.replace(' ','')
2653	at1,at2 = pair.split('-')
2654	ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms
2655	ncoef = Isotopes[at2][Isotope[at2]]['SL'][0]Natoms[at2]/sumatms
2656	if at1 != at2:
2657	ncoef *= 2.
2658	Ncoeff += [ncoef,]
2659	pairMin = [atPairs[pair] for pair in Pairs if pair in atPairs]
2660	maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2661	fname = Name+'.dat'
2662	fl = open(fname,'w')
2663	fl.write(' %% Hand edit the following as needed\n')
2664	fl.write('TITLE :: '+Name+'\n')
2665	fl.write('MATERIAL :: '+Meta['material']+'\n')
2666	fl.write('PHASE :: '+Meta['phase']+'\n')
2667	fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2668	fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2669	if RMCPdict.get('useGPU',False):
2670	fl.write('GPU_ACCELERATOR :: 0\n')
2671	minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2672	minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2673	maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2674	fl.write('MINIMUM_DISTANCES :: %s Angstrom\n'%minHD)
2675	maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2676	fl.write('MAXIMUM_MOVES :: %s Angstrom\n'%maxMv)
2677	fl.write('R_SPACING :: 0.0200 Angstrom\n')
2678	fl.write('PRINT_PERIOD :: 100\n')
2679	fl.write('TIME_LIMIT :: %.2f MINUTES\n'%Times[0])
2680	fl.write('SAVE_PERIOD :: %.2f MINUTES\n'%Times[1])
2681	fl.write('\n')
2682	fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2683	fl.write('\n')
2684	fl.write('FLAGS ::\n')
2685	fl.write(' > NO_MOVEOUT\n')
2686	fl.write(' > NO_SAVE_CONFIGURATIONS\n')
2687	fl.write(' > NO_RESOLUTION_CONVOLUTION\n')
2688	fl.write('\n')
2689	fl.write('INPUT_CONFIGURATION_FORMAT :: rmc6f\n')
2690	fl.write('SAVE_CONFIGURATION_FORMAT :: rmc6f\n')
2691	fl.write('IGNORE_HISTORY_FILE ::\n')
2692	fl.write('\n')
2693	fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n')
2694	fl.write('DISTANCE_WINDOW ::\n')
2695	fl.write(' > MNDIST :: %s\n'%minD)
2696	fl.write(' > MXDIST :: %s\n'%maxD)
2697	if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2698	fl.write('\n')
2699	fl.write('POTENTIALS ::\n')
2700	fl.write(' > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2701	fl.write(' > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2702	if len(RMCPdict['Potentials']['Stretch']):
2703	fl.write(' > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2704	for bond in RMCPdict['Potentials']['Stretch']:
2705	fl.write(' > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))
2706	if len(RMCPdict['Potentials']['Angles']):
2707	fl.write(' > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2708	for angle in RMCPdict['Potentials']['Angles']:
2709	fl.write(' > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2710	(angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2711	if RMCPdict['useBVS']:
2712	fl.write('BVS ::\n')
2713	fl.write(' > ATOM :: '+' '.join(Atseq)+'\n')
2714	fl.write(' > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2715	oxid = []
2716	for val in RMCPdict['Oxid']:
2717	if len(val) == 3:
2718	oxid.append(val[0][1:])
2719	else:
2720	oxid.append(val[0][2:])
2721	fl.write(' > OXID :: %s\n'%' '.join(oxid))
2722	fl.write(' > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2723	fl.write(' > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2724	fl.write(' > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2725	fl.write(' > SAVE :: 100000\n')
2726	fl.write(' > UPDATE :: 100000\n')
2727	if len(RMCPdict['Swap']):
2728	fl.write('\n')
2729	fl.write('SWAP_MULTI ::\n')
2730	for swap in RMCPdict['Swap']:
2731	try:
2732	at1 = Atseq.index(swap[0])
2733	at2 = Atseq.index(swap[1])
2734	except ValueError:
2735	break
2736	fl.write(' > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2737
2738	if len(RMCPdict['FxCN']):
2739	fl.write('FIXED_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['FxCN']))
2740	for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2741	try:
2742	at1 = Atseq.index(fxcn[0])
2743	at2 = Atseq.index(fxcn[1])
2744	except ValueError:
2745	break
2746	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]))
2747	if len(RMCPdict['AveCN']):
2748	fl.write('AVERAGE_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['AveCN']))
2749	for iav,avcn in enumerate(RMCPdict['AveCN']):
2750	try:
2751	at1 = Atseq.index(avcn[0])
2752	at2 = Atseq.index(avcn[1])
2753	except ValueError:
2754	break
2755	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]))
2756	for File in Files:
2757	if Files[File][0] and Files[File][0] != 'Select':
2758	if 'Xray' in File and 'F(Q)' in File:
2759	fqdata = open(Files[File][0],'r')
2760	lines = int(fqdata.readline()[:-1])
2761	fl.write('\n')
2762	fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2763	fl.write(' > FILENAME :: %s\n'%Files[File][0])
2764	fl.write(' > DATA_TYPE :: %s\n'%Files[File][2])
2765	fl.write(' > FIT_TYPE :: %s\n'%Files[File][2])
2766	if 'Xray' not in File:
2767	fl.write(' > START_POINT :: 1\n')
2768	fl.write(' > END_POINT :: 3000\n')
2769	fl.write(' > WEIGHT :: %.4f\n'%Files[File][1])
2770	fl.write(' > CONSTANT_OFFSET 0.000\n')
2771	fl.write(' > NO_FITTED_OFFSET\n')
2772	if RMCPdict['FitScale']:
2773	fl.write(' > FITTED_SCALE\n')
2774	else:
2775	fl.write(' > NO_FITTED_SCALE\n')
2776	if Files[File][3] !='RMC':
2777	fl.write(' > %s\n'%Files[File][3])
2778	if 'reciprocal' in File:
2779	fl.write(' > CONVOLVE ::\n')
2780	if 'Xray' in File:
2781	fl.write(' > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2782	fl.write(' > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2783	fl.write(' > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2784	fl.write(' > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2785	fl.write('\n')
2786	fl.write('BRAGG ::\n')
2787	fl.write(' > BRAGG_SHAPE :: %s\n'%gsasType)
2788	fl.write(' > RECALCUATE\n')
2789	fl.write(' > DMIN :: %.2f\n'%(dMin-0.02))
2790	fl.write(' > WEIGHT :: %10.3f\n'%BraggWt)
2791	fl.write(' > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n')
2792	fl.write('\n')
2793	fl.write('END ::\n')
2794	fl.close()
2795	return fname
2796
2797	# def FindBonds(Phase,RMCPdict):
2798	# generalData = Phase['General']
2799	# cx,ct,cs,cia = generalData['AtomPtrs']
2800	# atomData = Phase['Atoms']
2801	# Res = 'RMC'
2802	# if 'macro' in generalData['Type']:
2803	# Res = atomData[0][ct-3]
2804	# AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2805	# Pairs = RMCPdict['Pairs'] #dict!
2806	# BondList = []
2807	# notNames = []
2808	# for FrstName in AtDict:
2809	# nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2810	# Atyp1 = AtDict[FrstName]
2811	# if 'Va' in Atyp1:
2812	# continue
2813	# for nbr in nbrs:
2814	# Atyp2 = AtDict[nbr[0]]
2815	# if 'Va' in Atyp2:
2816	# continue
2817	# try:
2818	# bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2819	# except KeyError:
2820	# bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2821	# if any(bndData):
2822	# if bndData[0] <= nbr[1] <= bndData[1]:
2823	# bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2824	# revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2825	# if bondStr not in BondList and revbondStr not in BondList:
2826	# BondList.append(bondStr)
2827	# notNames.append(FrstName)
2828	# return Res,BondList
2829
2830	# def FindAngles(Phase,RMCPdict):
2831	# generalData = Phase['General']
2832	# Cell = generalData['Cell'][1:7]
2833	# Amat = G2lat.cell2AB(Cell)[0]
2834	# cx,ct,cs,cia = generalData['AtomPtrs']
2835	# atomData = Phase['Atoms']
2836	# AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2837	# AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2838	# Angles = RMCPdict['Angles']
2839	# AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2840	# AngleList = []
2841	# for MidName in AtDict:
2842	# nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2843	# if len(nbrs) < 2: #need 2 neighbors to make an angle
2844	# continue
2845	# Atyp2 = AtDict[MidName]
2846	# for i,nbr1 in enumerate(nbrs):
2847	# Atyp1 = AtDict[nbr1[0]]
2848	# for j,nbr3 in enumerate(nbrs[i+1:]):
2849	# Atyp3 = AtDict[nbr3[0]]
2850	# IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2851	# try:
2852	# angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2853	# except KeyError:
2854	# try:
2855	# angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2856	# except KeyError:
2857	# continue
2858	# XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2859	# calAngle = G2mth.getRestAngle(XYZ,Amat)
2860	# if angData[0] <= calAngle <= angData[1]:
2861	# angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2862	# revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2863	# if angStr not in AngleList and revangStr not in AngleList:
2864	# AngleList.append(angStr)
2865	# return AngleList
2866
2867	# def GetSqConvolution(XY,d):
2868
2869	# n = XY.shape[1]
2870	# snew = np.zeros(n)
2871	# dq = np.zeros(n)
2872	# sold = XY[1]
2873	# q = XY[0]
2874	# dq[1:] = np.diff(q)
2875	# dq[0] = dq[1]
2876
2877	# for j in range(n):
2878	# for i in range(n):
2879	# b = abs(q[i]-q[j])
2880	# t = q[i]+q[j]
2881	# if j == i:
2882	# snew[j] += q[i]sold[i](d-np.sin(td)/t)dq[i]
2883	# else:
2884	# snew[j] += q[i]sold[i](np.sin(bd)/b-np.sin(td)/t)*dq[i]
2885	# snew[j] /= np.pi*q[j]
2886
2887	# snew[0] = snew[1]
2888	# return snew
2889
2890	# def GetMaxSphere(pdbName):
2891	# try:
2892	# pFil = open(pdbName,'r')
2893	# except FileNotFoundError:
2894	# return None
2895	# while True:
2896	# line = pFil.readline()
2897	# if 'Boundary' in line:
2898	# line = line.split()[3:]
2899	# G = np.array([float(item) for item in line])
2900	# G = np.reshape(G,(3,3))**2
2901	# G = nl.inv(G)
2902	# pFil.close()
2903	# break
2904	# dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
2905	# return min(dspaces)
2906
2907	def findfullrmc():
2908	'''Find where fullrmc is installed. Tries the following:
2909
2910	1. Returns the Config var 'fullrmc_exec', if defined. No check
2911	is done that the interpreter has fullrmc
2912	2. The current Python interpreter if fullrmc can be imported
2913	and fullrmc is version 5+
2914	3. The path is checked for a fullrmc image as named by Bachir
2915
2916	:returns: the full path to a python executable that is assumed to
2917	have fullrmc installed or None, if it was not found.
2918	'''
2919	is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
2920	if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
2921	GSASIIpath.GetConfigValue('fullrmc_exec')):
2922	return GSASIIpath.GetConfigValue('fullrmc_exec')
2923	try:
2924	import fullrmc
2925	if int(fullrmc.__version__.split('.')[0]) >= 5:
2926	return sys.executable
2927	except:
2928	pass
2929	pathlist = os.environ["PATH"].split(os.pathsep)
2930	for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
2931	os.path.split(sys.executable)[0]):
2932	if p not in pathlist: pathlist.insert(0,p)
2933	import glob
2934	for p in pathlist:
2935	if sys.platform == "darwin":
2936	lookfor = "fullrmcmacOSi386-64bit"
2937	elif sys.platform == "win32":
2938	lookfor = "fullrmc*.exe"
2939	else:
2940	lookfor = "fullrmc*"
2941	fl = glob.glob(lookfor)
2942	if len(fl) > 0:
2943	return os.path.abspath(sorted(fl)[0])
2944
2945	def findPDFfit():
2946	'''Find where PDFfit2 is installed. Does the following:
2947	:returns: the full path to a python executable that is assumed to
2948	have PDFfit2 installed or None, if it was not found.
2949	'''
2950	try:
2951	import diffpy.pdffit2
2952	return sys.executable
2953	except:
2954	return None
2955
2956	def MakePDFfitAtomsFile(Phase,RMCPdict):
2957	'''Make the PDFfit atoms file
2958	'''
2959	General = Phase['General']
2960	fName = General['Name']+'-PDFfit.stru'
2961	fatm = open(fName.replace(' ','_'),'w')
2962	fatm.write('title structure of '+General['Name']+'\n')
2963	fatm.write('format pdffit\n')
2964	fatm.write('scale 1.000000\n') #fixed
2965	if RMCPdict['shape'] == 'sphere':
2966	sharp = '%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0])
2967	else:
2968	sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
2969	fatm.write('sharp '+sharp)
2970	shape = ''
2971	if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.:
2972	shape = ' sphere, %10.6f\n'%RMCPdict['spdiameter'][0]
2973	elif RMCPdict['stepcut'] > 0.:
2974	shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut']
2975	if shape:
2976	fatm.write('shape '+shape)
2977	fatm.write('spcgr %s\n'%General['SGData']['SpGrp'].replace(' ',''))
2978	cell = General['Cell'][1:7]
2979	fatm.write('cell %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%(
2980	cell[0],cell[1],cell[2],cell[3],cell[4],cell[5]))
2981	fatm.write('dcell '+5*' 0.000000,'+' 0.000000\n')
2982	Atoms = Phase['Atoms']
2983	fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms)))
2984	fatm.write('atoms\n')
2985	cx,ct,cs,cia = General['AtomPtrs']
2986	for atom in Atoms:
2987	fatm.write('%4s%18.8f%18.8f%18.8f%13.4f\n'%(atom[ct][:2].ljust(2),atom[cx],atom[cx+1],atom[cx+2],atom[cx+3]))
2988	fatm.write(' '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.))
2989	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4]))
2990	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,))
2991	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7]))
2992	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.))
2993	fatm.close()
2994
2995	def MakePDFfitRunFile(Phase,RMCPdict):
2996	'''Make the PDFfit python run file
2997	'''
2998
2999	def GetCellConstr(SGData):
3000	if SGData['SGLaue'] in ['m3', 'm3m']:
3001	return [1,1,1,0,0,0]
3002	elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3003	return [1,1,2,0,0,0]
3004	elif SGData['SGLaue'] in ['3R','3mR']:
3005	return [1,1,1,2,2,2]
3006	elif SGData['SGLaue'] == 'mmm':
3007	return [1,2,3,0,0,0]
3008	elif SGData['SGLaue'] == '2/m':
3009	if SGData['SGUniq'] == 'a':
3010	return [1,2,3,4,0,0]
3011	elif SGData['SGUniq'] == 'b':
3012	return [1,2,3,0,4,0]
3013	elif SGData['SGUniq'] == 'c':
3014	return [1,2,3,0,0,4]
3015	else:
3016	return [1,2,3,4,5,6]
3017
3018	General = Phase['General']
3019	rundata = '''
3020	#!/usr/bin/env python
3021	# -- coding: utf-8 --
3022	from diffpy.pdffit2 import PdfFit
3023	pf = PdfFit()
3024	'''
3025	Nd = 0
3026	Np = 0
3027	for file in RMCPdict['files']:
3028	if 'Neutron' in file:
3029	Nd += 1
3030	dType = 'Ndata'
3031	else:
3032	Nd += 1
3033	dType = 'Xdata'
3034	rundata += "pf.read_data(%s, '%s', 30.0, %.4f)\n"%(dType[0],RMCPdict['files'][file][0],RMCPdict[dType]['qdamp'][0])
3035	rundata += 'pf.pdfrange(%d, %6.2f, %6.2f\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1])
3036	rundata += 'pf.setdata(%d)\n'%Nd
3037	for item in ['dscale','qdamp','qbroad']:
3038	rundata += "pf.setvar('%s', %.2f)\n"%(item,RMCPdict[dType][item][0])
3039	if RMCPdict[dType][item][1]:
3040	Np += 1
3041	rundata += 'pf.constrain("%s","@%d")\n'%(item,Np)
3042	rundata += "pf.read_struct(%s)\n"%(General['Name']+'-PDFfit.stru')
3043	for item in ['delta1','delta2','sratio']:
3044	if RMCPdict[item][1]:
3045	Np += 1
3046	rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np)
3047	if 'sphere' in RMCPdict['shape'][0] and RMCPdict['spdiameter'][1]:
3048	Np += 1
3049	rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np
3050
3051	if RMCPdict['cellref']:
3052	cellconst = GetCellConstr(RMCPdict['SGData'])
3053	for ic in range(6):
3054	if cellconst[ic]:
3055	rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic])
3056	#Atom constraints here -------------------------------------------------------
3057
3058
3059
3060	# Refine & Save results ---------------------------------------------------------------
3061	rundata += 'pf.refine()\n'
3062	fName = General['Name'].replace(' ','_')+'-PDFfit'
3063	Nd = 0
3064	for file in RMCPdict['files']:
3065	Nd += 1
3066	rundata += 'pf.save_pdf(%d, %s)\n'%(Nd,fName+file[0]+'.fgr')
3067
3068	rundata += 'pf.save_struct(1, %s)\n'%(fName+'.rstr')
3069	rundata += 'pf.save_res(%s)\n'%(fName+'.res')
3070
3071
3072	print(rundata)
3073
3074
3075
3076	# rfile = open(fName+.py','w')
3077	# rfile.writelines(rundata)
3078	# rfile.close()
3079
3080
3081
3082	def UpdatePDFfit(Phase,RMCPdict):
3083
3084	General = Phase['General']
3085	fName = General['Name']+'-PDFfit.rstr'
3086	rstr = open(fName.replace(' ','_'),'r')
3087	lines = rstr.readlines()
3088	rstr.close()
3089	header = [line[:-1].split(' ',1) for line in lines[:7]]
3090	resdict = dict(header)
3091	for item in ['scale','sharp','cell']:
3092	resdict[item] = [float(val) for val in resdict[item].split(',')]
3093	General['Cell'][1:7] = resdict['cell']
3094	for inam,name in enumerate(['delta2','delta1','sratio']):
3095	RMCPdict[name][0] = resdict['sharp'][inam]
3096	if 'shape' in resdict and 'sphere' in resdict['shape']:
3097	RMCPdict['spdiameter'][0] = resdict['shape'][-1]
3098	cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)
3099	Atoms = Phase['Atoms']
3100	atmBeg = 0
3101	for line in lines:
3102	atmBeg += 1
3103	if 'atoms' in line:
3104	break
3105	for atom in Atoms:
3106	atstr = lines[atmBeg][:-1].split()
3107	Uiistr = lines[atmBeg+2][:-1].split()
3108	Uijstr = lines[atmBeg+4][:-1].split()
3109	atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3110	atom[ci] = 'A'
3111	atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3112	atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3113	atmBeg += 6
3114
3115	def MakefullrmcRun(pName,Phase,RMCPdict):
3116	'''Creates a script to run fullrmc. Returns the name of the file that was
3117	created.
3118	'''
3119	BondList = {}
3120	for k in RMCPdict['Pairs']:
3121	if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3122	BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3123	AngleList = []
3124	for angle in RMCPdict['Angles']:
3125	if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3126	continue
3127	for i in (0,1,2):
3128	angle[i] = angle[i].strip()
3129	AngleList.append(angle)
3130	# rmin = RMCPdict['min Contact']
3131	cell = Phase['General']['Cell'][1:7]
3132	SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3133	cx,ct,cs,cia = Phase['General']['AtomPtrs']
3134	atomsList = []
3135	for atom in Phase['Atoms']:
3136	el = ''.join([i for i in atom[ct] if i.isalpha()])
3137	atomsList.append([el] + atom[cx:cx+4])
3138	projDir,projName = os.path.split(os.path.abspath(pName))
3139	scrname = pName+'-fullrmc.py'
3140	restart = '%s_restart.pdb'%pName
3141	Files = RMCPdict['files']
3142	rundata = ''
3143	rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3144	rundata += '# created in '+__file__+" v"+filversion.split()[1]
3145	rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3146	rundata += '''
3147	# fullrmc imports (all that are potentially useful)
3148	import os,glob
3149	import time
3150	import pickle
3151	import numpy as np
3152	from fullrmc.Core import Collection
3153	from fullrmc.Engine import Engine
3154	import fullrmc.Constraints.PairDistributionConstraints as fPDF
3155	from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3156	from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3157	from fullrmc.Constraints.BondConstraints import BondConstraint
3158	from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3159	from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3160	from fullrmc.Generators.Swaps import SwapPositionsGenerator
3161	# utility routines
3162	def writeHeader(ENGINE,statFP):
3163	'header for stats file'
3164	statFP.write('generated-steps, total-error, ')
3165	for c in ENGINE.constraints:
3166	statFP.write(c.constraintName)
3167	statFP.write(', ')
3168	statFP.write('\\n')
3169	statFP.flush()
3170
3171	def writeCurrentStatus(ENGINE,statFP,plotF):
3172	'line in stats file & current constraint plots'
3173	statFP.write(str(ENGINE.generated))
3174	statFP.write(', ')
3175	statFP.write(str(ENGINE.totalStandardError))
3176	statFP.write(', ')
3177	for c in ENGINE.constraints:
3178	statFP.write(str(c.standardError))
3179	statFP.write(', ')
3180	statFP.write('\\n')
3181	statFP.flush()
3182	mpl.use('agg')
3183	fp = open(plotF,'wb')
3184	for c in ENGINE.constraints:
3185	p = c.plot(show=False)
3186	p[0].canvas.draw()
3187	image = p[0].canvas.buffer_rgba()
3188	pickle.dump(c.constraintName,fp)
3189	pickle.dump(np.array(image),fp)
3190	fp.close()
3191
3192	def calcRmax(ENGINE):
3193	'from Bachir, works for non-othorhombic'
3194	a,b,c = ENGINE.basisVectors
3195	lens = []
3196	ts = np.linalg.norm(np.cross(a,b))/2
3197	lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3198	ts = np.linalg.norm(np.cross(b,c))/2
3199	lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3200	ts = np.linalg.norm(np.cross(a,c))/2
3201	lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3202	return min(lens)
3203	'''
3204	rundata += '''
3205	### When True, erases an existing engine to provide a fresh start
3206	FRESH_START = {:}
3207	dirName = "{:}"
3208	prefix = "{:}"
3209	project = prefix + "-fullrmc"
3210	time0 = time.time()
3211	'''.format(RMCPdict['ReStart'][0],projDir,projName)
3212
3213	rundata += '# setup structure\n'
3214	rundata += 'cell = ' + str(cell) + '\n'
3215	rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3216	rundata += 'atomList = ' + str(atomsList).replace('],','],\n ') + '\n'
3217	rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3218
3219	rundata += '\n# initialize engine\n'
3220	rundata += '''
3221	engineFileName = os.path.join(dirName, project + '.rmc')
3222	projectStats = os.path.join(dirName, project + '.stats')
3223	projectPlots = os.path.join(dirName, project + '.plots')
3224	pdbFile = os.path.join(dirName, project + '_restart.pdb')
3225	# check Engine exists if so (and not FRESH_START) load it
3226	# otherwise build it
3227	ENGINE = Engine(path=None)
3228	if not ENGINE.is_engine(engineFileName) or FRESH_START:
3229	## create structure
3230	ENGINE = Engine(path=engineFileName, freshStart=True)
3231	ENGINE.build_crystal_set_pdb(symOps = SymOpList,
3232	atoms = atomList,
3233	unitcellBC = cell,
3234	supercell = supercell)
3235	'''
3236	import atmdata
3237	# rundata += ' # conversion factors (may be needed)\n'
3238	# rundata += ' sumCiBi2 = 0.\n'
3239	# for elem in Phase['General']['AtomTypes']:
3240	# rundata += ' Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3241	# rundata += ' sumCiBi2 += (Ci{})*2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3242	rundata += ' rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3243	# settings that require a new Engine
3244	for File in Files:
3245	filDat = RMCPdict['files'][File]
3246	if not os.path.exists(filDat[0]): continue
3247	sfwt = 'neutronCohb'
3248	if 'Xray' in File:
3249	sfwt = 'atomicNumber'
3250	if 'G(r)' in File:
3251	rundata += ' GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3252	if filDat[3] == 0:
3253	#rundata += ''' # read and xform G(r) as defined in RMCProfile
3254	# see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3255	#rundata += ' GR[1] = 4 np.pi * GR[0] * rho0 / sumCiBi2\n'
3256	#rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3257	rundata += ' # G(r) as defined in RMCProfile\n'
3258	rundata += ' GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3259	elif filDat[3] == 1:
3260	rundata += ' # This is G(r) as defined in PDFFIT\n'
3261	rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3262	elif filDat[3] == 2:
3263	rundata += ' # This is g(r)\n'
3264	rundata += ' GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3265	else:
3266	raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3267	rundata += ' ENGINE.add_constraints([GofR])\n'
3268	rundata += ' GofR.set_limits((None, calcRmax(ENGINE)))\n'
3269	elif '(Q)' in File:
3270	rundata += ' SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3271	if filDat[3] == 0:
3272	rundata += ' # F(Q) as defined in RMCProfile\n'
3273	#rundata += ' SOQ[1] *= 1 / sumCiBi2\n'
3274	if filDat[4]:
3275	rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3276	rundata += ' SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3277	elif filDat[3] == 1:
3278	rundata += ' # S(Q) as defined in PDFFIT\n'
3279	rundata += ' SOQ[1] -= 1\n'
3280	if filDat[4]:
3281	rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3282	rundata += ' SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3283	else:
3284	raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3285	rundata += ' ENGINE.add_constraints([SofQ])\n'
3286	else:
3287	print('What is this?')
3288	rundata += ' ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
3289	if BondList:
3290	rundata += ''' B_CONSTRAINT = BondConstraint()
3291	ENGINE.add_constraints(B_CONSTRAINT)
3292	B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3293	'''
3294	for pair in BondList:
3295	e1,e2 = pair.split('-')
3296	d1,d2 = BondList[pair]
3297	if d1 == 0: continue
3298	if d2 == 0:
3299	print('Ignoring min distance without maximum')
3300	#rundata += ' ("element","{}","{}",{}),\n'.format(
3301	# e1.strip(),e2.strip(),d1)
3302	else:
3303	rundata += ' ("element","{}","{}",{},{}),\n'.format(
3304	e1.strip(),e2.strip(),d1,d2)
3305	rundata += ' ])\n'
3306	if AngleList:
3307	rundata += ''' A_CONSTRAINT = BondsAngleConstraint()
3308	ENGINE.add_constraints(A_CONSTRAINT)
3309	A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3310	'''
3311	for item in AngleList:
3312	rundata += (' '+
3313	'("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3314	rundata += ' ])\n'
3315	rundata += '''
3316	for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3317	ENGINE.save()
3318	else:
3319	ENGINE = ENGINE.load(path=engineFileName)
3320
3321	ENGINE.set_log_file(os.path.join(dirName,prefix))
3322	'''
3323	if RMCPdict['Swaps']:
3324	rundata += '\n#set up for site swaps\n'
3325	rundata += 'aN = ENGINE.allNames\n'
3326	rundata += 'SwapGen = {}\n'
3327	for swap in RMCPdict['Swaps']:
3328	rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3329	rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3330	rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3331	rundata += ' for swaps in SwapGen:\n'
3332	rundata += ' AB = swaps.split("-")\n'
3333	rundata += ' ENGINE.set_groups_as_atoms()\n'
3334	rundata += ' for g in ENGINE.groups:\n'
3335	rundata += ' if aN[g.indexes[0]]==AB[0]:\n'
3336	rundata += ' g.set_move_generator(SwapGen[swaps][0])\n'
3337	rundata += ' elif aN[g.indexes[0]]==AB[1]:\n'
3338	rundata += ' g.set_move_generator(SwapGen[swaps][1])\n'
3339	rundata += ' sProb = SwapGen[swaps][2]\n'
3340	rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3341	# rundata += 'wtDict = {}\n'
3342	# for File in Files:
3343	# filDat = RMCPdict['files'][File]
3344	# if not os.path.exists(filDat[0]): continue
3345	# if 'Xray' in File:
3346	# sfwt = 'atomicNumber'
3347	# else:
3348	# sfwt = 'neutronCohb'
3349	# if 'G(r)' in File:
3350	# typ = 'Pair'
3351	# elif '(Q)' in File:
3352	# typ = 'Struct'
3353	# rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3354	rundata += 'for c in ENGINE.constraints: # loop over predefined constraints\n'
3355	rundata += ' if type(c) is fPDF.PairDistributionConstraint:\n'
3356	# rundata += ' c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3357	rundata += ' c.set_limits((None,calcRmax(ENGINE)))\n'
3358	if RMCPdict['FitScale']:
3359	rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3360	# rundata += ' c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3361	if RMCPdict['FitScale']:
3362	rundata += ' elif type(c) is ReducedStructureFactorConstraint:\n'
3363	rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3364	# torsions difficult to implement, must be internal to cell & named with
3365	# fullrmc atom names
3366	# if len(RMCPdict['Torsions']): # Torsions currently commented out in GUI
3367	# rundata += 'for c in ENGINE.constraints: # look for Dihedral Angle Constraints\n'
3368	# rundata += ' if type(c) is DihedralAngleConstraint:\n'
3369	# rundata += ' c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3370	# rundata += ' c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3371	# for torsion in RMCPdict['Torsions']:
3372	# rundata += ' %s\n'%str(tuple(torsion))
3373	# rundata += ' ]})\n'
3374	rundata += '''
3375	if FRESH_START:
3376	# initialize engine with one step to get starting config energetics
3377	ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3378	statFP = open(projectStats,'w')
3379	writeHeader(ENGINE,statFP)
3380	writeCurrentStatus(ENGINE,statFP,projectPlots)
3381	else:
3382	statFP = open(projectStats,'a')
3383
3384	# setup runs for fullrmc
3385	'''
3386	rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3387	rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3388	rundata += ' ENGINE.set_groups_as_atoms()\n'
3389	rundata += ' expected = ENGINE.generated+steps\n'
3390
3391	rundata += ' ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3392	rundata += ' writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3393	rundata += ' if ENGINE.generated != expected: break # run was stopped\n'
3394	rundata += 'statFP.close()\n'
3395	rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3396	rfile = open(scrname,'w')
3397	rfile.writelines(rundata)
3398	rfile.close()
3399	return scrname
3400
3401	def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3402	bondDist = []
3403	Cell = general['Cell'][1:7]
3404	Supercell = RMCPdict['SuperCell']
3405	Trans = np.eye(3)*np.array(Supercell)
3406	Cell = G2lat.TransformCell(Cell,Trans)[:6]
3407	Amat,Bmat = G2lat.cell2AB(Cell)
3408	indices = (-1,0,1)
3409	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3410	for bonds in bondList:
3411	Oxyz = np.array(Atoms[bonds[0]][1:])
3412	Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])
3413	Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3414	Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3415	for dx in Dx.T:
3416	bondDist.append(np.min(dx))
3417	return np.array(bondDist)
3418
3419	def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3420	bondAngles = []
3421	Cell = general['Cell'][1:7]
3422	Supercell = RMCPdict['SuperCell']
3423	Trans = np.eye(3)*np.array(Supercell)
3424	Cell = G2lat.TransformCell(Cell,Trans)[:6]
3425	Amat,Bmat = G2lat.cell2AB(Cell)
3426	indices = (-1,0,1)
3427	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3428	for angle in angleList:
3429	Oxyz = np.array(Atoms[angle[0]][1:])
3430	TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])
3431	TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
3432	DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3433	DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3434	DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3435	DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3436	iDAx = np.argmin(DAx,axis=0)
3437	iDBx = np.argmin(DBx,axis=0)
3438	for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3439	DAv = DAxV[iA,i]/DAx[iA,i]
3440	DBv = DBxV[iB,i]/DBx[iB,i]
3441	bondAngles.append(npacosd(np.sum(DAv*DBv)))
3442	return np.array(bondAngles)
3443
3444
3445	#### Reflectometry calculations ################################################################################
3446	def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3447	G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3448
3449	class RandomDisplacementBounds(object):
3450	"""random displacement with bounds"""
3451	def __init__(self, xmin, xmax, stepsize=0.5):
3452	self.xmin = xmin
3453	self.xmax = xmax
3454	self.stepsize = stepsize
3455
3456	def __call__(self, x):
3457	"""take a random step but ensure the new position is within the bounds"""
3458	while True:
3459	# this could be done in a much more clever way, but it will work for example purposes
3460	steps = self.xmax-self.xmin
3461	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
3462	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3463	break
3464	return xnew
3465
3466	def GetModelParms():
3467	parmDict = {}
3468	varyList = []
3469	values = []
3470	bounds = []
3471	parmDict['dQ type'] = data['dQ type']
3472	parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig
3473	for parm in ['Scale','FltBack']:
3474	parmDict[parm] = data[parm][0]
3475	if data[parm][1]:
3476	varyList.append(parm)
3477	values.append(data[parm][0])
3478	bounds.append(Bounds[parm])
3479	parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3480	parmDict['nLayers'] = len(parmDict['Layer Seq'])
3481	for ilay,layer in enumerate(data['Layers']):
3482	name = layer['Name']
3483	cid = str(ilay)+';'
3484	parmDict[cid+'Name'] = name
3485	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3486	parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3487	if layer.get(parm,[0.,False])[1]:
3488	varyList.append(cid+parm)
3489	value = layer[parm][0]
3490	values.append(value)
3491	if value:
3492	bound = [value*Bfac,value/Bfac]
3493	else:
3494	bound = [0.,10.]
3495	bounds.append(bound)
3496	if name not in ['vacuum','unit scatter']:
3497	parmDict[cid+'rho'] = Substances[name]['Scatt density']
3498	parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3499	return parmDict,varyList,values,bounds
3500
3501	def SetModelParms():
3502	line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3503	if 'Scale' in varyList:
3504	data['Scale'][0] = parmDict['Scale']
3505	line += ' esd: %.4g'%(sigDict['Scale'])
3506	G2fil.G2Print (line)
3507	line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3508	if 'FltBack' in varyList:
3509	data['FltBack'][0] = parmDict['FltBack']
3510	line += ' esd: %15.3g'%(sigDict['FltBack'])
3511	G2fil.G2Print (line)
3512	for ilay,layer in enumerate(data['Layers']):
3513	name = layer['Name']
3514	G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3515	cid = str(ilay)+';'
3516	line = ' '
3517	line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3518	line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3519	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3520	if parm in layer:
3521	if parm == 'Rough':
3522	layer[parm][0] = abs(parmDict[cid+parm]) #make positive
3523	else:
3524	layer[parm][0] = parmDict[cid+parm]
3525	line += ' %s: %.3f'%(parm,layer[parm][0])
3526	if cid+parm in varyList:
3527	line += ' esd: %.3g'%(sigDict[cid+parm])
3528	G2fil.G2Print (line)
3529	G2fil.G2Print (line2)
3530
3531	def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3532	parmDict.update(zip(varyList,values))
3533	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3534	return M
3535
3536	def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3537	parmDict.update(zip(varyList,values))
3538	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3539	return np.sum(M**2)
3540
3541	def getREFD(Q,Qsig,parmDict):
3542	Ic = np.ones_like(Q)*parmDict['FltBack']
3543	Scale = parmDict['Scale']
3544	Nlayers = parmDict['nLayers']
3545	Res = parmDict['Res']
3546	depth = np.zeros(Nlayers)
3547	rho = np.zeros(Nlayers)
3548	irho = np.zeros(Nlayers)
3549	sigma = np.zeros(Nlayers)
3550	for ilay,lay in enumerate(parmDict['Layer Seq']):
3551	cid = str(lay)+';'
3552	depth[ilay] = parmDict[cid+'Thick']
3553	sigma[ilay] = parmDict[cid+'Rough']
3554	if parmDict[cid+'Name'] == u'unit scatter':
3555	rho[ilay] = parmDict[cid+'DenMul']
3556	irho[ilay] = parmDict[cid+'iDenMul']
3557	elif 'vacuum' != parmDict[cid+'Name']:
3558	rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3559	irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3560	if cid+'Mag SLD' in parmDict:
3561	rho[ilay] += parmDict[cid+'Mag SLD']
3562	if parmDict['dQ type'] == 'None':
3563	AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
3564	elif 'const' in parmDict['dQ type']:
3565	AB = SmearAbeles(0.5Q,QRes,depth,rho,irho,sigma[1:])
3566	else: #dQ/Q in data
3567	AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3568	Ic += AB*Scale
3569	return Ic
3570
3571	def estimateT0(takestep):
3572	Mmax = -1.e-10
3573	Mmin = 1.e10
3574	for i in range(100):
3575	x0 = takestep(values)
3576	M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3577	Mmin = min(M,Mmin)
3578	MMax = max(M,Mmax)
3579	return 1.5*(MMax-Mmin)
3580
3581	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3582	if data.get('2% weight'):
3583	wt = 1./(0.02Io)*2
3584	Qmin = Limits[1][0]
3585	Qmax = Limits[1][1]
3586	wtFactor = ProfDict['wtFactor']
3587	Bfac = data['Toler']
3588	Ibeg = np.searchsorted(Q,Qmin)
3589	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
3590	Ic[:] = 0
3591	Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3592	'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3593	parmDict,varyList,values,bounds = GetModelParms()
3594	Msg = 'Failed to converge'
3595	if varyList:
3596	if data['Minimizer'] == 'LMLS':
3597	result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3598	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3599	parmDict.update(zip(varyList,result[0]))
3600	chisq = np.sum(result[2]['fvec']**2)
3601	ncalc = result[2]['nfev']
3602	covM = result[1]
3603	newVals = result[0]
3604	elif data['Minimizer'] == 'Basin Hopping':
3605	xyrng = np.array(bounds).T
3606	take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3607	T0 = estimateT0(take_step)
3608	G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3609	result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3610	interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3611	'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3612	chisq = result.fun
3613	ncalc = result.nfev
3614	newVals = result.x
3615	covM = []
3616	elif data['Minimizer'] == 'MC/SA Anneal':
3617	xyrng = np.array(bounds).T
3618	result = G2mth.anneal(sumREFD, values,
3619	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3620	schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3621	boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3622	ranRange=0.20,autoRan=False,dlg=None)
3623	newVals = result[0]
3624	parmDict.update(zip(varyList,newVals))
3625	chisq = result[1]
3626	ncalc = result[3]
3627	covM = []
3628	G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3629	elif data['Minimizer'] == 'L-BFGS-B':
3630	result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol,
3631	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3632	parmDict.update(zip(varyList,result['x']))
3633	chisq = result.fun
3634	ncalc = result.nfev
3635	newVals = result.x
3636	covM = []
3637	else: #nothing varied
3638	M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3639	chisq = np.sum(M**2)
3640	ncalc = 0
3641	covM = []
3642	sig = []
3643	sigDict = {}
3644	result = []
3645	Rvals = {}
3646	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]Io[Ibeg:Ifin]2))100. #to %
3647	Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2
3648	Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3649	Ib[Ibeg:Ifin] = parmDict['FltBack']
3650	try:
3651	if not len(varyList):
3652	Msg += ' - nothing refined'
3653	raise ValueError
3654	Nans = np.isnan(newVals)
3655	if np.any(Nans):
3656	name = varyList[Nans.nonzero(True)[0]]
3657	Msg += ' Nan result for '+name+'!'
3658	raise ValueError
3659	Negs = np.less_equal(newVals,0.)
3660	if np.any(Negs):
3661	indx = Negs.nonzero()
3662	name = varyList[indx[0][0]]
3663	if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3664	Msg += ' negative coefficient for '+name+'!'
3665	raise ValueError
3666	if len(covM):
3667	sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3668	covMatrix = covM*Rvals['GOF']
3669	else:
3670	sig = np.zeros(len(varyList))
3671	covMatrix = []
3672	sigDict = dict(zip(varyList,sig))
3673	G2fil.G2Print (' Results of reflectometry data modelling fit:')
3674	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3675	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3676	SetModelParms()
3677	return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3678	except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans
3679	G2fil.G2Print (Msg)
3680	return False,0,0,0,0,0,0,Msg
3681
3682	def makeSLDprofile(data,Substances):
3683
3684	sq2 = np.sqrt(2.)
3685	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3686	Nlayers = len(laySeq)
3687	laySeq = np.array(laySeq,dtype=int)
3688	interfaces = np.zeros(Nlayers)
3689	rho = np.zeros(Nlayers)
3690	sigma = np.zeros(Nlayers)
3691	name = data['Layers'][0]['Name']
3692	thick = 0.
3693	for ilay,lay in enumerate(laySeq):
3694	layer = data['Layers'][lay]
3695	name = layer['Name']
3696	if 'Thick' in layer:
3697	thick += layer['Thick'][0]
3698	interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3699	if 'Rough' in layer:
3700	sigma[ilay] = max(0.001,layer['Rough'][0])
3701	if name != 'vacuum':
3702	if name == 'unit scatter':
3703	rho[ilay] = np.sqrt(layer['DenMul'][0]2+layer['iDenMul'][0]2)
3704	else:
3705	rrho = Substances[name]['Scatt density']
3706	irho = Substances[name]['XImag density']
3707	rho[ilay] = np.sqrt(rrho2+irho2)*layer['DenMul'][0]
3708	if 'Mag SLD' in layer:
3709	rho[ilay] += layer['Mag SLD'][0]
3710	name = data['Layers'][-1]['Name']
3711	x = np.linspace(-0.15thick,1.15thick,1000,endpoint=True)
3712	xr = np.flipud(x)
3713	interfaces[-1] = x[-1]
3714	y = np.ones_like(x)*rho[0]
3715	iBeg = 0
3716	for ilayer in range(Nlayers-1):
3717	delt = rho[ilayer+1]-rho[ilayer]
3718	iPos = np.searchsorted(x,interfaces[ilayer])
3719	y[iBeg:] += (delt/2.)sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2sigma[ilayer+1]))
3720	iBeg = iPos
3721	return x,xr,y
3722
3723	def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3724
3725	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3726	Qmin = Limits[1][0]
3727	Qmax = Limits[1][1]
3728	iBeg = np.searchsorted(Q,Qmin)
3729	iFin = np.searchsorted(Q,Qmax)+1 #include last point
3730	Ib[:] = data['FltBack'][0]
3731	Ic[:] = 0
3732	Scale = data['Scale'][0]
3733	if data['Layer Seq'] == []:
3734	return
3735	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3736	Nlayers = len(laySeq)
3737	depth = np.zeros(Nlayers)
3738	rho = np.zeros(Nlayers)
3739	irho = np.zeros(Nlayers)
3740	sigma = np.zeros(Nlayers)
3741	for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3742	layer = data['Layers'][lay]
3743	name = layer['Name']
3744	if 'Thick' in layer: #skips first & last layers
3745	depth[ilay] = layer['Thick'][0]
3746	if 'Rough' in layer: #skips first layer
3747	sigma[ilay] = layer['Rough'][0]
3748	if 'unit scatter' == name:
3749	rho[ilay] = layer['DenMul'][0]
3750	irho[ilay] = layer['iDenMul'][0]
3751	else:
3752	rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3753	irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3754	if 'Mag SLD' in layer:
3755	rho[ilay] += layer['Mag SLD'][0]
3756	if data['dQ type'] == 'None':
3757	AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
3758	elif 'const' in data['dQ type']:
3759	res = data['Resolution'][0]/(100.*sateln2)
3760	AB = SmearAbeles(0.5Q[iBeg:iFin],resQ[iBeg:iFin],depth,rho,irho,sigma[1:])
3761	else: #dQ/Q in data
3762	AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3763	Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3764
3765	def abeles(kz, depth, rho, irho=0, sigma=0):
3766	"""
3767	Optical matrix form of the reflectivity calculation.
3768	O.S. Heavens, Optical Properties of Thin Solid Films
3769
3770	Reflectometry as a function of kz for a set of slabs.
3771
3772	:param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
3773	This is :math:`\\tfrac12 Q_z`.
3774	:param depth: float[m] (Ang).
3775	thickness of each layer. The thickness of the incident medium
3776	and substrate are ignored.
3777	:param rho: float[n,k] (1e-6/Ang^2)
3778	Real scattering length density for each layer for each kz
3779	:param irho: float[n,k] (1e-6/Ang^2)
3780	Imaginary scattering length density for each layer for each kz
3781	Note: absorption cross section mu = 2 irho/lambda for neutrons
3782	:param sigma: float[m-1] (Ang)
3783	interfacial roughness. This is the roughness between a layer
3784	and the previous layer. The sigma array should have m-1 entries.
3785
3786	Slabs are ordered with the surface SLD at index 0 and substrate at
3787	index -1, or reversed if kz < 0.
3788	"""
3789	def calc(kz, depth, rho, irho, sigma):
3790	if len(kz) == 0: return kz
3791
3792	# Complex index of refraction is relative to the incident medium.
3793	# We can get the same effect using kz_rel^2 = kz^2 + 4pirho_o
3794	# in place of kz^2, and ignoring rho_o
3795	kz_sq = kz*2 + 4e-6np.pi*rho[:,0]
3796	k = kz
3797
3798	# According to Heavens, the initial matrix should be [ 1 F; F 1],
3799	# which we do by setting B=I and M0 to [1 F; F 1]. An extra matrix
3800	# multiply versus some coding convenience.
3801	B11 = 1
3802	B22 = 1
3803	B21 = 0
3804	B12 = 0
3805	for i in range(0, len(depth)-1):
3806	k_next = np.sqrt(kz_sq - 4e-6np.pi(rho[:,i+1] + 1j*irho[:,i+1]))
3807	F = (k - k_next) / (k + k_next)
3808	F = np.exp(-2kk_nextsigma[i]**2)
3809	#print "==== layer",i
3810	#print "kz:", kz
3811	#print "k:", k
3812	#print "k_next:",k_next
3813	#print "F:",F
3814	#print "rho:",rho[:,i+1]
3815	#print "irho:",irho[:,i+1]
3816	#print "d:",depth[i],"sigma:",sigma[i]
3817	M11 = np.exp(1jkdepth[i]) if i>0 else 1
3818	M22 = np.exp(-1jkdepth[i]) if i>0 else 1
3819	M21 = F*M11
3820	M12 = F*M22
3821	C1 = B11M11 + B21M12
3822	C2 = B11M21 + B21M22
3823	B11 = C1
3824	B21 = C2
3825	C1 = B12M11 + B22M12
3826	C2 = B12M21 + B22M22
3827	B12 = C1
3828	B22 = C2
3829	k = k_next
3830
3831	r = B12/B11
3832	return np.absolute(r)**2
3833
3834	if np.isscalar(kz): kz = np.asarray([kz], 'd')
3835
3836	m = len(depth)
3837
3838	# Make everything into arrays
3839	depth = np.asarray(depth,'d')
3840	rho = np.asarray(rho,'d')
3841	irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
3842	sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
3843
3844	# Repeat rho,irho columns as needed
3845	if len(rho.shape) == 1:
3846	rho = rho[None,:]
3847	irho = irho[None,:]
3848
3849	return calc(kz, depth, rho, irho, sigma)
3850
3851	def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
3852	y = abeles(kz, depth, rho, irho, sigma)
3853	s = dq/2.
3854	y += 0.1354(abeles(kz+2s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
3855	y += 0.24935(abeles(kz-5s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma))
3856	y += 0.4111(abeles(kz-4s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma))
3857	y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
3858	y += 0.80074(abeles(kz-2s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
3859	y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
3860	y *= 0.137023
3861	return y
3862
3863	def makeRefdFFT(Limits,Profile):
3864	G2fil.G2Print ('make fft')
3865	Q,Io = Profile[:2]
3866	Qmin = Limits[1][0]
3867	Qmax = Limits[1][1]
3868	iBeg = np.searchsorted(Q,Qmin)
3869	iFin = np.searchsorted(Q,Qmax)+1 #include last point
3870	Qf = np.linspace(0.,Q[iFin-1],5000)
3871	QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
3872	If = QI(Qf)Qf*4
3873	R = np.linspace(0.,5000.,5000)
3874	F = fft.rfft(If)
3875	return R,F
3876
3877
3878	#### Stacking fault simulation codes ################################################################################
3879	def GetStackParms(Layers):
3880
3881	Parms = []
3882	#cell parms
3883	if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
3884	Parms.append('cellA')
3885	Parms.append('cellC')
3886	else:
3887	Parms.append('cellA')
3888	Parms.append('cellB')
3889	Parms.append('cellC')
3890	if Layers['Laue'] != 'mmm':
3891	Parms.append('cellG')
3892	#Transition parms
3893	for iY in range(len(Layers['Layers'])):
3894	for iX in range(len(Layers['Layers'])):
3895	Parms.append('TransP;%d;%d'%(iY,iX))
3896	Parms.append('TransX;%d;%d'%(iY,iX))
3897	Parms.append('TransY;%d;%d'%(iY,iX))
3898	Parms.append('TransZ;%d;%d'%(iY,iX))
3899	return Parms
3900
3901	def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
3902	'''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
3903
3904	:param dict Layers: dict with following items
3905
3906	::
3907
3908	{'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
3909	'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
3910	'Layers':[],'Stacking':[],'Transitions':[]}
3911
3912	:param str ctrls: controls string to be written on DIFFaX controls.dif file
3913	:param float scale: scale factor
3914	:param dict background: background parameters
3915	:param list limits: min/max 2-theta to be calculated
3916	:param dict inst: instrument parameters dictionary
3917	:param list profile: powder pattern data
3918
3919	Note that parameters all updated in place
3920	'''
3921	import atmdata
3922	path = sys.path
3923	for name in path:
3924	if 'bin' in name:
3925	DIFFaX = name+'/DIFFaX.exe'
3926	G2fil.G2Print (' Execute '+DIFFaX)
3927	break
3928	# make form factor file that DIFFaX wants - atom types are GSASII style
3929	sf = open('data.sfc','w')
3930	sf.write('GSASII special form factor file for DIFFaX\n\n')
3931	atTypes = list(Layers['AtInfo'].keys())
3932	if 'H' not in atTypes:
3933	atTypes.insert(0,'H')
3934	for atType in atTypes:
3935	if atType == 'H':
3936	blen = -.3741
3937	else:
3938	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3939	Adat = atmdata.XrayFF[atType]
3940	text = '%4s'%(atType.ljust(4))
3941	for i in range(4):
3942	text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3943	text += '%11.6f%11.6f'%(Adat['fc'],blen)
3944	text += '%3d\n'%(Adat['Z'])
3945	sf.write(text)
3946	sf.close()
3947	#make DIFFaX control.dif file - future use GUI to set some of these flags
3948	cf = open('control.dif','w')
3949	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
3950	x0 = profile[0]
3951	iBeg = np.searchsorted(x0,limits[0])
3952	iFin = np.searchsorted(x0,limits[1])+1
3953	if iFin-iBeg > 20000:
3954	iFin = iBeg+20000
3955	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3956	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3957	cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3958	else:
3959	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3960	inst = {'Type':['XSC','XSC',]}
3961	cf.close()
3962	#make DIFFaX data file
3963	df = open('GSASII-DIFFaX.dat','w')
3964	df.write('INSTRUMENTAL\n')
3965	if 'X' in inst['Type'][0]:
3966	df.write('X-RAY\n')
3967	elif 'N' in inst['Type'][0]:
3968	df.write('NEUTRON\n')
3969	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
3970	df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3971	U = ateln2*inst['U'][1]/10000.
3972	V = ateln2*inst['V'][1]/10000.
3973	W = ateln2*inst['W'][1]/10000.
3974	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
3975	HW = np.sqrt(np.mean(HWHM))
3976	# df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3977	if 'Mean' in Layers['selInst']:
3978	df.write('GAUSSIAN %.6f TRIM\n'%(HW)) #fast option - might not really matter
3979	elif 'Gaussian' in Layers['selInst']:
3980	df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W)) #slow - make a GUI option?
3981	else:
3982	df.write('None\n')
3983	else:
3984	df.write('0.10\nNone\n')
3985	df.write('STRUCTURAL\n')
3986	a,b,c = Layers['Cell'][1:4]
3987	gam = Layers['Cell'][6]
3988	df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3989	laue = Layers['Laue']
3990	if laue == '2/m(ab)':
3991	laue = '2/m(1)'
3992	elif laue == '2/m(c)':
3993	laue = '2/m(2)'
3994	if 'unknown' in Layers['Laue']:
3995	df.write('%s %.3f\n'%(laue,Layers['Toler']))
3996	else:
3997	df.write('%s\n'%(laue))
3998	df.write('%d\n'%(len(Layers['Layers'])))
3999	if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
4000	df.write('%.1f %.1f\n'%(Layers['Width'][0][0]10000.,Layers['Width'][0][0]10000.)) #mum to A
4001	layerNames = []
4002	for layer in Layers['Layers']:
4003	layerNames.append(layer['Name'])
4004	for il,layer in enumerate(Layers['Layers']):
4005	if layer['SameAs']:
4006	df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
4007	continue
4008	df.write('LAYER %d\n'%(il+1))
4009	if '-1' in layer['Symm']:
4010	df.write('CENTROSYMMETRIC\n')
4011	else:
4012	df.write('NONE\n')
4013	for ia,atom in enumerate(layer['Atoms']):
4014	[name,atype,x,y,z,frac,Uiso] = atom
4015	if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
4016	frac /= 2.
4017	df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
4018	df.write('STACKING\n')
4019	df.write('%s\n'%(Layers['Stacking'][0]))
4020	if 'recursive' in Layers['Stacking'][0]:
4021	df.write('%s\n'%Layers['Stacking'][1])
4022	else:
4023	if 'list' in Layers['Stacking'][1]:
4024	Slen = len(Layers['Stacking'][2])
4025	iB = 0
4026	iF = 0
4027	while True:
4028	iF += 68
4029	if iF >= Slen:
4030	break
4031	iF = min(iF,Slen)
4032	df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
4033	iB = iF
4034	else:
4035	df.write('%s\n'%Layers['Stacking'][1])
4036	df.write('TRANSITIONS\n')
4037	for iY in range(len(Layers['Layers'])):
4038	sumPx = 0.
4039	for iX in range(len(Layers['Layers'])):
4040	p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
4041	p = round(p,3)
4042	df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
4043	sumPx += p
4044	if sumPx != 1.0: #this has to be picky since DIFFaX is.
4045	G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
4046	df.close()
4047	os.remove('data.sfc')
4048	os.remove('control.dif')
4049	os.remove('GSASII-DIFFaX.dat')
4050	return
4051	df.close()
4052	time0 = time.time()
4053	try:
4054	subp.call(DIFFaX)
4055	except OSError:
4056	G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
4057	G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
4058	if os.path.exists('GSASII-DIFFaX.spc'):
4059	Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
4060	iFin = iBeg+Xpat.shape[1]
4061	bakType,backDict,backVary = SetBackgroundParms(background)
4062	backDict['Lam1'] = G2mth.getWave(inst)
4063	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
4064	profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
4065	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
4066	rv = st.poisson(profile[3][iBeg:iFin])
4067	profile[1][iBeg:iFin] = rv.rvs()
4068	Z = np.ones_like(profile[3][iBeg:iFin])
4069	Z[1::2] *= -1
4070	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4071	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
4072	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
4073	#cleanup files..
4074	os.remove('GSASII-DIFFaX.spc')
4075	elif os.path.exists('GSASII-DIFFaX.sadp'):
4076	Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
4077	Sadp = np.reshape(Sadp,(256,-1))
4078	Layers['Sadp']['Img'] = Sadp
4079	os.remove('GSASII-DIFFaX.sadp')
4080	os.remove('data.sfc')
4081	os.remove('control.dif')
4082	os.remove('GSASII-DIFFaX.dat')
4083
4084	def SetPWDRscan(inst,limits,profile):
4085
4086	wave = G2mth.getMeanWave(inst)
4087	x0 = profile[0]
4088	iBeg = np.searchsorted(x0,limits[0])
4089	iFin = np.searchsorted(x0,limits[1])
4090	if iFin-iBeg > 20000:
4091	iFin = iBeg+20000
4092	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4093	pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
4094	return iFin-iBeg
4095
4096	def SetStackingSF(Layers,debug):
4097	# Load scattering factors into DIFFaX arrays
4098	import atmdata
4099	atTypes = Layers['AtInfo'].keys()
4100	aTypes = []
4101	for atype in atTypes:
4102	aTypes.append('%4s'%(atype.ljust(4)))
4103	SFdat = []
4104	for atType in atTypes:
4105	Adat = atmdata.XrayFF[atType]
4106	SF = np.zeros(9)
4107	SF[:8:2] = Adat['fa']
4108	SF[1:8:2] = Adat['fb']
4109	SF[8] = Adat['fc']
4110	SFdat.append(SF)
4111	SFdat = np.array(SFdat)
4112	pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
4113
4114	def SetStackingClay(Layers,Type):
4115	# Controls
4116	rand.seed()
4117	ranSeed = rand.randint(1,2**16-1)
4118	try:
4119	laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
4120	'6/m','6/mmm'].index(Layers['Laue'])+1
4121	except ValueError: #for 'unknown'
4122	laueId = -1
4123	if 'SADP' in Type:
4124	planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
4125	lmax = int(Layers['Sadp']['Lmax'])
4126	else:
4127	planeId = 0
4128	lmax = 0
4129	# Sequences
4130	StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
4131	try:
4132	StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
4133	except ValueError:
4134	StkParm = -1
4135	if StkParm == 2: #list
4136	StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
4137	Nstk = len(StkSeq)
4138	else:
4139	Nstk = 1
4140	StkSeq = [0,]
4141	if StkParm == -1:
4142	StkParm = int(Layers['Stacking'][1])
4143	Wdth = Layers['Width'][0]
4144	mult = 1
4145	controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
4146	LaueSym = Layers['Laue'].ljust(12)
4147	pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
4148	return laueId,controls
4149
4150	def SetCellAtoms(Layers):
4151	Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
4152	# atoms in layers
4153	atTypes = list(Layers['AtInfo'].keys())
4154	AtomXOU = []
4155	AtomTp = []
4156	LayerSymm = []
4157	LayerNum = []
4158	layerNames = []
4159	Natm = 0
4160	Nuniq = 0
4161	for layer in Layers['Layers']:
4162	layerNames.append(layer['Name'])
4163	for il,layer in enumerate(Layers['Layers']):
4164	if layer['SameAs']:
4165	LayerNum.append(layerNames.index(layer['SameAs'])+1)
4166	continue
4167	else:
4168	LayerNum.append(il+1)
4169	Nuniq += 1
4170	if '-1' in layer['Symm']:
4171	LayerSymm.append(1)
4172	else:
4173	LayerSymm.append(0)
4174	for ia,atom in enumerate(layer['Atoms']):
4175	[name,atype,x,y,z,frac,Uiso] = atom
4176	Natm += 1
4177	AtomTp.append('%4s'%(atype.ljust(4)))
4178	Ta = atTypes.index(atype)+1
4179	AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
4180	AtomXOU = np.array(AtomXOU)
4181	Nlayers = len(layerNames)
4182	pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
4183	return Nlayers
4184
4185	def SetStackingTrans(Layers,Nlayers):
4186	# Transitions
4187	TransX = []
4188	TransP = []
4189	for Ytrans in Layers['Transitions']:
4190	TransP.append([trans[0] for trans in Ytrans]) #get just the numbers
4191	TransX.append([trans[1:4] for trans in Ytrans]) #get just the numbers
4192	TransP = np.array(TransP,dtype='float').T
4193	TransX = np.array(TransX,dtype='float')
4194	# GSASIIpath.IPyBreak()
4195	pyx.pygettrans(Nlayers,TransP,TransX)
4196
4197	def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
4198	# Scattering factors
4199	SetStackingSF(Layers,debug)
4200	# Controls & sequences
4201	laueId,controls = SetStackingClay(Layers,'PWDR')
4202	# cell & atoms
4203	Nlayers = SetCellAtoms(Layers)
4204	Volume = Layers['Cell'][7]
4205	# Transitions
4206	SetStackingTrans(Layers,Nlayers)
4207	# PWDR scan
4208	Nsteps = SetPWDRscan(inst,limits,profile)
4209	# result as Spec
4210	x0 = profile[0]
4211	profile[3] = np.zeros(len(profile[0]))
4212	profile[4] = np.zeros(len(profile[0]))
4213	profile[5] = np.zeros(len(profile[0]))
4214	iBeg = np.searchsorted(x0,limits[0])
4215	iFin = np.searchsorted(x0,limits[1])+1
4216	if iFin-iBeg > 20000:
4217	iFin = iBeg+20000
4218	Nspec = 20001
4219	spec = np.zeros(Nspec,dtype='double')
4220	time0 = time.time()
4221	pyx.pygetspc(controls,Nspec,spec)
4222	G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
4223	time0 = time.time()
4224	U = ateln2*inst['U'][1]/10000.
4225	V = ateln2*inst['V'][1]/10000.
4226	W = ateln2*inst['W'][1]/10000.
4227	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
4228	HW = np.sqrt(np.mean(HWHM))
4229	BrdSpec = np.zeros(Nsteps)
4230	if 'Mean' in Layers['selInst']:
4231	pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
4232	elif 'Gaussian' in Layers['selInst']:
4233	pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
4234	else:
4235	BrdSpec = spec[:Nsteps]
4236	BrdSpec /= Volume
4237	iFin = iBeg+Nsteps
4238	bakType,backDict,backVary = SetBackgroundParms(background)
4239	backDict['Lam1'] = G2mth.getWave(inst)
4240	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
4241	profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
4242	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
4243	try:
4244	rv = st.poisson(profile[3][iBeg:iFin])
4245	profile[1][iBeg:iFin] = rv.rvs()
4246	except ValueError:
4247	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
4248	Z = np.ones_like(profile[3][iBeg:iFin])
4249	Z[1::2] *= -1
4250	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4251	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
4252	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
4253	G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
4254
4255	def CalcStackingSADP(Layers,debug):
4256
4257	# Scattering factors
4258	SetStackingSF(Layers,debug)
4259	# Controls & sequences
4260	laueId,controls = SetStackingClay(Layers,'SADP')
4261	# cell & atoms
4262	Nlayers = SetCellAtoms(Layers)
4263	# Transitions
4264	SetStackingTrans(Layers,Nlayers)
4265	# result as Sadp
4266	Nspec = 20001
4267	spec = np.zeros(Nspec,dtype='double')
4268	time0 = time.time()
4269	hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
4270	Sapd = np.zeros((256,256))
4271	iB = 0
4272	for i in range(hkLim):
4273	iF = iB+Nblk
4274	p1 = 127+int(i*Incr)
4275	p2 = 128-int(i*Incr)
4276	if Nblk == 128:
4277	if i:
4278	Sapd[128:,p1] = spec[iB:iF]
4279	Sapd[:128,p1] = spec[iF:iB:-1]
4280	Sapd[128:,p2] = spec[iB:iF]
4281	Sapd[:128,p2] = spec[iF:iB:-1]
4282	else:
4283	if i:
4284	Sapd[:,p1] = spec[iB:iF]
4285	Sapd[:,p2] = spec[iB:iF]
4286	iB += Nblk
4287	Layers['Sadp']['Img'] = Sapd
4288	G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
4289
4290	#### Maximum Entropy Method - Dysnomia ###############################################################################
4291	def makePRFfile(data,MEMtype):
4292	''' makes Dysnomia .prf control file from Dysnomia GUI controls
4293
4294	:param dict data: GSAS-II phase data
4295	:param int MEMtype: 1 for neutron data with negative scattering lengths
4296	0 otherwise
4297	:returns str: name of Dysnomia control file
4298	'''
4299
4300	generalData = data['General']
4301	pName = generalData['Name'].replace(' ','_')
4302	DysData = data['Dysnomia']
4303	prfName = pName+'.prf'
4304	prf = open(prfName,'w')
4305	prf.write('$PREFERENCES\n')
4306	prf.write(pName+'.mem\n') #or .fos?
4307	prf.write(pName+'.out\n')
4308	prf.write(pName+'.pgrid\n')
4309	prf.write(pName+'.fba\n')
4310	prf.write(pName+'_eps.raw\n')
4311	prf.write('%d\n'%MEMtype)
4312	if DysData['DenStart'] == 'uniform':
4313	prf.write('0\n')
4314	else:
4315	prf.write('1\n')
4316	if DysData['Optimize'] == 'ZSPA':
4317	prf.write('0\n')
4318	else:
4319	prf.write('1\n')
4320	prf.write('1\n')
4321	if DysData['Lagrange'][0] == 'user':
4322	prf.write('0\n')
4323	else:
4324	prf.write('1\n')
4325	prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
4326	prf.write('%.3f\n'%DysData['Lagrange'][2])
4327	prf.write('%.2f\n'%DysData['E_factor'])
4328	prf.write('1\n')
4329	prf.write('0\n')
4330	prf.write('%d\n'%DysData['Ncyc'])
4331	prf.write('1\n')
4332	prf.write('1 0 0 0 0 0 0 0\n')
4333	if DysData['prior'] == 'uniform':
4334	prf.write('0\n')
4335	else:
4336	prf.write('1\n')