function rect_wg_nodal_fem
    input();
    mesh();
%    outmesh();
    pltmesh();
    sysmatrix();
    eigsolv();
    outeig();
%    outpot(2);
    pltpot(2);
%    outfield(2);
    pltfield(2);
end

function input()
    global nx ny sizex sizey;

    disp('**** INPUT ****');

    nx = 20;            % x方向分割数
    ny = 10;            % y方向分割数
    sizex = 58.1e-3;    % 導波管幅(a) [m]
    sizey = 29.1e-3;    % 導波管高(b) [m]
end

function mesh()
    global nelem nnode nx ny sizex sizey elem gn2x gn2y;

    disp('**** MESH ****');

    nelem = nx * ny * 2;             % 三角要素の数
    nnode = (nx + 1) * (ny + 1);     % 系全体の節点数

    % ローカル→グローバルへの変換
    %
    % k4                              k3
    % **********************************
    % *                            *****
    % *                        ****    *
    % *                    ****        *
    % *                ****            *
    % *            ****                *
    % *        ****                    *
    % *    ****                        *
    % *****                            *
    % **********************************
    % k1                              k2
    %
    % グローバル節点番号の計算
    k = 0;
    for i = 1:nx
        for j = 1:ny
            k1 = (ny + 1) * (i - 1) + j;      % 左下
            k2 = (ny + 1) * i + j;            % 右下
            k3 = (ny + 1) * i + (j + 1);      % 右上
            k4 = (ny + 1) * (i - 1) + (j + 1);% 左上

            k = k + 1;
            elem(k).node = [k1, k2, k3];
            k = k + 1;
            elem(k).node = [k3, k4, k1];
        end
    end

    % 節点の x,y 座標を格納する
    k = 0;
    for i = 1:(nx + 1)
        for j = 1:(ny + 1)
            k = k + 1;
            gn2x(k) = sizex * (i - 1) / nx;
            gn2y(k) = sizey * (j - 1) / ny;
        end
    end
end

function outmesh()
    global nelem elem gn2x gn2y;

    disp('**** OUTMESH ****');

    fileID = fopen('mesh.dat', 'w');

    for i = 1:nelem
        n = zeros(1, 3);
        x = zeros(1, 3);
        y = zeros(1, 3);

        % グローバル節点番号の取得
        for j = 1:3
            n(j) = elem(i).node(j);
        end

        % x,y 座標の取得
        for j = 1:3
            x(j) = gn2x(n(j));
            y(j) = gn2y(n(j));
        end

        % ファイルへの書き込み
        for j = 1:3
            fprintf(fileID, '%f %f\n', x(j), y(j));
        end

        fprintf(fileID, '%f %f\n\n', x(1), y(1)); % ポリゴンを閉じるための追加の書き込み
    end

    fclose(fileID);
end

function pltmesh()
    global nelem elem gn2x gn2y;
    global fig_mesh;

    disp('**** PLOT MESH ****');

    fig_mesh=figure;
    for i = 1:nelem
        n = zeros(1, 3);
        x = zeros(1, 3);
        y = zeros(1, 3);

        % グローバル節点番号の取得
        for j = 1:3
            n(j) = elem(i).node(j);
        end

        % x,y 座標の取得
        for j = 1:3
            x(j) = gn2x(n(j));
            y(j) = gn2y(n(j));
        end

        % 描画
        patch(x,y,'white');
    end
end

function [t, u] = elemmatrix(ie)
    global elem gn2x gn2y;

    t = zeros(3, 3);
    u = zeros(3, 3);

    n1 = elem(ie).node(1);
    n2 = elem(ie).node(2);
    n3 = elem(ie).node(3);

    x = [gn2x(n1), gn2x(n2), gn2x(n3)];
    y = [gn2y(n1), gn2y(n2), gn2y(n3)];

    b = [y(2)-y(3), y(3)-y(1), y(1)-y(2)];
    c = [x(3)-x(2), x(1)-x(3), x(2)-x(1)];
    area = (x(1)*b(1) + x(2)*b(2) + x(3)*b(3)) / 2.0;

    % エラー処理
    if area <= 0.0
        disp('Sub(elemmatrix):');
        disp('  area <= 0');
        disp(['  element number=', num2str(ie)]);
        error('Error: Area is less than or equal to zero.');
    end

    for i = 1:3
        for j = 1:3
            % T=K∇ 行列の計算
            t(i, j) = (b(i)*b(j) + c(i)*c(j)) / (area * 4.0);

            % U=K 行列の計算
            if i == j
                u(i, j) = area / 6.0;
            else
                u(i, j) = area / 12.0;
            end
        end
    end
