planewavegui.m

This is the procedure for lab 1. This is a two-week lab. Prelab should be done BEFORE going to the l

eps0=1/(v0^2*mu0);
omega=2*pi*freq;
eps=eps0*er;
mu=mu0*mr;
kq=omega^2.*mu.*(eps-j.*sigma/omega);
k=sqrt(kq);
beta=real(k);
alfa=-imag(k);
wl=(2*pi)./real(k);

% define ky based on the input choice
if ev==0
    ky=real(k(1))*sind(teta);
elseif ev==1
    ky=teta*k(1);
end

for q=1:Nl   % for all layers
    if real(er(q))<0&real(mr(q))<0
        if sigma(q)==0&(ky.^2>kq(q))
            kz(q)=sqrt(kq(q)-ky.^2);
        else
            kz(q)=-sqrt(kq(q)-ky.^2);
        end
    else
        kz(q)=sqrt(kq(q)-ky.^2);
    end
    yTE(q)=kz(q)./(omega.*mu(q));
    zTE(q)=1./yTE(q);
    yTM(q)=(omega.*(eps(q)-j*sigma(q)/omega))./kz(q);
    zTM(q)=1./yTM(q);
end

if LastLay==0   % the last layer is an indefinite dielectric/conductor
    gammaTEdx(Nl-1)=(zTE(end)-zTE(end-1))./(zTE(end)+zTE(end-1));
    gammaTMdx(Nl-1)=-(zTM(end)-zTM(end-1))./(zTM(end)+zTM(end-1));
elseif LastLay==-1   % the last layer is a PEC (electric gamma=-1)
    gammaTEdx(Nl-1)=-1;
    gammaTMdx(Nl-1)=1;
elseif LastLay==1   % the last layer is a PMC (electric gamma=1)
    gammaTEdx(Nl-1)=1;
    gammaTMdx(Nl-1)=-1;
end
for q=Nl-1:-1:2   % for all layers
    gammaTEsx(q)=gammaTEdx(q).*exp(-j.*2.*kz(q).*(-(Zlay(q+1)-Zlay(q))));
    gammaTMsx(q)=gammaTMdx(q).*exp(-j.*2.*kz(q).*(-(Zlay(q+1)-Zlay(q))));
    zTEL(q)=zTE(q).*(1+gammaTEsx(q))./(1-gammaTEsx(q));
    zTML(q)=zTM(q).*(1+gammaTMsx(q))./(1-gammaTMsx(q));
    gammaTEdx(q-1)=(zTEL(q)-zTE(q-1))./(zTEL(q)+zTE(q-1));
    gammaTMdx(q-1)=(zTML(q)-zTM(q-1))./(zTML(q)+zTM(q-1));
end

