Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: trunk/GSASIIpwd.py

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