end

function sysmatrix()
    global nelem nnode syst sysu elem;

    disp('**** SYSTEM MATRIX ****');

    syst = zeros(nnode, nnode);
    sysu = zeros(nnode, nnode);

    for i = 1:nnode
        for j = 1:nnode
            syst(i, j) = 0.0;
            sysu(i, j) = 0.0;
        end
    end

    for ie = 1:nelem
        [t, u] = elemmatrix(ie);
        n = elem(ie).node;
        
        for i = 1:3
            for j = 1:3
                syst(n(i), n(j)) = syst(n(i), n(j)) + t(i, j);
                sysu(n(i), n(j)) = sysu(n(i), n(j)) + u(i, j);
            end
        end
    end
end

function eigsolv()
    global syst sysu nnode;
    global eig_values eig_vectors;

    disp('**** SOLVE ****');

    % syst と sysu を使って一般化固有値問題を解く処理を実装する必要があります。
    % 以下は一般化固有値問題を解く簡単な例です。

    % 例えば、syst と sysu を使って一般化固有値問題を解く場合:
    num_eigs=nnode;
    [V, D] = eigs(syst, sysu, num_eigs, 'sm'); % num_eigs: 取得する固有値の数, 'sm': 最小の固有値を取得

    % 固有値と固有ベクトルの取得例:
    eig_values = diag(D); % 固有値を取得
    eig_vectors = V; % 固有ベクトルを取得

    % 得られた固有値や固有ベクトルを必要に応じて処理してください。
end

function outeig()
    global nnode;
    global eig_values eig_vectors;

    velocity_of_light = 2.998e8;

    disp('**** OUTEIG ****');

    % 固有値の出力
    fid_eigval = fopen('eigval.dat', 'w');
    for i = 1:nnode
        fprintf(fid_eigval, '%d %f\n', i, eig_values(i));
    end
    fclose(fid_eigval);

    % カットオフ周波数 (GHz) の出力
    fid_fcu = fopen('fcu.dat', 'w');
    for i = 1:nnode
        cutoff_freq = (velocity_of_light / (2.0 * pi)) * sqrt(abs(eig_values(i))) / 1.0e9;
        fprintf(fid_fcu, '%3d %6.3f\n', i, cutoff_freq);
    end
    fclose(fid_fcu);
end

function outpot(mode)
    global nelem elem gn2x gn2y;
    global eig_values eig_vectors;

    disp('**** OUTPOT ****');

    fid = fopen('pot.dat', 'w');
    for i = 1:nelem
        n = zeros(3, 1);
        x = zeros(3, 1);
        y = zeros(3, 1);
        
        % グローバル節点番号の取得と座標の取得
        for j = 1:3
            n(j) = elem(i).node(j);
            x(j) = gn2x(n(j));
            y(j) = gn2y(n(j));
        end
        
        % x, y, eig_vectors(mode, n)の書き込み
        for j = 1:3
            fprintf(fid, '%f %f %f\n', x(j), y(j), eig_vectors(mode, n(j)));
        end
        
        % ポリゴンを閉じるための書き込み
        fprintf(fid, '%f %f %f\n', x(1), y(1), eig_vectors(mode, n(1)));
        fprintf(fid, '\n');  % 改行
        fprintf(fid, '\n');  % 改行
    end
    fclose(fid);
end

function pltpot(mode)
    global nelem elem gn2x gn2y;
    global eig_values eig_vectors;

    disp('**** PLTPOT ****');

    figure;
    hold on;
    for i = 1:nelem
        n = zeros(3, 1);
        x = zeros(3, 1);
        y = zeros(3, 1);
        z = zeros(3, 1);
        
        % グローバル節点番号の取得と座標の取得
        for j = 1:3
            n(j) = elem(i).node(j);
            x(j) = gn2x(n(j));
            y(j) = gn2y(n(j));
            z(j) = eig_vectors(n(j),mode);
        end

        % x, y, zの描画
        fill3(x,y,z,'white');
    end
end

