%
% 平面波（直線偏波）のアニメーション
% 
% MATLAB R2015a
% 2016/7/25 Copyright (c) Takuichi Hirano
%
function plane_wave
    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;          % 周期

    %----------------
    % 描画・動画パラメータ
    canvas_scale=1;     % 描画領域サイズのスケール(1: 画面サイズいっぱい)
    %----------------

    k0=omega.*sqrt(eps0.*mu0);  % 自由空間の波数
    eta0=sqrt(mu0./eps0);       % 自由空間の界インピーダンス

    % グラフィックを描画
    scrsz = get(groot,'ScreenSize');    % 画面サイズの取得
    canvas_width=(scrsz(3))*canvas_scale;
    canvas_height=(scrsz(4))*canvas_scale;
    font_size=(canvas_width/800)*12;    % フォントサイズ
    fig=figure('Position',[(scrsz(3)/2-canvas_width/2) (scrsz(4)/2-canvas_height/2) canvas_width canvas_height]);
    % 電界用プロット点
    [X1,Y1] = meshgrid(0:1:0.1,-1.5:0.05:1.5);
    Z1 = 0.*(X1+Y1);
    % 磁界用プロット点
    [Z2,Y2] = meshgrid(0:1:0.1,-1.5:0.05:1.5);
    X2 = 0.*(X1+Y1);

    % Ex=Ey=0
    Ex=@(X,Y) 0.*(X1+Y1);
    Ey=@(X,Y) 0.*(X1+Y1);
    % Hz=Hy=0
    Hz=@(Y,Z) 0.*(Y2+Z2);
    Hy=@(Y,Z) 0.*(Y2+Z2);
    
    nframe=16;
    for i=1:nframe
        clf;        % 現在のFigureウィンドウをクリア
        t=(T./nframe).*(i-1);

        Ez=@(X,Y) real(0.*(X1+Y1)+exp(-ci.*k0.*Y1).*exp(ci.*omega.*t));

        % Figure
        U=Ex(X1,Y1);
        V=Ey(X1,Y1);
        W=Ez(X1,Y1);
        h = quiver3(X1,Y1,Z1,U,V,W,2);
        set(h,'Color','r');
        axis equal;
        xlim([-1,1]);
        ylim([-1.5,1.5]);
        zlim([-1,1]);

        Hx=@(Y,Z) real(0.*(Y2+Z2)+0.8.*exp(-ci.*k0.*Y2).*exp(ci.*omega.*t));

        U=Hx(Y2,Z2);
        V=Hy(Y2,Z2);
        W=Hz(Y2,Z2);
        hold on;
        h = quiver3(X2,Y2,Z2,U,V,W,2);
        set(h,'Color','b');
        hold off;

        set(gca,'FontSize',font_size);  % フォントサイズ設定
        drawnow
        frm(i)=getframe(gcf);
    end
    movie(fig,frm,10);                  % 動画の再生(最後の引数は繰り返し回数)
    
    % 動画のファイル保存
    v = VideoWriter('plane_wave.avi');
    v.FrameRate = 4;	% Default 30
    v.Quality = 75;     % Default 75
    open(v);
    nrep=3;            % 繰り返し回数
    for k=1:nrep
        for i=1:nframe
            writeVideo(v,frm(i));
        end
    end
    close(v);
end


