Context Navigation

source: trunk/GSASIIpwd.py @ 5018

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