import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])
# rcPramsを使ってTeXのフォントをデフォルトにする
plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定

np.set_printoptions(precision=5)

pi2 = 2.0*np.pi
eps = 1.0e-10
n = 1024
span = 20.0
diameter = 1.0
print('データ数 n:',n)

データ数 n: 1024

n_x = n
x_span = span
d_x = x_span/n_x
x_min = -x_span/2.0
x_max = -x_min-d_x
x = np.linspace(x_min, x_max, n_x)
print('xmin =',x[0],', xmax =',x[n_x-1],', d_x =',d_x)

xmin = -10.0 , xmax = 9.98046875 , d_x = 0.01953125

#x

d_kx = pi2/x_span
kx_min = -d_kx*n_x/2
kx_max = -kx_min-d_kx
kx = np.linspace(kx_min, kx_max, n_x)
print('kx_min =',kx[0],', kx_max =',kx[n_x-1],', d_kx =',d_kx)

kx_min = -160.8495438637974 , kx_max = 160.53538459843844 , d_kx = 0.3141592653589793

xr = diameter/2# 方形関数の範囲
a = 1.25# 方形関数の振幅
f = np.where((x>=-xr) & (x<xr), a, 0.0) # 方形関数値の計算
fs = np.fft.ifftshift(f) # データ逆シフト
gs = d_x*np.fft.fft(fs) # 高速フーリエ変換（FFT）
# gs = d_x*np.fft.ifft(fs) # iFFT
g = np.fft.fftshift(gs) # データのシフト
g = g/pi2

#df = pd.DataFrame({
#    '$x$':x,
#    '$f(x)$':f,
#    '$k_x$':kx,
#    '$g(k_x)$':g,
#})
#df.style.format({"$x$":"{:.2f}"})

nn = 100*n
xx_min = x_min
xx_max = -xx_min
xx = np.linspace(xx_min, xx_max, nn)
ff = np.where((xx>=-xr) & (xx<xr), a, 0.0) # 方形関数値の計算
xx,ff

(array([-10.     ,  -9.9998 ,  -9.99961, ...,   9.99961,   9.9998 ,
         10.     ]),
 array([0., 0., 0., ..., 0., 0., 0.]))

kkx_min = kx_min
kkx_max = -kkx_min
kkx = np.linspace(kkx_min, kkx_max, nn)
kkx = kkx + (kkx==0)*eps # ゼロのときだけepsを加える
u = xr*kkx
gg = 2.0*a*xr*np.sin(u)/u # 方形関数のフーリエ変換（解析）
gg = gg/pi2

u[nn//2-3:nn//2+3]

array([-0.00393, -0.00236, -0.00079,  0.00079,  0.00236,  0.00393])

f.dtype, ff.dtype

(dtype('float64'), dtype('float64'))

fig = plt.figure() # グラフ領域の作成
plt.stem(x,f,'o',label='FFT')
plt.plot(xx,ff,'-',color='tab:orange',label='sinc')
plt.title('Space domain')
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')
fig.savefig('p2_ex10_1.pdf')
plt.legend(ncol=1, loc='best', fancybox=False, frameon = True)

<matplotlib.legend.Legend at 0x31b16aa90>

g.dtype, gg.dtype

(dtype('complex128'), dtype('float64'))

fig = plt.figure() # グラフ領域の作成
plt.stem(kx,g.real,'o',label='FFT')
plt.plot(kkx,gg,'-',color='tab:orange',label='sinc')
plt.title('wavenumber domain')
plt.xlabel(r'$k_x$')
plt.ylabel(r'$g(k_x)$')
fig.savefig('p2_ex10_2.pdf')
plt.legend(ncol=1, loc='best', fancybox=False, frameon = True)

<matplotlib.legend.Legend at 0x31c8d6950>

高速フーリエ変換¶