fdtd_1_5.m

Sullivan的经典 matlab FDTD 教程 源码

%**************************************************************************
%   1D FDTD code wirh two different medium (one being space and another according to epsilon)
%   Source sinusoidal, lossy medium
%   Reference from the 1st Chapter
%   "Electromagnetic Simulation Using the FDTD Method" by Dennis M. Sullivan
%   Coding by Arun K Gurung, BUPT, 2005
%***************************************************************************   
clear
KE = 200;	%KE the no. of cells to be used
NSTEPS = 1000;
%Intialize
ex(1:KE+1)=0.0;
hy(1:KE+1)=0.0;

ddx = .01;%cell size equal to 1cm
dt = ddx/(2*3e8);%time step S=.5
wave_freq = 700e6;

epsz = 8.85419e-12% permittivity of free space in Farads/meter
epsilon = 4; % dielctric constant of the another medium
sigma = .02% conductivity of the medium in Steradian/meter
ca(1:KE/2-1) = 1.;
cb(1:KE/2-1) = .5;
eaf = dt*sigma/(2*epsz*epsilon)
ca(KE/2:KE) = (1.-eaf)/(1.+eaf);
cb(KE/2:KE) = .5/(epsilon*(1.+eaf));

t0 = 40.0;	%Center of the incident pulse
spread = 12;	%Width of the incident pulse

%Left side boundary coundition initialize, no right side boundary condidtion
ex_low_m1 = 0;
ex_low_m2 = 0;

%Main FDTD Loop
for n=1:NSTEPS
    %Calculate the Ex field
	ex(2:KE) = ca(2:KE).*ex(2:KE)+cb(2:KE).*(hy(1:KE-1)-hy(2:KE));
     
    %Absorbing Boundary Conditions
    ex(1) = ex_low_m2;
    ex_low_m2 = ex_low_m1;
    ex_low_m1 = ex(2);
    
    %Put a sinusoidal source at cell 5
	pulse = sin(2*pi*wave_freq*dt*n);
	ex(5) = ex(5)+pulse;

	%Calculate the Hy field
	hy(1:KE) = hy(1:KE)+.5*(ex(1:KE)-ex(2:KE+1));
    
	if mod(n,2)==0
        ntime=num2str(n);
    
        subplot(2,1,1),plot(ex,'r'),axis([0 200 -1 +1]) ;
        title(['time ',ntime,' units']);ylabel('Ex');
        text(20,.7,'free space');text(100,.7,'medium with permittivity 4, sigma .02');
        text(100,-.7,'| medium \rightarrow');

        subplot(2,1,2),plot(hy,'b');axis([0 200 -1 +1]);
        xlabel('FDTD cells');ylabel('Hy');
        text(20,.7,'free space');text(100,.7,'medium with permittivity 4, sigma .02');
        text(100,-.7,'| medium \rightarrow');
        
        mov(n) = getframe;
    end
end