Context Navigation

source: trunk/GSASIIpwd.py @ 4588

Last change on this file since 4588 was 4588, checked in by toby, 5 years ago
use G2VarObj for param limits; add more info to seq. ref. done dialog; show Frozen in show LS parameters
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 169.2 KB

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

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

Download in other formats: