main routine of exterm5.x system
Python source code of exterm5.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 -*-
# exterm5.py
# T. Ida, Sep. 8, 2021
# Modified by T. Ida, March 1, 2023
# Modified by T. Ida, June 23, 2023 (version 4.3)
# Modified by T. Ida, September 21, 2023 (version 5.0)
# Modified by T. Ida, November 20, 2023 (version 5.1)
# Modified by T. Ida, December 4, 2023 (version 5.2)
# Modified by T. Ida, December 11, 2023 (version 5.3)
# Modified by T. Ida, December 29, 2023 (version 5.4)
# Modified by T. Ida, Jan. 12, 2024 (version 5.5)
import numpy as np # import "numpy" module as "np"
import scipy.special as sp # import "scipy.special" as "sp"
import sys # import "sys" module to get command-line argument "argv"
import os # import "os" module to use "os.path.exists()"
import configparser # import "configparser"
import common5 as cmn # import "common5.py" as cmn
import xray5 as xray # import "xray5.py" as xray
import axial5 as axial # import "axial5.py" as axial
import equatorial1 as equat1 # import "equatorial1.py" as equat1
import equatorial2 as equat2 # import "equatorial2.py" as equat2
import equatorial3 as equat3 # import "equatorial2.py" as equat2
import trnspr1 as trnspr1 # import 'trnspr1.py' as trnspr1
import trnspr2 as trnspr2 # import 'trnspr2.py' as trnspr2
import trnspr3 as trnspr3 # import 'trnspr3.py' as trnspr3
import cntmn5 as cntmn # import 'cntmn.py' for decontamination
#################
# 1. PREPARATION
#################
###################################
# 1-1 Treat Command-Line Arguments
###################################
argvs = sys.argv # list of command-line arguments
argc = len(argvs) # number of arguments + 1
# argvs[1] is input file name
# argvs[2] is output file name
# DPATH = os.path.dirname(sys.argv[0]) # directory of the EXE file
DPATH = ''
if (argc == 1): # in case no arguments
argvs = argvs + [os.path.join(DPATH, '00raw.csv')]
argc = 2
if (argc == 2): # in case 1 argument
argvs = argvs + [os.path.join(DPATH, '01dct.csv')]
argc = 3
###################################################
# 1-2 Read Configuration File "dct.cfg"
###################################################
# Read dct.cfg
dct_config = configparser.ConfigParser()
CFG_FILE = os.path.join(DPATH, 'dct.cfg')
if os.path.exists(CFG_FILE):
dct_config.read(CFG_FILE, encoding = 'utf-8')
else:
sys.stderr.write(CFG_FILE + " is not found.\n")
sys.exit(2)
################
# CONFIGURATION
################
# Which Way to Treat Error... (not supported)
error_type = '0'
if (dct_config.has_option('Control','error_type')):
error_type = dct_config.get('Control','error_type')
# Line Numbers of Header Part in Data File
headerlines = 0
if (dct_config.has_option('Control','headerLines')):
headerLines = dct_config.get('Control','headerLines')
# Skip Spectroscopic Treatment
skipXray = '0'
if (dct_config.has_option('Control','skipXray')):
skipXray = dct_config.get('Control','skipXray')
# Skip Axial Treatment
skipAxial = '0'
if (dct_config.has_option('Control','skipAxial')):
skipAxial = dct_config.get('Control','skipAxial')
# Skip Equatorial Treatment
skipEquat = '0'
if (dct_config.has_option('Control','skipEquat')):
skipEquat = dct_config.get('Control','skipEquat')
# Skip Transparency Treatment
skipTrans = '0'
if (dct_config.has_option('Control','skipTrans')):
skipTrans = dct_config.get('Control','skipTrans')
# Skip Contamination Treatment
skipContamination = '1'
if (dct_config.has_option('Control','skipContamination')):
skipContamination = dct_config.get('Control','skipContamination')
# Horizontal Error (deg.)
hrz_error = 0
if (dct_config.has_option('Control','hrz_errror')):
hrz_error = dct_config.getfloat('Control','hrz_error')
#**********************************************
# Configuration for Axial-Divergence Aberration
#**********************************************
# Axial-divergence angles (deg.)
degAxial = 1.25
if (dct_config.has_option('Control','degAxial')):
degAxial = dct_config.getfloat('Control','degAxial')
degAxial2 = 1.18
if (dct_config.has_option('Control','degAxial2')):
degAxial2 = dct_config.getfloat('Control','degAxial2')
# Shape parameter for axial component model
alphaAxial = 0.8
if (dct_config.has_option('Control','alphaAxial')):
alphaAxial = dct_config.getfloat('Control','alphaAxial')
# Margin for Fourier-based Deconvolutional Treatment
marginAxial = 10.0
if (dct_config.has_option('Control','marginAxial')):
marginAxial = dct_config.getfloat('Control','marginAxial')
# Select Mode of Axial-Divergence Treatment
modeAxial = '1'
if (dct_config.has_option('Control','modeAxial')):
modeAxial = dct_config.get('Control','modeAxial')
#*****************************************
# Configuration for Equatorial Aberration
#*****************************************
# Goniometer Radius (mm)
gonioR = 150
if (dct_config.has_option('Control','gonioR')):
gonioR = dct_config.getfloat('Control','gonioR')
# Specimen Width (mm)
specW = 20.0
if (dct_config.has_option('Control','specW')):
specW = dct_config.getfloat('Control','specW')
W_2R = specW / (2*gonioR)
# Divergence Slit Open Angle (deg.)
degDS = 1.25
if (dct_config.has_option('Control','degDS')):
degDS = dct_config.getfloat('Control','degDS')
# View Angle of Si-strip X-ray Detector (deg.)
degLPSD = 4.89
if (dct_config.has_option('Control','degLPSD')):
degLPSD = dct_config.getfloat('Control','degLPSD')
# Specimen Shape, disk or rect (rectangular)
specS = 'rect'
if (dct_config.has_option('Control','specS')):
specS = dct_config.get('Control','specS')
# X-ray beam width (mm)
beamW = 10.0
if (dct_config.has_option('Control','beamW')):
beamW = dct_config.getfloat('Control','beamW')
# Shape parameter for equatorial component model
alphaEquat = 0.5
if (dct_config.has_option('Control','alphaEquat')):
alphaEquat = dct_config.getfloat('Control','alphaEquat')
# Mode of Equatorial Treatment
modeEquat = '1'
if (dct_config.has_option('Control','modeEquat')):
modeEquat = dct_config.get('Control','modeEquat')
# Margin for Fourier-based Deconvolutional Treatment
marginEquat = 10.0
if (dct_config.has_option('Control','marginEquat')):
marginEquat = dct_config.getfloat('Control','marginEquat')
#*****************************************
# Configuration for Transparency Aberration
#*****************************************
# Specimen Thickness
specT = 0.5 # default value is 0.5 mm
if (dct_config.has_option('Control','specT')):
specT = dct_config.getfloat('Control','specT')
# Penetration depth
muInv = 0.1
if (dct_config.has_option('Control','gonioR')):
muInv = dct_config.getfloat('Control','muInv')
# Linear attenuation coefficient
mu = 1/muInv
muR = mu * gonioR
# Margin for Fourier-based Deconvolutional Treatment
marginTrans = 10.0
if (dct_config.has_option('Control','marginTrans')):
marginTrans = dct_config.getfloat('Control','marginTrans')
# Select Mode of Transparency Treatment
modeTrans = '1'
if (dct_config.has_option('Control','modeTrans')):
modeTrans = dct_config.get('Control','modeTrans')
# Shape Parameter for Transparency Treatment
alphaTrans = 0.7
if (dct_config.has_option('Control','alphaTrans')):
alphaTrans = dct_config.getfloat('Control','alphaTrans')
###################################################
# 1-3 Read angle, intensity [ & error ] data file
###################################################
# 1st argument argvs[1] should be data file name
if not os.path.exists(argvs[1]):
sys.stderr.write(argvs[1] + " is not found.\n")
sys.exit(3)
# Read 1st line of data file
f = open(argvs[1], 'r')
line = f.readline()
nSkip = int(headerLines)
# Read data file
source_data = np.loadtxt(argvs[1],dtype='float',delimiter=',',skiprows=nSkip)
# Extract angle, intensity (, error) data
source_angle = source_data[:,0] # 0th column is angle
source_intensity = source_data[:,1] # 1st column is intensity
nSource = len(source_angle)
# If number of columns are equal to or larger than 3,
if (source_data.shape[1] >= 3):
source_error = source_data[:,2] # 2nd column is error
else:
# Error estimation
angular_error = hrz_error # estimated angular error in [deg.]
source_error = np.array(source_intensity) # reserve array
dest_intensity = np.array(source_intensity) # reserve array
dest_error = np.array(source_intensity) # reserve array
# Forward differentiation for the first point
dx = source_angle[1] - source_angle[0]
y1 = source_intensity[0]
y2 = source_intensity[1]
err = (y2 - y1) / dx * angular_error
source_error[0] = np.sqrt(np.maximum(1,y1) + np.square(err))
# Center differentiation for the intermediate points
for index in range(1,nSource-1):
dx = source_angle[index+1] - source_angle[index-1]
y0 = source_intensity[index-1]
y2 = source_intensity[index+1]
err = (y2 - y0) / dx * angular_error
source_error[index] = np.sqrt(np.maximum(1,y1) + np.square(err))
# Backward differentiation for the last point
dx = source_angle[nSource-1] - source_angle[nSource-2]
y0 = source_intensity[nSource-2]
y1 = source_intensity[nSource-1]
err = (y1 - y0) / dx * angular_error
source_error[nSource-1] = np.sqrt(np.maximum(1,y1) + np.square(err))
# Save Copy of Source Data
source_intensity0 = source_intensity.copy() # make a copy
source_error0 = source_error.copy() # make a copy
##############
# Source Data
##############
if (error_type == '0'):
source = [source_angle, source_intensity]
if ((error_type == '1') or (error_type == '2')):
source = [source_angle, source_intensity, source_error]
if (error_type == '3'):
source = [source_angle, source_intensity, source_error, source_error]
# Make a copy
source0 = source
##########################################
# 2. TREAT SPECTROSCOPIC PROFILE OF X-RAY
##########################################
if (skipXray != '1'):
#*****************************
# Treat Spectroscopic Profile
#*****************************
# Deconvolution-Convolution Treatment about Spectroscopic Profile
print("'xray5.py' starts...")
dest = xray.treat(source, DPATH)
print("'xray5.py' finished.")
else:
dest = source
#######################################
# 3. TREAT AXIAL-DIVERGENCE ABERRATION
#######################################
if (skipAxial != '1'):
source = dest # Copy "dest" to "source"
#********************************************************
# 3-1 Calculate Cumulants of Axial-Divergence Aberration
#********************************************************
s_angle = np.array(source[0])
print("'axial5.py' starts...")
print(" degAxial = ",degAxial)
print(" degAxial2 = ",degAxial2)
kappas = axial.calc_kappas(s_angle,degAxial,degAxial2)
#----------------------------------------
# Optional Output for 'skipAxial' == '2'
#----------------------------------------
if (skipAxial == '2'):
CUMULANT_FILE = 'cumulants_axial.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)
#***********************************************
# 3-2 Treat Lower-Angle Part of Axial-Divergnce
#***********************************************
cfg = [alphaAxial,marginAxial,modeAxial,skipAxial,error_type]
dest = axial.treat(-1,source,kappas,cfg)
#*************************************************
# 3-3 Treat Higher-Angle Part of Axial-Divergence
#*************************************************
source = dest # Copy "dest" to "source"
dest = axial.treat(1,source,kappas,cfg)
print("'axial5.py' finished.")
else:
# Copy "source" to "dest"
dest = source
#################################
# 4. TREAT EQUATORIAL ABERRATION
#################################
if ((skipEquat != '1') and (degDS > 0.0)):
source = dest # Copy "dest" to "source"
#**************************************************
# 4-1 Calculate Cumulants of Equatorial Aberration
#**************************************************
s_angle = np.array(source[0])
cfg = [degDS,degLPSD,specW,gonioR,marginEquat]
cfg += [alphaEquat,modeEquat,skipEquat,error_type]
if (modeEquat == '1'):
print("'equatorial1.py' starts...")
equat = equat1
elif (modeEquat == '2'):
print("'equatorial2.py' starts...")
equat = equat2
elif (modeEquat == '3'):
print("'equatorial3.py' starts...")
equat = equat3
print(" specS = ",specS)
print(" degDS = ",degDS)
print(" degLPSD = ",degLPSD)
print(" specW = ",specW)
print(" gonioR = ",gonioR)
if (specS == 'rect'):
kappas = equat.calc_kappas(s_angle,cfg)
elif ((specS == 'disc') or (specS == 'disk')):
nSize = s_angle.size
nGC = 5 # number of sampling points for Gauss-Chebychev
roots,weights,mu = sp.roots_chebyu(nGC,mu=True)
kappas = np.zeros((5,nSize),dtype=np.float64)
for iGC in range(nGC):
yGC1 = beamW * roots[iGC] * 0.5
wGC1 = np.sqrt(specW**2 - 4*yGC1**2)
kappa1 = equat.calc_kappas(s_angle,cfg)
kappas += np.array(kappa1)*weights[iGC]/mu
#---------------------------------------------
# Optional Output for 'skipEquatorial' == '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 (modeEquat == '1'):
dest = equat1.treat(-1,source,kappas,cfg)
dest = equat1.treat(1,source,kappas,cfg)
print("'equatorial1.py' finished.")
elif (modeEquat == '2'):
dest = equat2.treat(-1,source,kappas,cfg)
dest = equat2.treat(1,source,kappas,cfg)
print("'equatorial2.py' finished.")
elif (modeEquat == '3'):
dest = equat3.treat(-1,source,kappas,cfg)
dest = equat3.treat(1,source,kappas,cfg)
print("'equatorial3.py' finished.")
else:
dest = source # Copy "source" to "dest"
##########################################
# 5. TREAT SAMPLE TRANSPARENCY ABERRATION
##########################################
if ((skipTrans != '1') and (muR > 0.0)):
source = dest # Copy "dest" to "source"
if (modeTrans == '1'):
trnspr = trnspr1
print("'trnspr1.py' start...")
elif (modeTrans == '2'):
trnspr = trnspr2
print("'trnspr2.py' start...")
elif (modeTrans == '3'):
trnspr = trnspr3
print("'trnspr3.py' start...")
#***********************************************************
# 5-1 Calculate Cumulants of Sample Transparency Aberration
#***********************************************************
s_angle = np.array(source[0])
cfg=[muInv,gonioR,specT,specW,degDS,alphaTrans,marginTrans,skipTrans,error_type]
print(" muInv = ",muInv)
print(" gonioR = ",gonioR)
print(" specT = ",specT)
print(" specW = ",specW)
print(" degDS = ",degDS)
if (specS == 'rect'):
kappas = trnspr.calc_kappas(s_angle,cfg)
elif ((specS == 'disk')or(specS =='disc')):
nSize = s_angle.size
nGC = 5 # number of sampling points for Gauss-Chebychev
roots,weights,mu = sp.roots_chebyu(nGC,mu=True)
kappas = np.zeros((5,nSize),dtype=np.float64)
for iGC in range(nGC):
yGC1 = beamW * roots[iGC] * 0.5
wGC1 = np.sqrt(specW**2 - 4*yGC1**2)
kappa1 = trnspr.calc_kappas(s_angle,cfg)
kappas += np.array(kappa1)*weights[iGC]/mu
#---------------------------------------
# Optional Output for 'SkipTrans' == '2'
#-----------------------------------------------
if (skipTrans == '2'):
# Cumulants of instrumental function
CUMULANT_FILE = 'cumulants_transparency.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)
#-----------------------------------------------
# Optional Output for 'SkipTransparency' == '2'
#-----------------------------------------------
if (modeTrans == '1'):
# Cumulants of symmetrized function
# cfg=[muInv,gonioR,specT,specW,degDS,alphaTrans,
# marginTrans,skipTrans,error_type]
kappas = trnspr1.calc_sym_kappas(s_angle, cfg)
CUMULANT_FILE = 'cum_untruncated_trans.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 (modeTrans == '1'):
dest = trnspr1.treat(-1,source,kappas,cfg)
source = dest
dest = trnspr1.treat(1,source,kappas,cfg)
print("'trnspr1.py' finished.")
elif (modeTrans == '2'):
dest = trnspr2.treat(source,kappas,cfg)
print("'trnspr2.py' finished.")
elif (modeTrans == '3'):
#**************************************************************
# 5-2 Treat Exponential Part of Sample Transparency Aberration
#**************************************************************
dest = trnspr3.treat(-1,source,kappas,cfg)
#**************************************************************
# 5-3 Treat Rectangular Part of Sample Transparency Aberration
#**************************************************************
source = dest # Copy "dest" to "source"
dest = trnspr3.treat(1,source,kappas,cfg)
print("'trnspr3.py' finished.")
else:
dest = source # Copy "source" to "dest"
#####################
# 6. DECONTAMINATION
#####################
if (skipContamination != '1'):
#*****************************************
# Treat Sub-peaks Caused by Contamination
#*****************************************
# Copy "dest" to "source"
source = dest
dest = cntmn.treat(source, DPATH)
else:
dest = source
################################
# Save Data to Text (.csv) File
################################
save = np.array(dest).T # transpose matrix
# print ("save = ",save)
if (error_type == '0'):
np.savetxt(argvs[2], save, delimiter=",",fmt=["%.3f","%.2f"])
elif ((error_type == '1') or (error_type == '2')):
np.savetxt(argvs[2], save, delimiter=",",fmt=["%.3f","%.2f","%.4e"])
elif (error_type == '3'):
np.savetxt(argvs[2], save, delimiter=",",fmt=["%.3f","%.2f","%.4e","%.4e"])