miniflex7.py

treats “MiniFlex deformation” additional deformation of peak profile observed for data collected with Rigaku MiniFlex 600-C, but not observed for instruments Bruker D2 Phaser or Malvern-Panalytical Aeris.

# -*- coding: utf-8 -*-
######################################################
# miniflex7.py 
# to treat MiniFlex data
# originally coded by T. Ida, Jan. 9, 2025
# updated by Takashi Ida, May 16, 2026 (ver.7.3)
######################################################
import numpy as np # import "numpy" modeule as "np"
import scipy.interpolate as interpolate # import "interpolate" from "scipy"
import scipy.special as special # import "special" from "scipy"
import common5 as cmn # import "common5.py" as "cmn"
import warnings
warnings.filterwarnings('ignore')

#
##############################################
# Define Peak Profile Model for Deconvolution
##############################################
def f_a(x,gam):
    EPS = 1.0e-9
    abs_gam = np.maximum(np.abs(gam),EPS)
#    ans = np.where(x < abs_gam,np.exp(x/abs_gam-1.0)/abs_gam,0.0)
    ans = np.where(x < 0,np.exp(x/abs_gam)/abs_gam,0.0)
    return ans
def f_2a(x,gam,r,dx):
    ans = (1-r)*f_a(x,gam)
    ans += r*f_a(x-dx,gam)
    return ans
def fDeconv(x,gam,r,dx):
    # x: NumPy array
    # gam: decay width
    # r: subpeak relative intensity
    # dx: subpeak relative location
    return f_2a(x,gam,r,dx)

####################################
# miniflex7.py" -> "miniflex.treat"
# TREATMENT FOR MINIFLEX DATA
####################################
def treat(source,cfg):
    # source: Source Data (Python List)
    #    [s_angle, s_intensity],
    #    [s_angle, s_intensity, s_error1], or
    #    [s_angle, s_intensity, s_error1, s_error2]
    # cfg={'error_type':error_type,'dropCenter':dropCenter
    #    'gamMiniFlex':gamMiniFlex,'rhoMiniFlex':rhoMiniFlex,
    #    'delMiniFlex':delMiniFlex,'marginMiniFlex':marginMiniFlex,
    #    'skipMiniFlex':skipMiniFlex}
    #******************************
    # Configuration from "dct.cfg"
    #******************************
    gam = float(cfg['gamMiniFlex']) 
    # -0.035 [deg.]
    r = float(cfg['rhoMiniFlex'])
    # 0.10
    dx = float(cfg['delMiniFlex']) 
    # -0.137 [deg.]
    margin = float(cfg['marginMiniFlex'])
    skipMiniFlex = cfg['skipMiniFlex']
    error_type = cfg['error_type']
    #*****************************************************
    # Step 1-0: Parse Python List "source" to numpy.array
    #*****************************************************
    s_angle = np.array(source[0])
    s_int = np.array(source[1])
    if ((error_type == 1) or (error_type == 3)):
        s_error1 = np.array(source[2])
    if ((error_type == 2) or (error_type == 3)):
        s_error2 = np.array(source[2])
    nData = s_angle.size # number of data points
    #*****************************************************
    # Step 1-1: Scale Transform for MiniFlex
    # f(x)=exp(x-1)/ for x<1
    # s_angle (2Theta in deg.) => source_x
    # s_int => source_y
    #*****************************************************
    #-----------------------------
    # d(chi)/d(2Theta) (dchi_d2T) 
    #-----------------------------
    dchi_d2T = np.tan(s_angle*np.pi/360)
    # dchi_d2T = np.tan(s_angle*np.pi/360)/2 (?)
    # dchi_d2T = np.tan(s_angle*np.pi/360) * 360/np.pi (?)
    #------------------------------------------------------------
    # Integration & creation of abscissa values source_chi
    # Note: chi should be the integral of d(chi)/d(2Theta)
    # by 2Theta (d2Theta).  (-2*ln(cos(Theta)), in this case)
    #------------------------------------------------------------
    source_chi = np.cumsum(dchi_d2T) - dchi_d2T[0]
    source_chi *= (s_angle[-1]-s_angle[0])/(nData-1)
    # source_chi = -2*np.log(np.cos(s_angle*np.pi/360))
    #------------------------------------------------------------
    # Ordinate Correction Factor
    #------------------------------------------------------------
    source_corr = cmn.g_corr(s_angle/2)/dchi_d2T
    source_eta = s_int * source_corr
    #*******************************************
    # Step 1-2: Preparation for Fourier treatment
    #*******************************************
    index = cmn.indices(source_chi, margin)
    nSource = index[0] # number of source data (=nData)
    nValid = index[1] # index of last valid data
    i2 = index[2] # index for first division point in marginal range
    i3 = index[3] # index for second division point in marginal range
    nDest = index[4] # total number of equidistant data
    iMin = index[5] # index at minimum separation in source data
    dx_min = index[6] # minimum separation found in source data
    dchi = index[7] # interval of equidistant data
    chi0 = source_chi[0] # (float)
    chiV = source_chi[-1]
    nDest_2 = (nDest/2).astype('int')
    chi = np.linspace(chi0,chiV,num=nValid,endpoint=False)
    #**************************************************
    # Step 1-3:
    # Interpolation of data onto transformed scale chi
    # (source_chi, source_eta) => (chi, eta)
    #**************************************************
    f3 = interpolate.interp1d(source_chi,source_eta,kind="cubic",fill_value='extrapolate')
    # chi[nValid-1]=np.minimum(chi[nValid-1],source_chi[nData-1])
    eta = f3(chi)
    # eta = np.interp(chi,source_chi,source_eta)
    #**********************************************************
    # Step 1-4:
    # Redimension & Extrapolation for marginal (padding) range
    #**********************************************************
    chi = np.linspace(chi0,chi0+nDest*dchi,num=nDest,endpoint=False)
    eta = np.resize(eta, nDest)
    cmn.pad_margin(eta, index) # => "cmn_common.py"
    #***********************************************
    # Step 1-5: Fourier transform of intensity data
    # (eta) => (ft_eta)
    #***********************************************
    ft_eta = np.fft.rfft(eta) # Fourier transform for real data
    # nFT = ft_eta.size
    # dxi = 1.0/(dchi*nDest)
    # Abscissa of Fourier transform
    # xi = np.linspace(0,nFT*dxi,num=nDest,endpoint=False)
    # nDest_2 = (nDest/2).astype('int')
    # xi[nDest_2:nDest] = xi[nDest_2:nDest] - nDest * dxi
    #********************************************************
    # Step 1-6: Fourier transform of approximate model profile
    # => (ft_deconv), (ft_conv)
    #********************************************************
    chi = np.linspace(0,nDest*dchi,num=nDest,endpoint=False)
    deconv = np.empty(nDest)
    deconv[1:nDest_2] = fDeconv(chi[1:nDest_2],gam,r,dx)
    deconv[nDest_2:nDest] = fDeconv(chi[nDest_2:nDest]-nDest*dchi,gam,r,dx)
    deconv[0] = 1.0/dchi - np.sum(deconv[1:nDest])
    #---------------------------------------------
    # Optional Output for 'skipMiniFlex' == '2'
    #---------------------------------------------
    if (skipMiniFlex == 2)or(skipMiniFlex == '2'):
        file_name = 'deconv_miniflex.csv'
        save_deconv_miniflex = np.array([chi,deconv]).T
        fmt_list = ["%.5e","%.5e"]
        np.savetxt(file_name, save_deconv_miniflex, delimiter=",",fmt=fmt_list)
    ft_deconv = np.fft.rfft(deconv)
    ft_conv = np.abs(ft_deconv)
    #***************************************************
    # Step 1-7: Fourier transformed instrumental function
    #***************************************************
    # np.seterr(invalid='ignore')
    np.seterr(invalid='warn')
    ft_inst = ft_deconv / ft_conv # deconvolutional
    #***********************************************
    # Step 1-8: Deconvolution/convolution treatment
    #***********************************************
    ft_eta = ft_eta / ft_inst # deconvolutional
    #*************************************
    # Step 1-9: Inverse Fourier transform
    #*************************************
    zeta = np.real(np.fft.irfft(ft_eta)) # inverse Fourier transform
    #******************************
    # Step 1-10: Pick up intensity
    #******************************
    chi = np.linspace(chi0,chi0+nDest*dchi,num=nDest,endpoint = False)
    # f3 = interpolate.interp1d(chi,zeta,kind="cubic")
    # dest_zeta = f3(source_chi)
    dest_zeta = np.interp(source_chi,chi,zeta) # linear interpolation
    #******************************
    # Step 1-11: Rescale transform
    # dest_y => d_int
    #******************************
    d_int = dest_zeta/source_corr
    dest = source
    dest[1] = d_int
    #**********************************************************
    # Step 1-12: Fourier treatment of error data (not validated)
    #**********************************************************
    #if (error_type != 0):
    #    dest = cmn.TreatErrors(error_type,source,source_chi,index,source_corr,ft_inst)
    # print('d_int.shape = ',d_int.shape)
    #
    return dest