equatorial3.py – treatment of equatorial aberration (mode 3)
Python code “equatorial3.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 -*-
# equatorial3.py
# Coded by T. Ida, September 14, 2023
# Updated by T. Ida, September 20, 2023
# Modified by T. Ida, Jan. 11, 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
# by 5x5 sampling-point Gauss-Legendre quadrature
# 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]
# 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.
nSamp = 5 # number of sampling points
theta = theta.repeat(nSamp).reshape(nSize,nSamp)
# shape of theta is (nSize,nSamp)
# 5-point Gauss-Legendre abscissa and weights
x_j = np.array([0.046910,0.230765,0.500000,0.769235,0.953090])
w_j = np.array([0.118463,0.239314,0.284444,0.239314,0.118463])
# Set values of psi
psi_j = -psi/2 + x_j * psi
# shape of psi_j is (5)
psi_j = np.tile(psi_j, nSize).reshape(nSize,nSamp)
# shape of psi_j is (nSize,5)
# 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,nSamp)
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,nSamp)
# Set matrix elements ...
phi0 = np.tile(phiMin,nSamp).reshape(nSize,nSamp,nSamp)
phi1 = np.tile(phiMax,nSamp).reshape(nSize,nSamp,nSamp)
phix = np.tile(x_j.repeat(nSamp),nSize).reshape(nSize,nSamp,nSamp)
phi_ij = phi0 + phix * (phi1 - phi0)
theta = np.tile(theta,nSamp).reshape(nSize,nSamp,nSamp)
# shape of theta is (nSize,nSamp,nSamp)
psi_ji = np.tile(psi_j,nSamp).reshape(nSize,nSamp,nSamp)
# shape of psi_j is (nSize,nSamp,nSamp)
# print("phi_ij = ",phi_ij)
# print("psi_ji = ",psi_ji)
deltaTwot_ij = deltaTwot(theta, phi_ij, psi_ji)
# print("deltaTwot_ij = ",deltaTwot_ij)
# results should be ...
# / d2T(phi0,psi0) d2T(phi0,psi1) d2T(phi0,psi2) d2T(phi0,psi3) d2T(phi0,psi4) \
# | d2T(phi1,psi0) d2T(phi1,psi1) d2T(phi1,psi2) d2T(phi1,psi3) d2T(phi1,psi4) |
# | d2T(phi2,psi0) d2T(phi2,psi1) d2T(phi2,psi2) d2T(phi2,psi3) d2T(phi2,psi4) |
# | d2T(phi3,psi0) d2T(phi3,psi1) d2T(phi3,psi2) d2T(phi3,psi3) d2T(phi3,psi4) |
# \ d2T(phi4,psi0) d2T(phi4,psi1) d2T(phi4,psi2) d2T(phi4,psi3) d2T(phi4,psi4) /
# ...
fPhi = (phi1-phi0)/phi # (nSize,nSamp,nSamp)
deltaTwot0_ij = fPhi
deltaTwot1_ij = deltaTwot0_ij*deltaTwot_ij
deltaTwot2_ij = deltaTwot1_ij*deltaTwot_ij
deltaTwot3_ij = deltaTwot2_ij*deltaTwot_ij
deltaTwot4_ij = deltaTwot3_ij*deltaTwot_ij
s0 = np.dot(np.dot(deltaTwot0_ij, w_j.T), w_j.T)
s1 = np.dot(np.dot(deltaTwot1_ij, w_j.T), w_j.T)
s2 = np.dot(np.dot(deltaTwot2_ij, w_j.T), w_j.T)
s3 = np.dot(np.dot(deltaTwot3_ij, w_j.T), w_j.T)
s4 = np.dot(np.dot(deltaTwot4_ij, w_j.T), 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 + 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]
##############################################
# Define Peak Profile Model for Deconvolution
# for sign == 1/-1,
##############################################
def fDeconv(alpha, sign, chi):
# alpha: shape parameter (float)
# sign: -1 or 1 (integer)
# chi: (NumPy array)
s_chi=sign*chi
ans = np.where(s_chi<0.0,0.0,np.where(s_chi<1.0,alpha*s_chi**(alpha-1),0))
return ans
###################################################################
# Define Fourier Transform of Peak Profile Model for Deconvolution
# for sign == 1/-1,
###################################################################
def fFtDeconv(alpha, sign, xi):
# alpha # (float)
# sign = 1, -1
gam = special.gammainc(alpha,1)
gam *= special.gamma(alpha)
ans = alpha * gam / np.power(sign * 2 * np.pi * 1j * xi, alpha)
return ans
################################################
# "equatorial3.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]
# + [alphaEquat,modeEquat,skipEquat,error_type]
#******************************
# Configuration from "dct.cfg"
#******************************
marginEquat=cfg[4] # Margin for FFT
alphaEquat=cfg[5] # Shape Parameter
modeEquat=cfg[6]
skipEquat=cfg[7]
error_type=cfg[8] # error type
#******************************************************
# Step 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)
nData = s_intensity.size
cum1 = kappas[1]
cum3 = np.sign(kappas[3])*np.fabs(kappas[3])**3
# print("kappa3 =",kappa3)
s1 = alphaEquat/(1+alphaEquat) # 1st-order cumulant of model function
s2 = alphaEquat/(2+alphaEquat)
s3 = alphaEquat/(3+alphaEquat)
k1 = s1
k3 = s3-3*s2+s1+2*s1**3
#-----------------------------
# d(chi)/d(2Theta) (dchi_d2T)
#-----------------------------
def f1(sign,cum1,cum3,k1,k3):
d = k1*cum3/(3*k3*cum1)-cum1**2/(12*k1**2)
# print(d)
ans = 1/( sign*cum1/(2*k1) + np.sqrt(d) )
return ans
dchi_d2T = f1(sign,cum1,cum3,k1,k3)
# Integration & creation of abscissa values source_x
# Note: chi should be the integral of d(chi)/d(2Theta) by 2Theta
# (non-equidistantly located on "chi3"-scale)
source_x = np.cumsum(dchi_d2T) - dchi_d2T[0]
source_x *= (s_angle[nData-1]-s_angle[0])/(nData-1)
# Ordinate Correction Factor
# Note: "source_y * d(chi)" should be proportional
# to "s_intensity * d(2Theta)"
# Ordinate Scale Transform Conbined with Geometrical Correction
# (Lorentz, polarization & powder diffractometry corrections)
source_corr = cmn.g_corr(s_angle / 2) / dchi_d2T
#------------------------------------------------------------
# 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
# print("source_y = ",source_y)
if (sign==-1): # intensity correction should be once...
source_y *= source_corr2
# print("source_y = ",source_y)
# 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 2: Preparation for Fourier treatment
# Determine Marginal (Padding) Range, Total Points & Step Interval
#******************************************************************
# margin: defines marginal (padding) range
margin = marginEquat
# 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]
dchi = index[7] # dchi: Step interval of equidistant data
nDest_2 = int(nDest/2) # half of the total number
# 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)
chi = np.linspace(0,nValid*dchi,nValid,endpoint=False)
chi[nValid-1]=np.minimum(source_x[nSource-1],chi[nValid-1])
# print("source_x = ",source_x)
# print("chi = ",chi)
#*********************************************************
# Step 3:
# Cubic spline interpolation of Data to Transformed Scale
# (source_x, source_y) => (chi, eta)
#*********************************************************
# Cubic spline interpolation
f3 = interpolate.interp1d(source_x, source_y, kind="cubic")
eta = f3(chi)
#********************************************
# Step 4:
# Extrapolation for marginal (padding) range
# (chi, eta) => (chi, eta)
#********************************************
# Fill marginal range with values
chi = np.linspace(0,nDest*dchi,nDest,endpoint = False)
eta = np.resize(eta,nDest)
cmn.pad_margin(eta,index) # => dct_common.pad_margin()
#***********************************************
# Step 5: Fourier Transform of Intensity Data
# (eta) => (ft_eta)
#***********************************************
# Fourier transform of intensity
ft_eta = np.fft.fft(eta)
# Step interval of Fourier transform
dxi = 1.0 / (dchi * nDest)
# Prepare Abscissa Data Using compliment expression
# Abscissa of Fourier transform
xi = np.linspace(0, nDest * dxi, num=nDest, endpoint=False)
xi[nDest_2:nDest] = xi[nDest_2:nDest] - nDest * dxi
#*********************************************************
# Step 6: Fourier Transform of Approximate Model Profile
# Function To Be Deconvolved (Equatorial Model)
# => (ft_deconv), (ft_conv)
#*********************************************************
# Function to Be Deconvolved (+/-)
if (sign == -1):
chi = np.linspace(-nDest*dchi,nDest,endpoint=False)
else:
chi = np.linspace(0,nDest*dChi,nDest,endpoint=False)
deconv = fDeconv(alphaEquat,sign,chi)
deconv[0] = 1.0/dchi - np.sum(deconv[1:nDest])
print("chi = ",chi)
print("deconv = ",deconv)
ft_deconv = np.fft.rfft(deconv)
# Fourier Transform of Function To Be Deconvolved (+/-)
# ft_deconv = fFtDeconv(alphaEquat, sign, xi)
# Fourier transform of function to be convolved
ft_conv = np.abs(ft_deconv) # (complex absolute value)
# Fourier transform of instrumental function
ft_inst = ft_deconv / ft_conv
#**********************************************************
# Step 7: Here is the Core of DCT,
# Deconvolution/Convolution Treatment on Fourier transform
# (ft_eta)/(ft_inst) => (ft_zeta)
#**********************************************************
# Deconvolution and Convolution Treatment
ft_zeta = ft_eta/ft_inst
#*************************************
# Step 8: Inverse Fourier Transform
# (ft_zeta) => (zeta)
#*************************************
# Inverse Fourier transform of DCT data
zeta = np.real(np.fft.ifft(ft_zeta))
#***********************************************
# Step 9: 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(0, nDest*dchi, nDest, endpoint = False)
# Cubic spline interpolation, (chi,zeta) => (source_x, dest_y)
f3 = interpolate.interp1d(chi, zeta, kind="cubic")
dest_y = f3(source_x)
dest_y = np.maximum(0,dest_y)
#**********************************************
# Step 10: Rescale transform
# (source_x, dest3_y) => (s_angle, d_intensity)
#**********************************************
d_intensity = dest_y / source_corr
dest = source
dest[1] = d_intensity
# if (sign==-1):
# print("*** Treatment of Equatorial Aberration ***")
# print("s_intensity = ",s_intensity)
# print("d_intensity = ",d_intensity)
#*********************************************************
# Step 11 Fourier Treatment of Error Data (not validated)
#*********************************************************
dest = cmn.TreatErrors(error_type, source, source_x, index, source_corr, ft_inst)
return dest