if TE_TM==1   % TE wave

    clear J;
    J=find(z>=0);
    % primary field (only on plane yz, with z >= 0)
    EpTE(J,:)=A0.*exp(-j*kz(1)*(-z(J)).')*exp(-j*ky*y);
    ApEsTEDx=A0.*exp(-j*ky*y);

    EsTE_progr_only(J,:)=EpTE(J,:);

    clear J;
    J=find(z<0);
    EpTE(J,:)=zeros(length(J),length(y));

    % reflected field for layer 1
    clear J;
    J=find(z>=0);
    znow=z(J);

    % TE case
    EsTE(J,:)=exp(-j*kz(1)*(znow).')*ApEsTEDx*gammaTEdx(1);

    % save only regressive field in layer 1
    EsTE_regr_only(J,:)=EsTE(J,:);

    for q=2:Nl-1   % progressive and regressive fields for layers 2 --> N-1
        clear J;
        J=find(z<Zlay(q)&z>=Zlay(q+1));
        znow=z(J);

        % TE case
        ApEsTESx=ApEsTEDx.*(1+gammaTEdx(q-1))./(1+gammaTEsx(q));
        EsTE(J,:)=exp(-j*kz(q)*(-(znow.'-Zlay(q))))*ApEsTESx;

        % save only progressive field in all intermediate layers
        EsTE_progr_only(J,:)=EsTE(J,:);

        ApEsTEDx=ApEsTESx.*exp(-j*kz(q)*(-(Zlay(q+1)-Zlay(q))));
        AmEsTEDx=ApEsTEDx.*gammaTEdx(q);
        EsTE(J,:)=EsTE(J,:)+exp(-j*kz(q)*((znow.'-Zlay(q+1))))*AmEsTEDx;

        % save only regressive field in all intermediate layers
        EsTE_regr_only(J,:)=exp(-j*kz(q)*((znow.'-Zlay(q+1))))*AmEsTEDx;
    end

    % only progressive field in the last layer
    clear J;
    J=find(z<Zlay(end));
    znow=z(J);

    if LastLay==1|LastLay==-1
        EsTE(J,:)=0;
    else
        % TE case
        EsTE(J,:)=exp(-j*kz(end)*(-(znow.'-Zlay(end))))*ApEsTEDx*(1+gammaTEdx(end));

        % save only progressive field in the last layer
        EsTE_progr_only(J,:)=EsTE(J,:);
    end
elseif TE_TM==2   % TM wave
    clear J;
    J=find(z>=0);
    % primary field (only on plane yz, with z >= 0)
    EpTE(J,:)=A0.*exp(-j*kz(1)*(-z(J)).')*exp(-j*ky*y);
    ApEsTEDx=A0.*exp(-j*ky*y);

    EsTE_progr_only(J,:)=EpTE(J,:);

    clear J;
    J=find(z<0);
    EpTE(J,:)=zeros(length(J),length(y));

    % reflected field for layer 1
    clear J;
    J=find(z>=0);
    znow=z(J);

    % TE case
    EsTE(J,:)=exp(-j*kz(1)*(znow).')*ApEsTEDx*gammaTMdx(1);

    % save only regressive field in layer 1
    EsTE_regr_only(J,:)=EsTE(J,:);

    for q=2:Nl-1   % progressive and regressive fields for layers 2 --> N-1
        clear J;
        J=find(z<Zlay(q)&z>=Zlay(q+1));
        znow=z(J);

        % TE case
        ApEsTESx=ApEsTEDx.*(1+gammaTMdx(q-1))./(1+gammaTMsx(q));
        EsTE(J,:)=exp(-j*kz(q)*(-(znow.'-Zlay(q))))*ApEsTESx;

        % save only progressive field in all intermediate layers
        EsTE_progr_only(J,:)=EsTE(J,:);

        ApEsTEDx=ApEsTESx.*exp(-j*kz(q)*(-(Zlay(q+1)-Zlay(q))));
        AmEsTEDx=ApEsTEDx.*gammaTMdx(q);
        EsTE(J,:)=EsTE(J,:)+exp(-j*kz(q)*((znow.'-Zlay(q+1))))*AmEsTEDx;

        % save only regressive field in all intermediate layers
        EsTE_regr_only(J,:)=exp(-j*kz(q)*((znow.'-Zlay(q+1))))*AmEsTEDx;
    end

    % only progressive field in the last layer
    clear J;
    J=find(z<Zlay(end));
    znow=z(J);

    if LastLay==1|LastLay==-1
        EsTE(J,:)=0;
    else
        % TE case
        EsTE(J,:)=exp(-j*kz(end)*(-(znow.'-Zlay(end))))*ApEsTEDx*(1+gammaTMdx(end));

        % save only progressive field in the last layer
        EsTE_progr_only(J,:)=EsTE(J,:);
    end
end

if what==3
    % sum fields
    EtTE=EpTE+EsTE;

    if RoA==1   % plot real values
        if anim==1   % animate fields (only vertical cut)

            FieldT=EtTE;
            FieldP=EsTE_progr_only;
            FieldR=EsTE_regr_only;

            axes(fH(1))
            imagesc(y,z,real(FieldT))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Total field E_x')
            elseif TE_TM==2
                title('Total field H_x')
            end

            axes(fH(2))
            imagesc(y,z,real(FieldP))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Progressive field E_x')
            elseif TE_TM==2
                title('Regressive field H_x')
            end

            axes(fH(3))
            imagesc(y,z,real(FieldR))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Regressive field E_x')
            elseif TE_TM==2
                title('Regressive field H_x')
            end

            % movie
            tv=0;
            dT=1./(32.*freq);
            risp=1;
            while risp==1
                tv=tv+dT;
                axes(fH(4))
                % plot the field for a fixed value of y (center of the axis) REAL
                plot(z,real(EtTE(:,round(ns/2)).*exp(j.*omega.*tv)),'LineWidth',1.5)
                axis([z(1) z(end) -2*A0 2*A0])
                hold on;
                for q=2:Nl
                    h=line([Zlay(q) Zlay(q)],[-2*A0 2*A0]);
                    set(h,'Color','k','LineWidth',1)
                end
                xlabel('z [m]')
                if TE_TM==1
                    ylabel('real(E) [V/m]')
                    title('E_x field for y = 0 m')
                elseif TE_TM==2
                    ylabel('real(H) [A/m]')
                    title('H_x field for y = 0 m')
                end
                set(gca,'XDir','Reverse')
                hold off;

            end
        else   % static fields
            FieldT=EtTE;
            FieldP=EsTE_progr_only;
            FieldR=EsTE_regr_only;

            axes(fH(1))
            imagesc(y,z,real(FieldT))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Total field E_x')
            elseif TE_TM==2
                title('Total field H_x')
            end

            axes(fH(2))
            imagesc(y,z,real(FieldP))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Progressive field E_x')
            elseif TE_TM==2
                title('Regressive field H_x')
            end

            axes(fH(3))
            imagesc(y,z,real(FieldR))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Regressive field E_x')
            elseif TE_TM==2
                title('Regressive field H_x')
            end

            axes(fH(4))
            % plot the field for a fixed value of y (center of the axis) REAL
            plot(z,real(EtTE(:,round(ns/2))),'LineWidth',1.5)
            axis([z(1) z(end) -2*A0 2*A0])
            hold on;
            for q=2:Nl
                h=line([Zlay(q) Zlay(q)],[-2*A0 2*A0]);
                set(h,'Color','k','LineWidth',1)
            end
            xlabel('z [m]')
            if TE_TM==1
                ylabel('real(E) [V/m]')
                title('E_x field for y = 0 m')
            elseif TE_TM==2
                ylabel('real(H) [A/m]')
                title('H_x field for y = 0 m')
            end
            set(gca,'XDir','Reverse')
            hold off;
        end
    else   % plot absolute values
        if anim==1   % animate fields (only vertical cut)
            FieldT=EtTE;
            FieldP=EsTE_progr_only;
            FieldR=EsTE_regr_only;

            axes(fH(1))
            imagesc(y,z,abs(FieldT))
            xlabel('y [m]')
            ylabel('z [m]')
            set(gca,'YDir','normal')
            axis square
            hold on;
            for q=2:Nl
                h=line([y(1) y(end)],[Zlay(q) Zlay(q)]);
                set(h,'Color','k','LineWidth',1)
            end
            colorbar
            hold off;
            if TE_TM==1
                title('Total field E_x')
            elseif TE_TM==2
                title('Total field H_x')
            end

            axes(fH(2))
            imagesc(y,z,abs(FieldP))
            xlabel('y [m]')