Context Navigation

source: trunk/GSASIIpwd.py @ 4588

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

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

Download in other formats: