円形導波管
モジュールのインポート
import numpy as np
from scipy.special import jv, yv # ベッセル関数
from scipy.special import jn_zeros, jnp_zeros, jnjnp_zeros, jnyn_zeros # ベッセル関数の零点
from scipy import optimize
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 uvector2(aa):
norm = np.linalg.norm(aa, axis=1, keepdims=True)
with np.errstate(divide="ignore", invalid="ignore"):
uv = np.where(norm==0, np.zeros(aa.shape[1]), aa / norm)
return norm, uv
円形導波管の遮断波数を求めるため，零点の計算は，
def circular_waveguide_xi2(imode,m, n): # 特性方程式の根
if imode == 1: # TE
xi = jnp_zeros(m,n)
elif imode == -1: # TM
xi = jn_zeros(m,n)
return xi[n-1]
円形導波管の伝搬モードの位相定数，遮断モードの減衰定数は，
def beta_alpha(sk, kc): # 位相定数，減衰定数
beta, wl_g, alpha = 0.0, 0.0, 0.0
kz2 = sk**2-kc**2
beta = np.where(sk>kc, np.sqrt(kz2), 0.0)
alpha = np.where(sk<=kc, np.sqrt(-kz2), 0.0)
return beta, alpha
円形導波管のモード関数（電界および磁界の $x$ 成分，$y$ 成分）の計算は，
# imode = 1 (TE), -1 (TM)
# ic = 1(sin-mode), 2(cos-mode)
def circular_waveguide_mode(rho,phi, a,imode,m,n_xi,xi,ic): # 円形導波管のモード関数
xia = xi/a
if m==0:
em = 1.0
if imode==1:# TE
ic = 2
elif imode==-1:# TM
ic = 1
else:
em = 2.0
if imode==1:# TE
besjmx = jv(m,xi)
amn = np.sqrt(em/np.pi/(xi*xi-m*m)) / np.abs(besjmx)
elif imode==-1:# TM
besjmx = jv(m+1,xi)
amn = np.sqrt(em/np.pi) /xi /np.abs(besjmx)
phim = m*phi
sp = np.sin(phim)
cp = np.cos(phim)
x = xia * rho
besjmm = np.where(np.abs(rho)>a, 0.0, jv(m-1,x))
besjmp = np.where(np.abs(rho)>a, 0.0, jv(m+1,x))
bmn = amn*xi/2.0/a
e_rho = (besjmm + imode*besjmp) * bmn
e_phi = (-besjmm + imode*besjmp) * bmn
if ic==1:
e_rho = e_rho * cp
e_phi = e_phi * sp
elif ic==2:
e_rho = e_rho * sp
e_phi = -e_phi * cp
h_rho = -e_phi
h_phi = e_rho
return e_rho, e_phi
位置ベクトルの極座標$(\rho, \phi)$から直角座標$(x,y)$への変換は，
def rhophi_xy(e_rho,e_phi, phi):
sp = np.sin(phi)
cp = np.cos(phi)
e_x = e_rho*cp - e_phi*sp
e_y = e_rho*sp + e_phi*cp
return e_x, e_y
電界のモード関数の極座標成分 $e_\rho$，$e_\phi$ から電界の直角座標成分 $e_x$，$e_y$，磁界の直角座標成分 $h_x$，$h_y$ への変換は，
def rhophi_xy_eh(e_rho,e_phi, phi): #
sp = np.sin(phi)
cp = np.cos(phi)
e_x = e_rho*cp - e_phi*sp
e_y = e_rho*sp + e_phi*cp
h_x = -e_y
h_y = e_x
return e_x, e_y, h_x, h_y
2次元データからカット面データを切り取るため，
def two_one_dim(aa): # 2次元データからカット面データの切り取り revised
co_h = aa[::-1,n_phi//2]
co_h = np.append(co_h, aa[:,0])
co_e = aa[::-1,n_phi*3//4]
co_e = np.append(co_e, aa[:,n_phi//4])
co_45 = aa[::-1,n_phi*5//8]
co_45 = np.append(co_45, aa[:,n_phi//8])
with np.errstate(divide="ignore", invalid="ignore"):
co_e_db = 20.0*np.log10(np.abs(co_e))
co_h_db = 20.0*np.log10(np.abs(co_h))
co_45_db = 20.0*np.log10(np.abs(co_45))
return co_e, co_h, co_45, co_e_db, co_h_db, co_45_db
円形導波管のモードの遮断特性
まず，TE$_{1n}$モードに対しては，ベッセル関数の導関数 $J_1'(x)$ の零点，つまり$J_1'(x) =0$を満たす解を求め，
n = 5
chi_TE = jnp_zeros(1,n) # TE-mode
chi_TE
array([ 1.84118378, 5.33144277, 8.53631637, 11.7060049 , 14.86358863])
一方，TM$_{1n}$モードに対しては，ベッセル関数 $J_1(x)$ の零点，つまり$J_1(x) =0$を満たす解を求め，
chi_TM = jn_zeros(1,n) # TM-mode
chi_TM
array([ 3.83170597, 7.01558667, 10.17346814, 13.32369194, 16.47063005])
TE$_{1n}$およびTM$_{1n}$モードを遮断周波数の低い順に並べて，
n2 = n*2
chi = np.empty(n2)
chi[0::2] = chi_TE[:]
chi[1::2] = chi_TM[:]
chi
array([ 1.84118378, 3.83170597, 5.33144277, 7.01558667, 8.53631637,
10.17346814, 11.7060049 , 13.32369194, 14.86358863, 16.47063005])
モードの次数などを配列に設定して，
mm = np.ones(n2, dtype=np.int8) # モードの次数 m
i = np.arange(n, dtype=np.int8)
nn = np.empty(n2, dtype=np.int8) # モードの次数 n
nn[0::2] = i[:]+1
nn[1::2] = i[:]+1
mode = np.empty(n2, dtype=object)
mode[0::2] = 'TE'
mode[1::2] = 'TM'
imode = np.empty(n2, dtype=np.int8)
imode[mode=='TE'] = 1
imode[mode=='TM'] = -1
mm, nn, mode, imode
(array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int8),
array([1, 1, 2, 2, 3, 3, 4, 4, 5, 5], dtype=int8),
array(['TE', 'TM', 'TE', 'TM', 'TE', 'TM', 'TE', 'TM', 'TE', 'TM'],
dtype=object),
array([ 1, -1, 1, -1, 1, -1, 1, -1, 1, -1], dtype=int8))
X-bandで用いられる円形導波管の直径を，
D = 23.4 # X-bandで用いられる円形導波管の直径
a = D/2 # 円形の半径 [mm]
遮断波数は，
k_c = chi/a # 遮断波数
k_c
array([0.15736614, 0.32749624, 0.45567887, 0.59962279, 0.72959969,
0.86952719, 1.00051324, 1.13877709, 1.27039219, 1.40774616])
遮断波長は，
wl_c = 2*np.pi/k_c # 遮断波長
wl_c
array([39.92717557, 19.1855191 , 13.78862556, 10.47856317, 8.61182563,
7.2259791 , 6.27996218, 5.51748483, 4.94586267, 4.46329423])
遮断周波数は，
CCC = 299.79
f_c = CCC/wl_c # 遮断周波数
f_c
array([ 7.50841991, 15.62584772, 21.74183342, 28.6098385 , 34.81143404,
41.48780338, 47.73754862, 54.33454 , 60.61429932, 67.16787746])
フォートマットを指定して出力すると，
print("a =",a,'[mm]')
print(f"{'No.':>7s}",f"{'mode':>7s}",f"{'chi':>10s}",f"{'k_c':>10s}",f"{'wl_c[mm]':>10s}",f"{'f_c[GHz]':>10s}")
for i in range(0,n2):
print(f'{i+1:6.0f} {mode[i]:2s} {mm[i]:2d} {nn[i]:2d} \
{chi[i]:10.6f} {k_c[i]:10.6f} {wl_c[i]:10.6f} {f_c[i]:10.6f}')
a = 11.7 [mm]
No. mode chi k_c wl_c[mm] f_c[GHz]
1 TE 1 1 1.841184 0.157366 39.927176 7.508420
2 TM 1 1 3.831706 0.327496 19.185519 15.625848
3 TE 1 2 5.331443 0.455679 13.788626 21.741833
4 TM 1 2 7.015587 0.599623 10.478563 28.609838
5 TE 1 3 8.536316 0.729600 8.611826 34.811434
6 TM 1 3 10.173468 0.869527 7.225979 41.487803
7 TE 1 4 11.706005 1.000513 6.279962 47.737549
8 TM 1 4 13.323692 1.138777 5.517485 54.334540
9 TE 1 5 14.863589 1.270392 4.945863 60.614299
10 TM 1 5 16.470630 1.407746 4.463294 67.167877
分散特性のデータを計算して，
f_stop = 50.0
nf = 31
fff = np.empty((n2,nf))
beta_f = np.empty((n2,nf))
alpha_f = np.empty((n2,nf))
for i in range(n2):
ff = np.array([f_c[i]+(f_stop-f_c[i])*(j/(nf-1))**4 for j in range(nf)]) # リスト内包表記
fff[i,:] = ff
omega = 2*np.pi*fff
wwl = CCC/ff # 自由空間波長[mm]
ssk = 2*np.pi/wwl # 自由空間波数
beta_f[i,:] = np.sqrt(ssk**2-k_c[i]**2)
k_0 = np.linspace(0.0,1.0,101)
omega_0 = k_0*CCC
プロットすると，
fig = plt.figure(figsize=(8, 6))
plt.text( 0.35, 10, f"$a=${a:.1f} mm", fontsize=16)
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # 補助目盛りをつける
plt.xlabel(r"$\beta$") # x軸のラベル設定
plt.ylabel(r"$\omega$") # y軸のラベル設定
plt.xlim(0, 0.7) # y軸範囲の設定
plt.ylim(0, 200.0) # y軸範囲の設定
for i in range(4):
plt.plot(beta_f[i,:], omega[i,:], label=mode[i]+str(mm[i])+str(nn[i]))
plt.plot(k_0, omega_0, '--', label='light')
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon=True)
fig.savefig('p3_tlt_s13_1.pdf')
plt.show()
fig = plt.figure(figsize=(8, 6))
plt.text( 31.0, 0.9, f"$a=${a:.1f} mm", fontsize=16)
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlabel(r"$f$ [GHz]") # x軸のラベル設定
plt.ylabel(r"$\beta/k_0$") # y軸のラベル設定
plt.xlim(5, 40.0) # y軸範囲の設定
plt.ylim(0, 1.0) # y軸範囲の設定
for i in range(5):
plt.plot(fff[i,:], beta_f[i,:]/ssk, label=mode[i]+str(mm[i])+str(nn[i]))
plt.legend(ncol=1, loc='upper left', fancybox=False, frameon=True)
fig.savefig('p3_tlt_s13_2.pdf')
plt.show()
モードの分布
計算するモードを設定して，
imode, m, n = 1, 1, 1
if imode==1:
s_mode = 'TE_'+ str(m) +'_' + str(n)
t_mode = 'TE '+ str(m) +',' + str(n)
else:
s_mode = 'TM_'+ str(m) +'_' + str(n)
t_mode = 'TM '+ str(m) +',' + str(n)
xi = circular_waveguide_xi2(imode, m, n)
s_mode, t_mode, xi
('TE_1_1', 'TE 1,1', 1.8411837813406593)
モードの分布を計算して，
qq = 180.0/np.pi
rho_min, rho_max, n_rho = 0.001, 1.0, 61
phi_min, phi_max, n_phi = 0.0, 2.0*np.pi, 361
rho = np.linspace(rho_min, rho_max, n_rho)
phi = np.linspace(phi_min, phi_max, n_phi)
pphi, rrho = np.meshgrid(phi,rho, indexing='xy') # xy座標の並び (デフォルト)
#ic = 1# x-pol.
ic = 2# y-pol.
a = 1.0
e_rho, e_phi = circular_waveguide_mode(rrho,pphi, a,imode,m,n,xi,ic)
e_x, e_y, h_x, h_y = rhophi_xy_eh(e_rho,e_phi, pphi)
emax = np.array([e_x.max(), e_y.max(), e_rho.max(), e_phi.max()])
emaxx = emax.max()
e_x = e_x/emaxx
e_y = e_y/emaxx
e_rho = e_rho/emaxx
e_phi = e_phi/emaxx
e_x.max(), e_y.max()
カラーマップ名をリストに予め設定して，
Diverging_cmaps = ['PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu', 'RdYlBu',\
'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic'] # 0-11
電界のモード関数の極座標成分 $e_\rho$，$e_\phi$，電界の直角座標成分 $e_x$，$e_y$ をプロットして，
fig = plt.figure(figsize=(12, 8)) # グラフ領域の作成
v = np.linspace(-1.0, 1.0, 101)
mycolor1 = Diverging_cmaps[1]# *0（緑，ピンク），1（緑，紫）
mycolor2 = Diverging_cmaps[5]# *3（紫，茶）5, 8,10,11（赤，青）
irho, iphi = 1, 1
#irho, iphi = 0, 0
if irho == 1:
ax = plt.subplot(2,2,1, polar=True)
ax.set_rgrids([0.2, 0.4, 0.6, 0.8, 1], angle=325., fontname="Arial", fontsize=9)
ax = plt.contour(pphi, rrho, e_rho, colors='gray', linewidths=0.4)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
ax = plt.contourf(pphi, rrho, e_rho, levels=v, cmap=mycolor1)
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="vertical")
ax.set_label(r"Relative amplitude $e_\rho$", fontname="Arial", fontsize=9)
if iphi == 1:
ax = plt.subplot(2,2,2, polar=True)
ax.set_rgrids([0.2, 0.4, 0.6, 0.8, 1], angle=325., fontname="Arial", fontsize=9)
ax = plt.contour(pphi, rrho, e_phi, colors='gray', linewidths=0.4)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
ax = plt.contourf(pphi, rrho, e_phi, levels=v, cmap=mycolor1)
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="vertical")
ax.set_label(r"Relative amplitude $e_\phi$", fontname="Arial", fontsize=9)
if irho == 0 and iphi ==0:
ax = plt.subplot(1,2,1, polar=True)
else:
ax = plt.subplot(2,2,3, polar=True)
ax.set_rgrids([0.2, 0.4, 0.6, 0.8, 1], angle=325., fontname="Arial", fontsize=9)
ax = plt.contour(pphi, rrho, e_x, colors='gray', linewidths=0.4)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
ax = plt.contourf(pphi, rrho, e_x, levels=v, cmap=mycolor2)
if irho == 0 and iphi ==0:
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="horizontal")
else:
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="vertical")
ax.set_label(r"Relative amplitude $e_x$", fontname="Arial", fontsize=9)
if irho == 0 and iphi ==0:
ax = plt.subplot(1,2,2, polar=True)
else:
ax = plt.subplot(2,2,4, polar=True)
ax.set_rgrids([0.2, 0.4, 0.6, 0.8, 1], angle=325., fontname="Arial", fontsize=9)
ax = plt.contour(pphi, rrho, e_y, colors='gray', linewidths=0.4)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
ax = plt.contourf(pphi, rrho, e_y, levels=v, cmap=mycolor2)
if irho == 0 and iphi ==0:
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="horizontal")
else:
ax = plt.colorbar(pad=0.1, shrink=1.0, orientation="vertical")
ax.set_label(r"Relative amplitude $e_y$", fontname="Arial", fontsize=9)
plt.tight_layout()
fig.savefig('p3_tlt_s13_3_'+s_mode+'.pdf')
plt.show()
#fig = plt.figure(figsize=(12,4))
fig, ax = plt.subplots(1,2,figsize=(12, 4)) # グラフ領域の作成
v = np.linspace(-1.0, 1.0, 101)
mycolor1 = Diverging_cmaps[1]# *0（緑，ピンク），1（緑，紫）
mycolor2 = Diverging_cmaps[5]# *3（紫，茶）5, 8,10,11（赤，青）
ax[0] = plt.subplot(1,2,1, polar=True)
CS1 = ax[0].contour(pphi, rrho, h_x, colors='gray', linewidths=0.4)# 等高線表示
ax[0].clabel(CS1, fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
CS2 = ax[0].contourf(pphi, rrho, h_x, levels=v, cmap=mycolor2)
#ax = fig.colorbar(CS2, ax=ax[0], pad=0.1, shrink=1.0, orientation="horizontal")
cbar = fig.colorbar(CS2, ax=ax[0], pad=0.1, shrink=1.0, orientation="vertical")
cbar.set_label(r"Relative amplitude $h_x$", fontname="Arial", fontsize=12)
ax[1] = plt.subplot(1,2,2, polar=True)
CS3 = ax[1].contour(pphi, rrho, h_y, colors='gray', linewidths=0.4)# 等高線表示
ax[1].clabel(CS3, fmt='%1.1f', inline=1, fontsize=9)# 等高線の値を表示
CS4 = ax[1].contourf(pphi, rrho, h_y, levels=v, cmap=mycolor2)
#ax = fig.colorbar(CS4, ax=ax[1], pad=0.1, shrink=1.0, orientation="horizontal")
cbar = fig.colorbar(CS4, ax=ax[1], pad=0.1, shrink=1.0, orientation="vertical")
cbar.set_label(r"Relative amplitude $h_y$", fontname="Arial", fontsize=12)
[i.set_rgrids([0.2, 0.4, 0.6, 0.8, 1], angle=325., fontname="Arial", fontsize=9) for i in ax]
fig.tight_layout()
fig.savefig('p3_tlt_s13_5_'+s_mode+'.pdf')
plt.show()
rho2 = -rho[::-1]
rho2 = np.append(rho2, rho)
co_e, co_h, co_45, co_e_db, co_h_db, co_45_db = two_one_dim(e_x)
xx_e, xx_h, xx_45, xx_e_db, xx_h_db, xx_45_db = two_one_dim(e_y)
カット面での電界の直角座標成分 $e_x$，$e_y$ のデシベル値をプロットして，
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 4)) # グラフ領域の作成，nrows=縦に並べる数，ncols=横に並べる数
ax[0].set_title("$E_x$ ("+t_mode+")")
ax[0].set_xlabel(r"Normalized radius $\rho/a$", fontsize=10) # x軸のラベル設定
[ax[j].set_ylim(-40.0, 10.0) for j in range(2)]
[ax[j].set_ylabel("Relative amplitude [dB]", fontsize=10) for j in range(2)]
ax[0].plot(rho2, co_e_db, label="E-plane")
ax[0].plot(rho2, co_h_db, label="H-plane")
ax[0].plot(rho2, co_45_db, label="45-plane")
ax[1].set_title("$E_y$ ("+t_mode+")")
ax[1].set_xlabel(r"Normalized radius $\rho/a$", fontsize=10) # x軸のラベル設定
ax[1].plot(rho2, xx_e_db, label="E-plane")
ax[1].plot(rho2, xx_h_db, label="H-plane")
ax[1].plot(rho2, xx_45_db, label="45-plane")
[ax[j].grid(color = "gray", linestyle="--") for j in range(2)]
[ax[j].minorticks_on() for j in range(2)]
[ax[j].set_xlim(-1.0, 1.0) for j in range(2)]
[ax[j].set_xlabel(r"Normalized radius $\rho/a$", fontsize=10) for j in range(2)]
[ax[j].legend(ncol=1, loc='lower right', fancybox=False, frameon = True, fontsize=8) for j in range(2)]
fig.tight_layout()
fig.savefig('p3_tlt_s13_7_'+s_mode+'.pdf')
plt.show()
# 直角座標系(x,y)
x_min, x_max, n_x = -1.0, 1.0, 101
y_min, y_max, n_y = -1.0, 1.0, 101
x = np.linspace(x_min, x_max, n_x)
y = np.linspace(y_min, y_max, n_y)
xx, yy = np.meshgrid(x,y, indexing='xy')# xy座標の並び (デフォルト）
rrrho = np.sqrt(xx**2+yy**2)
ppphi = np.arctan2(yy,xx)
e_rho, e_phi = circular_waveguide_mode(rrrho,ppphi, a,imode,m,n,xi,ic)
e_x, e_y = rhophi_xy(e_rho,e_phi, ppphi)
emax = np.array([e_x.max(), e_y.max(), e_rho.max(), e_phi.max()])
emaxx = emax.max()
e_x = e_x/emaxx
e_y = e_y/emaxx
x_c, y_c = np.cos(phi), np.sin(phi)
電界のモード関数より電気力線と等高線図を重ねて描くと，
fig = plt.figure(figsize=(7, 7)) # グラフ領域の作成
v = np.linspace(0.0, 1.0, 21)
plt.title(t_mode)
e_amp = np.sqrt(e_x**2+e_y**2)
lw = 4*e_amp / e_amp.max()
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=0.7, linewidth=lw)
plt.contourf(x, y, e_amp, levels=v, cmap="Blues")# YlOrRd, OrRd
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(x_c, y_c, 'grey')
fig.tight_layout()
fig.savefig('p3_tlt_s13_9_'+s_mode+'.pdf')
plt.show()
前のページに戻る
「目次」のページに戻る