mur matlab.txt

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

clear,clc
tic

IE=100;
JE=100;
ic=IE/2;
jc=JE/2;
X=1:IE;
Y=1:JE;

n0=50;
spread=12;

Ez=zeros(IE,JE);
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:end-1,2:end-1)=Ez(2:end-1,2:end-1)+0.5*(Hy(2:end-1,2:end)-...
            Hy(2:end-1,1:end-1)-Hx(2:end,2:end-1)+Hx(1:end-1,2:end-1));
    
        pulse=exp(-(n-n0)^2/spread^2);
        Ez(ic,jc)=pulse+Ez(ic,jc);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Mur ABCs   ------> 4 Edges
        % Left Boundary
        Ez(2:end-1,1)=Ez_prev(2:end-1,2)-1/3*(Ez(2:end-1,2)-Ez_prev(2:end-1,1))...
            -1/6*(Hx(2:end,1)+Hx(2:end,2)-Hx(1:end-1,1)-Hx(1:end-1,2));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Right Boundary
        Ez(2:end-1,end)=Ez_prev(2:end-1,end-1)-1/3*(Ez(2:end-1,end-1)-Ez_prev(2:end-1,end))...
            -1/6*(Hx(2:end,end-1)+Hx(2:end,end)-Hx(1:end-1,end-1)-Hx(1:end-1,end));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Top Boundary
        Ez(end,2:end-1)=Ez_prev(end-1,2:end-1)-1/3*(Ez(end-1,2:end-1)-Ez_prev(end,2:end-1))...
            -1/6*(Hy(end-1,1:end-1)+Hy(end,1:end-1)-Hy(end-1,2:end)-Hy(end,2:end));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Bottom Boundary
        Ez(1,2:end-1)=Ez_prev(2,2:end-1)-1/3*(Ez(2,2:end-1)-Ez_prev(1,2:end-1))...
            -1/6*(Hy(1,1:end-1)+Hy(2,1:end-1)-Hy(1,2:end)-Hy(2,2:end));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 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(end,1)=Ez_prev(end-1,2)+cF*(Ez(end-1,2)-Ez_prev(end,1));
        % Right_Bottom Corner
        Ez(1,end)=Ez_prev(2,end-1)+cF*(Ez(2,end-1)-Ez_prev(1,end));
        % Right_Top Corner
        Ez(end,end)=Ez_prev(end-1,end-1)...
            +cF*(Ez(end-1,end-1)-Ez_prev(end,end));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

toc