２項アレーのアレーファクタ
モジュールのインポートや初期設定など
import numpy as np
from scipy.special import comb
from scipy.interpolate import interp1d
from scipy import optimize
import matplotlib.pyplot as plt
np.set_printoptions(precision=4)
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'])
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
ユーザ関数
任意の数の配列をフォーマット付き出力するため，
def table(*args):
n1, n2 = np.shape(args)
for i in range(n2):
[print(f'{args[j][i]:8.3f}',end=' ') for j in range(n1)]
print("")
return
2項アレー素子の励振振幅は，
def binominal_array(n):
k = np.linspace(0,n-1,n, dtype='int')
cnk = comb(n-1, k, exact=False)
return cnk
アレーファクタをフーリエ級数展開によって表すと，
def afactor_fourier(n,an,u):
f = 0
if np.mod(n,2)==0: # 偶数
for i0 in range(n//2):
i = i0+1
f = f+an[i0]*np.cos((2*i-1)/2*u)
f = 2.0*f
else: # 奇数
for i0 in range((n-1)//2):
i = i0+1
f = f+an[i0+1]*np.cos(i*u)
f = an[0]+2*f
return f
パラーメータ $u$ は，
def theta_to_u(th,dd,th0):
akd = 2*np.pi*dd
u = akd * (np.cos(th)-np.cos(th0))
return u
ビーム幅を数値的に求めるための評価関数は，
def func_bw2(u,n,an,bw):
return np.abs(afactor_fourier(n,an,u))-bw
励振振幅係数
2項アレー素子の励振振幅を素子数20まで求めると（パスカルの三角形となることが確認できる），
nb=19
for i in range(2,nb+2):
an2 = binominal_array(i)
print(f"{i:2d}",end=':')
[print(end=' ') for k in range(nb+1-i)]
[print(f"{an2[j]:5.0f} ",end='') for j in range(i)]
print()
2: 1 1
3: 1 2 1
4: 1 3 3 1
5: 1 4 6 4 1
6: 1 5 10 10 5 1
7: 1 6 15 20 15 6 1
8: 1 7 21 35 35 21 7 1
9: 1 8 28 56 70 56 28 8 1
10: 1 9 36 84 126 126 84 36 9 1
11: 1 10 45 120 210 252 210 120 45 10 1
12: 1 11 55 165 330 462 462 330 165 55 11 1
13: 1 12 66 220 495 792 924 792 495 220 66 12 1
14: 1 13 78 286 715 1287 1716 1716 1287 715 286 78 13 1
15: 1 14 91 364 1001 2002 3003 3432 3003 2002 1001 364 91 14 1
16: 1 15 105 455 1365 3003 5005 6435 6435 5005 3003 1365 455 105 15 1
17: 1 16 120 560 1820 4368 8008 11440 12870 11440 8008 4368 1820 560 120 16 1
18: 1 17 136 680 2380 6188 12376 19448 24310 24310 19448 12376 6188 2380 680 136 17 1
19: 1 18 153 816 3060 8568 18564 31824 43758 48620 43758 31824 18564 8568 3060 816 153 18 1
20: 1 19 171 969 3876 11628 27132 50388 75582 92378 92378 75582 50388 27132 11628 3876 969 171 19 1
素子間隔 $d = 0.435\lambda$ の15素子アレーについて，
n = 15 # 36 # N elements
dd = 0.435 # Element spacing normalized by wave-length [mm]
この条件より励振振幅を求めるため，
cnk = binominal_array(n)
cnk_total = np.sum(cnk)
cn = cnk/cnk_total
if np.mod(n,2) == 0:
ii = n//2
else:
ii = (n+1)//2
bn = cn[:ii]
an = bn[::-1] # 配列要素の順番の入れ替え
求めた励振振幅を最大値で規格化してプロットし，スプライン補間によって滑らかに点をつないで図示すると，
ann = an/np.max(an)
ne = np.linspace(0,len(an)-1,len(an),dtype='int')
spline = interp1d(ne,ann, kind='cubic') # ３次スプライン補間
ne_dense = np.linspace(ne[0],ne[-1],101)
ann_dense = spline(ne_dense)
fig = plt.figure() # グラフ領域の作成
plt.plot(ne,an/np.max(an),'o',color='tab:blue')
plt.plot(-ne,an/np.max(an),'o',color='tab:blue')
plt.plot(ne_dense, ann_dense, color='tab:blue')
plt.plot(-ne_dense, ann_dense, color='tab:blue')
plt.xlabel(r'$i$-th elements')
plt.ylabel('Amplitude')
fig.tight_layout()
fig.savefig('array_elements_N'+str(n)+'_d'+str(dd)+'wl.pdf')
plt.show()
共相励振した放射パターンの計算
ビームの方向を60度にするよう励振位相を設定すると，
th0deg = 60 # Beam direction [deg]
akd = 2*np.pi*dd
th0 = th0deg*np.pi/180
n_th = 180001
th = np.linspace(0, np.pi, n_th)
thdeg = th*180/np.pi
u = akd * (np.cos(th)-np.cos(th0))
nu = np.arange(n_th)
uoverpi = u/np.pi
nu_th0 = np.int64(np.round((n_th-1)*th0deg/180.0))
放射パターンを求めて，
f = afactor_fourier(n,an,u)
p = 20.0*np.log10(np.abs(f))
放射パターンからセカント法によりビーム幅を求めると，
qq = 180/np.pi
bw = np.sqrt(0.5)
fbw = func_bw2(u,n,an,bw)
nu_zero = nu[fbw==0.0]
fbw2 = fbw[:-1]*fbw[1:]
nu_zero = np.append(nu_zero,nu[:-1][fbw2<0.0])
if len(nu_zero[nu_zero<nu_th0])==0:
nu_zero_m = -1
else:
nu_zero_m = nu_zero[nu_zero<nu_th0][-1]
if len(nu_zero[nu_zero>nu_th0])==0:
nu_zero_p = -1
else:
nu_zero_p = nu_zero[nu_zero>nu_th0][0]
th_bw3_p, th_bw3_m = 0.0, 0.0
if nu_zero_p!=-1:
up1 = np.min([u[nu_zero_p],u[nu_zero_p+1]])
up2 = np.max([u[nu_zero_p],u[nu_zero_p+1]])
u_zero_p = optimize.brentq(func_bw2, up1,up2, args=(n,an,bw))
th_bw3_p = np.arccos(u_zero_p/akd +np.cos(th0))
if nu_zero_m!=-1:
um1 = np.min([u[nu_zero_m],u[nu_zero_m+1]])
um2 = np.max([u[nu_zero_m],u[nu_zero_m+1]])
u_zero_m = optimize.brentq(func_bw2, um1,um2, args=(n,an,bw))
th_bw3_m = np.arccos(u_zero_m/akd +np.cos(th0))
if (nu_zero_p!=-1)&(nu_zero_m==-1):
if th0deg<90.0:
bw3_deg = 2*np.abs(th_bw3_p)*qq
else:
bw3_deg = 2*(180.0-np.abs(th_bw3_p)*qq)
elif (nu_zero_p==-1)&(nu_zero_m!=-1):
if th0deg<90.0:
bw3_deg = 2*np.abs(th_bw3_m)*qq
else:
bw3_deg = 2*(180.0-np.abs(th_bw3_m)*qq)
elif (nu_zero_p!=-1)&(nu_zero_m!=-1):
bw3_deg = np.abs(th_bw3_m-th_bw3_p)*qq # 3-dBビーム幅
else:
bw3_deg = 0.0
放射パターンを図示して，
fig = plt.figure() # グラフ領域の作成
plt.ion()
plt.subplot(2,1,1)
plt.plot([0,0],[-1,1],'--',color='gray',alpha=0.5)
plt.plot([-4,4],[0,0],'--',color='gray',alpha=0.5)
plt.plot(uoverpi,f)
plt.axis([-4,4,-1,1])
plt.xlabel(r'$u/\pi = 2d ( \cos \theta - \cos \theta_0)/\lambda$')
plt.ylabel('Array factor $f(u)$')
plt.text(0.8,-0.3,'$N=$'+str(n))
plt.text(0.8,-0.55,'Element spacing: $d = $'+str(dd)+r'$\lambda$')
plt.text(0.8,-0.8,r'Beam direction: $\theta_0 = $'+str(th0deg)+r'$^\circ$')
plt.subplot(2,1,2)
plt.plot(np.ones(2)*th0deg,[-50,0],'--',color='gray',alpha=0.5)
plt.plot([0,180],np.ones(2)*(-3),'--',color='gray',alpha=0.5)
if bw3_deg!=0.0:
plt.plot(np.ones(2)*th_bw3_p*qq,[-50,0],'--',color='gray',alpha=0.5)
plt.plot(np.ones(2)*th_bw3_m*qq,[-50,0],'--',color='gray',alpha=0.5)
plt.plot(thdeg,p)
plt.axis([0, 180, -40, 0])
plt.xlabel(r'Angle $\theta$ [deg]')
plt.ylabel(r'$20 \log_{10} |f(\theta)|$ [dB]')
plt.text(th0deg-7,2.5,r'$\theta_0 = $'+str(th0deg)+r'$^\circ$')
if bw3_deg!=0.0:
plt.text(120,-7,f'3-dB beamwidth $={bw3_deg:.1f}^\circ$')
plt.xticks(np.linspace(0,180,5))
fig.tight_layout()
fig.savefig('array_factor_N'+str(n)+'_beam'+str(th0deg)+'deg_d'+str(dd)+'wl.pdf')
plt.show()
print(f"{'u/pi[rad]':>10s} {'th[deg]':>10s}"\
f"{'f(u)':>10s} {'f(th)[dB]':>10s}")
for i in range(0,len(th),5000):
print(f"{uoverpi[i]:10.5f} {thdeg[i]:10.5f}"\
f"{f[i]:10.5f} {p[i]:10.5f}")
u/pi[rad] th[deg] f(u) f(th)[dB]
0.43500 0.00000 0.02845 -30.91776
0.43169 5.00000 0.03018 -30.40566
0.42178 10.00000 0.03587 -28.90546
0.40536 15.00000 0.04720 -26.52169
0.38253 20.00000 0.06748 -23.41678
0.35349 25.00000 0.10237 -19.79635
0.31844 30.00000 0.16045 -15.89316
0.27766 35.00000 0.25257 -11.95252
0.23146 40.00000 0.38812 -8.22066
0.18018 45.00000 0.56643 -4.93701
0.12423 50.00000 0.76472 -2.32995
0.06401 55.00000 0.93156 -0.61575
-0.00000 60.00000 1.00000 0.00000
-0.06732 65.00000 0.92457 -0.68121
-0.13744 70.00000 0.71976 -2.85625
-0.20983 75.00000 0.46088 -6.72823
-0.28393 80.00000 0.23666 -12.51736
-0.35917 85.00000 0.09466 -20.47695
-0.43500 90.00000 0.02845 -30.91776
-0.51083 95.00000 0.00613 -44.24757
-0.58607 100.00000 0.00089 -61.04008
-0.66017 105.00000 0.00008 -82.16488
-0.73256 110.00000 0.00000 -109.05998
-0.80268 115.00000 0.00000 -144.39057
-0.87000 120.00000 0.00000 -194.02859
-0.93401 125.00000 0.00000 -275.84947
-0.99423 130.00000 0.00000 -325.11240
-1.05018 135.00000 0.00000 -308.85413
-1.10146 140.00000 0.00000 -223.84070
-1.14766 145.00000 0.00000 -178.78328
-1.18844 150.00000 0.00000 -149.81760
-1.22349 155.00000 0.00000 -129.80341
-1.25253 160.00000 0.00000 -115.64354
-1.27536 165.00000 0.00001 -105.73169
-1.29178 170.00000 0.00001 -99.15742
-1.30169 175.00000 0.00002 -95.39565
-1.30500 180.00000 0.00002 -94.17050
前のページに戻る
「目次」のページに戻る