Implementation of deconvolutional treatment about sample transparency aberration in XRD data collected for finite-thickness and finite-width specimen [Ida, in preparation].
Python code “trnspr3.py” 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 -*-
###########################################
# trnspr3.py
# to treat Sample-Transparency Aberration
# by a two-step deconvolutioal method
# originally coded by T. Ida, Sep.13, 2023
# modified by T. Ida, Sep. 19, 2023
# modified by T. Ida, April 5, 2024
###########################################
import numpy as np # import "numpy" modeule as "np"
from scipy import interpolate # import "interpolate" from "scipy"
from scipy import special # import "special" from "scipy"
import common5 as cmn # import "common5.py" as "cmn"
########################################################
# calc_kappas()
# Calculate Cumulants of Sample-Transparency Aberration
########################################################
def calc_kappas(degTwot, cfg):
# degTwot: 2Theta in [deg.] (NumPy array)
# cfg=[muInv,gonioR,specT,specW,degDS,alphaTrans,
# marginTrans,skipTrans,error_type]
muInv = cfg[0]
gonioR = cfg[1]
specT = cfg[2]
specW = cfg[3]
degDS = cfg[4]
t_R = specT/gonioR
muR = gonioR/muInv
# R: goniometer radius in [mm], (float)
# degPhiDS: divergence slit open angle in[deg.], (float)
# mu: linear attenuation coefficient in [1/mm], (float)
# t: thickness of sample in [mm], (float)
# W: width of sample in [mm], (float)
theta = degTwot * np.pi/360 # theta in [rad.] (NumPy array)
cos_theta = np.cos(theta) # cos(theta) (NumPy array)
sin_theta = np.sin(theta) # sin(theta) (NumPy array)
tan_theta = sin_theta/cos_theta # tan(theta) (NumPy array)
R_phiDS = gonioR * degDS * np.pi/180 # (float)
Omega = R_phiDS/sin_theta # (NumPy array)
tau = 2*specT/tan_theta # (NumPy array)
g = cos_theta*sin_theta/(muR) # (NumPy array)
u = 2*specT*cos_theta/gonioR # (NumPy array)
Omega0 = (specW+Omega)/2-tau
Omega1 = (specW-Omega)/2
Omega2 = specW-tau
Omega3 = (specW+Omega)/2
u1 = Omega1*sin_theta/gonioR
u2 = specW*sin_theta/gonioR
u3 = Omega3*sin_theta/gonioR
# Find boundaries
case_a = (Omega+2*tau <= specW) # (a)
case_b = (np.maximum(Omega,2*tau-Omega) <= specW)
case_b &= (specW < Omega+2*tau)# (b)
case_c = (tau <= specW) & (specW < Omega) # (c)
case_cp = (Omega <= specW) & (specW < 2*tau-Omega)
case_d = (specW <= np.minimum(Omega,tau))
# Coefficients for case (b)
fB0 = Omega3/Omega # (k)(u)
fB1 = tau/(u*Omega) # (k+1)(u)
fB2 = -Omega1/Omega # (k)(u1)
fB3 = -Omega1/(u1*Omega) # (k+1)(u1)
# Coefficients for case (c)
fC0 = 1 # (k)(u)
fC1 = tau/(specW*u) # (k+1)(u)
# Coefficients for case (cp)
fCP0 = Omega3/Omega # (k)(u3)
fCP1 = Omega3/(Omega*u3) # (k+1)(u3)
fCP2 = -Omega1/Omega # (k)(u1)
fCP3 = -Omega1/(Omega*u1) # (k+1)(u1)
# Coefficients for case (d)
fD0 = 1 # (k)(u2)
fD1 = 1/u2 # (k+1)(u2)
# Case (a) and (b)
exp_u_g = np.exp(-u/g) # exp(-u/g) (NumPy array)
s0i = 1 - exp_u_g # exp(-u/g) (NumPy array)
s1i = -(-u)*exp_u_g - g*s0i
s2i = -(-u)**2*exp_u_g - g*2*s1i
s3i = -(-u)**3*exp_u_g - g*3*s2i
s4i = -(-u)**4*exp_u_g - g*4*s3i
s5i = -(-u)**5*exp_u_g - g*5*s4i
# Case (c)
exp_u1_g = np.exp(-u1/g) # exp(-u1/g) (NumPy array)
s01 = 1 - exp_u1_g # exp(-u1/g) (NumPy array)
s11 = -(-u1)*exp_u1_g - g*s01
s21 = -(-u1)**2*exp_u1_g - g*2*s11
s31 = -(-u1)**3*exp_u1_g - g*3*s21
s41 = -(-u1)**4*exp_u1_g - g*4*s31
s51 = -(-u1)**5*exp_u1_g - g*5*s41
# Case (d)
exp_u2_g = np.exp(-u2/g) # exp(-u2/g) (NumPy array)
s02 = 1 - exp_u2_g # exp(-u2/g) (NumPy array)
s12 = -(-u2)*exp_u2_g - g*s02
s22 = -(-u2)**2*exp_u2_g - g*2*s12
s32 = -(-u2)**3*exp_u2_g - g*3*s22
s42 = -(-u2)**4*exp_u2_g - g*4*s32
s52 = -(-u2)**5*exp_u2_g - g*5*s42
# Case (cp)
exp_u3_g = np.exp(-u3/g) # exp(-u3/g) (NumPy array)
s03 = 1 - exp_u3_g # exp(-u3/g) (NumPy array)
s13 = -(-u3)*exp_u3_g - g*s02
s23 = -(-u3)**2*exp_u3_g - g*2*s12
s33 = -(-u3)**3*exp_u3_g - g*3*s22
s43 = -(-u3)**4*exp_u3_g - g*4*s32
s53 = -(-u3)**5*exp_u3_g - g*5*s42
# Coefficients for s_{k}(u)
f0 = np.where(case_a,1,np.where(case_b,fB0,np.where(case_c,fC0,0)))
# Coefficients for s_{k+1}(u)
f1 = np.where(case_b,fB1,np.where(case_c,fC1,0))
# Coefficients for s_{k}(u1)
f2 = np.where(case_b,fB2,np.where(case_cp,fCP2,0))
# Coefficients for s_{k+1}(u1)
f3 = np.where(case_b,fB3,np.where(case_cp,fCP3,0))
# Coefficients for s_{k}(u3)
f4 = np.where(case_cp,fCP0,0)
# Coefficients for s_{k+1}(u3)
f5 = np.where(case_cp,fCP1,0)
# Coefficients for s_{k}(u2)
f6 = np.where(case_d,fD0,0)
# Coefficients for s_{k+1}(u2)
f7 = np.where(case_d,fD1,0)
# Calculate Raw (Non-Central) Moments up to 4th order
#
s0 = f0*s0i+f1*s1i+f2*s01+f3*s11+f4*s03+f5*s13+f6*s02+f7*s12
s1 = f0*s1i+f1*s2i+f2*s11+f3*s21+f4*s13+f5*s22+f6*s12+f7*s22
s2 = f0*s2i+f1*s3i+f2*s21+f3*s31+f4*s23+f5*s32+f6*s22+f7*s32
s3 = f0*s3i+f1*s4i+f2*s31+f3*s41+f4*s33+f5*s42+f6*s32+f7*s42
s4 = f0*s4i+f1*s5i+f2*s41+f3*s51+f4*s43+f5*s52+f6*s42+f7*s52
#
m1 = s1 / s0
m2 = s2 / s0
m3 = s3 / s0
m4 = s4 / s0
# Calculate Cumulants up to 4th order
kappa1 = m1
kappa2 = m2 - m1**2
kappa3 = m3 - 3*m2*m1 + 2*m1**3
kappa4 = m4 - 4*m3*m1 - 3*m2**2 + 12*m2*m1**2 - 6*m1**4
# Convert to Reduced Cumulants in [deg.]
kappa1 = 180.0/np.pi*kappa1
kappa2 = 180.0/np.pi*np.sqrt(kappa2)
kappa3 = 180.0/np.pi*np.sign(kappa3)*np.abs(kappa3)**(1/3)
kappa4 = 180.0/np.pi*np.sign(kappa4)*np.abs(kappa4)**(1/4)
kappas = np.array([s0, kappa1,kappa2,kappa3,kappa4])
return kappas
# returns kappas: (Python List of NumPy array)
# [s0,kappa1,kappa2,kappa3,kappa4]
############################
# Component model functions
############################
# fTexp(), Truncated Exponential Function
def fTexp(x):
return np.where(x<0.0,np.exp(x),0.0)
# fFtTexp(), Fourier Transform of fTexp()
def fFtTexp(xi):
return 1 / (1 + 2*np.pi*1j*xi)
# fRect(), Rectangular Function
def fRect(x):
return np.where((-1.0<x)&(x<0.0),1.0,0.0)
# fFtRect(), Fourier Transform of fRect()
def fFtRect(xi):
farZero = np.abs(xi) > 1.0E-6
numer=1.0-np.exp(-2*np.pi*1j*xi) # numerator
denom=2*np.pi*1j*xi # denominator
ones=np.ones_like(numer)
np.seterr(invalid='ignore')
ans = np.divide(numer,denom,out=ones,where=farZero)
np.seterr(divide='warn')
return ans
##############################################
# "trnspr5.treat"
# TREATMENT OF SAMPLE TRANSPARENCY ABERRATION
##############################################
def treat(sign,source,kappas,cfg):
# sign: -1 or 1
# (sign==-1: exponential)(sign==1: rectangular)
# source: Source Data (Python List)
# [s_angle, s_intensity],
# [s_angle, s_intensity, s_error1], or
# [s_angle, s_intensity, s_error1, s_error2]
# cfg=[muInv,gonioR,specT,specW,degDS,alphaTrans,marginTrans,skipTrans,error_type]
debug = 0
if (debug == 1):
if (sign==-1):
print("*** Treatment of Sample Transparency Aberration (exponential) ***")
else:
print("*** Treatment of Sample Transparency Aberration (rectangular) ***")
#******************************
# Configuration from "dct.cfg"
#******************************
marginTrans = cfg[6]
modeTrans = 3
skipTrans = cfg[7]
error_type=cfg[8]
#***************************************************
# Step 0: Parse Python List "source" to numpy.array
#***************************************************
s_angle = np.array(source[0])
s_int = np.array(source[1])
if ((error_type == '1') or (error_type == '3')):
s_error1 = np.array(source[2])
if ((error_type == '2') or (error_type == '3')):
s_error2 = np.array(source[2])
if (debug == 1):
print("s_int[0] = ",s_int[0])
print("s_int[-1] = ",s_int[-1])
#*****************************************************
# Step 1: Scale Transform for Exponential
# Rectangulsr Components
# s_angle (2Theta in deg.) => source_chi
# s_int => source_eta
#*****************************************************
nData = s_angle.size # number of data points
kappa0 = kappas[0]
kappa1 = kappas[1]
kappa3 = np.sign(kappas[3])*np.abs(kappas[3])**3
# 1st & 3rd Cumulants of Truncated Exponential Component
k1Texp = -1.0
k3Texp = -2.0
# 1st & 3rd Cumulants of Rectangular Component
k1Rect = -0.5
k3Rect = 0.0
dchi_d2T_exp = np.power(k3Texp/kappa3, 1.0/3)
margin = marginTrans
iStart = 0
if (sign==-1): # (truncated exponential)
print('margin = ',margin)
dchi_d2T = dchi_d2T_exp
source_chi = np.cumsum(dchi_d2T)
source_chi -= dchi_d2T[0] # !
source_chi *= (s_angle[nData-1] - s_angle[0])/(nData-1)
# print('chi(E) = ',source_chi)
else: # (rectangular)
denom = kappa1 - k1Texp/dchi_d2T_exp # (maybe almost 0)
vEps = 1.0E-6 # (small number)
vInf = 1.0E+7 # (large number)
k1Rect = -0.5
dchi_d2T = np.minimum(vInf,np.maximum(vEps,k1Rect/denom))
# cum_dchi_d2T = np.cumsum(dchi_d2T)
# print('cum_dchi_d2T = ',cum_dchi_d2T)
# source_chi = np.flip(-cum_dchi_d2T) # flip
# source_chi -= dchi_d2T[0] # !
source_chi = np.cumsum(dchi_d2T)
source_chi *= (s_angle[nData-1] - s_angle[0])/(nData-1)
diff_chi = np.diff(source_chi)
# print('chi(R) = ',source_chi)
# print('diff_chi = ',diff_chi)
# non0_diff_chi = np.where(diff_chi > 0)
# print('non0_diff_chi = ',non0_diff_chi)
# source_chi_size = source_chi.size
# print('source_chi.size() = ',source_chi_size)
# expand = np.linspace(1.,1.000001,num=source_chi_size,endpoint=False)
# source_chi *= expand
if (debug==1):
print("source_chi[0] = ",source_chi[0])
print("source_chi[-1] = ",source_chi[-1],"(size: ",source_chi.size,")")
source_corr = cmn.g_corr(s_angle/2)/dchi_d2T
# cmn.g_corr(): geometrical correction
source_eta = s_int * source_corr
if (sign==-1):
source_corr2 = 1 / kappa0
source_eta *= source_corr2
if ((error_type == '1') or (error_type == '2')):
source[2] *= source_corr2
elif (error_type == '3'):
source[3] *= source_corr2
if (debug==1):
print("max(s_int[0:500]) is at angle = ",s_angle[np.argmax(s_int[0:500])])
print("d(chi)/d(2T)[0] = ",dchi_d2T[0])
print("d(chi)/d(2T)[-1] = ",dchi_d2T[-1],"(size: ",dchi_d2T.size,")")
print("source_chi[0] = ",source_chi[0])
print("source_chi[-1] = ",source_chi[-1],"(size: ",source_chi.size,")")
print("source_corr[0] = ",source_corr[0])
print("source_corr[-1] = ",source_corr[-1],"(size: ",source_corr.size,")")
print("source_eta[0] = ",source_eta[0])
print("source_eta[-1] = ",source_eta[-1],"(size: ",source_eta.size,")")
#*******************************************
# Step 2: Preparation for Fourier treatment
#*******************************************
# print('source_chi = ',source_chi)
index = cmn.indices(source_chi, margin)
nSource = index[0] # number of source data (=nData)
nValid = index[1] # index of last valid data
i1 = index[2] # index for first division point in marginal range
i2 = index[3] # index for second division point in marginal range
nDest = index[4] # total number of equidistant data
iMin = index[5] # index of minimum separation in source data
dchi = index[6] # interval of equidistant data
# dx_min = index[7] # minimum separation found in source data (not used)
if (debug==1):
print("margin = ",margin)
print("index = ",index)
print("dchi = ",dchi)
# chi0 = source_x[0] # (float)
chi=np.linspace(start=0,stop=nValid*dchi,num=nValid,endpoint=False)
#*********************************************************
# Step 3:
# Cubic spline interpolation of data to transformed scale
# (source_x, source_y) => (chi, eta)
#*********************************************************
f3 = interpolate.interp1d(source_chi,source_eta,kind="cubic",fill_value='extrapolate')
chi[nValid-1]=np.minimum(chi[nValid-1],source_chi[nData-1])
eta = f3(chi)
if (debug == 1):
print("eta[0] = ",eta[0],"eta[1] = ",eta[1])
print("eta[-1] (before padding) = ",eta[-1],"(size=",eta.size,")")
#**********************************************************
# Step 4:
# Redimension & Extrapolation for marginal (padding) range
#**********************************************************
chi = np.linspace(0, nDest*dchi, nDest, endpoint=False)
eta = np.resize(eta, nDest)
cmn.pad_margin(eta, index) # => "cmn_common.py"
if (debug == 1):
print("chi[0]=",chi[0],"chi[-1]=",chi[-1])
print("eta[-1] (after padding) = ",eta[-1],"(size=",eta.size,")")
# print("xi = ",xi,"(",xi.size,")")
# print("eta = ",eta,"(",eta.size,")")
# print("index = ",index)
#***********************************************
# Step 5: Fourier transform of intensity data
# (eta) => (ft_eta)
#***********************************************
ft_eta = np.fft.rfft(eta)
nFT = ft_eta.size
if (debug == 1):
print("ft_eta[0] = ",ft_eta[0],", ft_eta[1] = ",ft_eta[1])
print("ft_eta[-1] = ",ft_eta[-1],"(size: ",nFT,")")
#********************************************************
# Step 6: Fourier transform of approximate model profile
# => (ft_deconv), (ft_conv)
#********************************************************
if (modeTrans=='0'):
dxi=1.0/(nDest*dchi)
xi=np.linspace(0,nFT*dxi,nFT,endpoint="False")
if (sign==-1):
ft_deconv = fFtTexp(xi)
ft_conv = np.abs(ft_deconv)
elif (sign==1):
ft_deconv = fFtRect(xi)
ft_conv = np.abs(ft_deconv)
else: # if (modeTrans==1):
chi = np.linspace(-nDest*dchi, 0, nDest, endpoint=False)
if (sign==-1): # (exponential)
deconv = fTexp(chi) # see Line 118
else: # (rectangular)
deconv = fRect(chi)
deconv[0] = 1.0/dchi-np.sum(deconv[1:nDest]) # normalization
ft_deconv = np.fft.rfft(deconv) # Fourier transform
ft_conv = np.array(np.abs(ft_deconv),dtype=np.complex)
if (debug==1):
print("modeTrans = ",modeTrans)
print("ft_deconv[0] = ",ft_deconv[0],", ft_deconv[1] = ",ft_deconv[1])
print("ft_deconv[-1] = ",ft_deconv[-1],"(size:",ft_deconv.size,")")
print("ft_conv[0] = ",ft_conv[0],", ft_conv[1] = ",ft_conv[1])
print("ft_conv[-1] = ",ft_conv[-1],"(size:",ft_conv.size,")")
#***************************************************
# Step 7: Fourier transformed instrumental function
#***************************************************
# np.seterr(invalid='ignore')
np.seterr(invalid='warn')
ft_inst = ft_deconv / ft_conv # deconvolutional
if (debug==1):
print("ft_inst[0] = ",ft_inst[0],", ft_inst[1] = ",ft_inst[1])
print("ft_inst[-1] = ",ft_inst[-1],"(size:",ft_inst.size,")")
#***********************************************
# Step 8: Deconvolution/convolution treatment
#***********************************************
ft_zeta = ft_eta / ft_inst # deconvolutional
if (debug==1):
print("ft_zeta[0] = ",ft_zeta[0],", ft_zeta[1] = ",ft_zeta[1])
print("ft_zeta[-1] = ",ft_zeta[-1])
np.seterr(invalid='warn')
#*************************************
# Step 9: Inverse Fourier transform
#*************************************
zeta = np.real(np.fft.irfft(ft_zeta)) # inverse Fourier transform
if (debug==1):
print("zeta[0] = ",zeta[0])
print("zeta[-1] = ",zeta[-1],"(size: ",zeta.size,")")
# print("ft_deconv = ",ft_deconv,"(",ft_deconv.size,")")
# print("ft_conv = ",ft_conv)
# print("ft_conv/ft_deconv =",ft_conv/ft_deconv)
# print("ft_eta = ",ft_eta)
#******************************
# Step 10: Pick up intensity
#******************************
chi = np.linspace(0, nDest*dchi, nDest, endpoint = False)
f3 = interpolate.interp1d(chi, zeta, kind="cubic")
dest_eta = f3(source_chi)
#******************************
# Step 11: Rescale transform
# dest_y => d_int
#******************************
d_int = dest_eta / source_corr
d_int = np.maximum(0,d_int)
dest = source
dest[1] = d_int
if (debug==1):
print("d_int[0] = ",d_int[0],"d_int[1] = ",d_int[1])
print("d_int[-1] = ",d_int[-1],"(size:",d_int.size,")")
print("max(d_int[0:500]) is at angle = ",s_angle[np.argmax(d_int[0:500])])
print("")
#**********************************************************
# Step 12: Fourier treatment of error data (not validated)
#**********************************************************
dest = cmn.TreatErrors(error_type, source, source_chi, index, source_corr, ft_inst)
return dest