平面波展開によるフィールド変換

モジュールのインポート

import numpy as np
import matplotlib.pyplot as plt
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'])

ユーザ関数

等高線図のための関数を， 
def hd_contour2(sl,sm,sn, xx,yy,zz, xmin, xmax, ymin, ymax, v, s, scmap):
    plt.subplot(sl,sm,sn)
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.contourf(xx, yy, zz, v, alpha=.4, cmap=scmap)
    ax1 = plt.colorbar(aspect=40, pad=0.08, shrink=1.0, orientation="horizontal") # 水平カラーバー
    ax1.set_label(s, fontname="Arial", fontsize=10)
    cs = plt.contour(xx, yy, zz, v, colors='gray') # 等高線の表示
    plt.xlim(xmin, xmax) # x軸範囲の設定
    plt.ylim(ymin, ymax) # y軸範囲の設定
    plt.clabel(cs, fmt='%1.1f',inline=True, fontsize=8)# 等高線レベルの表示
    return
def hd_contour(sl,sm,sn, xx,yy,zz, xmin, xmax, ymin, ymax, v, s, scmap):
    ax = plt.subplot(sl,sm,sn)
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    cs1 = plt.contourf(xx,yy,zz, v, alpha=0.4, cmap=scmap)
    ax1 = plt.colorbar(cs1, pad=0.1, shrink=0.9, orientation="horizontal") # カラーバーを表示
    ax1.set_label(s, fontname="Arial", fontsize=10)
    cs = ax.contour(xx, yy, zz, v, colors='gray') # 等高線表示
    ax.set_xlim(xmin, xmax) # x軸範囲の設定
    ax.set_ylim(ymin, ymax) # y軸範囲の設定
    ax.clabel(cs, fmt='%1.1f', inline=True, fontsize=8) # 等高線の値を表示
    return
dBの計算のために， 
def dbdeg(z): # 電磁界の複素係数から絶対値，デジベル(dB)，偏角(deg)の計算
    q = 180.0/np.pi
    amp = np.abs(z)
    db = 20.0*np.log10(amp)
    rad = np.angle(z)
    rad = np.unwrap(rad, discont=0.95*np.pi)
    deg = rad*q
    return amp,db,deg

高速フーリエ変換

開口面分布を設定して， 
q = 180.0/np.pi
pi2 = 2.0*np.pi
eps = 1.0e-10
span, n = 30.0, 128# FFTの計算条件
diameter = 3.0 # 開口径（正方形の1辺の長さ，あるいは円形の直径）
amp = 1.0 # 方形関数の振幅
ica = 11 # ******** 開口面振幅分布の選定コード ***************
icp = 1 # ********* 開口面位相分布の設定コード ***************
tp = 0.2*np.pi # ** tパラメータ ***************************
f_z = 1.0 # ******* 開口面(z=0)からの波長で規格化した距離 *****
print('データ数 n:',n)
n2 = int(n/2)
データ数 n: 128
波長で規格化した座標を $(x,y)$ とし，空間領域のサンプル点を設定して， 
n_x, n_y = n, n
x_span, y_span = span, span
d_x, d_y = x_span/n_x, y_span/n_y
x_min, y_min = -x_span/2.0, -y_span/2.0
x_max, y_max = -x_min-d_x, -y_min-d_y
x = np.linspace(x_min, x_max, n_x)
y = np.linspace(y_min, y_max, n_y)
print('xmin =',x[0],', xmax =',x[n_x-1],', d_x =',d_x)
print('ymin =',y[0],', ymax =',y[n_y-1],', d_y =',d_y)
xxx, yyy = np.meshgrid(x,y)
xxx = xxx + (xxx==0)*eps
yyy = yyy + (yyy==0)*eps
xmin = -15.0 , xmax = 14.765625 , d_x = 0.234375
ymin = -15.0 , ymax = 14.765625 , d_y = 0.234375
波数 $k$ で規格化した波数ベクトル成分を $(k_x,k_y)$ とし，波数領域のサンプル点を設定して， 
d_kx, d_ky = 1.0/x_span, 1.0/y_span
kx_min, ky_min = -d_kx*n_x/2, -d_ky*n_y/2
kx_max, ky_max = -kx_min-d_kx, -ky_min-d_ky
kx = np.linspace(kx_min, kx_max, n_x)
ky = np.linspace(ky_min, ky_max, n_y)
print('kx_min =',kx[0],', kx_max =',kx[n_x-1],', d_kx =',d_kx)
print('ky_min =',ky[0],', ky_max =',ky[n_y-1],', d_ky =',d_ky)
kkkx, kkky = np.meshgrid(kx,ky)
kkkz2 = 1.0-kkkx**2-kkky**2
kx_min = -2.1333333333333333 , kx_max = 2.1 , d_kx = 0.03333333333333333
ky_min = -2.1333333333333333 , ky_max = 2.1 , d_ky = 0.03333333333333333
kkkz2[kkkz2>0]
array([0.01111111, 0.02555556, 0.03777778, ..., 0.03777778, 0.02555556,
       0.01111111])
