exterm7.py

Main routine of exterm7.3 system for deconvolutional treatment of XRD data

Source code:

# -*- coding: utf-8 -*-

# exterm7.py
# Coded by Takashi Ida on September 8, 2021
# Updated by Takashi Ida on January 9, 2025 (version 6.0)
# Updated by Takashi Ida on September 24, 2025 (version 6.7)
# Updated by Takashi Ida on Jun 4, 2026 (version 7.3)

import numpy as np # import "numpy" module as "np"
import scipy.special as sp # import "scipy.special" as "sp"
import scipy.interpolate as interpolate
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 json #

import common5 as cmn # import "common5.py" as cmn
import xray6 as xray # import "xray5.py" as xray
import axial6 as axial # import "axial6.py" as axial
import equat_trnspr7 as eq_tr # import "equat-transpr.py" as eq_tr
import cntmn5 as cntmn  # import 'cntmn5.py' for decontamination
import miniflex7 as miniflex

#################
# 1. PREPARATION
#################

###################################
# 1-1 Treat Command-Line Arguments
###################################
#*******************************************************
# Treatment of command-line argument
# argvs[0] : execute file name
# argvs[1] : input data file name
# argvs[2] : output data file name
# argvs[3] : configuration file name
# argvs[4] : equatorial-transparency cumulant file name
#*******************************************************
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
    # argvs[3] is configuration file name
    # argvs[4] is output cumulant file name 
# DPATH = os.path.dirname(sys.argv[0]) # directory of the EXE file
DPATH = ''
if (argc == 1):    # in case no argument
    argvs = argvs + [os.path.join(DPATH, '00raw.csv')]
    argc = 2
if (argc == 2): # in case less than 2 argument
    argvs = argvs + [os.path.join(DPATH, '01dct.csv')]
    argc = 3
if (argc == 3): # in case less than 3 arguments
    argvs = argvs + [os.path.join(DPATH, 'dct.cfg')]
    argc = 4
if (argc == 4): # in case less than 3 arguments
    argvs = argvs + [os.path.join(DPATH, 'cum_eq_tr.csv')]
    argc = 5
print('**********************************')
print('Input file name: ',argvs[1])
print('Output file name: ',argvs[2])
print('Main configuration filename: ',argvs[3])
print('Cumulant file name: ',argvs[4])
print('**********************************')
cum_file_name = argvs[4]
###################################################
# 1-2 Read Configuration File "dct.cfg" 
###################################################
# Read dct.cfg
dct_config = configparser.ConfigParser()
CFG_FILE = argvs[3]
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
################
# Delimiter to read data
in_delim = ','
if (dct_config.has_option('Control','in_delim')):
    in_delim = dct_config.get('Control','in_delim')
# Delimiter to write data
out_delim = ','
if (dct_config.has_option('Control','out_delim')):
    out_delim = dct_config.get('Control','out_delim')

# Which Way to Treat Error... (not supported)
error_type = 0
if (dct_config.has_option('Control','error_type')):
    error_type = dct_config.getint('Control','error_type')
# Line Numbers of Header Part in Data File
headerlines = 0
if (dct_config.has_option('Control','headerLines')):
    headerLines = dct_config.getint('Control','headerLines')
# Skip Spectroscopic Treatment
skipXray = 0
if (dct_config.has_option('Control','skipXray')):
    skipXray = dct_config.getint('Control','skipXray')
# Skip Axial Treatment
skipAxial = 0
if (dct_config.has_option('Control','skipAxial')):
    skipAxial = dct_config.getint('Control','skipAxial')
# Skip Equatorial Treatment
skipEquat = 2
if (dct_config.has_option('Control','skipEquat')):
    skipEquat = dct_config.getint('Control','skipEquat')
# Skip Transparency Treatment
skipTrans = 2
if (dct_config.has_option('Control','skipTrans')):
    skipTrans = dct_config.getint('Control','skipTrans')
# Skip MiniFlex Treatment
skipMiniFlex = 0
if (dct_config.has_option('Control','skipMiniFlex')):
    skipMiniFlex = dct_config.getint('Control','skipMiniFlex')
