WE-FDTD TMz(TE mode) Mur ABC. This script calculates S11 and S21 of a dielectric resonator in a hollow waveguide.

main_we_mur_wg2d_problem01_2.m

###########################################################
#  TE波(TMz) 導波管誘電体共振器 (WE-FDTD Murの吸収境界条件)
#
#  coded by ryujimiya April 2013
###########################################################
# ウィンドウをすべて閉じる
close all;
# メモリを解放する
clear all;

#------------------------------------------------------------------
# settings
#------------------------------------------------------------------
# 計算回数
#loopCnt = 8000;
loopCnt = 10000;
# クーラン数
#courantNumber = 0.5;
courantNumber = 0.5;
# ガウシアンパルスの種類
typeGaussian = 0;  # 素のガウシアンパルス
#typeGaussian = 1;  # 微分ガウシアンパルス
#typeGaussian = 2; # 正弦波変調
# 正弦波変調の時の搬送波規格化周波数
normalizedFreqCarrier = 1.5;
#normalizedFreqCarrier = 2.5;

#------------------------------------------------------------------
# FDTD calculation
#------------------------------------------------------------------
# TMz
# 直線導波管 入射波の時間変化を取得
printf("calculating straight waveguide ...\n");
source "we_mur_wg2d_straight.m";
[null, null, null, null, ...
    null, ...
    null, timeAry, ...
    EzTimePort1Inc, null] = we_fdtd_2d_TMz_straight(loopCnt, courantNumber, typeGaussian, normalizedFreqCarrier);
printf("... done.\n");

# 誘電体装荷共振器
printf("calculating discontinuities ...\n");
source "we_mur_wg2d_dielectric_filter.m";
[ndivX, ndivY, distanceX, distanceY, ...
    qzz, ...
    Ez, timeAry, ...
    EzTimePort1, EzTimePort2] = we_fdtd_2d_TMz_dielectric_filter(loopCnt, courantNumber, typeGaussian, normalizedFreqCarrier);
printf("... done.\n");

#------------------------------------------------------------------
# FFT function
#------------------------------------------------------------------
function [f, fd, pd] = fft_field(tAryAll, waveAryAll, n1, n2, showFlg)
    tAry = tAryAll(n1:n2);
    waveAry = waveAryAll(n1:n2);
    dt = tAry(2) - tAry(1);

    # フーリエ変換
    dataCnt = length(waveAry);
    fd = fft(waveAry) / dataCnt;
    # フーリエ振幅
    spec = abs(fd);
    # パワースペクトルの計算
    #  .*(要素同士のかけ算)とconj(共役複素数) 
    pd = fd .* conj(fd);

    # 時間長さ
    tl = dt * dataCnt;
    # 周波数刻み
    df = 1.0 / tl;
    # 周波数刻みの値をを行列で作成
    f = (0:(dataCnt - 1)) * df;

    if (showFlg)
        # 波形表示
        subplot(2, 1, 1); # 上側
        plot(tAry, waveAry);
        # スペクトルを表示
        subplot(2, 1, 2); # 下側
        # パワースペクトルとフーリエ振幅)の表示
        plot(f, spec, "-;fd;", f, pd, "-;pd;");
    endif
endfunction

