mtlab pml.txt

摩尔边界的源码在计算电磁学中有很大应用

% CHx_Py
CHx_Py_III=(1-eaf_yH)./(1+eaf_yH);
CHx_Py_III=repmat(CHx_Py_III,1,Lx);
CHx_Py_IV=flipud(CHx_Py_III);
% CEzHx_Py
CEzHx_Py_III=0.5/mu_r./(1+eaf_yH);
CEzHx_Py_III=repmat(CEzHx_Py_III,1,Lx);
CEzHx_Py_IV=flipud(CEzHx_Py_III);
% CEzBy_Py
CEzBy_Py=0.5/mu_r;
% CByHy_New_Py
CByHy_New_Py_III=1+eaf_yE;
CByHy_New_Py_III=repmat(CByHy_New_Py_III,1,Lx+1);
CByHy_New_Py_IV=flipud(CByHy_New_Py_III);
% CByHy_Old_Py
CByHy_Old_Py_III=1-eaf_yE;
CByHy_Old_Py_III=repmat(CByHy_Old_Py_III,1,Lx+1);
CByHy_Old_Py_IV=flipud(CByHy_Old_Py_III);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Clear the variables no longer needed
clear Npml Lx Ly c0 dx epsz sigma_M
clear eaf_xE eaf_xH eaf_yE eaf_yH sigma_xE 
clear sigma_xH sigma_yE sigma_yH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FDTD loop
figure(1)
if SeeContour==1        % Plot the FDTD lattice
    line([1,KX],[Tpml,Tpml],'color','red');
    line([1,KX],[Bpml,Bpml],'color','red');
    line([Lpml,Lpml],[1,KY],'color','red');
    line([Rpml,Rpml],[1,KY],'color','red');
    hold on
end
    
