Context Navigation

source: trunk/GSASIIpwd.py @ 4527

Last change on this file since 4527 was 4527, checked in by toby, 3 years ago
limits bug in hist.EditSimulated?; add default for do_refinements(); fullrmc bond angles
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 169.4 KB

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

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

Download in other formats: