treats equatorial aberration and sample transparency effects
# -*- coding: utf-8 -*-
######################################################
# equat_trnspr7.py
# to treat Equatorial-Transparency Aberration
# by a two-step deconvolutioal method
# originally coded by T. Ida, Jan. 9, 2025
# updated by Takashi Ida, Aprril 2, 2025 (ver.6.2)
# updated by Takashi Ida, June 20, 2025 (ver.6.3)
# updated by Takashi Ida, June 26, 2025 (ver.6.4)
# updated by Takashi Ida, August 21, 2025 (ver.6.5)
# updated by Takashi Ida, December 23, 2025 (ver.6.9)
# updated by Takashi Ida, May 18, 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"
############################################################
# calc_kappas()
# Calculate Cumulants of Equatorial-Transparency Aberration
############################################################
def sign(x):
return np.sign(x)
def abs(x):
return np.abs(x)
def expc(x): # 1-exp(-x)=x-x**2/2+x**3/6-...
return np.where(abs(x)<1E-8,x*(1-x/2*(1-x/3*(1-x/4))),1-np.exp(-x))
#
def calc_kappas(degTwot, cfg):
# degTwot: 2Theta in [deg.] (NumPy array)
# cfg={'gonioR':gonioR,'xrayD',xrayD,'degDS':degDS,'degLPSD':degLPSD,\
# '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}
#
# Goniometer radius
gonioR = float(cfg['gonioR'])
# Displacement of X-ray source (not used)
xrayD = float(cfg['xrayD'])
# Divergence slit
# print(" cfg['degDS'] = ",cfg['degDS'])
degDS = float(cfg['degDS'])
# Detector view angle
degSSD = float(cfg['degSSD'])
# Penetration depth of powder & holder
muInv,muInv2 = float(cfg['muInv']),float(cfg['muInv2'])
# Specimen thickness, width, displacement...
specT,specW,specD = float(cfg['specT']),float(cfg['specW']),float(cfg['specD'])
# Sampling numbers for numerical integral...
nU,nPhi,nPsi = int(cfg['nZ']),int(cfg['nPhi']),int(cfg['nPsi'])
# Coarse step to calculate cumulants of
# Equatorial,Sample Transparency, and Sample Displacement
# Aberrations
stepEquat = float(cfg['stepEquat'])
# Parameters for twisted (x-ray or divergence slit)
degDS2 = float(cfg['degDS2'])
relY = float(cfg['relY'])
# Parameter for inhomogeneous sensitivity of SSXD
dropCenter = float(cfg['dropCenter'])
#
# print(' gonioR,specT,specW,specD = ',gonioR,specT,specW,specD)
# print(' degDS,degSSD,muInv,muInv2 = ',degDS,degSSD,muInv,muInv2)
print(' stepEquat,nU,nPhi,nPsi = ',stepEquat,nU,nPhi,nPsi)
#
t_R,w_2R = specT/gonioR,0.5*specW/gonioR
mu,muR,muT = 1/muInv,gonioR/muInv,specT/muInv
mu2,muR2,muT2 = 1/muInv2,gonioR/muInv2,specT/muInv2
phiDS = degDS*np.pi/180
psiSSD = degSSD*np.pi/360
phiDS2 = degDS2*np.pi/180
nSize = degTwot.size
#
ones_size = np.ones(nSize) # (nSize,)
ones_Psi = np.ones(nPsi) # (nPsi,)
ones_U = np.ones(nU) # (nU,)
ones_Psi_U = np.ones((nPsi,nU))
ones_Phi = np.ones(nPhi) # (nPhi,)
ones_U_Phi = np.ones((nU,nPhi))
ones_Psi_U_Phi = np.ones((nPsi,nU,nPhi))
#
theta = degTwot*np.pi/360 # theta in [rad.] (nSize,)
cosT = np.cos(theta) # cos(theta) (nSize,)
sinT = np.sin(theta) # sin(theta) (nSize,)
tanT = sinT/cosT # tan(theta) (nSize,)
theta_Psi_U_Phi = np.outer(theta,ones_Psi_U_Phi).reshape(nSize,nPsi,nU,nPhi)
sinT_Psi_U_Phi = np.outer(sinT,ones_Psi_U_Phi).reshape(nSize,nPsi,nU,nPhi)
################################
# PREPARE INTEGRATION ABOUT Psi
################################
# Mid-point method for numerical integral...
# Sampling points:
# x_Psi=[-0.95,-0.85,...,-0.05,0.05,...,0.75,0.95] for nPsi=20
x_Psi = np.linspace(-1+1/nPsi,1-1/nPsi,num=nPsi,endpoint=True)
# Weights for each sampling point:
# w_Psi = np.full_like(x_Psi,1/nPsi)
# Inhomogeneous sensitivity of detector strips
# Weights for each sampling point:
height = dropCenter
tangent = height**2/2
w_Psi = np.maximum(0.0,height-tangent*(x_Psi+1))
w_Psi += np.maximum(0.0,height-tangent*(1-x_Psi))
w_Psi *= 1.0/nPsi # ????
# Lower and Upper limits for psi
psi_L,psi_U = -psiSSD/2,psiSSD/2 # (float)
print(' height,tangent = ',height,tangent)
# print(' x_Psi = ',x_Psi)
# print(' w_Psi = ',w_Psi)
# psi = (psi_L*(1-x_Psi)+psi_U*(1+x_Psi)) # ? (nPsi,)
psi = (psi_L*(1-x_Psi)+psi_U*(1+x_Psi))/2 # (nPsi,)
# Goniometer angle (nSize,nPsi) thetaG = theta + psi
thetaG = np.outer(theta,ones_Psi).reshape(nSize,nPsi)
thetaG -= np.outer(ones_size,psi).reshape(nSize,nPsi)
sinTG,cosTG = np.sin(thetaG),np.cos(thetaG) # (nSize,nPsi)
##################################
# PREPARE INTEGRATION ABOUT Z
# u = exp(2 mu Z / sin(theta))
##################################
# Mid-point method along Z
x_U = np.linspace(-1+1/nU,1-1/nU,num=nU,endpoint=True) # (nU,)
w_U = np.full_like(x_U,1/nU) # (nU,)
exp_z_L = np.exp(-2*muT/sinT) # (nSize,)
u_L = exp_z_L # Lower Limit of u (nSize,)
u_U = np.ones(nSize,) # Upper Limit of u (nSize,)
ones_U = np.ones(nU,dtype='float') # (nU,)
# print('-------------------------')
# print(' Shape of u_L : ',u_L.shape)
# print(' Shape of u_U : ',u_U.shape)
# print(' Shape of ones_U : ',ones_U.shape)
# print(' Shape of x_U : ',x_U.shape)
# print('-------------------------')
u_L_Psi = np.outer(u_L,ones_Psi).reshape(nSize,nPsi)
u_U_Psi = np.outer(u_U,ones_Psi).reshape(nSize,nPsi)
u = np.outer(u_L_Psi,(1-x_U))
u += np.outer(u_U_Psi,(1+x_U))
u /= 2 # (nSize,nPsi,nU)
# u = u.reshape(nSize,nPsi,nU)
# u = np.outer(u,ones_Phi).reshape(nSize,nPsi,nU,nPhi)
###################################
# Z-coordinate of reflection point
# z = ΔS + (sinΘ / 2μ) ln(u)
###################################
# sinT_array = np.outer(sinT,ones_Psi_U_Phi).reshape(nSize,nPsi,nU,nPhi)
sinT_array = np.outer(sinT,ones_Psi_U).reshape(nSize,nPsi,nU)
u = u.reshape(nSize,nPsi,nU)
z = sinT_array * np.log(u) # (nSize,nPsi,nU)
z /= np.full_like(z,2*mu).reshape(nSize,nPsi,nU)
z += specD # add specimen displacement (nSize,nPsi,nU)
################################
# PREPARE INTEGRATION ABOUT phi
################################
# Mid-point method is used ...
x_Phi = np.linspace(-1+1/nPhi,1-1/nPhi,num=nPhi,endpoint=True)
# sampling location (nPhi,)
# w_Phi = np.full_like(x_Phi,0.5/nPhi) # (nPhi,)
w_Phi = np.full_like(x_Phi,1/nPhi) # (nPhi,)
# Lower limit of effective phi... (nSize,nPsi,nU)
# print('specD = ',specD)
# print('specT = ',specT)
tanT_L = sinTG - (specD-specT)/gonioR # (nSize,nPsi)
tanT_L /= cosTG - w_2R # (nSize,nPsi)
tanT_L = tanT_L.reshape(nSize,nPsi)
# print('thetaG.shape = ',thetaG.shape)
# print('tanT_L.shape = ',tanT_L.shape)
# print('thetaG.dtype = ',thetaG.dtype)
# print('tanT_L.dtype = ',tanT_L.dtype)
phi_L1 = thetaG-np.arctan(tanT_L) # (nSize,nPsi)
# print('type(phiDS) : ',type(phiDS))
# print('type(phiDS2) : ',type(phiDS2))
# print('type(degDS2) : ',type(degDS2))
phi_L2 = np.full_like(phi_L1,-phiDS/2+phiDS2) # (nSize,nPsi)
phi_L = np.maximum(phi_L1,phi_L2) # (nSize,nPsi)
phi_L = np.outer(phi_L,ones_U).reshape(nSize,nPsi,nU)
# Upper limit of effective phi...
tanT_U = sinTG - specD/gonioR
tanT_U /= cosTG + w_2R
phi_U1 = thetaG-np.arctan(tanT_U) # (nSize,nPsi)
phi_U2 = np.full_like(phi_U1,phiDS/2+phiDS2) # (nSize,nPsi)
phi_U = np.minimum(phi_U1,phi_U2) # (nSize,nPsi)
phi_U = np.outer(phi_U,ones_U).reshape(nSize,nPsi,nU)
# Equatorial deviation angle phi (nSize,nPsi,nU,nPhi)
phi = np.outer(phi_L,1-x_Phi)/2 # (nSize,nPsi,nU,nPhi)
phi += np.outer(phi_U,1+x_Phi)/2 # (nSize,nPsi,nU,nPhi)
phi = phi.reshape(nSize,nPsi,nU,nPhi)
#####
# Incident glancing angle ...
theta_I = np.outer(thetaG,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
theta_I -= phi
sinT_I = np.sin(theta_I) # (nSize,nPsi,nU,nPhi)
cosT_I = np.cos(theta_I) # (nSize,nPsi,nU,nPhi)
tanT_I = sinT_I/cosT_I # (nSize,nPsi,nU,nPhi)
#####
# X-coordinate of incidence point...
specD_gonioR = np.full_like(sinTG,specD/gonioR)
x_I = np.outer(sinTG-specD_gonioR,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
x_I = x_I/tanT_I
x_I -= np.outer(cosTG,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
x_I *= gonioR # (nSize,nPsi,nU,nPhi)
#####
# Z-coordinate of reflection point...
z_R = np.outer(z,ones_Phi).reshape(nSize,nPsi,nU,nPhi)
#####
# X-coordinate of reflection point...
x_R = x_I+(specD-z_R)/tanT_I # (nSize,nPsi,nU,nPhi)
#####
# X & Z-coordinates of detection point...
tan2psi = np.tan(2*psi) # -> (nPsi,)
tan2psi = np.outer(ones_size,tan2psi) # (nSize,nPsi)
x_D = gonioR*(cosTG-sinTG*tan2psi) # (nSize,nPsi)
x_D = np.outer(x_D,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
z_D = gonioR*(sinTG+cosTG*tan2psi) # (nSize,nPsi)
z_D = np.outer(z_D,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
#####
# Emission glancing angle...
tanT_E = (z_D-z_R)/(x_D-x_R) # (nSize,nPsi,nU,nPhi)
theta_E = np.arctan(tanT_E) # (nSize,nPsi,nU,nPhi)
sinT_E = np.sin(theta_E) # (nSize,nPsi,nU,nPhi)
cosT_E = np.cos(theta_E) # (nSize,nPsi,nU,nPhi)
#####
# X-coordinate of Emission point...
x_E = x_R + (specD-z_R)/tanT_E
#####
# Deviation angle about 2Theta...
del2T = np.outer(2*theta,ones_Psi_U_Phi).reshape(nSize,nPsi,nU,nPhi)
del2T -= theta_I + theta_E # (nSize,nPsi,nU,nPhi)
#
#####################
# RELATIVE INTENSITY
#####################
# zeros = np.zeros((nSize,nPsi,nU,nPhi),dtype='float') # !
zeros = np.full_like(del2T,1E-6) # !
#
phiDS_array=np.full_like(x_I,phiDS)
#####
# Incident path length in holder
l2_I = ((-specW/2-x_I)/cosT_I).reshape(nSize,nPsi,nU,nPhi)
l2_I = np.where(-specW/2 < x_I, zeros,l2_I)
l2_I = np.maximum(0.0,l2_I)
#####
# Incident path length in sample powder
l1_I = ((x_R+specW/2)/cosT_I).reshape(nSize,nPsi,nU,nPhi)
l1_I = np.where(-specW/2 < x_I,(specD-z_R)/sinT_I,l1_I)
l1_I = np.maximum(0.0,l1_I)
#####
# Emission path length in sample powder
l1_E = ((specW/2-x_R)/cosT_E).reshape(nSize,nPsi,nU,nPhi)
l1_E = np.where(x_E<specW/2,(specD-z_R)/sinT_E,l1_E)
l1_E = np.maximum(0.0,l1_E)
# Emission path length in holder
l2_E = ((x_E-specW/2)/cosT_E).reshape(nSize,nPsi,nU,nPhi)
l2_E = np.where(x_E < specW/2,zeros,l2_E)
l2_E = np.maximum(0.0,l2_E)
######
# Relative intensity
thetaG_U_Phi = np.outer(thetaG,ones_U_Phi).reshape(nSize,nPsi,nU,nPhi)
isValid = np.where(theta_E > phiDS/2,True,False)
isValid = isValid & np.where(thetaG_U_Phi > phiDS/2,True,False)
# isValid = isValid & np.where(-specW/2 < x_I,True,False)
isValid = isValid & np.where(x_I < specW/2,True,False)
isValid = isValid & np.where(-specW/2 < x_E,True,False)
uu = np.outer(u,ones_Phi).reshape(nSize,nPsi,nU,nPhi)
phi_W = phi_U - phi_L # effective width of phi [rad.] (nSize,nPsi,nU)
phi_W = np.outer(phi_W,ones_Phi).reshape(nSize,nPsi,nU,nPhi)
# isValid = isValid & np.where(phiDS/2>0,True,False)
# exponent = -mu*(l1_I+l1_E)-mu2*(l2_I+l2_E) # (nSize,nPsi,nU)
exponent = -mu*(l1_I+l1_E)-mu2*(l2_I+l2_E)-np.log(uu)
#
#### ##############
# For debugging...
# debug = True
###################
debug = False
# debug = True
if (debug):
print('---------------------')
print(' degDS = ',degDS)
print(' width for degPhi = ',(np.full_like(phi_W,180/np.pi)*phi_W)[-1,0,0,0])
print(' lower limit for degPhi = ',(np.full_like(phi_L,180/np.pi)*phi_L)[-1,0,0])
print(' upper limit for degPhi = ',(np.full_like(phi_U,180/np.pi)*phi_U)[-1,0,0])
print(' at deg2T = ',degTwot[-1])
print('---------------------')
nPsi2,nU2,nPhi2 = int(nPsi/2),int(nU/2),int(nPhi/2)
print('nPsi2,nU2,nPhi2 = ',nPsi2,nU2,nPhi2)
print("----------------------")
print(" exponent[0,-1,-1,-1] = ",exponent[0,-1,-1,-1])
print(" exponent[1,-1,-1,-1] = ",exponent[1,-1,-1,-1])
print(" exponent[2,-1,-1,-1] = ",exponent[2,-1,-1,-1])
print(" exponent[3,-1,-1,-1] = ",exponent[3,-1,-1,-1])
print(" mu = ",mu)
print(" x_I[0,-1,-1,-1] = ",x_I.reshape(nSize,nPsi,nU,nPhi)[0,-1,-1,-1])
print(" x_I[1,-1,-1,-1] = ",x_E.reshape(nSize,nPsi,nU,nPhi)[1,-1,-1,-1])
print(" x_I[2,-1,-1,-1] = ",x_I.reshape(nSize,nPsi,nU,nPhi)[2,-1,-1,-1])
print(" x_I[3,-1,-1,-1] = ",x_E.reshape(nSize,nPsi,nU,nPhi)[3,-1,-1,-1])
print(" l1_I[0,-1,-1,-1] = ",l1_I.reshape(nSize,nPsi,nU,nPhi)[0,-1,-1,-1])
print(" l1_E[0,-1,-1,-1] = ",l1_E.reshape(nSize,nPsi,nU,nPhi)[0,-1,-1,-1])
print(" l2_I[0,-1,-1,-1] = ",l2_I.reshape(nSize,nPsi,nU,nPhi)[0,-1,-1,-1])
print(" l2_E[0,-1,-1,-1] = ",l2_E.reshape(nSize,nPsi,nU,nPhi)[0,-1,-1,-1])
print(" l1_I[1,-1,-1,-1] = ",l1_I.reshape(nSize,nPsi,nU,nPhi)[1,-1,-1,-1])
print(" l1_E[1,-1,-1,-1] = ",l1_E.reshape(nSize,nPsi,nU,nPhi)[1,-1,-1,-1])
print(" l2_I[1,-1,-1,-1] = ",l2_I.reshape(nSize,nPsi,nU,nPhi)[1,-1,-1,-1])
print(" l2_E[1,-1,-1,-1] = ",l2_E.reshape(nSize,nPsi,nU,nPhi)[1,-1,-1,-1])
print(" l1_I[2,-1,-1,-1] = ",l1_I.reshape(nSize,nPsi,nU,nPhi)[2,-1,-1,-1])
print(" l1_E[2,-1,-1,-1] = ",l1_E.reshape(nSize,nPsi,nU,nPhi)[2,-1,-1,-1])
print(" l2_I[2,-1,-1,-1] = ",l2_I.reshape(nSize,nPsi,nU,nPhi)[2,-1,-1,-1])
print(" l2_E[2,-1,-1,-1] = ",l2_E.reshape(nSize,nPsi,nU,nPhi)[2,-1,-1,-1])
print(" -2*muT = ",-2*muT)
print(" sinT = ",sinT[-1])
print(" u_L = ",u_L[-1])
print(" u_U = ",u_U[-1])
print(" Shape of u(3) : ",u.shape)
print(" uu = ",uu.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
# print(" z = ",z.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" z_R = ",z_R.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" x_I = ",x_I.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" x_R = ",x_R.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" x_E = ",x_E.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" sinT = ",sinT[-1])
print(" sinT_I = ",sinT_I.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" cosT_I = ",cosT_I.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
print(" tanT_I = ",tanT_I.reshape(nSize,nPsi,nU,nPhi)[-1,nPsi2,nU2,nPhi2])
# Calculation of g_ijk
exponent = np.minimum(exponent,100)
g = np.where(isValid, phi_W*np.exp(exponent)/phiDS,0.0)
g = g.reshape(nSize,nPsi,nU,nPhi)
# if (debug):
# print(" g (1) = ",g)
# g *= 1-np.exp(-2*muT/sinT_Psi_U_Phi)
# expc(x) = 1 - exp(-x)
# NOTE: If x is small value,
# the expression 1 - exp(-x) will cause
# a serious loss of significance !!!
g *= expc(2*muT/sinT_Psi_U_Phi)
g = g.reshape(nSize,nPsi,nU,nPhi)
# g *= 2
if (debug and (nSize < 15)):
print('specD_array[-1] = ',specD_array[-1])
print('z_R[-1]= ',z_R[-1])
print('(specD_array-z_R)[-1] = ',(specD_array-z_R)[-1])
print('l1_I[-1] = ',l1_I[-1])
print('x_I[-1] = ',x_I[-1])
print('g[-1] = ',g[-1])
##############
# INTEGRATION
##############
w_Phi = np.array(w_Phi).T
w_U = np.array(w_U).T
w_Psi = np.array(w_Psi).T
s0 = g @ w_Phi @ w_U @ w_Psi # (nSize,)
s1 = (del2T*g) @ w_Phi @ w_U @ w_Psi # (nSize,)
s2 = (del2T**2*g) @ w_Phi @ w_U @ w_Psi # (nSize,)
s3 = (del2T**3*g) @ w_Phi @ w_U @ w_Psi # (nSize,)
s4 = (del2T**4*g) @ w_Phi @ w_U @ w_Psi # (nSize,)
# Moments ...
m1 = s1 / s0
m2 = s2 / s0
m3 = s3 / s0
m4 = s4 / s0
# Calculate Cumulants up to 4th order
kappa1 = m1
kappa2 = m2 - m1**2
kappa3 = m3 - 3*m2*m1 + 2*sign(m1)*abs(m1)**3
kappa4 = m4 - 4*m3*m1 - 3*m2**2
kappa4 += 12*m2*sign(m1)*abs(m1)**2
kappa4 += -6*sign(m1)*abs(m1)**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])
debug = False
return kappas
# return kappas: (List of NumPy array)
# [s0,kappa1,kappa2,kappa3,kappa4]
##############################################
# Define Peak Profile Model for Deconvolution
# for sign == 1/-1,
##############################################
def fDeconv(alpha, sign, x):
# alpha # (float)
# sign = -1, 1
# x: NumPy array
def fFunc(gam_alpha,abs_x):
# np.seterr(divide="ignore")
# x**(alpha-1) * e(-x) / Gamma(alpha)
abs_x = np.where(abs_x>0,abs_x,np.full_like(abs_x,1E-8))
ans = abs_x**(alpha-1)*np.exp(-abs_x)/gam_alpha
# np.seterr(divide="warn")
return ans
gam_alpha = special.gamma(alpha) # Complete Gamma function
ans = np.where(x*sign>0,fFunc(gam_alpha,np.abs(x)),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
ans = np.power(1 + sign * 2*np.pi*1j * xi,-alpha)
return ans
##################################################
# "equat_trnspr6.py" -> "eq_tr6.treat"
# TREATMENT OF EQUATORIAL-TRANSPARENCY ABERRATION
##################################################
def treat(sign,source,kappas,cfg):
# sign: -1 or 1 (assumed to be -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: Reduced cumulants(NumPy array)
# cfg={'gonioR':gonioR,'degDS':degDS,'degLPSD':degLPSD,\
# 'specW':specW,'specT':specT,'muInv':muInv,'muInv2':muInv2,\
# 'nZ':nZ,'nPhi':nPhi,'nPsi':nPsi,\
# 'alphaEquat':alphaEquat,'marginEquat':marginEquat,\
# 'skipEquat':skipEquat,'error_type':error_type}
#******************************
# Configuration from "dct.cfg"
#******************************
marginEquat = cfg['marginEquat']
alphaEquat = cfg['alphaEquat']
skip = cfg['skipEquat']
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 Inverted Gamma
# f(x)=|x|**(alp-1)*exp(x)/Gamma(alp) for x<0
# s_angle (2Theta in deg.) => source_x
# s_int => source_y
#*****************************************************
#-----------------------------
# d(chi)/d(2Theta) (dchi_d2T)
#-----------------------------
dchi_d2T = -(2*alphaEquat)**(1/3)/kappas[3] # (positive, if kappas[3] is negative)
#------------------------------------------------------------
# Integration & creation of abscissa values source_chi
# Note: chi should be the integral of d(chi)/d(2Theta)
# by 2Theta (d2Theta)
#------------------------------------------------------------
source_chi = np.cumsum(dchi_d2T) - dchi_d2T[0]
source_chi *= (s_angle[-1]-s_angle[0])/(nData-1)
#------------------------------------------------------------
# Ordinate Correction Factor
#------------------------------------------------------------
source_corr = cmn.g_corr(s_angle/2)/dchi_d2T
#------------------------------------------------------------
# Another correction for lost intensity,
#------------------------------------------------------------
source_corr2 = 1/kappas[0]
# Corrected and Transformed Intensity
source_eta = s_int * source_corr
# print("source_y = ",source_y)
if (sign==-1): # intensity correction should be once...
source_eta *= source_corr2
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 1-2: Preparation for Fourier treatment
#*******************************************
margin = marginEquat
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]
chi = np.linspace(chi0,chiV,num=nValid,endpoint=False)
#*********************************************************
# Step 1-3:
# Interpolation of data to transformed scale
# (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(-nDest*dchi,0,num=nDest,endpoint=False)
deconv = fDeconv(alphaEquat,-1,chi)
deconv[0] = 1.0/dchi - np.sum(deconv[1:nDest])
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)
#
##########################
# Step 2: Shift treatment
##########################
skipShift = False
if (skipShift):
return dest
#-------------------------------
# Shift expected in the results
#-------------------------------
cum1_R = kappas[1] + alphaEquat/dchi_d2T
deg2T = s_angle
s_int = d_int
nSize = deg2T.size
ddeg2T = (deg2T[-1]-deg2T[0])/(nSize-1)
# Number of Left margin needed...
# If first cumulant is negative, extrapolation is needed
nLeft = np.maximum(np.ceil(-cum1_R[0]/ddeg2T),0).astype(int)
print('nLeft=',nLeft)
# Number of Right margin needed...
# If last cumulant is positive, extrapolation is needed
nRight = np.maximum(np.floor(cum1_R[-1]/ddeg2T),0).astype(int)
# Number of Total data...
# print('cum1_R[0],cum1_R[-1]=',cum1_R[0],cum1_R[-1])
nExpand = nLeft+nSize+nRight # Expanded size of data
deg2T_E = np.empty(nExpand,dtype='float')
dct_E = np.empty(nExpand,dtype='float')
iRight = nLeft+nSize
# print('nLeft,nRight,nExpand,iRight = ',nLeft,nRight,nExpand,iRight)
# Left data...
deg2T_E_0 = deg2T[0]-cum1_R[0]-nLeft*ddeg2T # First value of Left data
deg2T_E_1 = deg2T[0]-cum1_R[0] # Last value of Left data
if (nLeft > 0):
deg2T_E[0:nLeft] = np.linspace(deg2T_E_0,deg2T_E_1,num=nLeft,endpoint=False)
dct_E[0:nLeft] = s_int[0]
# Middle data...
deg2T_E[nLeft:iRight] = deg2T[0:nSize]-cum1_R[0:nSize] # Core data !!!
dct_E[nLeft:iRight] = s_int
# Right data...
deg2T_E_0 = deg2T[-1]-cum1_R[-1]+ddeg2T # First value of Right data
deg2T_E_1 = deg2T_E_0+nRight*ddeg2T # Last value of Right data
if (nRight > 0):
deg2T_E[iRight:nExpand] = np.linspace(deg2T_E_0,deg2T_E_1,num=nRight,endpoint=False)
dct_E[iRight:nExpand] = s_int[-1]
# print('deg2T_E_0,deg2T_E_1,ddeg2T=',deg2T_E_0,deg2T_E_1,ddeg2T)
# deg2T_E[iRight:nExpand] = np.arange(deg2T_E_0,deg2T_E_1,ddeg2T)
d_int = np.interp(deg2T,deg2T_E,dct_E) # linear interpolation
dest[1] = d_int
print('error_type = ',error_type)
return dest