function outfield(mode)
    global nelem elem gn2x gn2y;
    global nx ny sizex sizey;
    global eig_values eig_vectors;

    disp('**** OUTFIELD ****');

    fid = fopen('field.dat', 'w');
    absmax = 0.0;
    for i = 1:nelem
        n1 = elem(i).node(1);
        n2 = elem(i).node(2);
        n3 = elem(i).node(3);

        x1 = gn2x(n1);
        y1 = gn2y(n1);
        z1 = eig_vectors(mode, n1);

        x2 = gn2x(n2);
        y2 = gn2y(n2);
        z2 = eig_vectors(mode, n2);

        x3 = gn2x(n3);
        y3 = gn2y(n3);
        z3 = eig_vectors(mode, n3);

        xc = (x1 + x2 + x3) / 3.0;
        yc = (y1 + y2 + y3) / 3.0;

        nvecx = -y2 * z1 + y3 * z1 + y1 * z2 - y3 * z2 - y1 * z3 + y2 * z3;
        nvecy = x2 * z1 - x3 * z1 - x1 * z2 + x3 * z2 + x1 * z3 - x2 * z3;
        nvecz = -x2 * y1 + x3 * y1 + x1 * y2 - x3 * y2 - x1 * y3 + x2 * y3;

        x = nvecx;
        y = nvecy;

        if sqrt(x^2 + y^2) > absmax
            absmax = sqrt(x^2 + y^2);
        end

        fprintf(fid, '%23.15e %23.15e %23.15e %23.15e\n', xc, yc, xc + y / absmax * (sizey / ny), yc + x / absmax * (sizey / ny));
        fprintf(fid, '\n');  % 改行
        fprintf(fid, '\n');  % 改行
    end
    fclose(fid);
end

function pltfield(mode)
    global nelem elem gn2x gn2y;
    global nx ny sizex sizey;
    global eig_values eig_vectors;
    global fig_mesh;

    disp('**** PLTFIELD ****');
    figure(fig_mesh);
    hold on;

    absmax = 0.0;
    for i = 1:nelem
        n1 = elem(i).node(1);
        n2 = elem(i).node(2);
        n3 = elem(i).node(3);

        x1 = gn2x(n1);
        y1 = gn2y(n1);
        z1 = eig_vectors(n1, mode);

        x2 = gn2x(n2);
        y2 = gn2y(n2);
        z2 = eig_vectors(n2, mode);

        x3 = gn2x(n3);
        y3 = gn2y(n3);
        z3 = eig_vectors(n3, mode);

        xc = (x1 + x2 + x3) / 3.0;
        yc = (y1 + y2 + y3) / 3.0;

        nvecx = -y2 * z1 + y3 * z1 + y1 * z2 - y3 * z2 - y1 * z3 + y2 * z3;
        nvecy = x2 * z1 - x3 * z1 - x1 * z2 + x3 * z2 + x1 * z3 - x2 * z3;
        nvecz = -x2 * y1 + x3 * y1 + x1 * y2 - x3 * y2 - x1 * y3 + x2 * y3;

        x = nvecx;
        y = nvecy;

        if sqrt(x^2 + y^2) > absmax
            absmax = sqrt(x^2 + y^2);
        end

        quiver(xc, yc, y * (sizey / ny), x * (sizex / nx));
    end

    for i = 1:nelem
        n1 = elem(i).node(1);
        n2 = elem(i).node(2);
        n3 = elem(i).node(3);

        x1 = gn2x(n1);
        y1 = gn2y(n1);
        z1 = eig_vectors(n1, mode);

        x2 = gn2x(n2);
        y2 = gn2y(n2);
        z2 = eig_vectors(n2, mode);

        x3 = gn2x(n3);
        y3 = gn2y(n3);
        z3 = eig_vectors(n3, mode);

        xc = (x1 + x2 + x3) / 3.0;
        yc = (y1 + y2 + y3) / 3.0;

        nvecx = -y2 * z1 + y3 * z1 + y1 * z2 - y3 * z2 - y1 * z3 + y2 * z3;
        nvecy = x2 * z1 - x3 * z1 - x1 * z2 + x3 * z2 + x1 * z3 - x2 * z3;
        nvecz = -x2 * y1 + x3 * y1 + x1 * y2 - x3 * y2 - x1 * y3 + x2 * y3;

        x = nvecx;
        y = nvecy;

        quiver(xc, yc, y / absmax * (sizey / ny), x / absmax * (sizex / nx),'red');
    end
end