Context Navigation

source: trunk/GSASIIpwd.py @ 4397

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

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

Download in other formats: