Treatment of spectroscopic profile of source x-ray
Python source code “xray5.py” for treatment of (effective) spectroscopic profile of source x-ray is shown below.
Click (tap) [Python] icon on the right-top edge of the code area to copy the code to clipboard.
# -*- coding: utf-8 -*-
###############################################
# xray5.py
# to treat Spectroscopic intensity distribution
# modified by T. Ida, Sep. 26, 2023
###############################################
import numpy as np # import "numpy" as "np"
from scipy import interpolate # import "interpolate" from "scipi"
from scipy import special # import "special" from "scipi"
import sys # import "sys" module to use "sys.stderr()" & "sys.exit()"
import os.path # import "os.path" module
import configparser # import "configparser"
import common5 as cmn # import "dct_common.py" as "cmn"
################################################
# Definition of Horizontal Scale Transformation
################################################
def scale_hrz(degTwoTheta):
x = np.log(np.sin(np.pi * degTwoTheta / 360))
return x
##################################################
# Definition of Ordinate Scale Transform (Factor)
##################################################
def scale_ord(degTwoTheta):
f = 2 * np.tan(np.pi * degTwoTheta / 360)
return f
########################################
# Definition of Lorentzian Peak Profile
########################################
def f_Lor(x, x0, w):
x_w = (x - x0) / w
y = 2 / (np.pi * w * (1.0 + 4 * x_w**2))
return y
################################################
# Definition of Exponential Background Function
################################################
def f_Exp(x, x0, w):
ans = np.exp(-(x - x0) / w) / w
return ans
#################################################
# Deutsch et al. (2004) model
# E11 = 8047.837 eV, W11 = 2.285 eV, I11 = 0.579
# E12 = 8045.367 eV, W12 = 3.358 eV, I12 = 0.080
# E21 = 8027.993 eV, W21 = 2.666 eV, I21 = 0.236
# E22 = 8026.504 eV, W22 = 3.571 eV, I22 = 0.105
##################################################
#################################
# Parse Spectroscopic Parameters
#################################
def DefFunc(config, section):
nFunc = config.getint(section,'nFunc')
paramSet = []
nCuAlpha = config.getint(section,'nFuncCuAlpha')
nCuBeta = config.getint(section,'nFuncCuBeta')
nEdge = config.getint(section,'nFuncEdge')
for i in range(0, nFunc):
key_t = ('t{:02d}').format(i)
type = config.get(section,key_t)
key_x = ('x{:02d}').format(i)
key_w = ('w{:02d}').format(i)
key_I = ('I{:02d}').format(i)
key_c = ('c{:02d}').format(i)
# Common Parameters
posit = config.getfloat(section,key_x)
width = config.getfloat(section,key_w)
integ = config.getfloat(section,key_I)
params = {} # Python Dictionary Type
params["type"] = type
params["posit"] = posit
params["width"] = width
params["integ"] = integ
# In Case of Truncated Lorentzian or Gamma ...
if ((type == 'Lorentz') or (type == 'Gamma')):
cutof = float(config[section][key_c])
# In case of Cu Kbeta...
if ((nCuAlpha <= i) & (i < nCuAlpha + nCuBeta)):
# If rBeta_Alpha is indicated...
if (config[section]['rBeta_Alpha']):
fBeta_Alpha = config.getfloat(section,'rBeta_Alpha')
integ *= fBeta_Alpha
# If shiftBeta is indicated...
sBeta = config.getfloat(section,'shiftBeta')
posit += sBeta
# If smearBeta is indicated...
if (config[section]['smearBeta']):
wBeta = config.getfloat(section,'smearBeta')
width += wBeta
params["posit"] = posit
params["width"] = width
params["integ"] = integ
params["cutof"] = cutof
paramSet = paramSet + [params]
return paramSet
#####################################
# Create Spectroscopic Profile Model
#####################################
def FormFunc(paramSet, x):
nP = x.size
y = np.zeros(nP) # intitialize to zero
nFunc = len(paramSet)
dx = x[1] - x[0]
for params in paramSet:
type = params['type']
posit = params['posit']
width = params['width']
integ = params['integ']
# In case of truncated Lorentzian...
if (type == 'Lorentz'):
cutof = params['cutof']
nC = max(0, int((cutof - x[0]) / dx))
y[nC:nP] += integ * f_Lor(x[nC:nP], posit, width)
# tricky? use relation: x[n1:n2][0] == x[n1]
# In case of exponential edge...
elif (type == 'Exponential'):
if (width > 0): # normal case
dx = x[1] - x[0]
nC = max(0, int((posit - x[0]) / dx))
SB = np.zeros(nP) # initialized to zero
SB[nC:nP] = np.exp(-(x[nC:nP] - posit) / width) / width
y[nC:nP] += integ * SB[nC:nP]
else: # inversed case (width < 0)
dx = x[1] - x[0]
nC = max(0, int((posit - x[0]) / dx))
SB = np.zeros(nP) # initialized to zero
SB[0:nC] = np.exp(-(x[0:nC] - posit) / width) / (-width)
y[0:nC] += integ * SB[0:nC]
elif (type == 'HyperbolicSecant'):
dx = x[1] - x[0]
nC = max(0, int((posit - x[0]) / dx))
SB = np.zeros(nP)
SB = 1 / np.cosh((x - posit) / width) / (np.pi * width)
y += integ * SB
elif (type == 'Gamma'):
cutof = params['cutof']
sqrtAlpha = (posit - cutof) / width
alpha = sqrtAlpha**2
gamma = width / sqrtAlpha
dx = x[1] - x[0]
nC = max(0,int((cutof - x[0]) / dx))
SB = np.zeros(nP)
xcg = (x - cutof)/gamma
SB[nC:nP] = np.abs(xcg[nC:nP])**(alpha-1)*np.exp(-xcg[nC:nP])
SB[nC:nP] = SB[nC:nP] / ( special.gamma(alpha) * gamma )
y[nC:nP] += integ * SB[nC:nP]
return y
###################################################
# "xray.treat"
# TREATMENT OF SPECTROSCOPIC PROFILE, MAIN BODY
###################################################
def treat(source, DPATH):
# source: Source Data (Python List)
# [s_angle, s_intensity],
# [s_angle, s_intensity, s_error1], or
# [s_angle, s_intensity, s_error1, s_error2]
# DPATH: Path of Configuration File
#************************************
# Read Configuration File "xray.cfg"
# xray configuration
#************************************
CFG_FILE = os.path.join(DPATH, 'xray.cfg')
xray_config = configparser.ConfigParser()
if os.path.exists(CFG_FILE):
xray_config.read(CFG_FILE, encoding = 'utf-8')
# xray_config.readfp(codecs.open(CFG_FILE, "r", "utf8"))
else:
sys.stderr.write('"' + CFG_FILE + '" is not found.\n')
sys.exit(2)
#************************************
# Parse Parameters for Deconvolution
#************************************
paramsDeconv = DefFunc(xray_config, 'xray_deconv')
#**********************************
# Parse Parameters for Convolution
#**********************************
paramsConv = DefFunc(xray_config, 'xray_conv')
#***********************************
# Read Configuration File "dct.cfg"
# dct configuration
#***********************************
CFG_FILE = os.path.join(DPATH, 'dct.cfg')
dct_config = configparser.ConfigParser()
if os.path.exists(CFG_FILE):
dct_config.read(CFG_FILE, encoding = 'utf-8')
# dct_config.readfp(codecs.open(CFG_FILE, "r", "utf8"))
else:
sys.stderr.write('"'+CFG_FILE+'" is not found.\n')
sys.exit(2)
error_type = dct_config.get('Control', 'error_type')
# Margin for Fourier-based Treatment
marginXray = dct_config.getfloat('Control','marginXray')
#*****************************************************
# Assign Items in Python List "source" to numpy.array
#*****************************************************
s_angle = np.array(source[0]) # angle data
s_intensity = np.array(source[1]) # intensity data
if (error_type == '1'):
s_error1 = np.array(source[2])
elif (error_type == '2'):
s_error2 = np.array(source[2])
elif (error_type == '3'):
s_error1 = np.array(source[2])
s_error2 = np.array(source[3])
d_intensity = np.array(s_intensity) # reserve array for output
s_angle = np.array(source[0])
s_intensity = np.array(source[1])
d_intensity = np.array(s_intensity)
#*******************************************************
# Step 1: Scale Transform
# s_angle (2Theta in deg.) => source_x (on chi3-scale)
# s_intensity => source_y (on chi3-scale)
#*******************************************************
# Transform horizontal scale
source_chiX = np.array(s_angle)
source_chiX = scale_hrz(s_angle) # horizontal
# Correction on sclae transform
source_corr = scale_ord(s_angle)
# Correction with Geometrical Factors
source_corr *= cmn.g_corr(s_angle / 2)
# Correction of Ordinate Values
source_etaX = s_intensity * source_corr # ordinate
#******************************************************************
# Step 2: Preparation for Fourier treatment
# Determine Marginal (Padding) Range, Total Points & Step Interval
#******************************************************************
margin = marginXray
# Value read from the file'dct.cfg'
# margin = 2.0, for example, but user can change the value .
index = cmn.indices(source_chiX, margin)
nSource = index[0] # total number of valid data
i1 = index[1] # index of 1st anker point on interpolated abscissa
i2 = index[2] # 2nd index
i3 = index[3] # 3rd index
nDest = index[4] # total number of interpolated data
nDest_2 = int(nDest/2) # half of the total number
# Determin Step Interval, dx
chi0 = source_chiX[0]
chiV = source_chiX[nSource-1]
dchi = (chiV - chi0) / (i1 - 1)
# Prepare Arrays for Interpolated Data in Valid Range
nValid = i1
chiD = chi0 + nDest * dchi
chi = np.linspace(chi0, chiV, nValid, endpoint = False)
#*********************************************************
# Step 3:
# Cubic spline interpolation of Data to Transformed Scale
# (source_chi, source_eta) => (chi, eta)
#*********************************************************
# Cubic spline interpolation
# eta = np.array(chi)
f3 = interpolate.interp1d(source_chiX, source_etaX, kind="cubic")
eta = f3(chi)
#********************************************
# Step 4:
# Extrapolation for marginal (padding) range
# (chi, eta) => (chi, eta)
#********************************************
# Fill marginal range with values
chi = np.linspace(chi0, chiD, nDest, endpoint = False)
eta = np.resize(eta, nDest)
iLin = np.linspace(0, i3-i2, i3-i2)
etaL = source_etaX[nSource - 1]
etaR = source_etaX[0]
eta[i1:i2] = etaL
eta[i2:i3] = etaL + (etaR - etaL) * iLin[0:i3-i2] / (i3-i2)
eta[i3:nDest] = etaR
#***********************************************
# Step 5: Fourier Transform of Intensity Data
# (eta) => (ft_eta)
#***********************************************
# Fourier transform of intensity data
ft_eta = np.fft.fft(eta)
# Step interval of Fourier transform
dxi = 1.0 / (dchi * nDest)
# Abscissa of Fourier transform
xi = np.linspace(0, nDest * dxi, nDest, endpoint=False)\
# Use compliment expression for Fourier transform
xi[nDest_2:nDest] = xi[nDest_2:nDest] - nDest * dxi
# Fourier transform of function to be deconvolved
deconv = np.zeros(nDest)
# Fourier transform of function to be convolved
conv = np.zeros(nDest)
#**************************************************
# Step 6: Preparation of Approximate Model Profile
# Function to be Deconvolved / Convolved
# => (deconv), (conv)
#**************************************************
# Set start point to be zero・・
wchi = nDest * dchi
chi = np.linspace(0, wchi, nDest, endpoint = False)
# chi[nDest_2:nDest] -= wchi
# Function To Be Deconvolved (Realistic Profile Model)
deconv[0:nDest_2] = FormFunc(paramsDeconv, chi[0:nDest_2])
deconv[nDest_2:nDest] = FormFunc(paramsDeconv, chi[nDest_2:nDest]-wchi)
# Function To Be Convolved
# (Idealistic Symmetric Single-Peak Profile Model)
conv[0:nDest_2] = FormFunc(paramsConv, chi[0:nDest_2])
conv[nDest_2:nDest] = FormFunc(paramsConv, chi[nDest_2:nDest]-wchi)
#**********************************************************
# Step 7: Fourier Transform of Approximate Model Profile
# Function To Be Deconvolved / Convolved
# (deconv), (conv) => (ft_deconv), (ft_conv), ft_inst
#**********************************************************
# Fourier transform of function to be deconvolved
ft_deconv = np.fft.fft(deconv)
# Fourier transform of function to be convolved
ft_conv = np.fft.fft(conv)
# Fourier transform of instrumental function
ft_inst = np.array(ft_deconv) # (make a copy)
ft_inst = ft_deconv / ft_conv
#**********************************************************
# Step 8: Here is the Core of DCT,
# Deconvolution/Convolution Treatment on Fourier transform
# (ft_eta)/(ft_inst) => (ft_eta)
#**********************************************************
# Deconvolution and Convolution Treatment
ft_eta = ft_eta / ft_inst
#***********************************
# Step 9: Inverse Fourier Transform
# (ft_eta) => (zeta)
#***********************************
zeta = np.real(np.fft.ifft(ft_eta))
#**********************************************
# Step 10: Pick up Intensity from DCT data
# on "chi" scale by cubic spline interpolation
# (chi, zeta) => (source_x, dest_y)
#**********************************************
# Reset equidistant abscissa
chi = np.linspace(chi0, chiD, nDest, endpoint = False)
# Cubic spline interpolation
f3 = interpolate.interp1d(chi, zeta, kind="cubic")
dest_eta = f3(source_chiX)
#**********************************************
# Step 11: Rescale transform
# (source_x, dest_y) => (s_angle, d_intensity)
#**********************************************
d_intensity = dest_eta / source_corr
dest = source
dest[1] = d_intensity
#***********************************************************
# Step 1-12 Fourier Treatment of Error Data in "xray.treat"
#***********************************************************
dest = cmn.TreatErrors(error_type, source, source_chiX, index, source_corr, ft_inst)
return dest