Created
April 25, 2013 14:28
-
-
Save ryujimiya/5460127 to your computer and use it in GitHub Desktop.
WE-FDTD TMz(TE mode) Mur ABC. This script calculates field Ez of a dielectrc resonator in a hollow waveguide.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
########################################################### | |
# TE波(TMz) 導波管誘電体共振器 (WE-FDTD Murの吸収境界条件) | |
# | |
# coded by ryujimiya April 2013 | |
########################################################### | |
function [ndivX, ndivY, distanceX, distanceY, ... | |
qzz, ... | |
Ez, timeAry, ... | |
EzTimePort1, EzTimePort2] = we_fdtd_2d_TMz_dielectric_filter(loopCnt, courantNumber = 0.5, typeGaussian = 0, normalizedFreqCarrier = 1.5) | |
#--------------------------------------------- | |
# 定数 | |
#--------------------------------------------- | |
eps0 = 8.8541878e-12; | |
c0 = 2.99792458e8; | |
mu0 = 4.0e-7 * pi; | |
#--------------------------------------------- | |
# 領域 | |
#--------------------------------------------- | |
#ndivPerUnit = 4; | |
ndivPerUnit = 4; | |
ndivX1 = 50 * ndivPerUnit; | |
ndivX2 = 19 * ndivPerUnit; | |
ndivX3 = 50 * ndivPerUnit; | |
ndivX = ndivX1 + ndivX2 + ndivX3; | |
ndivY1 = 20 * ndivPerUnit; | |
ndivY2 = 10 * ndivPerUnit; | |
ndivY = ndivY1; | |
ndivYOfs = (ndivY1 - ndivY2) / 2; | |
distanceX = 3.0e-3 * ndivX; | |
distanceY = 3.0e-3 * ndivY; | |
#--------------------------------------------- | |
# 分割幅の逆数 | |
#--------------------------------------------- | |
dx = distanceX / ndivX; | |
dy = distanceY / ndivY; | |
cx = 1.0 / dx; | |
cy = 1.0 / dy; | |
#dt = 0.5 * 1.0 / (c0 * sqrt(cx * cx + cy * cy)); | |
dt = courantNumber * 1.0 / (c0 * sqrt(cx * cx + cy * cy)); | |
# X方向分割幅とY方向分割幅は同じとする | |
dl = dx; | |
cl = cx; | |
assert(abs(dx - dy) < 1.0e-12); | |
#--------------------------------------------- | |
# 初期化 | |
#--------------------------------------------- | |
Ez = zeros(ndivX + 1, ndivY + 1); | |
EzPrev = zeros(size(Ez)); | |
EzPrevPrev = zeros(size(Ez)); | |
# 比誘電率の逆数 | |
qzz = ones(ndivX + 1, ndivY + 1); | |
# Ez更新式の係数 | |
c1 = zeros(ndivX + 1, ndivY + 1); | |
c2 = zeros(ndivX + 1, ndivY + 1); | |
#c3 = zeros(ndivX + 1, ndivY + 1); | |
c4 = zeros(ndivX + 1, ndivY + 1); | |
# 吸収境界条件で使用する界 | |
prevEz1 = zeros((ndivY + 1), 1); # 境界1 + dx | |
prevEz2 = zeros((ndivY + 1), 1); # 境界2 - dx | |
# ポート上の1節点の界の時間変化 | |
timeAry = zeros(loopCnt, 1); | |
EzTimePort1 = zeros(loopCnt, 1); | |
EzTimePort2 = zeros(loopCnt, 1); | |
#--------------------------------------------- | |
# 励振源 | |
#--------------------------------------------- | |
# 周波数 | |
#freq = 5.0e+9; | |
#freq = 1.5 * c0 / (2.0 * distanceY); | |
#freq = 2.5 * c0 / (2.0 * distanceY); | |
freq = normalizedFreqCarrier * c0 / (2.0 * distanceY); | |
# 角周波数 | |
omega = 2.0 * pi * freq; | |
# 波長 | |
waveLength = c0 / freq; | |
# 波数 | |
k0 = 2.0 * pi / waveLength; | |
# 平面波インピーダンス | |
z0 = sqrt(mu0 / eps0); | |
printf("waveLength / dx: %f\n", waveLength / dx); | |
printf("waveLength / dy: %f\n", waveLength / dy); | |
#assert(false); | |
# 励振源の位置 | |
#srcPosX = 5; | |
srcPosX = 1 * ndivPerUnit; | |
# 付加励振源のY方向界分布 | |
srcProfileY = []; | |
# 付加励振源のY方向界分布 | |
srcProfileY = zeros((ndivY + 1), 1); | |
for i = 1: (ndivY + 1) | |
srcProfileY(i) = sin((pi/ndivY) * (i - 1)); # TE波(導波管TE10モード) | |
endfor | |
# 伝搬定数(X方向) | |
betaX = sqrt(k0 * k0 - (pi / distanceY) * (pi / distanceY)); # TE波(導波管TE10モード) | |
# 界インピーダンス | |
zf = (omega * mu0) / betaX; # TE波 | |
# ガウシアンパルス | |
#isGaussianPulse = false; | |
isGaussianPulse = true; | |
#typeGaussian = 0; # 素のガウシアンパルス | |
##typeGaussian = 2; # 正弦波変調 | |
if (typeGaussian == 0) | |
# 素のガウシアンパルス | |
#Tp = 30.0 * dt / sqrt(2.0); | |
#delay = 4.0 * sqrt(2.0) * Tp; | |
delay = 120.0 * dt; | |
Tp = delay / (4.0 * sqrt(2.0)); | |
elseif (typeGaussian == 1) | |
# 微分ガウシアンパルス | |
delay = 120.0 * dt; | |
Tp = delay / (4.0 * sqrt(2.0)); | |
elseif (typeGaussian == 2) | |
# 正弦波変調 | |
nCycle = 5; | |
delay = (1.0 / freq) * nCycle / 2; | |
Tp = delay / (2.0 * sqrt(2.0 * log(2.0))); | |
else | |
assert(false); | |
endif | |
nTp = Tp / dt; | |
ndelay = delay / dt; | |
# Murの吸収境界条件 | |
vpx = (omega / betaX); | |
cx_MurAbc = (vpx * dt - dx) / (vpx * dt + dx); | |
# ポート観測位置 | |
#port1RefPosXOfs = 2; # ポート1:励振源からのオフセット | |
#port2RefPosXOfs = 2; # ポート2:出力境界からのオフセット | |
port1RefPosXOfs = 1 * ndivPerUnit; # ポート1:励振源からのオフセット | |
port2RefPosXOfs = 1 * ndivPerUnit; # ポート2:出力境界からのオフセット | |
port1PosX = srcPosX + port1RefPosXOfs; | |
port1PosY = floor(ndivY / 2); | |
port2PosX = ndivX - port2RefPosXOfs; | |
port2PosY = floor(ndivY / 2); | |
assert(port1PosX < ndivX1); | |
assert(port2PosX > (ndivX - ndivX3)); | |
printf("waveguide TMz(TE) : isGaussianPulse %d typeGaussian %d normalizedFreqCarrier %f\n", isGaussianPulse, typeGaussian, normalizedFreqCarrier); | |
printf("nTp: %d ndelay: %d\n", nTp, ndelay); | |
#assert(false); | |
#--------------------------------------------- | |
# 媒質設定 | |
#--------------------------------------------- | |
# 境界条件 Ez = qzz Dz | |
# ポート1導波路 | |
# 上下の導体 | |
qzz(1:(ndivX1 + 1), 1) = 0.0; | |
qzz(1:(ndivX1 + 1), (ndivY + 1)) = 0.0; | |
# 不連続部 | |
# 上下の導体 | |
qzz((ndivX1 + 1):(ndivX - ndivX3 + 1), 1:(ndivYOfs + 1)) = 0.0; | |
qzz((ndivX1 + 1):(ndivX - ndivX3 + 1), (ndivY - ndivYOfs + 1):(ndivY + 1)) = 0.0; | |
# ポート2導波路 | |
# 上下の導体 | |
qzz((ndivX - ndivX3 + 1):(ndivX + 1), 1) = 0.0; | |
qzz((ndivX - ndivX3 + 1):(ndivX + 1), (ndivY + 1)) = 0.0; | |
#qzz(:, 1) = 0.0; | |
#qzz(:, (ndivY + 1)) = 0.0; | |
# 散乱オブジェクト | |
scatterEr = 2.62; | |
#scatterEr = 1.00; | |
scatterPosX = ndivX1 + 6 * ndivPerUnit + 1; | |
scatterPosY = ndivYOfs + 2 * ndivPerUnit + 1; | |
#scatterLenX = 7 * ndivPerUnit; | |
#scatterLenY = 6 * ndivPerUnit; | |
scatterLenX = 7 * ndivPerUnit - 1; | |
scatterLenY = 6 * ndivPerUnit; | |
# 左下 | |
qzz(scatterPosX, scatterPosY) = 1.0 / ((1.0 + 1.0 + 1.0 + scatterEr) / 4.0); | |
# 左上 | |
qzz(scatterPosX, (scatterPosY + scatterLenY)) = 1.0 / ((1.0 + 1.0 + 1.0 + scatterEr) / 4.0); | |
# 右下 | |
qzz((scatterPosX + scatterLenX), scatterPosY) = 1.0 / ((1.0 + 1.0 + 1.0 + scatterEr) / 4.0); | |
# 右上 | |
qzz((scatterPosX + scatterLenX), (scatterPosY + scatterLenY)) = 1.0 / ((1.0 + 1.0 + 1.0 + scatterEr) / 4.0); | |
# 左境界 | |
qzz(scatterPosX, (scatterPosY + 1):(scatterPosY + scatterLenY - 1)) = 1.0 / ((1.0 + scatterEr) / 2.0); | |
# 右境界 | |
qzz((scatterPosX + scatterLenX), (scatterPosY + 1):(scatterPosY + scatterLenY - 1)) = 1.0 / ((1.0 + scatterEr) / 2.0); | |
# 下境界 | |
qzz((scatterPosX + 1):(scatterPosX + scatterLenX - 1), scatterPosY) = 1.0 / ((1.0 + scatterEr) / 2.0); | |
# 上境界 | |
qzz((scatterPosX + 1):(scatterPosX + scatterLenX - 1), (scatterPosY + scatterLenY)) = 1.0 / ((1.0 + scatterEr) / 2.0); | |
# 内部 | |
qzz((scatterPosX + 1):(scatterPosX + scatterLenX - 1), (scatterPosY + 1):(scatterPosY + scatterLenY - 1)) = 1.0 / scatterEr; | |
#disp(qzz((scatterPosX):(scatterPosX + scatterLenX), (scatterPosY):(scatterPosY + scatterLenY))); | |
#assert(false); | |
# 係数 | |
c4(:, :) = dt * dt * qzz(:, :) * cl * cl / (mu0 * eps0); | |
c1(:, :) = 2.0 - 4.0 * c4(:, :); | |
c2(:, :) = -1.0; | |
#c3(:, :) = -0.5 * dt * qzz(:, :) / eps0; | |
#--------------------------------------------- | |
# Yeeアルゴリズム | |
#--------------------------------------------- | |
for n = 1:loopCnt | |
# 電界 | |
EzPrevPrev(:, :) = EzPrev(:, :); | |
EzPrev(:, :) = Ez(:, :); | |
Ez(2:ndivX, 2:ndivY) = ... | |
c1(2:ndivX, 2:ndivY) .* EzPrev(2:ndivX, 2:ndivY) ... | |
+ c2(2:ndivX, 2:ndivY) .* EzPrevPrev(2:ndivX, 2:ndivY) ... | |
+ c4(2:ndivX, 2:ndivY) .* ( ... | |
EzPrev(3:(ndivX + 1), 2:ndivY) + EzPrev(1:(ndivX - 1), 2:ndivY) ... | |
+ EzPrev(2:ndivX, 3:(ndivY + 1)) + EzPrev(2:ndivX, 1:(ndivY - 1)) ... | |
); | |
# 励振 | |
if (isGaussianPulse) | |
#for i = 2 : ndivY | |
for i = 1 : (ndivY + 1) | |
srcJz1 = 0.0; | |
srcJz2 = 0.0; | |
if ((n) <= floor(2.0 * ndelay + 1.0)) | |
if (typeGaussian == 0) | |
srcJz1 = srcProfileY(i) ... | |
* exp(-1.0 * ((n + 1) - ndelay) * ((n + 1) - ndelay)/(2.0 * nTp * nTp)); | |
srcJz2 = srcProfileY(i) ... | |
* exp(-1.0 * ((n - 1) - ndelay) * ((n - 1) - ndelay)/(2.0 * nTp * nTp)); | |
elseif (typeGaussian == 1) | |
srcJz1 = srcProfileY(i) ... | |
* (-sqrt(2.0) * ((n + 1) - ndelay) / nTp) * exp(-1.0 * ((n + 1) - ndelay) * ((n + 1) - ndelay)/(2.0 * nTp * nTp)); | |
srcJz2 = srcProfileY(i) ... | |
* (-sqrt(2.0) * ((n - 1) - ndelay) / nTp) * exp(-1.0 * ((n - 1) - ndelay) * ((n - 1) - ndelay)/(2.0 * nTp * nTp)); | |
elseif (typeGaussian == 2) | |
srcJz1 = srcProfileY(i) * sin(omega * ((n + 1) - ndelay) * dt) ... | |
* exp(-1.0 * ((n + 1) - ndelay) * ((n + 1) - ndelay)/(2.0 * nTp * nTp)); | |
srcJz2 = srcProfileY(i) * sin(omega * ((n - 1) - ndelay) * dt) ... | |
* exp(-1.0 * ((n - 1) - ndelay) * ((n - 1) - ndelay)/(2.0 * nTp * nTp)); | |
else | |
assert(false); | |
endif | |
else | |
srcJz1 = 0.0; | |
srcJz2 = 0.0; | |
endif | |
srcJz1 *= -1.0; | |
srcJz2 *= -1.0; | |
Ez(srcPosX, i) = Ez(srcPosX, i) - (0.5 * dt / eps0) * (srcJz1 - srcJz2); | |
endfor | |
else | |
#for i = 2 : ndivY | |
for i = 1 : (ndivY + 1) | |
srcJz1 = srcProfileY(i) * sin(omega * (n + 1) * dt); | |
srcJz2 = srcProfileY(i) * sin(omega * (n - 1) * dt); | |
srcJz1 *= -1.0; | |
srcJz2 *= -1.0; | |
Ez(srcPosX, i) = Ez(srcPosX, i) - (0.5 * dt / eps0) * (srcJz1 - srcJz2); | |
endfor | |
endif | |
# Murの吸収境界条件 | |
# 左の境界 | |
workEzb1 = getRowVec(Ez, 1); # 境界1 | |
workEz1 = getRowVec(Ez, 2); # 境界1 + dx | |
workEzb1(:) = prevEz1(:) + cx_MurAbc * (workEz1(:) - workEzb1(:)); | |
# 上下の導体 | |
workEzb1(1) = 0; | |
workEzb1(ndivY + 1) = 0; | |
prevEz1(:) = workEz1(:); | |
Ez = setRowVec(Ez, 1, workEzb1); # 境界1 | |
# 右の境界 | |
workEzb2 = getRowVec(Ez, ndivX + 1); # 境界2 | |
workEz2 = getRowVec(Ez, ndivX); # 境界2 - dx | |
workEzb2(:) = prevEz2(:) + cx_MurAbc * (workEz2(:) - workEzb2(:)); | |
# 上下の導体 | |
workEzb2(1) = 0; | |
workEzb2(ndivY + 1) = 0; | |
prevEz2(:) = workEz2(:); | |
Ez = setRowVec(Ez, ndivX + 1, workEzb2); # 境界2 | |
# ポート観測点の界を格納 | |
timeAry(n) = n * dt; | |
EzTimePort1(n) = Ez(port1PosX, port1PosY); | |
EzTimePort2(n) = Ez(port2PosX, port2PosY); | |
endfor; | |
endfunction | |
#------------------------------------------------------- | |
# rowベクトルを取得する | |
#------------------------------------------------------- | |
function vec = getRowVec(ary, rowIndex) | |
arySize = size(ary); | |
rowSize = arySize(1); | |
colSize = arySize(2); | |
vec = zeros(colSize, 1); | |
for i = 1 : colSize | |
vec(i) = ary(rowIndex, i); | |
endfor | |
endfunction | |
#------------------------------------------------------- | |
# rowベクトルを配列に格納する | |
#------------------------------------------------------- | |
function retAry = setRowVec(ary, rowIndex, vec) | |
arySize = size(ary); | |
rowSize = arySize(1); | |
colSize = arySize(2); | |
retAry = ary; | |
for i = 1 : colSize | |
retAry(rowIndex, i) = vec(i); | |
endfor | |
endfunction | |
#------------------------------------------------------- | |
# colベクトルを取得する | |
#------------------------------------------------------- | |
function vec = getColVec(ary, colIndex) | |
arySize = size(ary); | |
rowSize = arySize(1); | |
colSize = arySize(2); | |
vec = zeros(rowSize, 1); | |
for i = 1 : rowSize | |
vec(i) = ary(i, colIndex); | |
endfor | |
endfunction | |
#------------------------------------------------------- | |
# colベクトルを配列に格納する | |
#------------------------------------------------------- | |
function retAry = setColVec(ary, colIndex, vec) | |
arySize = size(ary); | |
rowSize = arySize(1); | |
colSize = arySize(2); | |
retAry = ary; | |
for i = 1 : rowSize | |
retAry(i, colIndex) = vec(i); | |
endfor | |
endfunction |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment