mur1order.m

极为实用的TM波一阶Mur吸收边界条件编程

clear,clc
IE=100;
JE=100;
ic=IE/2;
jc=JE/2;

n0=50;
spread=12;

Ez=zeros(IE+1,JE+1);
Hx=zeros(IE+1,JE);
Hy=zeros(IE,JE+1);

% Matrices for the Mur ABCs
Ez_prev=Ez;

N_Steps=10;      % Increment for NSTEPS
N_Temp=N_Steps;  % Temporary variable for time steps
T_Des=300;       % Destination value of time steps
n_start=1;        % Record the start value of n
n=n_start;

cF=(1-2*sqrt(2))/(1+2*sqrt(2));
while n<=T_Des-N_Steps
    for n=n_start:N_Temp
        % Update the Ez field
        Ez(2:IE,2:JE)=Ez(2:IE,2:JE)+0.5*(Hy(2:IE,2:JE)-...
            Hy(1:IE-1,2:JE)-Hx(2:IE,2:JE)+Hx(2:IE,1:JE-1));
    
        pulse=exp(-(n-n0)^2/spread^2);
        Ez(ic,jc)=pulse+Ez(ic,jc);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Mur ABCs   ------> 4 Edges
        % Left Boundary
        Ez(1,2:JE)=Ez_prev(2,2:JE)-1/3*(Ez(2,2:JE)-Ez_prev(1,2:JE));%normalized H 
        % Right Boundary
        Ez(IE+1,2:JE)=Ez_prev(IE,2:JE)-1/3*(Ez(IE,2:JE)-Ez_prev(IE+1,2:JE));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Top Boundary
        Ez(2:IE,JE+1)=Ez_prev(2:IE,JE)-1/3*(Ez(2:IE,JE)-Ez_prev(2:IE,JE+1));
        % Bottom Boundary
        Ez(2:IE,1)=Ez_prev(2:IE,2)    -1/3*(Ez(2:IE,2)-Ez_prev(2:IE,1));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Mur ABCs   ------> 4 Corners
        % Left_Bottom Corner
        Ez(1,1)=Ez_prev(2,2)+cF*(Ez(2,2)-Ez_prev(1,1));
        % Left_Top Corner
        Ez(1,JE+1)=Ez_prev(2,JE)+cF*(Ez(2,JE)-Ez_prev(1,JE+1));
        % Right_Bottom Corner
        Ez(IE+1,1)=Ez_prev(IE,2)+cF*(Ez(IE,2)-Ez_prev(IE+1,1));
        % Right_Top Corner
        Ez(IE+1,JE+1)=Ez_prev(IE,JE)+cF*(Ez(IE,JE)-Ez_prev(IE+1,JE+1));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        Ez_prev=Ez;
     
        Hx(1:IE+1,1:JE)=Hx(1:IE+1,1:JE)-0.5*(Ez(1:IE+1,2:JE+1)-Ez(1:IE+1,1:JE));
        Hy(1:IE,1:JE+1)=Hy(1:IE,1:JE+1)+0.5*(Ez(2:IE+1,1:JE+1)-Ez(1:IE,1:JE+1));
    end
    mesh(Ez);
    xlabel('FDTD cells');
    ylabel('Ez');
    title(['timestep=',num2str(n)]);
    %axis([0 IE 0 JE -1 1]);
    %view(2),axis equal
    %j=1:JE+1;
    %plot(j,Ez(ic,j))
    pause(0.1);
    N_Temp=N_Temp+N_Steps;
    n_start=n+1;
end