while n<=T_Des-N_Steps
    for n=n_start:NSTEPS
        
        % Update the Ez field
        % Zone A
        DzA(1:end-1,2:end)=CDz_A(1:end-1,2:end).*DzA(1:end-1,2:end)...
            +CHxDz_A(1:end-1,2:end).*(Hx(Tpml-1:end-1,2:Lpml)-Hx(Tpml:end,2:Lpml)...
            -Hy(Tpml:end-1,1:Lpml-1)+Hy(Tpml:end-1,2:Lpml));
        
        Ez(Tpml:end-1,2:Lpml)=CEz_A(1:end-1,2:end).*Ez(Tpml:end-1,2:Lpml)...
            +CDzEz_A(1:end-1,2:end).*(DzA(1:end-1,2:end)-DzA_prev(1:end-1,2:end));
        
        % Zone B
        DzB(2:end,2:end)=CDz_B(2:end,2:end).*DzB(2:end,2:end)...
            +CHxDz_B(2:end,2:end).*(Hx(1:Bpml-1,2:Lpml)-Hx(2:Bpml,2:Lpml)...
            -Hy(2:Bpml,1:Lpml-1)+Hy(2:Bpml,2:Lpml));
        
        Ez(2:Bpml,2:Lpml)=CEz_B(2:end,2:end).*Ez(2:Bpml,2:Lpml)...
            +CDzEz_B(2:end,2:end).*(DzB(2:end,2:end)-DzB_prev(2:end,2:end));

        % Zone C
        DzC(2:end,1:end-1)=CDz_C(2:end,1:end-1).*DzC(2:end,1:end-1)...
            +CHxDz_C(2:end,1:end-1).*(Hx(1:Bpml-1,Rpml:end-1)-Hx(2:Bpml,Rpml:end-1)...
            -Hy(2:Bpml,Rpml-1:end-1)+Hy(2:Bpml,Rpml:end));
        
        Ez(2:Bpml,Rpml:end-1)=CEz_C(2:end,1:end-1).*Ez(2:Bpml,Rpml:end-1)...
            +CDzEz_C(2:end,1:end-1).*(DzC(2:end,1:end-1)-DzC_prev(2:end,1:end-1));
        
        % Zone D
        DzD(1:end-1,1:end-1)=CDz_D(1:end-1,1:end-1).*DzD(1:end-1,1:end-1)...
            +CHxDz_D(1:end-1,1:end-1).*(Hx(Tpml-1:end-1,Rpml:end-1)-Hx(Tpml:end,Rpml:end-1)...
            -Hy(Tpml:end-1,Rpml-1:end-1)+Hy(Tpml:end-1,Rpml:end));
        
        Ez(Tpml:end-1,Rpml:end-1)=CEz_D(1:end-1,1:end-1).*Ez(Tpml:end-1,Rpml:end-1)...
            +CDzEz_D(1:end-1,1:end-1).*(DzD(1:end-1,1:end-1)-DzD_prev(1:end-1,1:end-1));
        
        % Zone I
        Ez(Bpml+1:Tpml-1,2:Lpml)=CEz_Px_I(:,2:end).*Ez(Bpml+1:Tpml-1,2:Lpml)...
            +CHxEz_Px_I(:,2:end).*(Hx(Bpml:Tpml-2,2:Lpml)-Hx(Bpml+1:Tpml-1,2:Lpml)...
            -Hy(Bpml+1:Tpml-1,1:Lpml-1)+Hy(Bpml+1:Tpml-1,2:Lpml));
        
        % Zone II
        Ez(Bpml+1:Tpml-1,Rpml:end-1)=CEz_Px_II(:,1:end-1).*Ez(Bpml+1:Tpml-1,Rpml:end-1)...
            +CHxEz_Px_II(:,1:end-1).*(Hx(Bpml:Tpml-2,Rpml:end-1)-Hx(Bpml+1:Tpml-1,Rpml:end-1)...
            -Hy(Bpml+1:Tpml-1,Rpml-1:end-1)+Hy(Bpml+1:Tpml-1,Rpml:end));
        
        % Zone III
        Ez(Tpml:end-1,Lpml+1:Rpml-1)=CEz_Py_III(1:end-1,:).*Ez(Tpml:end-1,Lpml+1:Rpml-1)...
            +CHxEz_Py_III(1:end-1,:).*(Hx(Tpml-1:end-1,Lpml+1:Rpml-1)-Hx(Tpml:end,Lpml+1:Rpml-1)...
            -Hy(Tpml:end-1,Lpml:Rpml-2)+Hy(Tpml:end-1,Lpml+1:Rpml-1));
        
        % Zone IV
        Ez(2:Bpml,Lpml+1:Rpml-1)=CEz_Py_IV(2:end,:).*Ez(2:Bpml,Lpml+1:Rpml-1)...
            +CHxEz_Py_IV(2:end,:).*(Hx(1:Bpml-1,Lpml+1:Rpml-1)-Hx(2:Bpml,Lpml+1:Rpml-1)...
            -Hy(2:Bpml,Lpml:Rpml-2)+Hy(2:Bpml,Lpml+1:Rpml-1));
        
        % Zone FDTD
        Ez(Bpml+1:Tpml-1,Lpml+1:Rpml-1)=Ez(Bpml+1:Tpml-1,Lpml+1:Rpml-1)...
            +0.5/ep_r*(Hx(Bpml:Tpml-2,Lpml+1:Rpml-1)-Hx(Bpml+1:Tpml-1,Lpml+1:Rpml-1))...
            -0.5/ep_r*(Hy(Bpml+1:Tpml-1,Lpml:Rpml-2)-Hy(Bpml+1:Tpml-1,Lpml+1:Rpml-1));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Source excitation
        pulse=6.477*exp(-0.5*(t0-n)^2/spread^2);
        %pulse=sin(2*pi*freq*n*dt);
        Ez(kcx,kcy)=pulse+Ez(kcx,kcy);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Update the Hx field
        % Zone A
        BxA=CBx_A.*BxA+CEzBx_A.*(Ez(Tpml:end-1,1:Lpml)-Ez(Tpml+1:end,1:Lpml));
        Hx(Tpml:end,1:Lpml)=Hx(Tpml:end,1:Lpml)+CBxHx_New_A.*BxA-CBxHx_Old_A.*BxA_prev;
        
        % Zone B
        BxB=CBx_B.*BxB+CEzBx_B.*(Ez(1:Bpml-1,1:Lpml)-Ez(2:Bpml,1:Lpml));
        Hx(1:Bpml-1,1:Lpml)=Hx(1:Bpml-1,1:Lpml)+CBxHx_New_B.*BxB-CBxHx_Old_B.*BxB_prev;
        
        % Zone C
        BxC=CBx_C.*BxC+CEzBx_C.*(Ez(1:Bpml-1,Rpml:end)-Ez(2:Bpml,Rpml:end));
        Hx(1:Bpml-1,Rpml:end)=Hx(1:Bpml-1,Rpml:end)+CBxHx_New_C.*BxC-CBxHx_Old_C.*BxC_prev;
        
        % Zone D
        BxD=CBx_D.*BxD+CEzBx_D.*(Ez(Tpml:end-1,Rpml:end)-Ez(Tpml+1:end,Rpml:end));
        Hx(Tpml:end,Rpml:end)=Hx(Tpml:end,Rpml:end)+CBxHx_New_D.*BxD-CBxHx_Old_D.*BxD_prev;
        
        % Zone I
        BX1=BX1+CEzBx_Px*(Ez(Bpml:Tpml-1,1:Lpml)-Ez(Bpml+1:Tpml,1:Lpml));
        Hx(Bpml:Tpml-1,1:Lpml)=Hx(Bpml:Tpml-1,1:Lpml)+CBxHx_New_Px_I.*BX1...
            -CBxHx_Old_Px_I.*BX1_prev;
        
        % Zone II
        BX2=BX2+CEzBx_Px*(Ez(Bpml:Tpml-1,Rpml:end)-Ez(Bpml+1:Tpml,Rpml:end));
        Hx(Bpml:Tpml-1,Rpml:end)=Hx(Bpml:Tpml-1,Rpml:end)+CBxHx_New_Px_II.*BX2...
            -CBxHx_Old_Px_II.*BX2_prev;
        
        % Zone III
        Hx(Tpml:end,Lpml+1:Rpml-1)=CHx_Py_III.*Hx(Tpml:end,Lpml+1:Rpml-1)...
            +CEzHx_Py_III.*(Ez(Tpml:end-1,Lpml+1:Rpml-1)-Ez(Tpml+1:end,Lpml+1:Rpml-1));
        
        % Zone IV
        Hx(1:Bpml-1,Lpml+1:Rpml-1)=CHx_Py_IV.*Hx(1:Bpml-1,Lpml+1:Rpml-1)...
            +CEzHx_Py_IV.*(Ez(1:Bpml-1,Lpml+1:Rpml-1)-Ez(2:Bpml,Lpml+1:Rpml-1));
        
        % Zone FDTD
        Hx(Bpml:Tpml-1,Lpml+1:Rpml-1)=Hx(Bpml:Tpml-1,Lpml+1:Rpml-1)...
            +0.5/mu_r*(Ez(Bpml:Tpml-1,Lpml+1:Rpml-1)-Ez(Bpml+1:Tpml,Lpml+1:Rpml-1));
        
        % Update the Hy field
        % Zone A
        ByA=CBy_A.*ByA-CEzBy_A.*(Ez(Tpml:end,1:Lpml-1)-Ez(Tpml:end,2:Lpml));
        Hy(Tpml:end,1:Lpml-1)=Hy(Tpml:end,1:Lpml-1)+CByHy_New_A.*ByA-CByHy_Old_A.*ByA_prev;
        
        % Zone B
        ByB=CBy_B.*ByB-CEzBy_B.*(Ez(1:Bpml,1:Lpml-1)-Ez(1:Bpml,2:Lpml));
        Hy(1:Bpml,1:Lpml-1)=Hy(1:Bpml,1:Lpml-1)+CByHy_New_B.*ByB-CByHy_Old_B.*ByB_prev;
        
        % Zone C
        ByC=CBy_C.*ByC-CEzBy_C.*(Ez(1:Bpml,Rpml:end-1)-Ez(1:Bpml,Rpml+1:end));
        Hy(1:Bpml,Rpml:end)=Hy(1:Bpml,Rpml:end)+CByHy_New_C.*ByC-CByHy_Old_C.*ByC_prev;
        
        % Zone D
        ByD=CBy_D.*ByD-CEzBy_D.*(Ez(Tpml:end,Rpml:end-1)-Ez(Tpml:end,Rpml+1:end));
        Hy(Tpml:end,Rpml:end)=Hy(Tpml:end,Rpml:end)+CByHy_New_D.*ByD-CByHy_Old_D.*ByD_prev;
        
        % Zone I
        Hy(Bpml+1:Tpml-1,1:Lpml-1)=CHy_Px_I.*Hy(Bpml+1:Tpml-1,1:Lpml-1)...
            -CEzHy_Px_I.*(Ez(Bpml+1:Tpml-1,1:Lpml-1)-Ez(Bpml+1:Tpml-1,2:Lpml));
        
        % Zone II
        Hy(Bpml+1:Tpml-1,Rpml:end)=CHy_Px_II.*Hy(Bpml+1:Tpml-1,Rpml:end)...
            -CEzHy_Px_II.*(Ez(Bpml+1:Tpml-1,Rpml:end-1)-Ez(Bpml+1:Tpml-1,Rpml+1:end));
        
        % Zone III
        BY3=BY3-CEzBy_Py*(Ez(Tpml:end,Lpml:Rpml-1)-Ez(Tpml:end,Lpml+1:Rpml));
        Hy(Tpml:end,Lpml:Rpml-1)=Hy(Tpml:end,Lpml:Rpml-1)+CByHy_New_Py_III.*BY3...
            -CByHy_Old_Py_III.*BY3_prev;
        
        % Zone IV
        BY4=BY4-CEzBy_Py*(Ez(1:Bpml,Lpml:Rpml-1)-Ez(1:Bpml,Lpml+1:Rpml));
        Hy(1:Bpml,Lpml:Rpml-1)=Hy(1:Bpml,Lpml:Rpml-1)+CByHy_New_Py_IV.*BY4...
            -CByHy_Old_Py_IV.*BY4_prev;
        
        % Zone FDTD
        Hy(Bpml+1:Tpml-1,Lpml:Rpml-1)=Hy(Bpml+1:Tpml-1,Lpml:Rpml-1)...
            -0.5/mu_r*(Ez(Bpml+1:Tpml-1,Lpml:Rpml-1)-Ez(Bpml+1:Tpml-1,Lpml+1:Rpml));    
        
        % Record the previous values
        DzA_prev=DzA;
        DzB_prev=DzB;
        DzC_prev=DzC;
        DzD_prev=DzD;
        
        BxA_prev=BxA;
        BxB_prev=BxB;
        BxC_prev=BxC;
        BxD_prev=BxD;
        
        ByA_prev=ByA;
        ByB_prev=ByB;
        ByC_prev=ByC;
        ByD_prev=ByD;
        
        BX1_prev=BX1;
        BX2_prev=BX2;
        
        BY3_prev=BY3;
        BY4_prev=BY4;
        
        % Sample the field values
        Ez_samp(n)=Ez(kcx,kcy);
        Hx_samp(n)=Hx(kcx,kcy);
        Hy_samp(n)=Hy(kcx,kcy);
    end
    % Now we plot the Ez field
    if SeeContour==1
        contour(Ez,0.02,'b');
        title(['T=',num2str(n)])
        axis([1 KX 1 KY])  
        axis square
        pause
    else
        surf(Hy)%
        xlabel('X')
        ylabel('Y')
        title(['T=',num2str(n)])
        axis([1 KX 1 KY -1 1])
        %axis([1 KY -0.2 0.2])
        %axis tight%square
        %view(2)
    end
    drawnow;    
    
    NSTEPS=NSTEPS+N_Steps;
    n_start=n+1;
end
hold off