Skip to content

Instantly share code, notes, and snippets.

@ryujimiya

ryujimiya/we_mur_wg2d_straight.m

Created April 25, 2013 14:27
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ryujimiya/5460117 to your computer and use it in GitHub Desktop.
Save ryujimiya/5460117 to your computer and use it in GitHub Desktop.
Download ZIP
WE-FDTD TMz(TE mode) Mur ABC. This script calculates field Ez of straight hollow waveguides.
Raw
we_mur_wg2d_straight.m
###########################################################
# TE波(TMz) 導波管直線 (WE-FDTD Murの吸収境界条件)
#
# coded by ryujimiya April 2013
###########################################################
function [ndivX, ndivY, distanceX, distanceY, ...
qzz, ...
Ez, timeAry, ...
EzTimePort1, EzTimePort2] = we_fdtd_2d_TMz_straight(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;
#ndivX = 100 * ndivPerUnit;
ndivX = 20 * ndivPerUnit;
ndivY = 20 * ndivPerUnit;
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);
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) = 0.0;
qzz(:, (ndivY + 1)) = 0.0;
# 係数
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