p3_ap_Woodward

ウッドワード・ローソン法

モジュールのインポートや初期設定など

import numpy as np  
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
np.set_printoptions(precision=5)

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の設定

計算条件

波源分布の長さ，ウッドワード・ローソン法における合成するビームの数を設定して，

dd = 10.0# Diameter normalized by wave-length
nbeam = 6# ビームの数
if np.mod(nbeam,2) == 0:# 偶数
    mm = nbeam//2 
else:# 奇数
    mm = (nbeam-1)//2
akover2 = np.pi*dd
th = np.linspace(-np.pi*0.5, np.pi*0.5, 360+1)
thdeg = th*180/np.pi
nth = len(th)
u = akover2 * np.sin(th)
uoverpi = u/np.pi
an = np.ones(mm)
a0 = 1.0
if np.mod(nbeam,2) == 0:# 偶数
    f0 = np.zeros(nth)
else:# 奇数
    f0 = a0*np.sinc(uoverpi)

開口面分布および放射パターン

ビームを合成して得られる放射パターン，そのときの波源分布を求めると，

xx = np.linspace(-1.0, 1.0, 100+1)# 波源位置座標
if np.mod(nbeam,2) == 0:# 偶数
    af = 0.0
else:# 奇数
    af = 1.0
fp = np.zeros([mm,nth])
fm = np.zeros([mm,nth])
f = f0
for i in range(mm):
    if np.mod(nbeam,2) == 0:# 偶数
        um = akover2*(2*(i+1)-1)/2/dd
    else:# 奇数
        um = akover2*(i+1)/dd
    umpi = um/np.pi
    ump = uoverpi+umpi
    umm = uoverpi-umpi
    fp[i,:] = an[i]*np.sinc(ump)
    fm[i,:] = an[i]*np.sinc(umm)
    f = f+fp[i,:]+fm[i,:]
    af = af+ 2.0*np.cos(um*xx)
np.max(af)
afn = af/np.max(af)
p = 20.0*np.log10(np.abs(f))# 放射パターン[dB]

波源分布の計算値の一部を出力すると，

print(f"{'2x/a':>10s} {'e(x)':>10s}")
for i in range(0,len(xx),10):
    print(f"{xx[i]:10.5f} {afn[i]:10.5f}")

      2x/a       e(x)
  -1.00000    0.00000
  -0.80000    0.16667
  -0.60000   -0.12109
  -0.40000   -0.16667
  -0.20000    0.51295
   0.00000    1.00000
   0.20000    0.51295
   0.40000   -0.16667
   0.60000   -0.12109
   0.80000    0.16667
   1.00000    0.00000

また，放射パターンの計算値の一部を出力して，

print(f"{'th[deg]':>10s} {'u/pi':>10s} {'f':>10s} {'p[dB]':>10s}")
for i in range(0,len(uoverpi),20):
    print(f"{thdeg[i]:10.5f} {uoverpi[i]:10.5f} {f[i]:10.5f} {p[i]:10.5f}")

   th[deg]       u/pi          f      p[dB]
 -90.00000  -10.00000   -0.01040  -39.66058
 -80.00000   -9.84808   -0.00955  -40.40017
 -70.00000   -9.39693    0.00379  -48.42470
 -60.00000   -8.66025    0.00689  -43.23515
 -50.00000   -7.66044   -0.00912  -40.80137
 -40.00000   -6.42788   -0.00645  -43.81433
 -30.00000   -5.00000    0.05577  -25.07219
 -20.00000   -3.42020    0.05430  -25.30426
 -10.00000   -1.73648    1.09896    0.81963
   0.00000    0.00000    1.10347    0.85524
  10.00000    1.73648    1.09896    0.81963
  20.00000    3.42020    0.05430  -25.30426
  30.00000    5.00000    0.05577  -25.07219
  40.00000    6.42788   -0.00645  -43.81433
  50.00000    7.66044   -0.00912  -40.80137
  60.00000    8.66025    0.00689  -43.23515
  70.00000    9.39693    0.00379  -48.42470
  80.00000    9.84808   -0.00955  -40.40017
  90.00000   10.00000   -0.01040  -39.66058

放射パターンのピーク値を求めると，

locs, _ = find_peaks(p) # ピーク値
upeak = u[locs]
fpeak = f[locs]
thpeak = thdeg[locs]
ppeak = p[locs]
print(f"{'th[deg]':>10s} {'u/pi':>10s} {'f':>10s} {'p[dB]':>10s}")
for i in range(len(locs)):
    print(f"{thpeak[i]:10.5f} {upeak[i]:10.5f} {fpeak[i]:10.5f} {ppeak[i]:10.5f}")

   th[deg]       u/pi          f      p[dB]
 -86.00000  -31.33940   -0.01042  -39.63994
 -64.00000  -28.23645    0.01314  -37.63000
 -53.00000  -25.08987   -0.01714  -35.32029
 -44.00000  -21.82334    0.02343  -32.60479
 -36.50000  -18.68691   -0.03441  -29.26598
 -29.50000  -15.46994    0.05664  -24.93755
 -23.00000  -12.27518   -0.11911  -18.48131
 -12.00000   -6.53174    1.17123    1.37281
   0.00000    0.00000    1.10347    0.85524
  12.00000    6.53174    1.17123    1.37281
  23.00000   12.27518   -0.11911  -18.48131
  29.50000   15.46994    0.05664  -24.93755
  36.50000   18.68691   -0.03441  -29.26598
  44.00000   21.82334    0.02343  -32.60479
  53.00000   25.08987   -0.01714  -35.32029
  64.00000   28.23645    0.01314  -37.63000
  86.00000   31.33940   -0.01042  -39.63994

これらの結果を図示すると，

fig = plt.figure() # グラフ領域の作成
plt.subplot(2,2,1)
plt.plot(uoverpi,f)
if np.mod(nbeam,2) == 1:# 奇数
    plt.plot(uoverpi,f0,'--',linewidth=1.5)
for i in range(mm):
    plt.plot(uoverpi,fp[i,:],'--',linewidth=1.5)
    plt.plot(uoverpi,fm[i,:],'--',linewidth=1.5)
plt.axis([-10,10,-0.5,1.5])
plt.xlabel(r'$u/\pi$ [rad]')
plt.ylabel('Array factor $f(u)$')

plt.subplot(2,1,2)
plt.plot(thdeg,p)
plt.plot(thpeak,ppeak,'o')
plt.axis([-90, 90, -40, 10])
plt.xlabel(r'Angle $\theta$ [deg]')
plt.ylabel(r'$20 \log_{10} |f(u)|$ [dB]')
plt.xticks(np.linspace(-90,90,7))

plt.subplot(2,2,2)
plt.plot(xx,afn,color='tab:red')
plt.axis([-1, 1, -1, 2])
plt.xlabel(r'$2x/a$')
plt.ylabel(r'$e(x)$')

fig.tight_layout()
fig.savefig('Woodward_Lawson_nbeam'+str(nbeam)+'.pdf')
plt.show()

ビーム数を3として計算し，図示すると次のようになる．

同様にして，ビーム数4の場合，

ビーム数5の場合，