%   TL_01 - 1D FDTD solver for a single-conductor short-short TL line
%   FDTDTL01 solves the terminated TL network (a single-conductor line) and 
%   implements the second-order central differences in time/space 
%   ECE539 - Spring 2010
clear all;
%------------------------------------------------------------------
D       = 300.0e-3;            %   Trace trace length, m
Nz      = 400; dz  = D/Nz;     %   Number of spatial cells and spatial step
zkv     = [0:dz:D];            %   Z-coordinates for the voltage nodes, Nz+1 total
zki     = [dz/2:dz:D-dz/2];    %   Z-coordinates for the current nodes, Nz total
%------------------------------------------------------------------
C    =  1.22e-010;           %   Line capacitance per unit length
L    =  3.02e-007;           %   Line inductance per unit length
R    =  100;
%------------------------------------------------------------------
cph     = 1/sqrt(L*C);                  %   phase speed on the main line
dt      = (1/cph)*dz;                   %   magic time step
T       = 7.3e-9;                        %   observation time
NT  = round(T/dt); t = [0: dt: NT*dt];  %   Number of time steps and time vector
%------------------------------------------------------------------
%   FDTD - marching on in time
% Initialization
Inext   = zeros(1, Nz+0);   %   fourth layer in the series of four  (n+3/2) 
Vnext   = zeros(1, Nz+1);   %   third layer in the series of four   (n+1)  
Ipast   = zeros(1, Nz+0);   %   second layer in the series of four  (n+1/2)  
Vpast   = zeros(1, Nz+1);   %   first layer in the series of four   (n+0)
indv    = 2:Nz;             %   batch process in space (voltage)
indi    = 1:Nz;             %   batch process in space (current)
delta   = dt/dz;            %   standard ratio
%   Initial condition (voltage)
sigma   = 0.02; Vpast   = exp(-(D/2 - zkv).^2/sigma^2); 
%   Main loop - "bootstrapping"
n = 1;
temp1 = (1 + dt*R/L/2);
temp2 = (1 - dt*R/L/2);
while n <= NT-1
    Vnext(indv) = Vpast(indv) +  delta*(-1/C)*(Ipast(indv+0) - Ipast(indv-1));     
    Inext(indi) = Ipast(indi) +  (-dt/L)*(((Vnext(indi+1) - Vnext(indi-0))/dz) + ((R)*Ipast(indi))));      
    Vpast          = Vnext;
    Ipast          = Inext;
    n              = n + 1;  
    %------------------------------------------------------------------
    string = strcat(num2str(1e9*t(n)), ' ns'); 
    x = [min(zkv) max(zkv)]; y = [0.5 0.5];
    plot(zkv, Vnext(1, :), 'r', 'LineWidth', 2); grid on;    
    line(x, y, 'LineWidth', 2); line(x, -y, 'LineWidth', 2);
    axis([min(zkv), max(zkv) -1 1]); title (strcat('Line voltage at t=',string));     
    xlabel('TL length, m'); ylabel('Voltage instant. profile'); drawnow;
    %------------------------------------------------------------------  
end