Context Navigation

source: trunk/GSASIIpwd.py @ 4951

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