fdtd22.m

1D FDTD simulation between free space and dielectric medium with absorbing boundary conditions,calcu

%1D FDTD simulation between free space and dielectric medium with ABC
%New formulation using flux density
%The Fourier Transform has been added

%Define the parameters
KE=200;
kk=1:200;
kc=KE/2;
epsz=8.85419e-12;
ddx=0.01;
dt=ddx/6e8;
%Specify the input pulse
t0=50;
spread=10;
T=0;
%Initialize to the free space
for k=1:KE
    ga(k)=1;
    gb(k)=0;
    dx(k)=0;
    ex(k)=0;
    hy(k)=0;
    ix(k)=0;
    mag(k)=0;
    %Real,imaginary,amplitude and phase of the Fourier Transform 
    for m=1:3
        real_pt(m,k)=0;
        imag_pt(m,k)=0;
        ampn(m,k)=0;
        phasen(m,k)=0;
    end        
end
%Fourier Transform of input pulse
real_in(1:3)=0;
imag_in(1:3)=0;
amp_in(1:3)=0;
phase_in(1:3)=0;
%Initialize to absorbing boundary conditions
exlow1=0;
exlow2=0;
exhigh1=0;
exhigh2=0;
%Three frequencys for the Fourier Transform
freq(1:3)=[100e6 200e6 500e6];
arg(1:3)=2*pi*freq(1:3)*dt;
%Input parameters of this problem
kstart=input('Dielectric starts at--');
epsilon=input('Epsilon--');
sigma=input('Conductivity--');
ga(kstart:KE)=1/(epsilon+sigma*dt/epsz);
gb(kstart:KE)=sigma*dt/epsz;
nsteps=input('NSTEPS--');
tt=1:nsteps;
px(1:nsteps)=0;
%Main part of the program
for n=1:nsteps
    T=T+1;
    %Calculate the Dx field
    dx(2:KE)=dx(2:KE)+0.5*(hy(1:KE-1)-hy(2:KE));
    %Initialize with a pulse
    pulse=exp(-0.5*((t0-T)/spread).^2);
    dx(5)=dx(5)+pulse;
    px(T)=pulse;
    %Calculate Ex from Dx
    ex(2:KE)=ga(2:KE).*(dx(2:KE)-ix(2:KE));
    ix(2:KE)=ix(2:KE)+gb(2:KE).*ex(2:KE);
    %Calculate the Fourier Transform of Ex
   
     for m=1:3
        for k=2:KE
            real_pt(m,k)=real_pt(m,k)+cos(arg(m).*T).*ex(k);
            imag_pt(m,k)=imag_pt(m,k)-sin(arg(m).*T).*ex(k);
            if amp_in(m)>0
            ampn(m,k)=(1/amp_in(m)).*sqrt(real_pt(m,k).^2+imag_pt(m,k).^2);
            else
            ampn(m,k)=0;
            end
            phasen(m,k)=atan2(imag_pt(m,k),real_pt(m,k))-phase_in(m);
        end
    end
    %Calculate the Fourier Transform of the input pulse
    if T<100
        for m=1:3
            real_in(m)=real_in(m)+cos(arg(m).*T).*ex(10);
            imag_in(m)=imag_in(m)-sin(arg(m).*T).*ex(10);
            amp_in(m)=sqrt(imag_in(m).^2+real_in(m).^2);
            phase_in(m)=atan2(imag_in(m),real_in(m));
        end
    end
    %Boundary conditions
    ex(1)=exlow2;
    exlow2=exlow1;
    exlow1=ex(2);
    ex(KE)=exhigh2;
    exhigh2=exhigh1;
    exhigh1=ex(KE-1);
    %Calculate the Hy field
    hy(1:KE-1)=hy(1:KE-1)+0.5*(ex(1:KE-1)-ex(2:KE));
   
 subplot(4,1,1);
 stem(tt,px);
 xlabel('T');ylabel('source pulse');
 title('1D FDTD simulation with ABC between free space and dielectric medium ');
 subplot(4,1,2);
 plot(kk,ex);
 xlabel('k');ylabel('Ex in time domain');
 axis([0 200 -1 1]);
 subplot(4,1,3);
 plot(kk,hy);
 xlabel('k');ylabel('Hy in time domain');
 axis([0 200 -1 1]);
 subplot(4,1,4);
 plot(kk,ampn(3,:));
 xlabel('k');ylabel('Amp3-500MHz');
 axis([0 200 0 2]);
 drawnow;
end