Individual peak profile fitting of Si XRD data.
# ippf_Si.py
# 個別ピークデータの抽出と描画，
# scipy.optimize.leastsq() メソッドによる非対称化 Voigt 曲線あてはめ
# Peak profile model
# defined as an asymmetrized function of the convolution of the Lorentzian function
# with a symmetric function defined by standard deviation and kurtosis
#
# Coded by T. Ida, Apr. 14, 2023
# Updated by T. Ida, May 10, 2023
# Updated by T. Ida, Apr. 28, 2024
#
import sys
import numpy as np
import scipy.special as sp
import scipy.optimize as opt
import matplotlib.pyplot as plt
import warnings
from scipy.optimize import OptimizeWarning
def main():
# 入力ファイル名の指定...
sRaw = '00raw.csv' # 処理前ファイル名
sDCT = '01dct.csv' # 処理後ファイル名
deg_hw = 1.0 # 参照ピーク位置 ± deg_hw の範囲のデータを抽出する
# 入力ファイルからのデータの読み込み...
source1 = np.loadtxt(sRaw,dtype='float',delimiter=',') # 処理前データの読み込み
x = source1[:,0] # x は各行 0 列目
y1 = source1[:,1] # y (raw) は各行 1 列目
source2 = np.loadtxt(sDCT,dtype='float',delimiter=',') # 処理後データの読み込み
y2 = source2[:,1] # y (DCT) は各行 1 列目
# 単純な誤差モデル（正当化できるかは不明...）
err = np.sqrt(np.maximum(1,y2))
dx = x[1] - x[0] # 標本点の間隔
nP = x.size # 標本点の数
# Peak data [2theta, d, I, h, k, l] for Si & LaB6 / Cu K_alpha
# (NIST SRM640d & SRM660b)
srm = 640 # SRM640d
if (srm == 640): # SRM 640 シリーズ Si の場合
peaks=[[28.441,3.1380,100,1,1,1],[47.300,1.9200,60,2,2,0]]
peaks+=[[56.119,1.6380,35,3,1,1],[69.126,1.3570,8,4,0,0]]
peaks+=[[76.371,1.2460,13,3,3,1],[88.024,1.1083,17,4,2,2]]
peaks+=[[94.946,1.0450,9,5,1,1],[106.700,0.9599,5,4,4,0]]
peaks+=[[114.082,0.9178,11,5,3,1],[127.532,0.8586,9,6,2,0]]
peaks+=[[136.877,0.8281,5,5,3,3]]
else: # SRM 660 シリーズ LaB6 の場合
peaks=[[21.358,4.157660,54,1,0,0],[30.385,2.939180,100,1,1,0]]
peaks+=[[37.442,2.399800,41,1,1,1],[43.507,2.077980,22,2,0,0]]
peaks+=[[48.958,1.858620,46,2,1,0],[53.989,1.696870,24,2,1,1]]
peaks+=[[63.219,1.469700,8,2,2,0],[67.548,1.385340,23,3,0,0]]
peaks+=[[71.746,1.314350,16,3,1,0],[75.844,1.253290,10,3,1,1]]
peaks+=[[79.870,1.200010,2,2,2,2],[83.846,1.152880,6,3,2,0]]
peaks+=[[87.792,1.110970,13,3,2,1],[95.672,1.039280,2,4,0,0]]
peaks+=[[99.643,1.008220,8,4,1,0],[103.661,0.979826,7,4,1,1]]
peaks+=[[107.749,0.953625,3,3,3,1],[111.934,0.929477,4,4,2,0]]
peaks+=[[116.245,0.907086,9,4,2,1],[120.723,0.886245,3,3,3,2]]
peaks+=[[130.409,0.848500,3,4,2,2],[135.800,0.831399,3,5,0,0]]
#
peaks = np.array(peaks) # NumPy 配列 (11行6列)または（22行6列）
nRow = peaks.shape[0] # 行数
nCol = peaks.shape[1] # 列数
ttRef = np.zeros(nRow) # 参照ピーク位置用の NumPy 配列の確保
ttOpt = np.zeros(nRow) # 最適化ピーク位置用の NumPy 配列の確保
indRef = np.arange(nRow,dtype='int') # 参照番号用の整数配列の確保
deltaTT = np.zeros(nRow) # Δ2Θ
err_ttOpt = np.zeros(nRow) # 最適化ピーク位置の推定誤差用の NumPy 配列の確保
gamOpt = np.zeros(nRow) # 非対称性パラメータ用の NumPy 配列の確保
err_gamOpt = np.zeros(nRow) # 非対称性パラメータ推定誤差用の NumPy 配列の確保
#
iRow = 0 # 参照ピークデータの行番号
iOpt = 0 # 有効なピークの番号
while (iRow < nRow):
ttRef[iOpt] = peaks[iRow][0] # 初期ピーク位置
# ピーク位置周辺データの抽出
jFirst = np.maximum(0,np.floor((peaks[iRow][0] - deg_hw - x[0]) / dx).astype(int))
jLast = np.minimum(nP,np.ceil((peaks[iRow][0] + deg_hw - x[0]) / dx).astype(int))
xsub, y1sub, y2sub = x[jFirst:jLast], y1[jFirst:jLast], y2[jFirst:jLast]
errsub = err[jFirst:jLast]
# はじめの推定値
bkg = 0.5*(y2[jFirst]+y2[jLast]) # 定数背景強度推定値
intn = (y2[jFirst:jLast].mean()-bkg)*dx*(jLast-jFirst) # ピーク積分強度推定値
xpeak = x[np.round(0.5*(jFirst+jLast)).astype(int)] # ピーク位置推定値
sigm = 0.1 # Gauss 型成分標準偏差初期値
w = 180/np.pi*0.00035*np.tan(xpeak*np.pi/360) # Lorentz 型成分半値半幅初期値
w_sigm = w/sigm # Lorentzian 半値半幅と Gaussian 標準偏差の比
gam = -0.01 # 裁断指数函数型成分の減衰幅初期値
gam_sigm = gam/sigm # 裁断指数函数型成分減衰幅と Gaussian 標準偏差の比
param = np.array([bkg,intn,xpeak,sigm,w_sigm,gam_sigm]) # 初期推定パラメータ (NumPy 配列）
param0 = np.copy(param) # 初期推定パラメータの複製・保持
#
# 最適化アルゴリズム
# method = 'curve_fit' # 'curve_fit','leastsq',or 'minimize(BFGS)',...
method = 'leastsq' # 'curve_fit','leastsq',or 'minimize(BFGS)',...
# method = 'minimize(BFGS)' # 'curve_fit','leastsq',or 'minimize(BFGS)',...
#
# 初期推定値についてのコンソール表示
sText = "{:d}".format(peaks[iRow][3].astype(int)) # index h
sText += "{:d}".format(peaks[iRow][4].astype(int)) # index k
sText += "{:d}".format(peaks[iRow][5].astype(int)) # index l
print("*****")
print('Optimization for ' + sText + '-reflection, '+method+' method')
print('Initial Guess...')
print(' Background : {:.0f}'.format(bkg))
print(' Integrated intensity : {:.0f}'.format(intn))
print(' Peak location : {:.4f}'.format(xpeak))
print(' Gaussian sigma : {:.4f}'.format(sigm))
print(' Lorentzian HWHM : {:.4f}'.format(w))
print(' Asymmetry: {:.4f}'.format(gam))
chi_sq = f_cost(param,xsub,y2sub,errsub) # コスト函数
print(' ----')
print(' chi**2 = ',chi_sq)
# 曲線当て嵌め
opt_results={'success':False,'popt':np.array([]),'psigma':np.array([]),'nfev':0} # 辞書型変数の初期化
warnings.simplefilter("error", OptimizeWarning)
try:
if (method == 'curve_fit'): # "scipy.optimize.curve_fit()" method
popt,pcov,infodict,mesg,ier = opt.curve_fit(f_model,xsub,y2sub,p0=param0,sigma=errsub, \
absolute_sigma=True,method='lm',full_output=True)
psigma = np.sqrt(np.diag(pcov))
opt_results={'success':True,'popt':popt,'psigma':psigma,'nfev':infodict['nfev']}
elif (method == 'leastsq'): # "scipy.optimize.leastsq()" method
popt,pcov,infodict,mesg,ier = opt.leastsq(f_residual,param0,args=(xsub,y2sub,errsub),maxfev=2800,full_output=True)
# maxfev: default value = (6+1)*200
psigma = np.sqrt(np.diag(pcov))
opt_results={'success':True,'popt':popt,'psigma':psigma,'nfev':infodict['nfev']}
elif (method == 'minimize(BFGS)'): # "scipy.optimize.minimize()" method
result = opt.minimize(f_cost,param0,args=(xsub,y2sub,errsub),method='BFGS')
popt = result.x
psigma = np.sqrt(np.diag(result.hess_inv))
nfev = result.nfev
opt_results={'success':True,'popt':popt,'psigma':psigma,'nfev':nfev}
except ValueError as e: # "ValueError" が起こった場合
sys.stderr.write("*** ValueError in "+method+" ***\n")
sys.stderr.write(str(e) + "\n")
opt_results={'success':False,'popt':np.array([]),'psigma':np.array([])}
except OptimizeWarning as e: # "OptimizeWarning" が起こった場合
sys.stderr.write("*** OptimizeWarning in "+method+" ***\n")
sys.stderr.write(str(e) + "\n")
opt_results={'success':False,'popt':np.array([]),'psigma':np.array([])}
except Exception as e:
sys.stderr.write("*** Exception in "+method+" ***\n")
sys.stderr.write(str(e) + "\n")
opt_results={'success':False,'popt':np.array([]),'psigma':np.array([])}
# print('opt_results = ',opt_results)
if (opt_results['success']): # 最適化に成功した場合 ...
# コンソール出力...
print("=====")
print('Optimization succeeded.')
popt = opt_results['popt']
psigma = opt_results['psigma']
# 最適化の結果 ...
print('Optimized ...')
print(' Background : {:.0f}'.format(popt[0]),end='')
print(' ± {:.0f}'.format(psigma[0]))
print(' Integrated intensity : {:.0f}'.format(popt[1]),end='')
print(' ± {:.0f}'.format(psigma[1]))
ttRef[iOpt] = peaks[iRow][0]
ttOpt[iOpt] = popt[2] # Peak 2Theta optimized
err_ttOpt[iOpt] = psigma[2] # Estimated error in peak 2Theta
print(' Peak location : {:.4f}'.format(popt[2]),end='')
print(' ± {:.4f}'.format(psigma[2]))
sigm_opt = popt[3]
sigm_sigma = psigma[3]
print(' Gaussian sigma : {:.4f}'.format(sigm_opt),end='')
print(' ± {:.4f}'.format(sigm_sigma))
w_opt = popt[4] * sigm_opt
w_sigma = w_opt * np.sqrt((psigma[4]/popt[4])**2 + (sigm_sigma/sigm_opt)**2)
print(' Lorentzian HWHM : {:.4f}'.format(w_opt),end='')
print(' ± {:.4f}'.format(w_sigma))
gamOpt[iOpt] = popt[5]*popt[3]
err_gamOpt[iOpt] = popt[5]*popt[3] * np.sqrt((psigma[5]/popt[5])**2 + (psigma[3]/popt[3])**2)
print(' Asymmetry: {:.4f}'.format(gamOpt[iOpt]),end='')
print(' ± {:.4f}'.format(err_gamOpt[iOpt]))
print(' ----')
chi_sq = f_cost(popt,xsub,y2sub,errsub) # コスト函数
print(' chi**2 = ',chi_sq)
print(' Number of function calls = ',opt_results['nfev'])
#
fit = f_model_AsymVoigt(xsub,popt) # fit curve
res = y2sub - fit # residual
# y1 と y2 の最大の値 ...
y_max = np.max(y1[jFirst:jLast])
y_max = np.maximum(y_max,np.max(y2[jFirst:jLast]))
# グラフにテキストを書き入れる位置の指定 ...
xText = x[jFirst]+(x[jLast]-x[jFirst])*0.05
yText = y_max * 0.9
#
###################################################
# Make a fitting figure with a difference plot
# 差プロットを付けたピーク形状フィッティング図の描画 (ピークシフト）
###################################################
plt.rcParams['font.family'] = 'Serif' # セリフ体のフォントを使う
plt.rcParams['mathtext.fontset'] = 'stix' # 数式には stix フォントを使う
# plt.rcParams['font.size'] = 9 # フォントサイズは 9 pt
fig = plt.figure(figsize=(6.4,3.6)) # 幅 6.4 インチ，高さ 3.6 インチ
gs = fig.add_gridspec(nrows=2,ncols=1,hspace=0.0,height_ratios=(7,3))
# ２つのパネル を上下 7:3 で配置
ax0 = fig.add_subplot(gs[0,0]) # 0 行目 (フィッティング図）
ax1 = fig.add_subplot(gs[1,0]) # 1 行目 (差プロット）
#######################################
# ピーク形状フィッティング図（散布図と折れ線グラフ）
#######################################
# 逆畳込的処理前の強度図形，折れ線グラフ.
ax0.plot(xsub,y1sub,label= 'raw') # 折れ線グラフ，凡例ラベル設定
# 逆畳込的処理後の強度図形，散布図（マーカーでプロット）
ax0.scatter(xsub,y2sub,label= 'DCT',marker='o',color='orange') # 散布図，凡例ラベル設定
ax0.xaxis.set_label_position('top')
ax0.xaxis.tick_top()
ax0.minorticks_on() # 補助目盛表示
ax0.tick_params(length=6) # 目盛線の長さ 6 pt (約 2.1 mm)
ax0.plot(xsub,fit,label='fit',color='green')
ax0.text(xText,yText,sText+'\n'+method+' method')
ax0.axvline(x=peaks[iRow][0],linewidth=0.5) # 垂直線描画
# ax0.set_xlabel('$2\Theta (^\circ)$') # x 軸ラベル
ax0.set_ylabel('Intensity (counts)') # y 軸ラベル
ax0.set_xlim(x[jFirst],x[jLast]) # x 軸描画範囲
ax0.set_ylim(0,y_max/0.9) # y 軸描画範囲
handle, label = ax0.get_legend_handles_labels() # 凡例のハンドルとラベル
# print('handle = ',handle)
# print('label = ',label)
handle=[handle[0],handle[2],handle[1]] # 凡例の順序を変える
label=[label[0],label[2],label[1]] # 凡例ラベルの順序を変える
# print('handle = ',handle)
# print('label = ',label)
ax0.legend(handle,label)
#####################
# 差プロット（折れ線グラフ）
#####################
ax1.plot(xsub,res,label = 'residual',color='orange')
ax1.axvline(x=peaks[iRow][0],linewidth=0.5) # 垂直線描画
ax1.axhline(y=0,linewidth=0.5) # 水平線描画
ax1.set_xlim(x[jFirst],x[jLast])
res_y_max = y_max/0.9*3/7*0.5
ax1.set_ylim(-res_y_max,res_y_max)
ax1.set_ylabel('Difference (counts)',labelpad=20)
ax1.tick_params(labelleft=False)
ax1.set_xlabel('$2\Theta (^\circ)$')
ax1.minorticks_on() # 補助目盛表示
ax0.tick_params(length=6) # 目盛線の長さ 6 pt (約 2.1 mm)
plt.subplots_adjust(hspace=0.0)
ax1.legend() # 凡例描画
#
file_name = 'fig{0:03}.png'.format(iOpt+1)
plt.savefig(file_name,dpi=600)
plt.show()
#
iOpt += 1 # 最適化成功インデックスを 1 増やす
else: # 最適化に失敗した場合
print('Optimization failed for ' + sText + '-reflection')
iRow += 1
nOpt = iOpt # 最適化に成功したピークの数
print("Number of optimized/total peaks : ",nOpt,"/",nRow)
indRef = np.resize(indRef,nOpt)
ttRef = np.resize(ttRef,nOpt)
ttOpt = np.resize(ttOpt,nOpt)
deltaTT = ttOpt - ttRef
err_ttOpt = np.resize(err_ttOpt,nOpt)
gamOpt = np.resize(gamOpt,nOpt)
err_gamOpt = np.resize(err_gamOpt,nOpt)
#
########################################
# 素朴ピークシフトモデルによるピークシフトのモデル化
########################################
print('')
print("*****")
print('Optimization for peak shift by a naive model')
# 初期推定値についてのコンソール表示
param0 = np.array([0.01,0.01])
print("*****")
print('Initial Guess...')
print(' Δ2Θ_0 (º) : {:.4f}'.format(param0[0]))
print(' Specimen displacement (mm) : {:.4f}'.format(param0[1]))
res = f_res_naive(param0,ttRef,deltaTT,err_ttOpt) # 残差
chi_sq = np.sum(res**2)
print(' chi**2 = ',chi_sq)
print("Optimization starts...")
# print('ttRef = ',ttRef)
# print('deltaTT = ',deltaTT)
# print('err_ttOpt = ',err_ttOpt)
opt_results = {'success':False,'popt':np.array([]),'psigma': np.array([])}
try:
# popt,pcov = opt.curve_fit(f_model_naive,ttRef,deltaTT,err_ttOpt,p0=param0,sigma=err_gamOpt,method='lm')
popt,pcov,infodict,mesg,ier = opt.leastsq(f_res_naive,param0,args=(ttRef,deltaTT,err_ttOpt),maxfev=2800,full_output=True)
psigma = np.sqrt(np.diag(pcov))
opt_results = {'success':True,'popt':popt,'psigma':psigma}
except ValueError as e:
print("Optimization failed")
sys.stderr.write("*** ValueError in "+method+" ***\n")
sys.stderr.write(str(e) + "\n")
print('')
opt_results = {'success':False,'popt':np.array([]),'psigma': np.array([])}
except OptimizeWarning as e:
print("Optimization failed")
sys.stderr.write("*** OptimizeWarning in "+method+" ***\n")
sys.stderr.write(str(e) + "\n")
print('')
opt_results = {'success': False,'popt':np.array([]),'psigma': np.array([])}
if (opt_results['success']):
print("Optimization finished.")
popt = opt_results['popt']
psigma = opt_results['psigma']
print(' Δ2Θ_0 = {:.4f}º ± {:.4f}º'.format(popt[0],psigma[0]))
print(' Specimen displacement: {:.4f} ± {:.4f} mm'.format(popt[1],psigma[1]))
res = f_res_naive(popt,ttRef,deltaTT,err_ttOpt) # 残差
chi_sq = np.sum(res**2)
print(' chi**2 = ',chi_sq)
#
fig = plt.figure(figsize=(6.4,3.6)) # 寸法を指定して figure を作成
gs = fig.add_gridspec(nrows=2,ncols=1,hspace=0.0,height_ratios=(7,3))
# ２つのパネル を上下 7:3 で配置
ax0 = fig.add_subplot(gs[0,0]) # 0 行目（上パネル） (フィッティング図）
ax1 = fig.add_subplot(gs[1,0]) # 1 行目（下パネル）(差プロット）
#######################################################################
# （上パネル）ピークシフト Δ2Θ の回折角依存性とフィッティング図（散布図と折れ線グラフ） の描画
#######################################################################
ax0.xaxis.tick_top() # 上パネルでは x軸目盛を上辺に
ax0.minorticks_on() # 補助目盛表示
ax0.set_xlabel(r'$2\Theta\ (^\circ)$') # x軸ラベル: 2Θ (º)
ax0.set_ylabel(r'$\Delta 2\Theta\ (^\circ)$') # y軸ラベル: Δ2Θ (º)
# ピークシフトの回折角依存性 Δ2Θ，散布図（誤差棒つきマーカーでプロット）
error_bar_label = r'$\Delta 2 \Theta \,(^\circ)$' # ピークシフト，マーカー凡例ラベル
ax0.errorbar(ttRef,deltaTT,yerr=err_ttOpt,fmt='o',capsize=5,label=error_bar_label)
# 当て嵌め曲線（折れ線グラフ）
dttFit_x = np.linspace(ttRef[0],ttRef[-1],400)
dttFit_y = f_model_naive(dttFit_x,popt[0],popt[1])
fit_curve_label = 'Naive model' # 当て嵌め曲線凡例ラベル
ax0.plot(dttFit_x,dttFit_y,label=fit_curve_label)
ax0.axhline(y=0.0,linewidth=0.5) # 水平線描画
# ピークシフトと回折ピーク角の最大値と最小値 ...
x_min = np.min(ttRef)
x_max = np.max(ttRef)
y_min = np.minimum(np.min(deltaTT),np.min(dttFit_y))
y_max = np.maximum(np.max(deltaTT),np.max(dttFit_y))
# x軸と y軸の描画範囲...
ax0_x_min = np.maximum(0.0,x_min-(x_max-x_min)*0.1)
ax0_x_max = np.minimum(180.0,x_max+(x_max-x_min)*0.1)
ax0_y_min = y_min-(y_max-y_min)*0.1
ax0_y_max = y_max+(y_max-y_min)*0.1
ax0.set_xlim(ax0_x_min,ax0_x_max) # x 軸描画範囲
ax0.set_ylim(ax0_y_min,ax0_y_max) # y 軸描画範囲
handle, label = ax0.get_legend_handles_labels() # 凡例のハンドルとラベル
# print('handle,label = ',handle,lsigma)
handle=[handle[1],handle[0]] # 凡例の順序を変える
label=[label[1],label[0]] # 凡例ラベルの順序を変える
ax0.legend(handle,label) # 上パネル凡例表示
#######################################
# (下パネル）ピークシフトモデルについての差プロット
#######################################
ttRes = deltaTT - f_model_naive(ttRef,popt[0],popt[1])
ax1.plot(ttRef,ttRes,label='residual')
ax1.minorticks_on()
ax1.tick_params(length=6)
ax1.set_xlabel(r'$2\Theta\ (^\circ)$') # x軸ラベル 2Θ (º)
ax1.set_ylabel(r'${\rm Difference}\ (^\circ)$') # y軸ラベル Difference (º)
ax1.set_xlim(ax0_x_min,ax0_x_max) # x軸描画範囲は ax0 と揃える
ax1_y_max = (ax0_y_max-ax0_y_min)*3.0/7.0/2 # ax0 と尺度を揃える
ax1.set_ylim(-ax1_y_max,ax1_y_max) # y軸描画範囲
ax1.axhline(y=0,linewidth=0.5) # 水平線描画
ax1.legend() # 下パネル凡例表示
#
plt.savefig('fig{0:03}.png'.format(nOpt+1),dpi=600)
plt.show()
delTT = ttOpt - ttRef
################################
# 非対称性パラメータの角度依存性の表示
################################
print('')
print("*****")
print("Weighted average of asymmetry parameter...")
numer = np.sum(gamOpt/err_gamOpt**2) # numerator
denom = np.sum(1.0/err_gamOpt**2) # denominator
avg_gamOpt = numer / denom
print(" Average of asymmetry (º) : {0:4f}".format(avg_gamOpt))
print(" -----")
###################################################
# Show a figure with average level
# 平均値表示を付けたピーク形状フィッティング図の描画
# (非対称パラメータ）
###################################################
# 非対称性パラメータの回折角依存性 散布図（誤差棒つきマーカーでプロット）
fig = plt.figure(figsize=(6.4,3.6)) # 幅 6.4 インチ，高さ 3.6 インチ
gs = fig.add_gridspec(nrows=1,ncols=1)
ax0 = fig.add_subplot(gs[0,0]) # 0 行目 (フィッティング図）
error_bar_label = r'$\gamma \ (^\circ)$' # ピークシフト，マーカー凡例ラベル
ax0.errorbar(ttRef,gamOpt,yerr=err_gamOpt,fmt='o',capsize=5,label=error_bar_label)
# 平均線描画
avg_x = np.array([ttRef[0],ttRef[-1]])
avg_y = np.array([avg_gamOpt,avg_gamOpt])
avg_line_label = 'Average' # 平均線凡例ラベル
ax0.plot(avg_x,avg_y,label=avg_line_label)
# ゼロ線描画
ax0.axhline(y=0.0,linewidth=0.5)
# 回折ピーク角と非対称パラメータの最大値と最小値 ...
x_min = np.min(ttRef)
x_max = np.max(ttRef)
y_min = np.min(gamOpt)
y_max = np.max(gamOpt)
ax0_x_min = np.maximum(0.0,x_min-(x_max-x_min)*0.1) # ± 10 % の余白をつける
ax0_x_max = np.minimum(180.0,x_max+(x_max-x_min)*0.1)
ax0_y_min = np.minimum(0,y_min-(y_max-y_min)*0.1)
ax0_y_max = np.maximum(0,y_max+(y_max-y_min)*0.1)
ax0.minorticks_on() # 補助目盛表示
ax0.set_xlabel(r'$2\Theta\ (^\circ)$') # x軸ラベル: 2Θ (º)
ax0.set_ylabel(r'$\gamma\ (^\circ)$') # y軸ラベル: Δ2Θ (º)
ax0.set_xlim(ax0_x_min,ax0_x_max) # x 軸描画範囲
ax0.set_ylim(ax0_y_min,ax0_y_max) # y 軸描画範囲
ax0.xaxis.set_ticks_position('both')
ax0.yaxis.set_ticks_position('both')
handle, label = ax0.get_legend_handles_labels() # 凡例のハンドルとラベル
handle=[handle[1],handle[0]] # 凡例の順序を変える
label=[label[1],label[0]] # 凡例ラベルの順序を変える
ax0.legend(handle,label) # 上パネル凡例表示
plt.savefig('fig{0:03}.png'.format(nOpt+2),dpi=600)
plt.show()
#****************************************************
# Model function for scipy.optimize.curve_fit()
#****************************************************
def f_model(data_x,*param): # モデル函数
# param: parameters to be optimized (NumPy array)
# data_x: x data (NumPy arrays)
return f_model_AsymVoigt(data_x, param)
#**********************************************
# Residual function for scipy.optimize.leastsq
#**********************************************
def f_residual(param,*data): # 残差函数
# param: parameters to be optimized (NumPy array)
data_x,data_y,data_e = data
# data_x, data_y, data_e: x, y, error data (NumPy arrays)
fit_y = f_model_AsymVoigt(data_x, tuple(param))
residual = (data_y - fit_y) / data_e
return residual
#***********************************************
# Cost function for scipy.optimize.minimize
#***********************************************
def f_cost(param,data_x,data_y,data_e): # コスト函数
residual = f_residual(param,data_x,data_y,data_e)
chisq = np.sum(residual**2) # 重み付き平方和
# print('chisq = ',chisq)
return chisq
#*******************************
# Asymmetrized Voigtian profile
#*******************************
def asym_voigt_profile(x, sigm, w, gam):
ans = np.where(gam != 0, asym_sub(x,sigm,w,gam), sp.voigt_profile(x,sigm,w))
return ans
def asym_sub(x, sigm, w, gam):
n_x = x.size
n_z = 20
dz = 1.0 / n_z
# (1/n_z) * \sum_{i=1}^{n_z} fVoigt(x+gam*ln(z),sigm,w) * dz
z = dz/2 + np.linspace(0,1.0,num=n_z,endpoint=False,dtype='float')
a = np.repeat(x,n_z)
a = np.reshape(a,(n_x,n_z))
b = np.tile(z,n_x)
b = np.reshape(b,(n_x,n_z))
y = np.resize(gam*np.log(b),(n_x,n_z))
x_y = a + gam*np.log(b)
c = sp.voigt_profile(x_y,sigm,w) * dz
ans = np.sum(c, axis=1)
return ans
######################################
# Asymmetrized Voigtian profile model
######################################
def f_model_AsymVoigt(x, args):
b,I,x0,sigm,w_sigm,gam_sigm = args
# b: background
# I: intensity
# x0: peak location
# sigm: sigma of non-Lorentzian component
# w_sigm: Lorentzian HWHM / sigm
# gam_sigm: asymmetry parameter / sigm
return b + I * asym_voigt_profile(x - x0, sigm, w_sigm*sigm, gam_sigm*sigm)
################################################
# Naive model for peak shift ピークシフトの素朴モデル
################################################
# def f_model_naive(degTT,*param):
# degDelTT0,sDis = param
def f_model_naive(degTT,degDelTT0,sDis):
# degTT: 2Theta (deg.) 回折ピーク位置
# degDelTT0: Delta 2 Theta_0 (deg.) オフセット誤差
# sDis: specimen displacement (mm) 試料位置ずれ
R = 150 # goniometer radius (mm)
return degDelTT0 + (360/np.pi)*sDis*np.cos(degTT*np.pi/360)/R
def f_res_naive(param,ttRef,deltaTT,errTT): # 残差函数
dTTobs = deltaTT
dTTcal = f_model_naive(ttRef,param[0],param[1])
residual = (dTTobs-dTTcal)/errTT # 重み付き残差
return residual
##################
# Asymmetry model
##################
def f_model_asymmetry(degTT,muInv):
R =150
return -muInv*np.sin(degTT*np.pi/180)/(2*R)*180/np.pi
def f_res_asymmetry(param,ttRef,gamOpt,err_gamOpt):
gamCal = f_model_asymmetry(ttRef,param[0])
residual = (gamOpt-gamCal)/err_gamOpt
return residual
# Test routine ...
# testx = np.arange(-1.0,1.0,0.01)
# testy = fLorSymFit(testx, 10, 100, 0.01, 0.1,0.1,4)
# fig, ax = plt.subplots()
# ax.plot(testx,testy)
# plt.show()
main()