fdtd_1d.m

1 D FDTD matlab code

clc;
close all;
clear all;


vp=3e8/sqrt(4.4);
wavelength=vp/10e9;%0.03
dx=wavelength/20; %1/10 of wavelength  magic step 10ps dt
M=ceil(0.254/dx)
dt=0.75*dx/vp; % dx/vp =10ps, dt=1/2 dx/vp;



tmax=5e-09;
M=101;
Xmax=0.254;
L=77.1e-09/0.254;C=30.9e-12/0.254;
%L=146e-09/0.254;C=14.6e-12/0.254;
%Ro=0;
Ro=0.139/0.254;
RL=50;
Rs=50;

del_x=Xmax/(M-1);
del_t=1e-12;
N=tmax/del_t;
t0=2.5e-9;
%c=2e-19;
c=2e-19;

Vx(1:M)=0;   % 1  , 2 ...    ,101
Ix(1:M-1)=0; %   1,   2...100

%Vs(1:N)=0;

p=1;
for n=0:tmax/N:tmax-(tmax/N)
    Vs(p,1)=n/1e-09;
    
    a=-power(n-t0,2); % (n-t0).^2
    Vs(p,2)=exp(a/c); % 


    %Vs(p,2)=exp(-(n-t0).^2/c) * cos(2*pi*2e9*((n-t0))); % 

    p=p+1;
end
figure(1)
plot(Vs(:,1),Vs(:,2));
grid on;

Rt=(del_t*Ro)/(2*L);
Lt=del_t/(del_x*L);
Ct=del_t/(del_x*C);

figure(2)
for n=1:1:N    
    for m=2:1:M-1
         Vxp(m) = Vx(m) - (Ct*(Ix(m)-Ix(m-1)));
    end  
    Vxp(1) = (Vs(n,2)-(Ix(1)*Rs));
    Vxp(M)=Ix(M-1)*RL;
    
    %Vxp(1)......Vxp(M)
    
    
    for m=1:1:M-1
        %Ixp(m) = Ix(m)-Lt*(Vxp(m+1)-Vxp(m));
        Ixp(m) = (((1-Rt)*Ix(m))-(Lt*(Vxp(m+1)-Vxp(m))))/(1+Rt);
    end

    %Ixp(1)......Ixp(M-1)
    
        
    Vx=Vxp;
    Ix=Ixp;
    
    V(n,1)=Vs(n,1);% time
    V(n,2)=Vs(n,2);% source pulse
    V(n,3)=Vxp(1);
    V(n,4)=Vxp(10);
    V(n,5)=Vxp(M);
end

%plot(V(:,1),V(:,2),V(:,1),V(:,3),V(:,1),V(:,4),V(:,1),V(:,5));

hold on;
plot(V(:,1),V(:,2),'r-');
plot(V(:,1),V(:,3),'b-');
plot(V(:,1),V(:,4),'g-');
plot(V(:,1),V(:,5),'k-');

grid on;

    
    V1_fd=fft(V(:,2),2^12);
 %   V1_fd=fft(Vs,2^12);
    f = 1e12 *(0: (2^11-1))/(2^12);
    figure(3);
plot(f,V1_fd(1:2^11),'r-');