one_dimensional_pml.m

matlab 完全匹配层 matlab FDTD

%%%%%%%%%%%%%%%%%%%%%%% In The Name of God %%%%%%%%%%%%%%%%%%%%%%
% 
% Numerical Method in Electromagnetic
% Dr. M.S Abrishamian
% 
% Mohsen Bahrami Panah        8600374                11/12/86
% 
% Problem 11 ---------- *** One Dimensional Perfectly Matched Layer ***

    close all
    clear all
    clc
    clf
    whitebg([1 0.9 1])
    set(gcf,'Color',[1,0.9,0.9])
	IMAX = 100;         %  number of grid cells in z-direction
	NMAX = 130;         %  total number of time steps
 	PI = acos(-1.0);    %  ? = 3.1415
	C = 2.997925E8;     %  SPEED OF LIGHT
	E0 = 8.854223E-12;  %  ?0=PERMITTIVITY OF FREE SPACE
	U0 = 1.256640E-6;   %  ?0 =PERMEABLITY OF FEEE SPACE
    ETTA = sqrt(U0/E0);
    ER = 1;             %  RELATIVE PERMITTIVITY
    

% __________________________ Source Design _______________________________    
    
    TO = 40.0;          % CENTER OF GAUSSIAN PULSE
	SIGMA = 6;          % BANDWITH OF GAUSSIAN PULSE
    Loc = 50;           % Source Location
    
% __________________________ Grid parameters ______________________________

    DZ = 0.01;          %space increment
    DT = DZ/(C);        %time step size

% _____________________________ PML _______________________________________

    m = 3;             % variation
    R = 1e-10;          % reflecion coefficient
    Dpml = 10;         % PML Cells
    d = Dpml*DZ;
    SIGmax = -(m+1)*log(R)/(2*d*ETTA);
    SIG(1:Dpml) = SIGmax*(((1:Dpml)-0.5)./Dpml).^m;
    SIGm(1:Dpml) = (ETTA^2)*SIGmax*((1:Dpml)./Dpml).^m;

% _________________________________________________________________________
%     Coefficients for space region
      CA(1:IMAX) = 1;
      CB(1:IMAX) = DT/(E0*DZ);
      DA(1:IMAX) = 1;
      DB(1:IMAX) = DT/(U0*DZ);
      
      
%     Coefficients for PML region 
      T = Dpml;
      for m = 2:Dpml+1
          CA(m) = (1-SIG(T)*DT/(2*E0))/(1+SIG(T)*DT/(2*E0));
          CB(m) = (DT/(E0*DZ))/(1+SIG(T)*DT/(2*E0));
          T = T-1;
      end
      T = 1;
      for m = (IMAX-Dpml+1):IMAX
          CA(m) = (1-SIG(T)*DT/(2*E0))/(1+SIG(T)*DT/(2*E0));
          CB(m) = (DT/(E0*DZ))/(1+SIG(T)*DT/(2*E0));
          T = T+1;
      end
      
      T = Dpml;
      for m = 1:Dpml
          DA(m) = (1-SIGm(T)*DT/(2*U0))/(1+SIGm(T)*DT/(2*U0));
          DB(m) = (DT/(U0*DZ))/(1+SIGm(T)*DT/(2*U0));
          T = T-1;
      end
      T = 1;
      for m = (IMAX-Dpml+1):IMAX
          DA(m) = (1-SIGm(T)*DT/(2*U0))/(1+SIGm(T)*DT/(2*U0));
          DB(m) = (DT/(U0*DZ))/(1+SIGm(T)*DT/(2*U0));
          T = T+1;
      end
      
      
%     Illustrating the regions
      Z =  IMAX-Dpml;
      Z = [Z Z];
      X = [Dpml Dpml];
      Y = [-1.1 1.1];

% _________________________ << INITIALIZATION >> __________________________

	Ex = zeros(1,IMAX);
	Hy = zeros(1,IMAX);
    S = 100;

% ______________________ << TIME STEPING STARTS >> ________________________

	for N = 1 : NMAX
        
%                        <<   Excitation   >>
        if N<=S
            source = exp(-(N-TO)^2/(2*SIGMA^2));
            Ex(1,Loc) = source+Ex(1,Loc);
%             Hy(1,Loc) = Hy(1,Loc)+source/ETTA;
        end
        
%                      << Electric Field Update >>
        for I = 1 : IMAX-1
            Ex(1,I) = DA(1,I)*Ex(1,I) - DB(1,I)*(Hy(1,I+1) - Hy(1,I));
        end    
        
%                      << Magnetic Field Update >>
        for I = 2 : IMAX
            Hy(1,I) = CA(1,I)*Hy(1,I) - CB(1,I)*(Ex(1,I) - Ex(1,I-1));   % HY FIELD CALCULATION
        end     
        
        % getting picture at different instance
        fill([0;0;Dpml;Dpml;0],[-1.1;1.1;1.1;-1.1;-1.1],'c')
        hold on
        fill([IMAX-Dpml;IMAX-Dpml;IMAX;IMAX;IMAX-Dpml],[-1.1;1.1;1.1;-1.1;-1.1],'c')
        plot(1:IMAX,Ex,'m-',Z,Y,'b',X,Y,'b','linewidth',2);
        hold off
        title('\bf\it One Dimensional PML','color','b')
        ylabel('\bf\it Amplitude','color','b')
        axis([0,IMAX,-1.1,1.1]);
        text(.75*IMAX,1,['\bf\it t = ',num2str(N)],'color','m','fontsize',10,'BackgroundColor','c')
        text(Dpml/2,-0.6,'\bf PML','rotation',90,'color','b')
        text(IMAX-Dpml/2,-0.6,'\bf PML','rotation',90,'color','b')
        pause(0.01)
    end
% ________________________ END OF THE PROGRAM ___________________________