import mpmath # Python の任意精度浮動小数点演算パッケージ
import numpy as np
from scipy.special import factorial
from scipy.signal import find_peaks
from scipy import optimize
import matplotlib.pyplot as plt
np.set_printoptions(precision=4)
CCC = 299.79
qq = 180.0/np.pi
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]:11.8f}',end=' ') for j in range(n1)]
print("")
return
階乗の計算は，
def hfactorial(n):
k = 1
for i in range(1,n+1):
k = k*i
k = np.array(k,dtype=np.float128)
return k
$ n!/m!$ の計算は，
def hfactorial3(n,m): # n!/m!
k = 1
if n>m:
for i in range(m+1,n+1):
k = k*i
k = np.array(k,dtype=np.float128)
elif n<m:
for i in range(n+1,m+1):
k = k*i
k = np.array(k,dtype=np.float128)
k = 1.0/k
return k
ドルフ・チェビシェフアレーの一般式よりアレー素子の励振振幅係数は，
def dolph_tschebyscheff_array(n,R_db):
R = 10**(R_db/20)
R2 = np.sqrt(R**2-1)
nm1 = n-1
x0 = ((R+R2)**(1/nm1)+(R-R2)**(1/nm1))/2
if n == 1:
epn = 2
else:
epn = 1
if np.mod(n,2) == 0:
M = n//2
an = np.zeros(M)
aa = np.zeros((M+1,M+1),dtype=np.float128)
for n in range(1,M+1):
for q in range(n,M+1):
aa[n,q] = hfactorial3(q+M-2,q+n-1)/hfactorial(q-n)/hfactorial(M-q) * (-1)**(M-q) * x0**(2*q-1) * (2*M-1)
an[n-1] = np.sum(aa[n,:])
else:
M = (n-1)//2
an = np.zeros(M+1)
aa = np.zeros((M+2,M+2),dtype=np.float128)
for n in range(1,M+2):
for q in range(n,M+2):
aa[n,q] = hfactorial3(q+M-2,q+n-2)/epn/hfactorial(q-n)/hfactorial(M-q+1) * (-1)**(M-q+1) * x0**(2*(q-1)) * 2*M
an[n-1] = np.sum(aa[n,:])
return an,aa
フーリエ級数展開は，
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$と角度$\theta$との関係は，
def theta_to_u(th,dd,th0):
akd = 2*np.pi*dd
u = akd * (np.cos(th)-np.cos(th0))
return u
im = 1 # ピーク値およびナル点のプロットの有無, =0:なし，=1:あり
n = 30 # 36 # N elements
dd = 0.435 # Element spacing normalized by wave-length [mm]
R_db = 25.0 # サイドローブレベル（>0） [dB]
if np.mod(n,2) == 0:
ii = n//2
else:
ii = (n+1)//2
an2, aa = dolph_tschebyscheff_array(n,R_db)
print('an2:',an2)
an_ = an2/an2[-1]# エッジ素子の振幅で規格化
print('an_:',an_)
if np.mod(n,2)==0:
an = an_/(2*np.sum(an_))# u=0での指向性関数で規格化（素子数が偶数の場合）
else:
an = an_/(2*np.sum(an_)-an_[0])# u=0での指向性関数で規格化（素子数が奇数の場合）
print('an:',an)
fig = plt.figure() # グラフ領域の作成
plt.ion()
plt.subplot(2,1,1)
if im == 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)
if im == 1:
plt.plot(uoverpi[locs],f[locs],'.')
plt.plot(uoverpi[locs_null],np.zeros(len(locs_null)),'.')
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)
if im == 1:
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)
if im == 1:
plt.plot(thdeg[locs],p[locs],'o')
plt.plot(thdeg[locs_null],np.ones(len(locs_null))*-40,'o')
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$')
if (peak_sl>-80.0)&(peak_sl<0.0):
plt.text(35,-10.0,f'Peak sidelobe level $={peak_sl:.2f}$ dB')
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_im'+str(im)+'.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}")