DENNIS SULLIVAN

Hi.
i have a problem in FDTD 3D programme in free space with dipole antena.exaclty in T=60 i dont get the same resulte with SULLIVAN book.
In T=30 T=40 and T=50 it's ok.

can you help me please.



%%%%%%%%%%%%%%%Fundamental constants%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
c=3e8; % speed of light in free space
mu0=4.0*pi*1.0e-7; % permeability of free space
eps0=1.0/(c*c*mu0); % permittivity of free space

%%%%%%%%%%Grid parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nx=40; %number of grid cells in x-direction
Ny=40; %number of grid cells in y-direction
Nz=40; %number of grid cells in z-direction
Dx=1.0e-2; %space increment of square lattice in x direction
%Dx=5.0e-3; %space increment of square lattice in x direction
Dy=Dx; %space increment of square lattice in y direction
Dt=Dx/(2*c); %time step
itmax=60; %total number of time steps

%%%%%%%%%%%%%%%source position%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

is=Nx/2;
js=Ny/2;
ks=Nz/2;

%%%%%%%%%%%Material parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

epsr=1.0;
mur=1.0;
sig=0.0;
ro=0.0;

%%%%%%%%%%%%Updating coefficients%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

caex(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*sig)/(2.0*eps0*epsr))/(1.0+(Dt*sig)/(2.0*eps0*epsr));
cbex(1:Nx,1:Ny+1,1:Nz)=(Dt/(eps0*epsr*Dx))/(1.0+(Dt*sig)/(2.0*eps0*epsr));

caey(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*sig)/(2.0*eps0*epsr))/(1.0+(Dt*sig)/(2.0*eps0*epsr));
cbey(1:Nx,1:Ny,1:Nz)=(Dt/(eps0*epsr*Dx))/(1.0+(Dt*sig)/(2.0*eps0*epsr));

caez(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*sig)/(2.0*eps0*epsr))/(1.0+(Dt*sig)/(2.0*eps0*epsr));
cbez(1:Nx,1:Ny,1:Nz)=(Dt/(eps0*epsr*Dx))/(1.0+(Dt*sig)/(2.0*eps0*epsr));

dahx(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*ro)/(2.0*mu0*mur))/(1.0+(Dt*ro)/(2.0*mu0*mur));
dbhx(1:Nx,1:Ny,1:Nz)=(Dt/(mu0*mur*Dx))/(1.0+(Dt*ro)/(2.0*mu0*mur));

dahy(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*ro)/(2.0*mu0*mur))/(1.0+(Dt*ro)/(2.0*mu0*mur));
dbhy(1:Nx,1:Ny,1:Nz)=(Dt/(mu0*mur*Dx))/(1.0+(Dt*ro)/(2.0*mu0*mur));

dahz(1:Nx,1:Ny,1:Nz)=(1.0-(Dt*ro)/(2.0*mu0*mur))/(1.0+(Dt*ro)/(2.0*mu0*mur));
dbhz(1:Nx,1:Ny,1:Nz)=(Dt/(mu0*mur*Dx))/(1.0+(Dt*ro)/(2.0*mu0*mur));

caez(is,js,12:29)=0.0;
caez(is,js,ks)=0.0;
cbez(is,js,12:29)=0.0;
cbez(is,js,ks)=0.0;

%%%%%%%%%%%%%%%%%%%%initializ Field%%%%%%%%%%%%%%%%%%%%%%%%%%%

ex=zeros(Nx,Ny+1,Nz+1);
ey=zeros(Nx+1,Ny,Nz+1);
ez=zeros(Nx+1,Ny+1,Nz);
hx=zeros(Nx+1,Ny,Nz);
hy=zeros(Nx,Ny+1,Nz);
hz=zeros(Nx,Ny,Nz+1);

%%%%%%%%%%%%%%%%%%%%PEC conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ex(1:Nx+1,1,1:Nz+1)=0.0;
ex(1:Nx+1,Ny+1,1:Nz+1)=0.0;
ex(1:Nx+1,1:Ny+1,1)=0.0;
ex(1:Nx+1,1:Ny+1,Nz+1)=0.0;
%ez(1:Nx+1,1:Ny+1,1)=0.0;
%ez(1:Nx+1,1:Ny+1,Nz+1)=0.0;
ey(1,1:Ny+1,1:Nz+1)=0.0;
ey(Nx+1,1:Ny+1,1:Nz+1)=0.0;
ey(1:Nx+1,1:Ny+1,1)=0.0;
ey(1:Nx+1,1:Ny+1,Nz+1)=0.0;