#------------------------------------------------------------------
# plpt s11 & s21 function
#------------------------------------------------------------------
function plot_s_parameters( ...
    waveguideWidth, ...
    tAry, waveAryPort1Inc, waveAryPort1, waveAryPort2)
    # 定数
    c0 = 2.99792458e8;

    # check
    #plot(tAry / (tAry(2)- tAry(1)), waveAryPort1, "-;port1;", ...
    #    tAry / (tAry(2)- tAry(1)), waveAryPort2, "-;port2;");

    samplingCnt = length(tAry);
    samplingCntPort1Inc = length(waveAryPort1Inc);
    assert(samplingCntPort1Inc <= samplingCnt);

    # 入射波の波形データ数を反射波、透過波に合わせる(0埋め)
    workWaveAryPort1Inc = zeros(samplingCnt, 1);
    workWaveAryPort1Inc(1:samplingCntPort1Inc) = waveAryPort1Inc(1:samplingCntPort1Inc);
    waveAryPort1Inc = workWaveAryPort1Inc;
    workWaveAryPort1Inc = [];
    # ポート1の波形データから入射波を引く
    waveAryPort1Reflect = zeros(samplingCnt, 1);
    waveAryPort1Reflect(:) = waveAryPort1(:) - waveAryPort1Inc(:);
    
    # ポート1は入射波と反射波を分離し、0埋めする
    # ポート2はすべて使用
    workWaveAry = waveAryPort1Inc;
    [freq1, fd1, pd1] = fft_field(tAry, workWaveAry, 1, length(workWaveAry), false);
    workWaveAry = waveAryPort1Reflect;
    [freq1r, fd1r, pd1r] = fft_field(tAry, workWaveAry, 1, length(workWaveAry), false);
    workWaveAry = waveAryPort2;
    [freq2, fd2, pd2] = fft_field(tAry, workWaveAry, 1, length(workWaveAry), false);

    # check
    assert(length(freq1) == length(freq1r));
    assert(length(freq1) == length(freq2));
    assert(length(freq1) == length(fd1));
    assert(length(freq1) == length(fd1r));
    assert(length(freq1) == length(fd2));
    assert(length(freq1) == length(pd1));
    assert(length(freq1) == length(pd1r));
    assert(length(freq1) == length(pd2));

    # 切り捨てる値を決める
    thVal = 0.0;
    maxPd = 0.0;
    for i = 1:length(pd1)
        if (maxPd < pd1(i))
            maxPd = pd1(i);
        endif
    endfor
    thVal = 0.01 * maxPd;
    assert(abs(maxPd) >= 1.0e-12);
    # 反射係数
    s11Freq = zeros(length(pd1r), 1);
    # 透過係数
    #s21Freq = abs(fd2) ./ abs(fd1);
    #s21Freq = sqrt(pd2) ./ sqrt(pd1);
    s21Freq = zeros(length(pd2), 1);
    for i = 1:length(pd1)
        if (pd1(i) >= thVal)
            s11Freq(i) = sqrt(pd1r(i) / pd1(i));
            s21Freq(i) = sqrt(pd2(i) / pd1(i));
        endif
    endfor
    # 周波数特性の表示
    #freqMax = freq1(length(freq1));
    #plot(freq1, s21Freq);
    #axis([0, freqMax / 2]); # 周波数領域は半分だけ表示
    nFreq = length(freq1) / 2;
    normalizedFreq = zeros(nFreq, 1);
    normalizedFreq(1:nFreq) = freq1(1:nFreq) * 2.0 * waveguideWidth / c0;
    s11FreqHalf(1:nFreq) = s11Freq(1:nFreq);
    s21FreqHalf(1:nFreq) = s21Freq(1:nFreq);

    subplot(2, 1, 1); # 上側
    plot(
        tAry / (tAry(2)- tAry(1)), waveAryPort1Inc, "-;port1(Inc);", ...
        tAry / (tAry(2)- tAry(1)), waveAryPort1, "-;port1;", ...
        tAry / (tAry(2)- tAry(1)), waveAryPort2, "-;port2;");
    subplot(2, 1, 2); # 下側
    plot(normalizedFreq, s11FreqHalf, "-;s11;", normalizedFreq, s21FreqHalf, "-;s21;");
    axis([1.0, (3.0 + 1.0e-12), 0.0, (1.4 + 1.e-12)]);
endfunction

#------------------------------------------------------------------
# plot s11 & s21
#------------------------------------------------------------------

workWaveAryPort1Inc = EzTimePort1Inc;
workWaveAryPort1 = EzTimePort1;
workWaveAryPort2 = EzTimePort2;
# 反射係数,透過係数の計算
plot_s_parameters(distanceY, timeAry, workWaveAryPort1Inc, workWaveAryPort1, workWaveAryPort2);

#------------------------------------------------------------------
# plot function
#------------------------------------------------------------------
function plot_field(fVals, ndiv1, ndiv2)
    x = [];
    y = [];
    z = [];
    for i = 1 : ndiv1
        x(i) = i;
    endfor
    for i = 1 : ndiv2
        y(i) = i;
    endfor

    for i = 1 : ndiv1
        for j = 1 : ndiv2
           z(j, i) = fVals(i, j); 
        endfor
    endfor

    figure;
    mesh(x, y, z);
endfunction

#------------------------------------------------------------------
# plot
#------------------------------------------------------------------
# qzz
plot_field((1.0 ./ qzz), ndivX + 1, ndivY + 1);
#TMz
plot_field(Ez, ndivX + 1, ndivY + 1);