%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1-D FEM 解析プログラム (T. Hirano)
% 動作確認環境: MATLAB R2019a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fem_1d
    % メインのfem_1d関数の下部にあるサブルーチンで使えるように定義
    % 綺麗な汎用性の高いプログラムを組み場合は隠蔽して
    % オブジェクト指向にした方がいいのだが。
    global n_zs nelem glob_edge_no pos_z;

    ci=complex(0,1);        % 虚数単位i
    c=2.99792458e8;         % 光速

    f=2.45e9;               % 周波数
    w=2*pi*f;               % 角周波数
    lambda0=c/f;            % 波長
    k0=2*pi/lambda0;        % 波数
    eps0=8.854e-12;         % 真空の誘電率
    mu0=4*pi*1e-7;          % 真空の透磁率
    eta0=sqrt(mu0/eps0);    % 真空の波動インピーダンス

    len_z=2*lambda0;        % 波長単位での解析領域の長さ
    dz=lambda0/40;          % メッシュのサイズ
    nz=ceil(len_z/dz);      % メッシュ（セル）数
    nelem=nz;               % FEMのエレメント（要素）数

    % 電界、磁界メモリの初期化: (nz+1)行、1列の初期値0の配列を確保
    ex=zeros(nz+1,1);
    hy=zeros(nz+1,1);

    % 媒質の設定: 各エレメントの電気定数を設定する
    % 本来は、各セルごとに設定するとメモリが無駄になるので媒質番号を設定し、
    % 各セルには媒質番号を参照させた方がよい。
    er=ones(nz+1,1);        % εr
    mur=ones(nz+1,1);       % μr
    sig=zeros(nz+1,1);      % σ
    
    % 構造のモデル化
    n_z1=round(0.37*lambda0/dz);    % ガラスの左端のエレメント番号
    n_z2=round(0.75*lambda0/dz);    % ガラスの右端のエレメント番号
    n_zs=round(1.50*lambda0/dz);    % 波源位置のエレメント番号
    for i=n_z1:n_z2
        er(i)=4;   % ガラス SiO2 を想定した比誘電率（上書き更新）
    end

    % テーブルを作成する
    % glob_edge_no(要素番号,ローカルエッジ番号)->グローバルエッジ番号
    glob_edge_no=zeros(nz+1);
    for i=1:nelem
        glob_edge_no(i,1)=i;
        glob_edge_no(i,2)=i+1;
    end

    % グローバルエッジ番号からそのz座標を返す関数の定義
    pos_z=@(no_g_edge) dz.*(no_g_edge-1);

    A=zeros(nz+1);
    B=zeros(nz+1,1);
    % System Matrix（系行列）を作成する
    % エレメントごとにループし、それを全体の系行列に埋め込む
    % この手法は工夫がないが、実用的には疎行列の性質を利用して反復解法を用いる
    for i=1:nelem
        elem_matrix(1,1)=(1./mur(i)).*emn(i,1,1)-k0^2.*er(i).*fmn(i,1,1);
        elem_matrix(1,2)=(1./mur(i)).*emn(i,1,2)-k0^2.*er(i).*fmn(i,1,2);
        elem_matrix(2,1)=(1./mur(i)).*emn(i,2,1)-k0^2.*er(i).*fmn(i,2,1);
        elem_matrix(2,2)=(1./mur(i)).*emn(i,2,2)-k0^2.*er(i).*fmn(i,2,2);

        elem_vector(1)=-ci.*k0.*eta0.*ecur(i)./2;
        elem_vector(2)=-ci.*k0.*eta0.*ecur(i)./2;
        % System Matrixに組み込む（...はMATLABで次の行に続くという意味）
        A(glob_edge_no(i,1),glob_edge_no(i,1))= ...
            A(glob_edge_no(i,1),glob_edge_no(i,1))+elem_matrix(1,1);
        A(glob_edge_no(i,1),glob_edge_no(i,2))= ...
            A(glob_edge_no(i,1),glob_edge_no(i,2))+elem_matrix(1,2);
        A(glob_edge_no(i,2),glob_edge_no(i,1))= ...
            A(glob_edge_no(i,2),glob_edge_no(i,1))+elem_matrix(2,1);
        A(glob_edge_no(i,2),glob_edge_no(i,2))= ...
            A(glob_edge_no(i,2),glob_edge_no(i,2))+elem_matrix(2,2);
        B(glob_edge_no(i,1))=B(glob_edge_no(i,1))+elem_vector(1);
        B(glob_edge_no(i,2))=B(glob_edge_no(i,2))+elem_vector(2);
    end
    % 吸収境界条件
    % [Boundary]の項が方程式に追加されることになるので、左側だと
    % (1/μr(1))*A1*(β1-α1*e_{1x}^{1}(z1))となり、
    % e_{1x}^{1}(z1)は左で値が1, α1=-j k, β1=0と置くので次のようになる
    A(glob_edge_no(1,1),glob_edge_no(1,1))= ...
        A(glob_edge_no(1,1),glob_edge_no(1,1))...
        +(ci*k0*sqrt(mur(1)*er(1))/mur(1));     % 左側の壁
    A(glob_edge_no(nelem,2),glob_edge_no(nelem,2))= ...
        A(glob_edge_no(nelem,2),glob_edge_no(nelem,2)) ...
        +(ci*k0*sqrt(mur(nelem)*er(nelem))/mur(nelem));     % 右側の壁

    % PEC境界条件 (行列・ベクトルサイズの変更)
    % 左
    % A(:,1)=[];  % 第1列を取り除く
    % A(1,:)=[];  % 第1行を取り除く
    % B(1)=[];    % 第1要素を取り除く
    % 右
    % A(:,nz)=[]; % 第nz列を取り除く
    % A(nz,:)=[]; % 第nz行を取り除く
    % B(nz)=[];   % 第nz要素を取り除く

    % 行列方程式を解く
    efield_vector=A\B;

    % 磁界を計算する
    hfield_vector=zeros(nz,1);
    for i=2:nz-1
        hfield_vector(i)=-(1/(ci*w*mur(i)*mu0)) ...
                         *(efield_vector(i+1)-efield_vector(i))/dz;
    end

    % グラフ出力
    pos_z_vector=dz*[0:nelem];
    plot(pos_z_vector,abs(efield_vector));
    rectangle('Position',[n_z1*dz,0,n_z2*dz,max(abs(efield_vector))], ...
              'EdgeColor','cyan');

    figure;
    plot(pos_z_vector(2:nz-1),abs(hfield_vector(2:nz-1)));
    rectangle('Position',[n_z1*dz,0,n_z2*dz,max(abs(hfield_vector))], ...
              'EdgeColor','cyan');

    % CSVファイルに出力
    csvwrite('pos_z.dat',pos_z_vector);
    csvwrite('e_field.dat',efield_vector);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% サブルーチン
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% エレメント番号を入力して、励振されていたら1を返す。
function y=ecur(elem)
    global n_zs nelem glob_edge_no pos_z;

    if(elem==n_zs)
        y=1;
    else
        y=0;
    end;
end

% 2つの基底関数の積の積分(E_{ij})
function y=emn(elem,m,n)
    global n_zs nelem glob_edge_no pos_z;

    if(m==n)
        y=1./(pos_z(glob_edge_no(elem,2))-pos_z(glob_edge_no(elem,1)));
    else
        y=1./(pos_z(glob_edge_no(elem,1))-pos_z(glob_edge_no(elem,2)));
    end;
end

% 2つの基底関数の微分の積の積分(F_{ij})
function y=fmn(elem,m,n)
    global n_zs nelem glob_edge_no pos_z;

    if(m==n)
        y=(pos_z(glob_edge_no(elem,2))-pos_z(glob_edge_no(elem,1)))./3;
    else
        y=(pos_z(glob_edge_no(elem,2))-pos_z(glob_edge_no(elem,1)))./6;
    end;
end