# Skip Contamination Treatment
skipContamination = 1
if (dct_config.has_option('Control','skipContamination')):
    skipContamination = dct_config.getint('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 Spectroscopic Profile 
#***************************************************
# Margin for Fourier-based Deconvolutional Treatment
#***************************************************
marginXray = 10.0
if (dct_config.has_option('Control','marginXray')):
    marginXray = dct_config.getfloat('Control','marginXray')
#**********************************************
# 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.getint('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 = 0.625
if (dct_config.has_option('Control','degDS')):
    degDS = dct_config.getfloat('Control','degDS')
print(" degDS = ",degDS)
# Divergence Slit Twist Angle (deg.)
degDS2 = 0.0
if (dct_config.has_option('Control','degDS2')):
    degDS2 = dct_config.getfloat('Control','degDS2')
# Relative Y-coorinate to Correct Twisted Incident beam
relY = 0.0
if (dct_config.has_option('Control','relY')):
    relY = dct_config.getfloat('Control','relY')
# View Angle of Si-strip X-ray Detector (deg.)
degSSD = 4.89
if (dct_config.has_option('Control','degSSD')):
    degSSD = dct_config.getfloat('Control','degSSD')
# 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.getint('Control','modeEquat')
# Margin for Fourier-based Deconvolutional Treatment
marginEquat = 10.0
if (dct_config.has_option('Control','marginEquat')):
    marginEquat = dct_config.getfloat('Control','marginEquat')
# Step Interval for Coarse Sampling Step (deg.) for
# Equatorial-Transparency-Displacement Aberration
stepEquat = 1.0
if (dct_config.has_option('Control','stepEquat')):
    stepEquat = dct_config.getfloat('Control','stepEquat')
# Specimen Thickness
specT = 0.5 # default value is 0.5 mm
if (dct_config.has_option('Control','specT')):
    specT = dct_config.getfloat('Control','specT')
# Specimen Displacement
specD = 0.0 # default value is 0.00 mm
if (dct_config.has_option('Control','specD')):
    specD = dct_config.getfloat('Control','specD')
# Penetration depth os sample
muInv = 0.1
if (dct_config.has_option('Control','muInv')):
    muInv = dct_config.getfloat('Control','muInv')
# Penetration depth os sample holder
muInv2 = 0.1
if (dct_config.has_option('Control','muInv2')):
    muInv2 = dct_config.getfloat('Control','muInv2')
# Linear attenuation coefficient of sample
mu = 1/muInv
muR = mu * gonioR
# Number of sampling points (Z) for Transparency Treatment
nZ = 40
if (dct_config.has_option('Control','nZ')):
    nZ = dct_config.getint('Control','nZ')
# Number of sampling points (psi) for Equatorial Treatment
nPsi= 10
if (dct_config.has_option('Control','nPsi')):
    nPsi = dct_config.getint('Control','nPsi')
# Number of sampling points (phi) for Equatorial Treatment
nPhi = 10
if (dct_config.has_option('Control','nPhi')):
    nPhi = dct_config.getint('Control','nPhi')
# X-ray source displacement (xrayD) (not supported)
xrayD = 0.0
if (dct_config.has_option('Control','xrayD')):
    xrayD = dct_config.getfloat('Control','xrayD')
# Singular angle
if (dct_config.has_option('Control','sing2T')):
    treat_singular = True
    sing2T = dct_config.getfloat('Control','sing2T')
else:
    treat_singular = False
# Singular angles list
sing2T_list = [6.0]
if (dct_config.has_option('Control','sing2T_list')):
    treat_singular_list = True
    sing2T_list = json.loads(dct_config.get('Control','sing2T_list'))
else:
    treat_singular_list = False
# Inhomogeneous sensitivity
dropCenter = 1.0
if (dct_config.has_option('Control','dropCenter')):
    dropCenter = dct_config.getfloat('Control','dropCenter')
#*************************
# Parameters for MiniFlex
#*************************
# gamMiniFlex = -0.070
gamMiniFlex = -0.035 # decay width (deg.)
if (dct_config.has_option('Control','gamMiniFlex')):
    gamMiniFlex = dct_config.getfloat('Control','gamMiniFlex')
# rhoMiniFlex:  relative intensity of sub-peak for MiniFlex
rhoMiniFlex = 0.10
if (dct_config.has_option('Control','rhoMiniFlex')):
    rhoMiniFlex = dct_config.getfloat('Control','rhoMiniFlex')
# delMiniFlex: relative location of sub-peak for MiniFlex
# delMiniFlex = -0.137 # relative position of sub-peak (deg.)
delMiniFlex = -0.0258
if (dct_config.has_option('Control','delMiniFlex')):
    delMiniFlex = dct_config.getfloat('Control','delMiniFlex')
marginMiniFlex = 10 # margin for deconvolutional treatment
if (dct_config.has_option('Control','marginMiniFlex')):
    marginMiniFlex = dct_config.getfloat('Control','marginMiniFlex')
skipMiniFlex = 0
if (dct_config.has_option('Control','skipMiniFlex')):
    skipMiniFlex = dct_config.getint('Control','skipMiniFlex')
modeMiniFlex = 1
if (dct_config.has_option('Control','modeMiniFlex')):
    modeMiniFlex = dct_config.getint('Control','modeMiniFlex')

###################################################
# 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
if (in_delim == 'None'):
    source_data = np.loadtxt(argvs[1],dtype='float',skiprows=nSkip)
else:
    source_data = np.loadtxt(argvs[1],dtype='float',delimiter=in_delim,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) or (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
##########################################
print(' skipXray = ',skipXray)
if (skipXray != 1) and (skipXray != '1'):
    #*****************************
    # Treat Spectroscopic Profile
    #*****************************
    # Deconvolution-Convolution Treatment about Spectroscopic Profile
    cfg={'marginXray':marginXray,'error_type':error_type}
    print("**************************************************")
    print("Treatment of Spectroscopic profile of source X-ray")
    print("'xray6.py' starts...")
    dest = xray.treat(source, DPATH, cfg)
    print("'xray6.py' finished.")
    print("**************************************************")
else:
    dest = source

#######################################
# 3. TREAT AXIAL-DIVERGENCE ABERRATION 
#######################################
print('skipAxial = ',skipAxial)
if (skipAxial != 1) and (skipAxial != '1'):
    source = dest    # Copy "dest" to "source" 
    #********************************************************
    # 3-1 Calculate Cumulants of Axial-Divergence Aberration
    #********************************************************
    s_angle = np.array(source[0])
    print("****************************************")
    print("Treatment of Axial-Divergence Abberation")
    print("'axial6.py' starts...")
    print(" degAxial = ",degAxial)
    print(" degAxial2 = ",degAxial2)
    kappas = axial.calc_kappas(s_angle,degAxial,degAxial2)
    #----------------------------------------
    # Optional Output for 'skipAxial' == '2'
    #----------------------------------------
    if (skipAxial == 2)or(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('cum_eq_tr.csv', save_cumulant, delimiter=",",fmt=fmt_list)
    #***********************************************
    # 3-2 Treat Lower-Angle Part of Axial-Divergnce 
    #***********************************************
    cfg={'alphaAxial':alphaAxial,'marginAxial':marginAxial,\
        'modeAxial':modeAxial,'skipAxial':skipAxial,\
        'error_type':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("'axial6.py' finished.")
    print("****************************************")
else:
    # Copy "source" to "dest"
    dest = source

################################################
# 4. TREAT EQUATORIAL & TRANSPARENCY ABERRATION 
################################################
print('skipEquat = ',skipEquat)
if ((skipEquat != 1) and (skipEquat != '1')):
    source = dest    # Copy "dest" to "source" 
    #************************************
    # 4-1 Calculate Cumulants of 
    # Equatorial-Transparency Aberration
    #************************************
    s_angle = np.array(source[0]) # source angle
    # Set configuration about equatorial-transparency aberration
    cfg={'gonioR':gonioR,'xrayD':xrayD,'degDS':degDS,'degSSD':degSSD,\
        'specW':specW,'specT':specT,'specD':specD,'muInv':muInv,'muInv2':muInv2,\
        'nZ':nZ,'nPhi':nPhi,'nPsi':nPsi,\
        'alphaEquat':alphaEquat,'marginEquat':marginEquat,\
        'skipEquat':skipEquat,'error_type':error_type,\
        'stepEquat':stepEquat,'degDS2':degDS2,'relY':relY,\
        'dropCenter':dropCenter}
    print("************************************************")
    print("Treatment of Equatorial-Transparency Aberration")
    print("'equat_trnspr7.py' cumulant evaluation starts...")
    print(" gonioR = ",gonioR)
    print(" degDS, degSSD = ",degDS,degSSD)
    print(" specS, specW ,specT, specD = ",specS,specW,specT,specD)
    print(" muInv, muInv2 = ",muInv,muInv2)
    print(" alphaEquat,marginEquat = ",alphaEquat,marginEquat)
    print(" nZ, nPhi, nPsi = ",nZ,nPhi,nPsi)
    print(" stepEquat,degDS2, relY = ",stepEquat,degDS2,relY)
    print(" dropCenter = ",dropCenter)
    if (treat_singular):
        print(" sing2T = ",sing2T)
    if (treat_singular_list):
        print(' sing2T_list = ',sing2T_list)
    # Calculate cumulants at coarse (1 deg.) step, 
    # then interpolate it 
    def calc_kappas_eco(s_angle,sing2T_list): # Economized Caluculation
        stepEquat = float(cfg['stepEquat'])
        s_angle_1_L = np.floor(s_angle[0]/stepEquat)*stepEquat
        s_angle_1_R = (np.ceil(s_angle[-1]/stepEquat)+1)*stepEquat # ?
        s_angle_1 = np.arange(s_angle_1_L,s_angle_1_R,stepEquat)
        # Calculate cumulants for Coarse Step (stepEquat)
        kappas_1 = eq_tr.calc_kappas(s_angle_1,cfg)
        s0_1,k1_1,k2_1 = kappas_1[0],kappas_1[1],kappas_1[2]
        k3_1,k4_1 = kappas_1[3],kappas_1[4]
        if (treat_singular):
            # Treat "singular angle" ???
            iSing = int(np.ceil((sing2T-s_angle_1_L)/stepEquat)-1)
            s0_1[iSing] = (s0_1[iSing-1] + s0_1[iSing+1])/2
            k1_1[iSing] = (k1_1[iSing-1] + k1_1[iSing+1])/2
            k2_1[iSing] = (k2_1[iSing-1] + k2_1[iSing+1])/2
            k3_1[iSing] = (k3_1[iSing-1] + k3_1[iSing+1])/2
            k4_1[iSing] = (k4_1[iSing-1] + k4_1[iSing+1])/2
            print(f's0_1[{iSing-1:d}] = {s0_1[iSing-1]:f}')
            print(f's0_1[{iSing+1:d}] = {s0_1[iSing+1]:f}')
            print(f's0_1[{iSing:d}] = {s0_1[iSing]:f}')
        if (treat_singular_list):
            # sing2T_list = map(float,sing2T_list.split(','))
            print(' sing2T_list  = ',sing2T_list)
            sing2T_array = np.array(sing2T_list)
            print(' type of sing2T_array : ',type(sing2T_array))
            print(' shape of sing2T_array : ',sing2T_array.shape)
            sing2T_array = np.array(sing2T_array)
            print(' sing2T_array = ',sing2T_array)
            # Decending sort...
            sing2T_array = np.sort(sing2T_array)[::-1]
            for item in sing2T_array:
                iSing = int(np.ceil((item-s_angle_1_L)/stepEquat)-1)
                s0_1[iSing] = s0_1[iSing+1]
                k1_1[iSing] = k1_1[iSing+1]
                k2_1[iSing] = k2_1[iSing+1]
                k3_1[iSing] = k3_1[iSing+1]
                k4_1[iSing] = k4_1[iSing+1]
        s0 = np.interp(s_angle,s_angle_1,s0_1)
        k1 = np.interp(s_angle,s_angle_1,k1_1)
        k2 = np.interp(s_angle,s_angle_1,k2_1)
        k3 = np.interp(s_angle,s_angle_1,k3_1)
        k4 = np.interp(s_angle,s_angle_1,k4_1)
        kappas = np.array([s0,k1,k2,k3,k4])
        return kappas
    if (specS == 'rect'):
        # Call a rootine to calculate cumulants ...
        # kappas = eq_tr.calc_kappas(s_angle,cfg)
        kappas = calc_kappas_eco(s_angle,sing2T_list)
    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)
        # Initialize object for cumulants ...
        kappas = np.empty((5,nSize),dtype=np.float64)
        for iGC in range(nGC):
            yGC1 = beamW * roots[iGC] * 0.5
            wGC1 = np.sqrt(specW**2 - 4*yGC1**2)
            # Call a rootine to calculate cumulants ...
            kappa1 = calc_kappas_eco(s_angle)
            kappas += np.array(kappa1)*weights[iGC]/mu
    #
    print("'equat_trnspr7.py' cumulant evaluation finished.")
    print("************************************************")
    #
    print("'equat_trnspr7.py' deconvolutional treatment starts...")
    #---------------------------------------------
    # Optional Output for 'skipEquat' == '2'
    #---------------------------------------------
    if (skipEquat == 2)or(skipEquat == '2'):
        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(cum_file_name, save_cumulant, delimiter=",",fmt=fmt_list)
    cfg={'gonioR':gonioR,'degDS':degDS,'degSSD':degSSD,\
        'specW':specW,'specT':specT,'muInv':muInv,'muInv2':muInv2,\
        'nZ':nZ,'nPhi':nPhi,'nPsi':nPsi,\
        'alphaEquat':alphaEquat,'marginEquat':marginEquat,\
        'stepEquat':stepEquat,\
        'skipEquat':skipEquat,'error_type':error_type,\
        'dropCenter':dropCenter}
    if (modeEquat == 1):
        #*******************************************************
        # 4-2 Treat Equatorial-Transparency Aberration
        #*******************************************************
        dest = eq_tr.treat(-1,source,kappas,cfg)
        #*************************************************
        # 4-3 Treat Higher-Angle Part of Axial-Divergence 
        #*************************************************
        # dest = eq_tr6.treat(1,source,kappas,cfg)
        print("'equat_trnspr7.py' finished.")
        source = dest    # Copy "dest" to "source" 
else:
    dest = source # Copy "source" to "dest"

################################
# 5. TREAT MINIFLEX DEFORMATION
################################
print('skipMiniFlex = ',skipMiniFlex)
if ((skipMiniFlex != 1) and (skipMiniFlex != '1')):
    source = dest    # Copy "dest" to "source" 
    s_angle = np.array(source[0]) # source angle
    # Set configuration about equatorial-transparency aberration
    cfg={'error_type':error_type,
        'skipMiniFlex':skipMiniFlex,
        'dropCenter':dropCenter,'gamMiniFlex':gamMiniFlex,
        'rhoMiniFlex':rhoMiniFlex,'delMiniFlex':delMiniFlex,
        'modeMiniFlex':modeMiniFlex,'error_type':error_type,
        'marginMiniFlex':marginMiniFlex}
    print("************************************************")
    print("Treatment of MiniFlex Deformation")
    print("'miniflex7.py' starts...")
    # print(" dropCenter = ",dropCenter)
    print(" gamMiniFlex = ",gamMiniFlex)
    print(" rhoMiniFlex = ",rhoMiniFlex)
    print(" delMiniFlex = ",delMiniFlex)
    print(" marginMiniFlex = ",marginMiniFlex)
    if (skipMiniFlex != 1) and (skipMiniFlex != '1'):
        dest = miniflex.treat(source,cfg)
    dest = source # Copy "source" to "dest"

#####################
# 6. DECONTAMINATION 
#####################
if (skipContamination != '1')and(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) or (error_type == '0'):
    if (out_delim == 'None'):
        np.savetxt(argvs[2], save, delimiter=' ',fmt=["%.3f","%.2f"])
    else:
        np.savetxt(argvs[2], save, delimiter=out_delim,fmt=["%.3f","%.2f"])
elif ((error_type == '1') or (error_type == '2')):
    if (delim == 'None'):
        np.savetxt(argvs[2], save, delimiter=' ',fmt=["%.3f","%.2f","%.4e"])
    else:
        np.savetxt(argvs[2], save, delimiter=out_delim,fmt=["%.3f","%.2f","%.4e"])
elif (error_type == '3'):
    if (delim == 'None'):
        np.savetxt(argvs[2], save, delimiter=' ',fmt=["%.3f","%.2f","%.4e","%.4e"])
    else:
        np.savetxt(argvs[2], save, delimiter=out_delim,fmt=["%.3f","%.2f","%.4e","%.4e"])
print("Save data finished.")