kzc = (1+0j)*kkkz2
kzc = kzc**0.5
kzc.real[kkkz2>0.0]
array([0.10540926, 0.15986105, 0.19436506, ..., 0.19436506, 0.15986105,
       0.10540926])
ct = np.where(kkkz2>0.0, kzc.real, 0.0)
st = np.sqrt(1.0-ct*ct)
st = np.where(kkkx<0, -st, st)
th = np.empty(n)
th[:] = np.arcsin(st[n2,:])*q
近傍界（開口面分布）は， 
xr, yr, a = diameter/2.0, diameter/2.0, diameter/2.0# 開口の大きさ（半径）
rrr = np.sqrt(xxx**2+yyy**2)

if ica == 1: # 方形開口（一様分布）の計算
    f = np.where((xxx>=-xr) & (xxx<xr) & (yyy>=-yr) & (yyy<yr), amp, 0.0)
elif ica == 2: # 方形開口（コサインテーパ分布）の計算
    f = np.cos(np.pi*xxx/diameter) * np.cos(np.pi*yyy/diameter)
    f = np.where((xxx>=-xr) & (xxx<xr) & (yyy>=-yr) & (yyy<yr), f, 0.0)
elif ica == 11: # 円形開口（一様分布）の計算
    f = np.where((rrr>=0) & (rrr<=a), amp, 0.0)
elif ica == 12: # 円形開口（放物線テーパ分布）の計算
    f = np.where((rrr>=0) & (rrr<=a), 1.0-(rrr/a)**2, 0.0)

if icp == 1:# 球面波面を表す位相分布
    ph = np.where((rrr>=0) & (rrr<=a), -tp*(rrr/a)**2, 0.0)
f = f*np.exp(1j*ph)

pt = np.sum(np.abs(f)**2)*d_x*d_y
print("pt:",pt)
f = f/np.sqrt(pt)
pt = np.sum(np.abs(f)**2)*d_x*d_y
print("pt:",pt)
pt: 7.086181640625
pt: 0.9999999999999998
iFFTによる近傍界（開口面）ー遠方界（スペクトラム）の計算にあたって， ゼロ周波数の項を配列の中心にシフトする． ただし，2次元入力の場合，第1象限と第3象限，第2象限と第4象限を入れ替える． 
fs = np.fft.ifftshift(f) # データ逆シフト
gs = np.fft.ifft2(fs) # iFFT
g = np.fft.fftshift(gs) # データのシフト
g = g*d_x*d_y*n*n
b = g/pi2/pi2
pt2 = np.sum(np.abs(b)**2)*pi2
print("pt2:",pt2)
print("b[n2,n2]:",b[n2,n2])
b0 = b/b[n2,n2]
ice = 0
if ice == 0:# ユニバーサルパターン
    g0 = b0*1.0
elif ice == 1:# 平面波展開の無限遠での漸近展開
    g0 = b0*ct
elif ice == 2:# 開口面法の平面開口近似
    g0 = b0*(1.0+ct)/2.0
pt2: 3.6282976237349427
b[n2,n2]: (0.06304252764115055-0.02056687051563702j)
最大利得（100% gain）は， 
g100 = 4.0*np.pi*diameter**2
g100_db = 10.0*np.log10(g100)
print("100% gain=",g100,g100_db,"[dB]")
100% gain= 113.09733552923255 20.53452373461421 [dB]
利得は， 
gain = 4.0*np.pi*np.abs(b[n2,n2]*pi2*pi2)**2
gain_db = 10.0*np.log10(gain)
print("gain=",gain,gain_db,"[dB]")
gain= 86.12336870402707 19.351210088541556 [dB]
FFTによるスペクトラム（遠方界）から近傍界への計算は， 
g = g*np.exp(-1j*pi2*f_z*ct)
gs = np.fft.ifftshift(g) # データ逆シフト
fs = np.fft.fft2(gs)# FFT
f_nf = np.fft.fftshift(fs) # データのシフト
f_nf = f_nf/(n*n*d_x*d_y)
f_ph = np.angle(f)*q
g0_ph = np.angle(g0)*q
f_nf_ph = np.angle(f_nf)*q
設定した開口面分布をプロットすると， 
fig = plt.figure(figsize = (10, 8))# Figure
v = np.linspace(0.0, 1.0, 11) # 等高線
hd_contour2(1,2,1, xxx,yyy, np.abs(f),  -2.0*xr, 2.0*xr, -2.0*yr, 2.0*yr, v, 'Amplitude', "Blues")
plt.xlabel(r"$x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$y$", fontsize=14) # y軸のラベル
v = np.linspace(-180.0, 180.0, 13) # 等高線
hd_contour2(1,2,2, xxx,yyy, f_ph, -2.0*xr, 2.0*xr, -2.0*yr, 2.0*yr, v, 'Phase [deg]', "pink")
plt.xlabel(r"$x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$y$", fontsize=14) # y軸のラベル
plt.tight_layout()
fig.savefig('p3_ap_s1_1.pdf')




