Context Navigation

source: trunk/GSASIIpwd.py @ 4952

Last change on this file since 4952 was 4952, checked in by toby, 2 years ago
switch to native fullrmc input modes rather than convert files
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 178.2 KB

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