%
% 平面波(P波)入射の反射と透過のアニメーション
% 
% MATLAB R2015a
% 2016/9/14 Copyright (c) Takuichi Hirano
%
function p_wave_inc_ref_trans
    ci=complex(0,1);    % 虚数単位i
    eps0=8.854e-12;     % 真空の誘電率
    mu0=4.*pi.*1e-7;    % 真空の透磁率
    c=2.99792458e8;     % 光速

    lambda0=1.;         % 波長
    freq=c./lambda0;    % 周波数
    omega=2.*pi.*freq;  % 角周波数
    T=1./freq;          % 周期

    %----------------
    % 解析入力パラメータ
    er1=2.0;            % 媒質1の比誘電率
    mur1=1.0;           % 媒質1の比透磁率

    er2=1.0;            % 媒質2の比誘電率
    mur2=1.0;           % 媒質2の比透磁率
    sigma2=0.0;        % 媒質2の導電率
    er2c=er2+sigma2./(ci.*omega.*eps0); % 媒質2の複素比誘電率

    theta_i=35.2644*(pi./180.); % 入射角 (rad)
    %----------------
    % 描画・動画パラメータ
    xzmin=-1;           % 描画範囲の左・下 (波長単位)
    xzmax= 1;           % 描画範囲の右・上 (波長単位)

    n_grid_samp_field=15;   % quiver描画の格子間隔

    canvas_scale=1;     % 描画領域サイズのスケール(1: 画面サイズいっぱい)
    nframe=16;          % 1周期のフレーム数
    nrepeat=10;         % 動画の繰り返し回数

    avi_frate=2;        % フレーム レート (1 秒あたりに表示するフレームの数) 
    avi_quality=70;     % 画質設定 (0 から 100 までのパーセント値)
    avi_nrep=3;         % AVIファイル出力時の動画繰り返し回数
    %----------------

    % 臨界角の計算
    theta_c=asin(sqrt(er2/er1));
    % ブリュースター角の計算
    theta_b=asin(sqrt((er2/er1)/(1+(er2/er1))));

    k1=omega.*sqrt(er1.*eps0.*mur1.*mu0);
    k2=omega.*sqrt(er2c.*eps0.*mur2.*mu0);

    eta1=sqrt((mur1.*mu0)./(er1.*eps0));
    eta2=sqrt((mur2.*mu0)./(er2c.*eps0));

    theta_r=theta_i;                                    % 反射の法則
    theta_t=@(theta_i) asin((k1./k2).*sin(theta_i));    % スネルの法則
    % スネルの法則
    % sqrt(er1)*sin(theta_i)=sqrt(er2)*sin(theta_t)
    % 次のsin(theta_t), cos(theta_t)は全反射の場合のために計算しておく
    sin_theta_t=sqrt(er1./er2c).*sin(theta_i);
    cos_theta_t=sqrt(1.-(er1./er2c).*(sin(theta_i).^2));
    if imag(cos_theta_t)>0
        cos_theta_t=-cos_theta_t;
    end

    Rp=(eta1.*cos(theta_i)-eta2.*cos_theta_t) ...
        ./(eta1.*cos(theta_i)+eta2.*cos_theta_t);     % 反射係数
    Tp=(2.*eta2.*cos(theta_i)) ...
        ./(eta1.*cos(theta_i)+eta2.*cos_theta_t);     % 透過係数

    emax=(1+abs(Rp));           % 反対方向電界定在波の最大値
    hmax=(1+abs(Rp))/eta1;      % 反対方向磁界定在波の最大値

    % E and H fields
    hiy=@(X,Z) (1./eta1).*(exp(-ci.*k1.*(sin(theta_i).*X-cos(theta_i).*Z))) ...
        .*heaviside(Z);
    hry=@(X,Z) (Rp./eta1).*(exp(-ci.*k1.*(sin(theta_r).*X+cos(theta_r).*Z))) ...
        .*heaviside(Z);
    hty=@(X,Z) (Tp./eta2).*(exp(-ci.*k2.*(sin_theta_t.*X-cos_theta_t.*Z))) ...
        .*heaviside(-Z);

    eix=@(X,Z) exp(-ci.*k1.*(sin(theta_i).*X-cos(theta_i).*Z)) ...
        .*( cos(theta_i)).*heaviside(Z);
    eiz=@(X,Z) exp(-ci.*k1.*(sin(theta_i).*X-cos(theta_i).*Z)) ...
        .*( sin(theta_i)).*heaviside(Z);
    erx=@(X,Z) Rp.*exp(-ci.*k1.*(sin(theta_r).*X+cos(theta_r).*Z)) ...
        .*(-cos(theta_r)).*heaviside(Z);
    erz=@(X,Z) Rp.*exp(-ci.*k1.*(sin(theta_r).*X+cos(theta_r).*Z)) ...
        .*( sin(theta_r)).*heaviside(Z);
    etx=@(X,Z) Tp.*exp(-ci.*k2.*(sin_theta_t.*X-cos_theta_t.*Z)) ...
        .*( cos_theta_t).*heaviside(-Z);
    etz=@(X,Z) Tp.*exp(-ci.*k2.*(sin_theta_t.*X-cos_theta_t.*Z)) ...
        .*( sin_theta_t).*heaviside(-Z);

    % グラフィックを描画
    scrsz = get(groot,'ScreenSize');    % 画面サイズの取得
    canvas_width=(scrsz(3))*canvas_scale;
    canvas_height=(scrsz(4))*canvas_scale;
    font_size1=(canvas_width/800)*12;   % フォントサイズ1
    font_size2=(canvas_width/800)*10;   % フォントサイズ2
    fig=figure('Position',[(scrsz(3)/2-canvas_width/2) (scrsz(4)/2-canvas_height/2) canvas_width canvas_height]);

    [X,Z]=meshgrid(xzmin:xzmax./100:xzmax);
    [X2,Z2]=meshgrid(xzmin:xzmax./n_grid_samp_field:xzmax);

    for i=1:nframe
        clf;        % 現在のFigureウィンドウをクリア
        t=(T./nframe).*(i-1);

        % 入射波+透過波
        hy_it=real((hiy(X,Z)+hty(X,Z)).*exp(ci.*omega.*t));
        ex_it=real((eix(X,Z)+etx(X,Z)).*exp(ci.*omega.*t));
        ez_it=real((eiz(X,Z)+etz(X,Z)).*exp(ci.*omega.*t));

        hy_it2=real((hiy(X2,Z2)+hty(X2,Z2)).*exp(ci.*omega.*t));
        ex_it2=real((eix(X2,Z2)+etx(X2,Z2)).*exp(ci.*omega.*t));
        ez_it2=real((eiz(X2,Z2)+etz(X2,Z2)).*exp(ci.*omega.*t));

        % 反射波+透過波
        hy_rt=real((hry(X,Z)+hty(X,Z)).*exp(ci.*omega.*t));
        ex_rt=real((erx(X,Z)+etx(X,Z)).*exp(ci.*omega.*t));
        ez_rt=real((erz(X,Z)+etz(X,Z)).*exp(ci.*omega.*t));

        hy_rt2=real((hry(X2,Z2)+hty(X2,Z2)).*exp(ci.*omega.*t));
        ex_rt2=real((erx(X2,Z2)+etx(X2,Z2)).*exp(ci.*omega.*t));
        ez_rt2=real((erz(X2,Z2)+etz(X2,Z2)).*exp(ci.*omega.*t));

        % 入射波+反射波+透過波
        hy_irt=real((hiy(X,Z)+hry(X,Z)+hty(X,Z)).*exp(ci.*omega.*t));
        ex_irt=real((eix(X,Z)+erx(X,Z)+etx(X,Z)).*exp(ci.*omega.*t));
        ez_irt=real((eiz(X,Z)+erz(X,Z)+etz(X,Z)).*exp(ci.*omega.*t));

        hy_irt2=real((hiy(X2,Z2)+hry(X2,Z2)+hty(X2,Z2)).*exp(ci.*omega.*t));
        ex_irt2=real((eix(X2,Z2)+erx(X2,Z2)+etx(X2,Z2)).*exp(ci.*omega.*t));
        ez_irt2=real((eiz(X2,Z2)+erz(X2,Z2)+etz(X2,Z2)).*exp(ci.*omega.*t));

        %--------
        % Order of figures
        % 1 2
        %--------
        %---- 1
        subplot(1,2,1)                  % 複数グラフをタイル状に描画
        pcolor(X,Z,hy_irt);             % グラフを描画
        shading flat;                   % グリッドメッシュを削除
        colormap('jet');                % 色の設定
        caxis([-hmax hmax]);            % カラーマップの範囲を指定
        zlim([-hmax hmax]);             % z軸の範囲を指定
        axis square;                    % グラフ縦横比の設定
        xlabel({'x (wavelength)'});     % xlabel を作成
        ylabel({'z (wavelength)'});     % ylabel を作成
        title({'H-field'});               % title を作成

        %-- 入射・反射・透過平面波の向きの矢印を描く
        hold on;                        % 重ね書きのため
        r=(xzmax-xzmin)/4;
        v=quiver(-r*sin(theta_i),r*cos(theta_i),r*sin(theta_i),-r*cos(theta_i),0);  % 入射平面波の向き
        v.Color='w';
        v.LineWidth=2;
        v.MaxHeadSize=0.5;
        v=quiver(0,0,r*sin(theta_r),r*cos(theta_r),0);  % 反射平面波の向き
        v.Color='w';
        v.LineWidth=2;
        v.MaxHeadSize=0.5;
        v=quiver(0,0,r*sin(real(theta_t(theta_i))),-r*cos(real(theta_t(theta_i))),0);  % 透過平面波の向き
        v.Color='w';
        v.LineWidth=2;
        v.MaxHeadSize=0.5;
        hold off;
        set(gca,'FontSize',font_size1); % フォントサイズ設定

        %---- 2
        subplot(1,2,2)                  % 複数グラフをタイル状に描画

        hold on;                        % 重ね書きのため
        scale_fact=((3/emax)/(xzmax-xzmin))/n_grid_samp_field;       % 矢印の長さの自動スケーリング
        h=quiver(X2,Z2,scale_fact*ex_irt2,scale_fact*ez_irt2);
        h.AutoScale='off';              % 矢印の長さの自動スケーリング
        h.Color='b';                    % 矢印の色
        h.MaxHeadSize=10;
        h.Clipping='on';
        hold off;

        axis square;                    % グラフ縦横比の設定
        axis([xzmin xzmax xzmin xzmax]);% データの範囲と描画範囲を一致させる
        xlabel({'x (wavelength)'});     % xlabel を作成
        ylabel({'z (wavelength)'});     % ylabel を作成
        title({'E-field'});             % title を作成

        % 説明用 textbox を描画
        tb1=annotation('textbox',...
            [0.1 0.84 0.35 0.125],...
            'String',{sprintf('Inc. angle (deg): %4.1f',theta_i.*(180./pi)), ...
                      sprintf('Trans. angle (deg): %4.1f, Complex: %.3f+j(%.3f)',theta_t(theta_i).*(180./pi),real(theta_t(theta_i)),imag(theta_t(theta_i))), ...
                      sprintf('Critical angle (deg): %4.1f',theta_c.*(180./pi)), ...
                      sprintf('Brewster angle (deg): %4.1f',theta_b.*(180./pi))});
        tb1.FontSize=font_size2;        % フォントサイズ設定
        tb1.EdgeColor='none';           % 枠の色
        %tb1.Color='w';                  % テキストの色
        %tb1.FontWeight='bold';          % 太字
        %tb.BackgroundColor='w';         % 背景色
        %tb.FaceAlpha=0.3;               % 透明度(透明:0～不透明:1)
        tb2=annotation('textbox',...
            [0.5 0.9 0.15 0.06],...
            'String',{sprintf('R: %.3f+j(%.3f)',real(Rp),imag(Rp)), ...
                      sprintf('T: %.3f+j(%.3f)',real(Tp),imag(Tp))});
        tb2.EdgeColor='none';           % 枠の色
        tb2.FontSize=font_size2;        % フォントサイズ設定
        if sigma2>0
            tb3=annotation('textbox',...
                [0.7 0.88 0.22 0.08],...
                'String',{sprintf('Skin depth (wl.): %.3f',sqrt(2/(omega*(mur2*mu0)*sigma2))), ...
                          sprintf('Sigma (skin depth=1 wl.) (S/m): %.3f',(2/(omega*(mur2*mu0)))/((1)^2))});
            tb3.EdgeColor='none';           % 枠の色
            tb3.FontSize=font_size2;        % フォントサイズ設定
        end
        tb4=annotation('textbox',...
            [0.1 0.05 0.15 0.1],...
            'String',{sprintf('er1   : %.2f',er1), ...
                      sprintf('mur1  : %.2f',mur1)});
        tb4.EdgeColor='none';           % 枠の色
        tb4.FontSize=font_size2;        % フォントサイズ設定
        tb5=annotation('textbox',...
            [0.3 0.05 0.15 0.1],...
            'String',{sprintf('er2   : %.2f',er2), ...
                      sprintf('mur2  : %.2f',mur2), ...
                      sprintf('sigma2: %.3f',sigma2)});
        tb5.EdgeColor='none';           % 枠の色
        tb5.FontSize=font_size2;        % フォントサイズ設定

        set(gca,'FontSize',font_size1); % フォントサイズ設定
        drawnow;                        % 描画する
        frm(i)=getframe(gcf);           % gcf: 現在の Figure ハンドル番号, getframe: ムービーフレームのキャプチャ
    end
    movie(fig,frm,nrepeat);             % 動画の再生(最後の引数は繰り返し回数)

    % 動画AVIファイルの生成
    v = VideoWriter('p_wave.avi');
    v.FrameRate = avi_frate;    % Default 30
    v.Quality = avi_quality;    % Default 75
    open(v);
    for k=1:avi_nrep
        for i=1:nframe
            writeVideo(v,frm(i));
        end
    end
    close(v);
end