ez(1,1:Ny+1,1:Nz+1)=0.0;
ez(Nx+1,1:Ny+1,1:Nz+1)=0.0;
ez(1:Nx+1,1,1:Nz+1)=0.0;
ez(1:Nx+1,Ny+1,1:Nz+1)=0.0;
%%%%%%%%%BEGIN TIME-STEPPING LOOP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seg=1e-10;
m=4*seg;
for t=1:itmax

%s(t)=sin(2*pi*1500*1e6*Dt*t);

%s(t)=(-4/sqrt(2*pi)).*(((t*Dt)-m)/seg).*exp(-(((t*Dt)-m).^2/(2*(seg.^2))));
t0=20;
spread=6.0;


s(t)=exp(-0.5*((t0-t).^2/((spread.^2))));

%%%%%%%%%%%%%%%%%%Update electric fields%%%%%%%%%%%%%%%%%%%%%%


ex(1:Nx,2:Ny,2:Nz)=caex(1:Nx,2:Ny,2:Nz).*ex(1:Nx,2 :Ny,2:Nz)+...
cbex(1:Nx,2:Ny,2:Nz).*(hz(1:Nx,2:Ny,2:Nz)-hz(1:Nx,1:Ny-1,2:Nz)+...
hy(1:Nx,2:Ny,1:Nz-1)-hy(1:Nx,2:Ny,2:Nz));

ey(2:Nx,1:Ny,2:Nz)=caey(2:Nx,1:Ny,2:Nz).*ey(2:Nx,1 :Ny,2:Nz)+...
cbey(2:Nx,1:Ny,2:Nz).*(hx(2:Nx,1:Ny,2:Nz)-hx(2:Nx,1:Ny,1:Nz-1)+...
hz(1:Nx-1,1:Ny,2:Nz)-hz(2:Nx,1:Ny,2:Nz));

ez(2:Nx-1,2:Ny-1,1:Nz)=caez(2:Nx-1,2:Ny-1,1:Nz).*ez(2:Nx-1,2:Ny-1,1:Nz)+...
cbez(2:Nx-1,2:Ny-1,1:Nz).*(hx(2:Nx-1,1:Ny-2,1:Nz)-hx(2:Nx-1,2:Ny-1,1:Nz)+...
hy(2:Nx-1,2:Ny-1,1:Nz)-hy(1:Nx-2,2:Ny-1,1:Nz));

ez(is,js,ks)=s(t);


hx(1:Nx,1:Ny,1:Nz)=dahx(1:Nx,1:Ny,1:Nz).*hx(1:Nx,1 :Ny,1:Nz)+...
dbhx(1:Nx,1:Ny,1:Nz).*(ey(1:Nx,1:Ny,2:Nz+1)-ey(1:Nx,1:Ny,1:Nz)+...
ez(1:Nx,1:Ny,1:Nz)-ez(1:Nx,2:Ny+1,1:Nz));

hy(1:Nx,1:Ny,1:Nz)=dahy(1:Nx,1:Ny,1:Nz).*hy(1:Nx,1 :Ny,1:Nz)+...
dbhy(1:Nx,1:Ny,1:Nz).*(ex(1:Nx,1:Ny,1:Nz)-ex(1:Nx,1:Ny,2:Nz+1)+...
ez(2:Nx+1,1:Ny,1:Nz)-ez(1:Nx,1:Ny,1:Nz));

hz(1:Nx,1:Ny,1:Nz)=dahz(1:Nx,1:Ny,1:Nz).*hz(1:Nx,1 :Ny,1:Nz)+...
dbhz(1:Nx,1:Ny,1:Nz).*(ex(1:Nx,2:Ny+1,1:Nz)-ex(1:Nx,1:Ny,1:Nz)+...
ey(1:Nx,1:Ny,1:Nz)-ey(2:Nx+1,1:Ny,1:Nz));
%%%%%%%%%%%%%%Update magnetic fields%%%%%%%%%%%%%%%%%%%%%%%%%%



%************************************************* **********************
% END TIME-STEPPING LOOP
%************************************************* **********************

end
surf(ez(1:Nx,1:Ny,ks))[img][/img][img][/img]



tanks a lot.