fdtd2de[1].m

fdtd一维到三维算法

%2D finite-difference time-domain (FDTD) program
%E polarization (TMz)

close all;
clear all;

Lx = 8; %domain length in meters along x axis
Ly = 8; %domain length in meters along y axis
ds = 0.05; %spatial step in meters
Niter = 300; %# of iterations to perform
fs = 300e6; %source frequency in Hz
bw = 2; %fractional bandwidth of Gaussian pulse
t0 = 1/(fs*bw);
Nx = round(Lx/ds); %# spatial samples on x axis
Ny = round(Ly/ds); %# spatial samples on x axis
dt = ds/(300e6*sqrt(2)); %"magic time step"
eps0 = 8.854e-12; %permittivity of free space
mu0 = pi*4e-7; %permeability of free space
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);

%scale factors for E and H
ae = ones(Nx,Ny)*dt/(ds*eps0);
am = ones(Nx,Ny)*dt/(ds*mu0);
as = ones(Nx,Ny);
epsr = ones(Nx,Ny);
mur= ones(Nx,Ny);
sigma = zeros(Nx,Ny);

%here we can simulate epsilon, sigma, and mu profiles
for i=1:Nx
   for j=1:Ny
		if (sqrt((x(i)-Lx/2)^2+(y(j)-Ly/2)^2)<1) sigma(i,j) = 1e6; end         
   end
end

ae = ae./epsr;
am = am./mur;
ae = ae./(1+dt*(sigma./epsr)/(2*eps0));
as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0));

%plot the dielectric constant profile
figure(1)
map = [falsecolor(256);ones(4,3)];
%map = 1-gray(256);
colormap(map);
subplot(2,2,1);
im = transpose(epsr)-1;
image(x,y,255*im/max([1 max(max(im))]));
grid on;
axis equal;
set(gca,'YDir','normal');
title('relative permittivity');
subplot(2,2,2);
im = transpose(mur)-1;
image(x,y,255*im/max([1 max(max(im))]));
grid on;
axis equal;
set(gca,'YDir','normal');
title('relative permeabiliity');
subplot(2,2,3);
im = transpose(sigma);
image(x,y,255*im/max([1 max(max(im))]));
grid on;
axis equal;
set(gca,'YDir','normal');
title('conductivity');

%initialize fields to zero
Hx = zeros(Nx,Ny);
Hy = zeros(Nx,Ny);
Ez = zeros(Nx,Ny);

for iter=1:Niter
	%line source
	for j=3:Ny-3 Ez(3,j) = Ez(3,j)+2*sin(pi*(j-3)/(Ny-6))*sin(2*pi*fs*dt*iter); end
   
   for j=2:Ny-1 %ABC for left-propagating waves
      Hx(1,j) = Hx(2,j);
      Hy(1,j) = Hy(2,j);
   end
   for i=2:Nx-1 %ABC for down-propagating waves
      Hx(i,1) = Hx(i,2);
      Hy(i,1) = Hy(i,2);
   end
   for i=2:Nx-1
      for j=2:Ny-1
      	Hx(i,j) = Hx(i,j)-am(i,j)*(Ez(i,j+1)-Ez(i,j));
         Hy(i,j) = Hy(i,j)+am(i,j)*(Ez(i+1,j)-Ez(i,j));
      end
   end
   for j=2:Ny-1 %ABC for right-propagating waves
      Ez(Nx,j) = Ez(Nx-1,j);
   end
   for i=2:Nx-1 %ABC for up-propagating waves
      Ez(i,Ny) = Ez(i,Ny-1);
   end
   for i=2:Nx-1
      for j=2:Ny-1
      	Ez(i,j) = as(i,j)*Ez(i,j)+ae(i,j)*(Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
      end
   end
   figure(1)
   subplot(2,2,4)
   set(gcf,'doublebuffer','on');
   image(x,y,255*(transpose((abs(Ez)))));
   axis equal
   set(gca,'YDir','normal');
   grid on
   title('abs(Ez)');
   pause(0)
   %uncomment the following line to save image sequence to disk
   %imwrite(uint8(256*transpose((abs(Ez(:,Ny:-1:1))))),map,sprintf('images/f%03d.png',iter),'png');
   iter
end

%point source
%Ez(3,Ny/2) = Ez(3,Ny/2)+20*(1-exp(-((iter-1)/10)^2))*sin(2*pi*fs*dt*iter);
%Gaussian pulse
%Ez(3,Ny/2) = Ez(3,Ny/2)+20*exp(-pi*((iter*dt-t0)*bw*fs)^2)*sin(2*pi*fs*(dt*iter-t0));   
%line source
%for j=3:Ny-3 Ez(3,j) = Ez(3,j)+2*sin(pi*(j-3)/(Ny-6))*sin(2*pi*fs*dt*iter); end

%an array of wires
%if (i==round(Nx/2)) if (mod(j,5)>3) sigma(i,j)=100; end; end
%point source at middle left of page with slow turn on
%Ez(3,Ny/2) = Ez(3,Ny/2)+10*(1-exp(-((iter-1)/10)^2))*sin(2*pi*fs*dt*iter); 
%magnitude of H or B field
%contourf(x,y,transpose(377*sqrt(Hx.^2+Hy.^2)),V);
%contourf(x,y,transpose(377*mur(i,j).*sqrt(Hx.^2+Hy.^2)),V);
%cylinder
%if (sqrt((x(i)-1)^2+(y(j)-2)^2)<0.5) epsr(i,j)=9; end
%knife edge
%if ((y(j)<Ly/2)&(abs(x(i)-Lx/2)<2*ds)) sigma(i,j)=1; end
%waveguide
%if ((abs(y(j)-Ly/2)<0.5)&(x(i)>0.5)&(x(i)<6)) epsr(i,j)=2; end
%metal waveguide
%if ((j==round(1*Ny/8))|(j==round(7*Ny/8))) sigma(i,j) = 100; end
%pinhole
%if ((abs(y(j)-Ly/2)>0.25)&(abs(x(i)-Lx/2)<2*ds)) sigma(i,j)=10; end