Python scipyで離散データの内挿 interpolate 積分 integrateを使ってみる

線形補間

# サンプルデータを取得
import urllib.request
with urllib.request.urlopen("https://gitlab.com/gifuescan/atima_rangeenergy/-/raw/1e7c7d1432c302d0705a5ea9fe40d696ec468e54/splines_emul_1.41/BeamZ1_TargetZ6.txt") as web_file:
    data = web_file.read()
    with open("BeamZ1_TargetZ6.txt", mode='wb') as local_file:
        local_file.write(data)

import numpy as np
data = np.loadtxt("BeamZ1_TargetZ6.txt")

内挿は簡単でデータ点を詰めて関数を取得する。kindは‘linear’, ‘nearest’などが使えるが、今回は'cubic'にした。

# 内挿する関数
from scipy import interpolate
f = interpolate.interp1d(data.T[0], data.T[1], kind='cubic')
print(f(0.5))

以下、内挿の関数やデータは全く使わない。積分(Cumulatively integrate)も同じくデータ点を詰めて初期値を与えるだけである。dEdx[MeV cm2/g]の逆数をエネルギー[MeV]で積分するとRange[g/cm2]が得られるので、

#dE/dxとEnergyからRangeを数値計算の積分で求める
from scipy import integrate
Energy = data.T[0]
dEdx = data.T[1]/1000.0
Range = integrate.cumtrapz(1.0/dEdx, Energy, initial=0)

描画をしておく

#描画
import matplotlib.pyplot as plt
plt.plot(Energy, data.T[2], 'r.', label="Original")
plt.plot(Energy, Range, 'b-', label="Integral")
plt.xscale("log")
plt.xlabel("Energy [MeV]")
plt.yscale("log")
plt.ylabel(r"Range [mg/cm2]")
plt.legend()     
plt.show()

import matplotlib.pyplot as plt
plt.plot(Energy[1:], (Range[1:]/data.T[2][1:]-1)*100, 'b-')
plt.xscale("log")
plt.xlabel("Energy [MeV]")
plt.ylim(-0.1,0.1)
plt.ylabel("Difference [%]")
plt.grid(True)
plt.show()