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 =1for i in range(1,n+1):
k = k*i
k = np.array(k,dtype=np.float128)return k
の計算は，
def hfactorial3(n,m):# n!/m!
k =1if 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))/2if n ==1:
epn =2else:
epn =1if 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 =0if 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
パラメータと角度との関係は，
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)
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}")