変換した波数スペクトラムをプロットして，

fig = plt.figure(figsize = (10, 8))# Figure
v = np.linspace(-40.0, 0.0, 9) # 等高線
hd_contour2(1,2,1, kkkx, kkky, 20.0*np.log10(np.abs(g0)), -1.0, 1.0, -1.0, 1.0, v, 'Relative power [dB]', "Blues")
plt.xlabel(r"$k_x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$k_y$", fontsize=14) # y軸のラベル
v = np.linspace(-180.0, 180.0, 13)# 等高線
hd_contour2(1,2,2, kkkx, kkky, g0_ph, -1.0, 1.0, -1.0, 1.0, v, 'Relative phase [deg]', "RdBu")
plt.xlabel(r"$k_x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$k_y$", fontsize=14) # y軸のラベル
plt.tight_layout()
fig.savefig('p3_ap_s1_2.pdf')




近傍界に変換して，

plt.figure(figsize = (10, 8))# Figure
v = np.linspace(0.0, 1.0, 11)# 等高線
hd_contour2(1,2,1, xxx,yyy, np.abs(f_nf),  -2.0*xr, 2.0*xr, -2.0*yr, 2.0*yr, v, 'Amplitude', "Blues")
plt.xlabel(r"$x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$y$", fontsize=14) # y軸のラベル
v = np.linspace(-180.0, 180.0, 13)# 等高線
hd_contour2(1,2,2, xxx,yyy, f_nf_ph, -2.0*xr, 2.0*xr, -2.0*yr, 2.0*yr, v, 'Phase [deg]', "RdBu")
plt.xlabel(r"$x$", fontsize=14) # x軸のラベル
plt.ylabel(r"$y$", fontsize=14) # y軸のラベル
plt.tight_layout()
fig.savefig('p3_ap_s1_3.pdf')




これらのカット面をまとめてプロットすると，

fig = plt.figure(figsize = (10, 8))# Figure

plt.subplot(2, 2, 1) # 左上側の図
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlim(x_min, x_max) # x軸範囲の設定
plt.xticks(fontsize=8)
plt.xlabel(r"$x$", fontsize=8) # x軸のラベル設定
#plt.ylim(y_min, y_max) # y軸範囲の設定
plt.yticks(fontsize=8)
plt.ylabel(r"$f(x,y)$", fontsize=8) # y軸のラベル設定
plt.plot(x, np.abs(f[n2][:]), label="$y=0 (z=0)$")
plt.plot(y, np.abs(f[:][n2]), '.', label="$x=0 (z=0)$")
plt.plot(x, np.abs(f_nf[n2][:]), label="$y=0 (z=d)$")
plt.plot(y, np.abs(f_nf[:][n2]), 'o', label="$x=0 (z=d)$")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True, fontsize=8)

plt.subplot(2, 2, 2) # 左上側の図
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlim(-1.0, 1.0) # x軸範囲の設定
plt.xticks(fontsize=8)
plt.xlabel(r"$k_x$", fontsize=8) # x軸のラベル設定
#plt.ylim(-1.0, 1.0) # y軸範囲の設定
plt.yticks(fontsize=8)
plt.ylabel(r"$g(k_x,k_y)$", fontsize=8) # y軸のラベル設定
plt.plot(kx, np.abs(g0[n2][:]), label="$k_y=0$")
plt.plot(ky, np.abs(g0[:][n2]), '.', label="$k_x=0$")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True, fontsize=8)

plt.subplot(2, 2, 3)
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlim(x_min, x_max) # x軸範囲の設定
plt.xticks(fontsize=8)
plt.xlabel(r"$x$", fontsize=8) # x軸のラベル設定
#plt.ylim(y_min, y_max) # y軸範囲の設定
plt.yticks(fontsize=8)
plt.ylabel(r"$f(x,y)$", fontsize=8) # y軸のラベル設定
plt.plot(x, f_ph[n2][:], label="$y=0 (z=0)$")
plt.plot(y, f_ph[:][n2], label="$x=0 (z=0)$")
plt.plot(x, f_nf_ph[n2][:], '.', label="$y=0 (z=d)$")
plt.plot(y, f_nf_ph[:][n2], '.', label="$x=0 (z=d)$")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True, fontsize=8)

plt.subplot(2, 2, 4)
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlim(-90.0, 90.0) # x軸範囲の設定
plt.xticks(fontsize=8)
plt.xlabel(r"$\theta$ [deg]", fontsize=8) # x軸のラベル設定
plt.ylim(-40.0, 0.0) # y軸範囲の設定
plt.yticks(fontsize=8)
plt.ylabel(r"$g(\theta)$ [dB]", fontsize=8) # y軸のラベル設定
plt.plot(th, 20.0*np.log10(np.abs(g0[n2][:])), label="$\phi=0^\circ$")
plt.plot(th, 20.0*np.log10(np.abs(g0[:][n2])), '.', label="$\phi=90^\circ$")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True, fontsize=8)
fig.savefig('p3_ap_s1_4.pdf')