equatorial1.py – treatment of equatorial aberration (mode 1)
Implementation of deconvolutional treatment about equatorial aberration in XRD data collected for finite-width specimen with a semiconductor strip (PIN photodiode array) X-ray detector [Ida, 2020, 2021].
Python code “equatorial1.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 -*-
# equatorial1.py
# Coded by T. Ida, corrected by N. Sakurai,
# Updated by T. Ida, December 11, 2023
# Updated by T. Ida, February 7, 2024
import numpy as np # import "numpy" module
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 Reduced Cumulants of Equatorial Aberration
# Output units are [1,deg.,deg.,deg.,deg.]
#######################################################
def calc_kappas(degTwot, cfg):
# degTwot: 2Theta in [deg.] (NumPy array)
# cfg=[degDS,degLPSD,specW,gonioR,marginEquat,modeEquat,skipEquat,error_type]
# cfg = [degDS,degLPSD,specW,gonioR,marginEquat]
# cfg += [alphaEquat,modeEquat,skipEquat,error_type]
degDS = cfg[0]
degPsi = cfg[1]/2
W_2R = 0.5*cfg[2]/cfg[3]
#*******************************************
# Exact Expression of Deviation Angle
# Caused by Equatorial Aberration for given
# values of (theta, phi, psi)
# Output is in [rad.]
#*******************************************
def deltaTwot(theta, phi, psi):
# theta: Half of 2Theta in [rad.] (NumPy array)
# phi: equatorial deviation angle of incident beam in [rad.] (NumPy array)
# psi: half of offset angle of detector strip in [rad.] (float number)
g = np.cos(theta-psi) - np.sin(theta-psi) / np.tan(theta-psi-phi)
g = g * np.cos(2 * psi)
g = g + np.cos(theta + psi)
deltaTwot = theta + phi + psi
deltaTwot -= np.arctan(np.sin(theta + psi) / g)
return deltaTwot
nSize = degTwot.size # number of values of degTwot
theta = degTwot * np.pi / 360 # shape of theta is (nSize) here
sinTheta = np.sin(theta) # shape of sinTheta is (nSize) here
phi = degDS * np.pi / 180 # float number
psi = degPsi * np.pi / 180 # float number
# Note: effective equatorial divergence angle should be
# reduced from divergence slit open angle Phi_DS
# at 2Theta angles lower than 2ThetaC.
theta = theta.repeat(4).reshape(nSize,4) # shape of theta is (nSize,4)
# 4-point Gauss-Legendre abscissa and weights
x_j = np.array([0.069431844203,0.330009478208,0.669990521792,0.930568155797])
w_j = np.array([0.173927422569,0.326072577431,0.326072577431,0.173927422569])
# Set values of psi
psi_j = -psi/2 + x_j * psi # shape of psi_j is (4)
psi_j = np.tile(psi_j, nSize).reshape(nSize,4) # shape of psi_j is (nSize,4)
# Set limit values of phi for each of (theta, psi)
phiMin = theta - psi_j - np.arctan(np.sin(theta-psi_j)/(np.cos(theta-psi_j)-W_2R))
phiMin = np.maximum(-phi/2, phiMin) # shape of phiMin is (nSize, 4)
phiMax = theta - psi_j - np.arctan(np.sin(theta-psi_j)/(np.cos(theta-psi_j)+W_2R))
phiMax = np.minimum(phi/2, phiMax) # shape of phiMax is (nSize, 4)
# Set matrix elements ...
phi0 = np.tile(phiMin,4).reshape(nSize,4,4)
phi1 = np.tile(phiMax-phiMin,4).reshape(nSize,4,4)
phix = np.tile(x_j.repeat(4),nSize).reshape(nSize,4,4)
phi_ij = phi0 + phix * phi1
theta = np.tile(theta, 4).reshape(nSize,4,4) # shape of theta is (nSize,4,4)
psi_ji = np.tile(psi_j, 4).reshape(nSize,4,4) # shape of psi_j is (nSize,4,4)
deltaTwot_ij = deltaTwot(theta, phi_ij, psi_ji)
# print("deltaTwot_ij = ",deltaTwot_ij)
# results should be ...
# / d2theta(phi1,psi1) d2theta(phi1,psi2) d2theta(phi1,psi3) d2theta(phi1,psi4) \
# | d2theta(phi2,psi1) d2theta(phi2,psi2) d2theta(phi2,psi3) d2theta(phi2,psi4) |
# | d2theta(phi3,psi1) d2theta(phi3,psi2) d2theta(phi3,psi3) d2theta(phi3,psi4) |
# \ d2theta(phi4,psi1) d2theta(phi4,psi2) d2theta(phi4,psi3) d2theta(phi4,psi4) /
# ...
deltaTwot0_ij = np.array([1.0]).repeat(nSize*4*4).reshape(nSize,4,4)
deltaTwot2_ij = deltaTwot_ij**2 # not matrix, but element operation
deltaTwot3_ij = deltaTwot2_ij * deltaTwot_ij # element operation
deltaTwot4_ij = deltaTwot2_ij**2 # element operation
s0 = np.dot(np.dot(phi1 * deltaTwot0_ij, w_j.T) / phi, w_j.T)
s1 = np.dot(np.dot(phi1 * deltaTwot_ij, w_j.T) / phi, w_j.T)
s2 = np.dot(np.dot(phi1 * deltaTwot2_ij, w_j.T) / phi, w_j.T)
s3 = np.dot(np.dot(phi1 * deltaTwot3_ij, w_j.T) / phi, w_j.T)
s4 = np.dot(np.dot(phi1 * deltaTwot4_ij, w_j.T) / phi, w_j.T)
# Calculate Cumulants up to 4th order
kappa1 = s1/s0
kappa2 = s2/s0 - (s1/s0)**2
kappa3 = s3/s0 - 3*s2*s1/s0**2 + 2*(s1/s0)**3
kappa4 = s4/s0 - 4*s3*s1/s0**2 - 3*(s2/s0)**2
kappa4 += 12*(s2/s0)*(s1/s0)**2 - 6*(s1/s0)**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]
##################################################
# FormFuncEquat(chi, rho)
# Define Equatorial Aberration Model Function for
# "chiS" and "chiT" scales
# (2nd-order Approximation)
##################################################
def FormFuncEquat(chi, rho):
# chi: NumPy array
# rho: float number
chiMin = -3 * ( 1 + rho )
if ( rho < 2 ):
chiMax = 0.75 * rho ** 2
else:
chiMax = -3 * ( 1 - rho )
if ( rho == 0 ):
y = np.where((-3 < chi)&(chi < 0), 1/(6*np.sqrt(-chi/3)), 0)
else:
D = rho**2 - chi / 0.75
sqrtD = np.sqrt(np.maximum(0, D))
phiL = np.maximum(-rho + sqrtD, rho - sqrtD )
phiU = np.minimum(2, rho + sqrtD )
y = np.where((D>0)&(chiMin<chi)&(chi<chiMax), np.log(phiU/phiL)/(6*rho), 0)
return y
################################################
# "equatorial1.treat"
# TREATMENT OF EQUATORIAL ABERRATION, MAIN BODY
################################################
def treat(sign, source, kappas, cfg):
# sign: -1 or 1
# source: Source Data (Python List)
# [s_angle, s_intensity],
# [s_angle, s_intensity, s_error1], or
# [s_angle, s_intensity, s_error1, s_error2]
# kappas: Cumulants (NumPy array)
# cfg: Configuration Parameters (Python list)
# cfg = [degDS,degLPSD,specW,gonioR,marginEquat]
# cfg += [alphaEquat,modeEquat,skipEquat,error_type]
#******************************
# Configuration from "dct.cfg"
#******************************
degDS = cfg[0] # divergence slit open angle (deg.)
degLPSD = cfg[1] # LPSD view angle (deg.)
specW = cfg[2] # specien width (mm)
gonioR = cfg[3] # goniometer radius (mm)
W_R = specW / gonioR
marginEquat = cfg[4] # margin for equatorial treatment
alphaEquat = cfg[5]
modeEquat = cfg[6] # mode for equatorial treatment
skipEquat = cfg[7] # skip process
error_type = cfg[8] # error type
#******************************************************
# Step 0-0: Parse Python List "source" to NumPy arrays
#******************************************************
s_angle = np.array(source[0]) # angle data (NumPy array)
s_intensity = np.array(source[1]) # intensity data (NumPy array)
kappa1 = kappas[1]
kappa3 = np.sign(kappas[3])*np.abs(kappas[3])**3
#----------------------------------------
# Optional Output for 'skipEquat' == '2'
#----------------------------------------
if (skipEquat == '2'):
CUMULANT_FILE = 'cumulants_equatorial.csv'
save_cumulant = np.insert(kappas, 0, s_angle, axis=0)
save_cumulant = np.array(save_cumulant).T # transpose matrix
fmt_list = ["%.3f","%.5e","%.5e","%.5e","%.5e","%.5e"]
np.savetxt(CUMULANT_FILE, save_cumulant, delimiter=",",fmt=fmt_list)
if (sign==-1):
###################################################################
# Stage 1 (-1): Correction of the 3rd cumulant of equatorial aberration
# (Correction about Asymmetric Deformation)
###################################################################
rho = (degLPSD / 2) / degDS # Note: "psi" is half of the view angle
k3 = -16/35*(1+21*rho**2/4) # 3rd cumulant of model profile on chi3-scale
dchi3_d2Theta = np.power(k3/kappa3, 1/3) # = d(chi3)/d(2Theta)
# Integration & creation of abscissa values source_x
# Note: chi3 should be the integral of d(chi3)/d(2Theta) by 2Theta
# (non-equidistantly located on "chi3"-scale)
source_x = np.cumsum( dchi3_d2Theta )
source_x *= (s_angle[1]-s_angle[0])
# Geometrical Correction
# (Lorentz, polarization & powder diffractometry corrections)
source_corr = cmn.g_corr(s_angle / 2)
# Ordinate Correction Factor
# Note: "source_y * d(chi3)" should be proportional
# to "s_intensity * d(2Theta)"
source_corr /= dchi3_d2Theta
#------------------------------------------------------------
# Another correction for loss of intensity,
# caused by spill-over of irradiated area from specimen face
#------------------------------------------------------------
source_corr2 = 1 / kappas[0]
# Corrected and Transformed Intensity
source_y = s_intensity * source_corr * source_corr2
# Error data are also corrected
if ((error_type == '1') or (error_type == '2')):
source[2] = source[2] * source_corr2
elif (error_type == '3'):
source[3] = source[3] * source_corr2
#******************************************************************
# Step 1-2: Preparation for Fourier treatment
# Determine Marginal (Padding) Range, Total Points & Step Interval
#******************************************************************
# marginModel: Standard deviation of model profile
marginModel = np.sqrt(0.8*(1+1.25*rho**2))
# margin: defines marginal (padding) range
margin = marginEquat * marginModel
# marginEquat is the value read from the file'dct.cfg'
# User can change the value in the current version.
index = cmn.indices(source_x, margin)
# Derive total number and 1st,2nd,3rd anker points
# cmn.indices(x, w_margin) is defined in "common5.py"
# x[..] = source_x : abscissa (NumPy array)
nSource = index[0] # total number of source data
nValid = index[1] # number of valid equidistant data
i1 = index[2] # index of 1st anker point
i2 = index[3] # index of 2nd anker point
nDest = index[4] # number of total equidistant data
iMin = index[5]
dx_min = index[6]
nDest_2 = int(nDest/2) # half of the total number
dchi3 = index[7] # dchi3: Step interval of equidistant data
rho = (degLPSD / 2) / degDS
# Prepare equidistant data for valid range
# nValid = int((source_x[nSource-1] - source_x[0]) / dchi3) + 1
chi0 = source_x[0] # 1st location on chi3-scale (zero)
wchi3 = nValid * dchi3
chi3 = np.linspace(chi0, chi0 + wchi3, nValid, endpoint = False)
#*********************************************************
# Step 1-3:
# Cubic spline interpolation of Data to Transformed Scale
# (source_x, source_y) => (chi3, eta3)
#*********************************************************
# Cubic spline interpolation
# print("source_x = ",source_x)
# print("chi3 = ",chi3)
chi3[-1]=source_x[-1] # adjust the last value
f3 = interpolate.interp1d(source_x, source_y, kind="cubic")
eta3 = f3(chi3)
#********************************************
# Step 1-4:
# Extrapolation for marginal (padding) range
# (chi3, eta3) => (chi3, eta3)
#********************************************
# Fill marginal range with values
wchi3 = nDest * dchi3
chi3 = np.linspace(chi0, chi0 + wchi3, nDest, endpoint = False)
eta3 = np.resize(eta3, nDest)
cmn.pad_margin(eta3, index) # => common5.pad_margin()
#***********************************************
# Step 1-5: Fourier Transform of Intensity Data
# (eta3) => (ft_eta3)
#***********************************************
# Fourier transform of intensity
ft_eta3 = np.fft.fft(eta3)
#**********************************************************
# Step 1-6: Fourier Transform of Approximate Model Profile
# Function To Be Deconvolved (Equatorial Model)
# => (ft_deconv3), (ft_conv3)
#**********************************************************
# Prepare Fourier transform of function to be deconvolved
# Set start point to be zero, and use compliment expression ・・
chi3 = np.linspace(0, wchi3, nDest, endpoint=False)
chi3[nDest_2:nDest] -= nDest * dchi3 # make compliment
deconv3 = np.zeros(nDest)
deconv3[1:nDest] = FormFuncEquat(chi3[1:nDest], rho)
# Adjustment to avoid singularity of model function at x = 0
# Note (deconv3[0] + ... + deconv3[nDest-1]) * dchi3 = 1)
deconv3[0] = 1.0 / dchi3 - np.sum(deconv3[1:nDest])
# Fourier Transform of Function To Be Deconvolved (Equatorial Model)
ft_deconv3 = np.fft.fft(deconv3)
# Fourier Transform of Function To Be Convolved
ft_conv3 = np.abs(ft_deconv3)
#******************************************************************
# Step 1-7: Fourier Transform about Artificial Shift Introduction
# Function To Be Deconvolved is Dirac delta function "delta(x+k1)"
# => (ft_inst3)
#******************************************************************
# Fourier transform of shift function to be deonvolved
# Note: Shift by -k1 = +1.0 is introduced
# It will cause (leave) Shift by Kappa1 on 2Theta scale
# Shift by 1 should be introduced to function to be deconvolved.
# Shift by -k1 = +1.0 on Fourier transform is introduced.
# Step interval of Fourier transform
dxi3 = 1.0 / (dchi3 * nDest)
xi3 = np.linspace(0, nDest*dxi3, nDest, endpoint=False)
xi3[nDest_2:nDest] -= nDest * dxi3 # ... ?
# xi[0] = 0, xi[1] = dxi, ...,
# ft_shift3 = np.exp(2 * np.pi * 1.0j * xi3) # Ida's erroneous (?) code
ft_shift3 = np.exp(-2 * np.pi * 1.0j * xi3) # corrected by N. Sakurai
# Fourier transform of instrumental function
# ( to be used for error estimation )
ft_inst3 = ft_deconv3 * ft_shift3/ ft_conv3
#**********************************************************
# Step 1-8: Here is the Core of DCT,
# Deconvolution/Convolution Treatment on Fourier transform
# (ft_eta3)/(ft_inst3) => (ft_zeta3)
#**********************************************************
# Deconvolution and Convolution Treatment
ft_zeta3 = ft_eta3 / ft_inst3
#*************************************
# Step 1-9: Inverse Fourier Transform
# (ft_zeta3) => (zeta3)
#*************************************
# Inverse Fourier transform of DCT data
zeta3 = np.real(np.fft.ifft(ft_zeta3))
#***********************************************
# Step 1-10: Pick up Intensity from DCT data
# on "chi3" scale by cubic spline interpolation
# (chi3, zeta3) => (source_x, dest3_y)
#***********************************************
# Reset equidistant abscissa
wchi3 = nDest * dchi3
chi3 = np.linspace(chi0, chi0 + wchi3, nDest, endpoint = False)
# Cubic spline interpolation, (chi,zeta) => (source_x, dest_y)
f3 = interpolate.interp1d(chi3, zeta3, kind="cubic")
dest3_y = f3(source_x)
#**********************************************
# Step 1-11: Rescale transform
# (source_x, dest3_y) => (s_angle, d3_intensity)
#**********************************************
d3_intensity = dest3_y / source_corr
dest = source
dest[1] = d3_intensity
#*******************************************************
# Step 1-12 Fourier Treatment of Error Data, on Stage 1
#*******************************************************
dest = cmn.TreatErrors(error_type, source, source_x, index, source_corr, ft_inst3)
else:
##########################################
# Stage 2: Correction of the 1st cumulant
# (Shift Correction)
##########################################
#******************************************************
# Step 2-0: Parse Python List "source" to numpy.array
# Attention: Missed in a previous version
#******************************************************
s_angle = np.array(source[0])
s_intensity = np.array(source[1])
#********************************************************
# Step 2-1: Scale Transform for the 1st cumulant, "chi1"
# (s_angle, s_intensity) => (source_x, source_y)
#********************************************************
# Abscissa values (Non-equidistantly separated)
# kappa1 = kappas[1] # 1st order cumulant of Delta2Theta in [rad.]
k1 = -1 # 1st order cumulant of the model function on chi1-scale
dchi1_d2Theta = k1 / kappa1
# NOTE: the values should be similar to (k3 / kappa3)**(1/3)
# but should slightly be different.
source_x = np.cumsum( dchi1_d2Theta ) # integration
source_x *= (s_angle[1]-s_angle[0])
# Ordinate correction factor
source_corr = cmn.g_corr(s_angle / 2) / dchi1_d2Theta
# Corrected intensity
source_y = s_intensity * source_corr
#******************************************************************
# Step 2-2: Preparation for Fourier treatment
# Determine Marginal (Padding) Range, Total Points & Step Interval
#******************************************************************
rho = (degLPSD / 2) / degDS # Note: "psi" is half of the view angle
marginModel = np.sqrt(0.8*(1+1.25*rho**2))
margin = marginEquat * marginModel
# marginEquat is the value read from the file'dct.cfg'
# User can change the value in the current version.
# NOTE: The use of the same value of margin may not be justified.
index = cmn.indices(source_x, margin)
# Derive total number and 1st,2nd,3rd anker points
# dct_cmn.indices(x, w_margin) is defined in "dct_cmn.py"
# x[..] = source_x : abscissa (NumPy array)
nSource = index[0] # number of source data points
nValid = index[1] # index of 1st anker point on interpolated abscissa
i2 = index[2] # index of 2nd anker point
i3 = index[3] # index of 3rd anker point
nDest = index[4] # total number of interpolated data
dchi1 = index[7] # Step Interval for "chi1" scale
nDest_2 = int(nDest/2) # half of the total number
# Prepare equidistant data for valid range
chi0 = source_x[0]
wchi1 = nValid * dchi1
chi1 = np.linspace(chi0, chi0 + wchi1, nValid, endpoint = False)
#*********************************************************
# Step 2-3:
# Cubic spline interpolation of Data on Transformed Scale
# (source_x, source_y) => (chi1, eta1)
#*********************************************************
# Cubic spline interpolation
# eta = np.array(chi)
chi1[-1]=source_x[-1] # adjust the last value
f1 = interpolate.interp1d(source_x, source_y, kind="cubic")
eta1 = f1(chi1) # (source_x, source_y) => (chi1, eta1)
#********************************************
# Step 2-4:
# Extrapolation for marginal (padding) range
# (chi1, eta1) => (chi1, eta1)
#********************************************
# Fill marginal range with extrapolated values
wchi1 = nDest * dchi1
chi1 = np.linspace(chi0, chi0 + wchi1, nDest, endpoint = False)
eta1 = np.resize(eta1, nDest)
cmn.pad_margin(eta1, index) # => "common5.py"
#***********************************************
# Step 2-5: Fourier Transform of Intensity Data
# (eta1) => (ft_eta1)
#***********************************************
# Fourier transform of intensity
ft_eta1 = np.fft.fft(eta1)
# Step interval of Fourier transform
dxi1 = 1.0 / (dchi1 * nDest)
# Values of abscissa
xi1 = np.linspace(0, nDest * dxi1, nDest, endpoint=False)
xi1[nDest_2:nDest] -= nDest * dxi1
#***************************************************
# Step 2-6: Fourier Transform of Shit function
# Function To Be Deconvolved (Dirac delta function)
# => (ft_inst1)
#***************************************************
# Fourier transform of function to be deonvolved
# Note: Shift by k1 = -1 has been introduced on Stage 1
# Shift by +1 should be introduced here.
# Shift by k1 = -1.0 on Fourier transform is given by
# ft_shift1 = np.exp(-2 * np.pi * 1.0j * xi1) # Ida's erroneous (?) code
ft_shift1 = np.exp(2 * np.pi * 1.0j * xi1) # corrected by N. Sakurai
ft_inst1 = ft_shift1
#************************************************
# Step 2-7: Shift Correction by Fourier method
# Deconvolutional Treatment on Fourier transform
# (ft_eta1)/(ft_inst1) => (ft_zeta1)
#************************************************
ft_zeta1 = ft_eta1 / ft_inst1
#*************************************
# Step 2-8: Inverse Fourier Transform
# (ft_zeta1) => (zeta1)
#*************************************
zeta1 = np.real(np.fft.ifft(ft_zeta1))
# Reset equidistant abscissa
wchi1 = nDest * dchi1
chi1 = np.linspace(chi0, chi0 + wchi1, nDest, endpoint = False)
#******************************************************
# Step 2-9: Pick up Intensity from shift treated data
# on "chi3" scale by cubic spline interpolation
# (chi1, zeta1) => (source_x, dest_y)
#*****************************************************
f1 = interpolate.interp1d(chi1, zeta1, kind="cubic")
dest_y = f1(source_x)
#***********************************************
# Step 2-10: Rescale transform
# (source_x, dest_y) => (2Theta, d_intensity)
#***********************************************
d_intensity = dest_y / source_corr
dest = source
dest[1] = d_intensity
#*******************************************************
# Step 2-11: Fourier Treatment of Error Data on Stage 2
#*******************************************************
dest = cmn.TreatErrors(error_type, source, source_x, index, source_corr, ft_inst1)
return dest