Context Navigation

source: trunk/GSASIIpwd.py @ 4432

Last change on this file since 4432 was 4432, checked in by vondreele, 5 years ago
change ifConstraint flag in Transform to False flag missing reflection list in RMCProfile setup use Limits to bound PWDR sum routine - allows summing of unequal scans
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 163.9 KB

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

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

Download in other formats: