fdtd_2_1.m

Sullivan的经典 matlab FDTD 教程 源码

%**************************************************************************
%   1D FDTD code wirh new formulation using flux density
%   two different medium (one being space and another according to epsilon)
%   Reference from the 2nd Chapter
%   "Electromagnetic Simulation Using the FDTD Method" by Dennis M. Sullivan
%   Created by Arun K Gurung, BUPT, 2005
%***************************************************************************   
clear

epsilonR = 4;   %Relative permittivity of the another medium
epsilonO = 8.8e-12; %Permittivity of the air in F/m
sigma = 1e-2;    %Conductivity of the medium  in S/m
nc = 200;	%No. of cells
NSTEPS = 400;
ddx = .01;   %Cell size
dt = ddx/6e8;   %Time step S=.5

%Intialize
ex(1:nc+1) = 0.0;
hy(1:nc+1) = 0.0;
dx(1:nc+1) = 0.0;
ix(1:nc+1) = 0.0;
ga(1:nc/2-1) = 1.;
ga(nc/2:nc+1) = 1.0/(epsilonR+(sigma*dt/epsilonO));
gb(1:nc/2-1) = 0.0;
gb(nc/2:nc+1) = sigma*dt/epsilonO;

real_pt(1:3,1:nc+1) = 0.0;  %Real part of Fourier Transform
imag_pt(1:3,1:nc+1) = 0.0;  %Imagignary part of Fourier Transform
ampn(1:3,1:nc+1) = 0.0;     %Amplitude of Fourier Transform
phasen(1:3,1:nc+1) = 0.0;   %Phase of Fourier Transform

real_in(1:3) = 0.;  %Fourier Transform of input pulse, real part
imag_in(1:3) = 0.;  %Fourier Transform of input pulse, imaginary part 
t0 = 40.0;	%Center of the incident pulse
spread = 20;	%Width of the incident pulse

%Parameters for the Fourier Transform
freq(1) = 100.e6;   %Three frequencies
freq(2) = 300.e6;
freq(3) = 500.e6;
arg(1:3) = 2*pi*freq(1:3)*dt;

%Boundary counditions
ex_low_m1 = 0;
ex_low_m2 = 0;
ex_high_m1 = 0;
ex_high_m2 = 0;

%Main FDTD Loop
for n=1:NSTEPS
    %Calculate Dx (flux density) field
    dx(2:nc) = dx(2:nc)+.5*(hy(1:nc-1)-hy(2:nc));
    
    %Put a Gaussian pulse at x=5
	pulse = exp(-.5*(((n-t0)/spread).^2.0));
	dx(5) = dx(5)+pulse;
    
    %Calculate the Ex field
	ex(2:nc) = ga(2:nc).*(dx(2:nc)-ix(2:nc));
    ix(2:nc) = ix(2:nc)+gb(2:nc).*ex(2:nc);
     
    %Calculate Fourier Transform of Ex
    for j=1:3
        real_pt(j,1:nc) = real_pt(j,1:nc)+cos(arg(j)*n).*ex(1:nc);
        imag_pt(j,1:nc) = imag_pt(j,1:nc)-sin(arg(j)*n).*ex(1:nc);
    end
    
    %Fourier Transform of the input pulse
    if ( n < 100)
        real_in(1:3) = real_in(1:3)+cos(arg(1:3)*n)*ex(10);
        imag_in(1:3) = imag_in(1:3)-sin(arg(1:3)*n)*ex(10);
    end 
    
    %Absorbing Boundary Conditions
    ex(1) = ex_low_m2;  %left side
    ex_low_m2 = ex_low_m1;
    ex_low_m1 = ex(2);
%     ex(1) = ex_low_m2;  %
%     ex_low_m2 = ex_low_m1;
%     ex_low_m1 = ex(2);
    
    %Calculate the Hy field
	hy(1:nc) = hy(1:nc)+.5*(ex(1:nc)-ex(2:nc+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 .01');
        text(100,-.7,'| medium \rightarrow');        
    end
     mov(n) = getframe;
end 
amp_in(1:3) = sqrt(real_in(1:3).^2+imag_in(1:3).^2);
phase_in(1:3) = atan2(imag_in(1:3),real_in(1:3));

for j=1:3
    ampn(j,1:nc) = (1./amp_in(j))*sqrt(real_pt(j,1:nc).^2+imag_pt(j,1:nc).^2);
    phasen(j,1:nc) = atan2(imag_pt(j,1:nc),real_pt(j,1:nc))-phase_in(j);
end

subplot(2,1,2),plot(ampn(3,:),'g');axis([0 200 0 +2]);
xlabel('FDTD cells');ylabel('Amp at 500Mhz');
text(20,1.7,'free space');text(100,1.7,'medium with permittivity 4, sigma .01');
text(100,.3,'| medium \rightarrow');