common.py – subroutines commonly used for exterm5.x system
Python code “common5.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 -*-
import numpy as np # import "numpy" module
from scipy import interpolate # import "interpolate" from "scipy"
import datetime # import "datetime" module
# "common5.py"
# Python codes commonly used for deconvolutional treatment
# Updated by T. Ida, September 20, 2023
#***********************************************
# Find Total Points and Indices of Anker Points
# from Non-equidistant Data and Margin Width
#***********************************************
def indices(x, margin):
# x[0:nSource]
# abscissa values of non-equidistant data (NumPy array)
# margin = 10.0 : margin width
nSource = x.size # (integer)
diff_x = np.diff(x) # (NumPy array) (nSource-1,)
iMin = np.argmin(diff_x) # (integer)
dx_min = diff_x[iMin] # minimum separation
valid = x[nSource-1] - x[0] # valid width
nValid = int(np.ceil(valid/dx_min)+1)
dx = valid/(nValid-1)
power2 = np.log2(2*(valid + margin)/dx)
nPower2 = int(np.ceil(power2)) # ceiling function
nDest = np.power(2,nPower2) # total number
i1 = int(np.round((2*nValid+nDest)/3))
i2 = int(np.round((nValid+2*nDest)/3))
return [nSource,nValid,i1,i2,nDest,iMin,dx_min,dx]
###################################################
# Padding for Marginal Range for Fourier Treatment
###################################################
def pad_margin(y, index):
# y: y-array (NumPy array)
# index: 7-item list (Python List)
i1 = index[1] # end of valid range
i2 = index[2] # 1st anker
i3 = index[3] # 2nd anker
i4 = index[4] # total number
iLin = np.linspace(start=0,stop=i3-i2,num=i3-i2)
yL = y[i1-1]
yR = y[0]
y[i1:i2] = yL
y[i2:i3] = yL + (yR-yL)*iLin[0:i3-i2]/(i3-i2)
y[i3:i4] = yR
#################################################
# Modification of Errors for Change in Intervals
#################################################
def mod_error(s_chi, d, dchi):
n = s_chi.size
e = np.array(d) # make a copy
e[0] = d[0] * np.sqrt((s_chi[1]-s_chi[0]) / dchi)
e[1:n-1] = d[1:n-1] * np.sqrt((s_chi[2:n]-s_chi[0:n-2]) / (2*dchi))
e[n-1] = d[n-1] * np.sqrt((s_chi[n-1]-s_chi[n-2]) / dchi)
return e
####################################################
# Re-Modification of Errors for Change in Intervals
####################################################
def remod_error(s_chi, e, dchi):
n = s_chi.size
d = np.array(e) # make a copy
d[0] = e[0] / np.sqrt((s_chi[1]-s_chi[0]) / dchi)
d[1:n-1] = e[1:n-1] / np.sqrt((s_chi[2:n]-s_chi[0:n-2]) / (2*dchi))
d[n-1] = e[n-1] / np.sqrt((s_chi[n-1]-s_chi[n-2]) / dchi)
return d
###########################################################
# Geometrical Correction Factor for
# Bragg-Brentano Diffractometer Without Diffracting Optics
###########################################################
def g_corr(degTheta):
radTheta = degTheta * np.pi / 180
# ans = 2 * np.sin(radTheta) * np.sin(2*radTheta)
ans = np.sin(radTheta) * np.sin(2*radTheta)
ans = ans / (1 + np.square(np.cos(2*radTheta)))
return ans
##############################
# Fourier Treatment of Errors
##############################
def TreatErrors(error_type,source,source_chi,index,source_corr,ft_inst):
# error_type: type of error models [0,1,2,3]
# source: source data (Python list of NumPy arrays)
# source_chi: scale-transformed abscissa
# index: 7-element Python list of anker points
# source_corr: correction factor
# ft_inst: Fourier transform of instrumental function
dest = source
nSource = index[0]
i1 = index[1] # 1st index (nValid)
i2 = index[2] # 2nd index
i3 = index[3] # 3rd index
i4 = index[4] # 4th index (nDest)
iMin = index[5]
dx_min = index[6]
dx = index[7]
nDest = i4
nDest_2 = int(nDest / 2)
# Determine Step Interval (dchi)
chi0 = source_chi[0]
chiV = source_chi[nSource-1]
dchi = dx
chiD = chi0 + nDest * dchi
# Linear interpolation function
iLin = np.linspace(0,i3-i2,i3-i2)
#********************
# Pick up error 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])
#******************************************
# Fourier Treatment of Error Data (type I)
#******************************************
if ((error_type == '1') or (error_type == '3')):
# Correct errors for scale transform and geometrical factors
source_d1 = s_error1 * source_corr
# Correction for change in interval
source_e = mod_error(source_chi, source_d1, dchi)
# Abscissa points for interpolation
chi = np.linspace(chi0, chiV, nValid, endpoint = False)
# Cubic spline interpolation
f3 = interpolate.interp1d(source_chi, source_e, kind="cubic")
eps = f3(chi)
epsL = source_e[nSource - 1]
epsR = source_e[0]
# Add marginal range
eps = np.resize(eps, nDest)
# Fill marginal range with values
eps[i1:i2] = epsL
eps[i2:i3] = epsL+(epsR-epsL)*iLin[0:i3-i2]/(i3-i2)
eps[i3:nDest] = epsR
# Inverse squared errors (1/eps**2)
inv_sq_eps = 1.0 / (eps * eps)
# Fourier transform of inverse squared errors (FT(1/eps**2))
ft_inv_sq_eps = np.fft.fft(inv_sq_eps)
# Instrumental function (inst)
inst = np.real(np.fft.ifft(ft_inst))
# Squared instrumental function (inst**2)
sq_inst = inst * inst
# Fourier transform of squared instrumental function (FT(inst**2))
ft_sq_inst = np.fft.fft(sq_inst)
# Mutual correlation in Fourier transform (!)
# print ("ft_inv_sq_eps.shape = ",ft_inv_sq_eps.shape)
# print ("ft_sq_inst.shape = ",ft_sq_inst.shape)
ft_inv_sq_eps = ft_inv_sq_eps * np.conj(ft_sq_inst)
# Inverse Fourier transform
delta = np.real(np.fft.ifft(ft_inv_sq_eps))
# Estimated error (model I)
delta = 1 / np.sqrt(np.abs(delta))
# Abscissa of Fourier transform
chi = np.linspace(chi0, chiD, nDest, endpoint = False)
# Cubic spline interpolation
f3 = interpolate.interp1d(chi, delta, kind="cubic")
dest_e = f3(source_chi)
# Recorrect error values for change in intervals
dest_d = remod_error(source_chi, dest_e, dchi)
# Recorrect error data for scale transform and geometrical factors
d_error1 = dest_d / source_corr
#*******************************************
# Fourier Treatment of Error Data (type II)
#*******************************************
if ((error_type == '2') or (error_type == '3')):
# Correct error data for scale transform & geometrical factors
source_d2 = s_error2 * source_corr
# Correct error values for change in intervals
source_e = mod_error(source_chi, source_d2, dchi)
# Abscissa of Fourier transform
chi = np.linspace(chi0, chiV, nValid, endpoint = False)
# Cubic spline interpolation
f3 = interpolate.interp1d(source_chi, source_e, kind="cubic")
eps = f3(chi)
epsL = source_e[nSource - 1]
epsR = source_e[0]
# Add marginal range
eps = np.resize(eps, nDest)
# Fill marginal range with values
eps[i1:i2] = epsL
eps[i2:i3] = epsL+(epsR-epsL)* iLin[0:i3-i2]/(i3-i2)
eps[i3:nDest] = epsR
# Inverse instrumental function
inv_inst = np.real(np.fft.ifft(1 / ft_inst))
# Squared inverse instrumental function
sq_inv_inst = inv_inst * inv_inst
# Fourier transform of squared inversse instrumental function
ft_sq_inv_inst = np.fft.fft(sq_inv_inst)
# Squared errors
sq_eps = eps * eps
# Fourier transform of squared errors
ft_sq_eps = np.fft.fft(sq_eps)
# Convolution in Fourier transform
ft_sq_eps = ft_sq_eps * ft_sq_inv_inst
# Invers Fourier transform
sq_eps = np.real(np.fft.ifft(ft_sq_eps))
# Estimated error in interpolated data
delta = np.sqrt(np.abs(sq_eps))
# Abscissa of Fourier transform
chi = np.linspace(chi0, chiD, nDest, endpoint = False)
# Cubic spline interpolation
f3 = interpolate.interp1d(chi, delta, kind="cubic")
dest_e = f3(source_chi)
# Recorrect error values for change in intervals
dest_d = remod_error(source_chi, dest_e, dchi)
# Recorrect error data for scale transform & geometrical factors
d_error2 = dest_d / source_corr
#********************
# Replace error data
#********************
if (error_type == '1'):
dest[2] = d_error1
elif (error_type == '2'):
dest[2] = d_error2
elif (error_type == '3'):
dest[2] = d_error1
dest[3] = d_error2
return dest
##############
# Report Time
##############
def ReportTime(s):
# Argument s: strings to be shown after date time
d = datetime.datetime.today() # current date time
# Show current date time
print('{0:02d}:{1:02d}:{2:02d}'.format(d.hour, d.minute, d.second),' ',s)
return 0