Context Navigation

source: trunk/GSASIIpwd.py @ 4266

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

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

Download in other formats: