%   TL_03 - 1D FDTD solver for a terminated TL line with a generator/load;
%   it implements the second-order central differences in time/space including BCs.  
%   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
%------------------------------------------------------------------
cph     = 1/sqrt(L*C);                  %   phase speed on the main line
dt      = 1.00*(1/cph)*dz;              %   magic time step
T       = 15e-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

%   Main loop - "bootstrapping"
n = 1;
%------------------------------------------------------------------
%   BCs
f = 1e9; Vg  = 1*sin(2*pi*f*t);         %   Generator
RG = 50; RL = 1000;                       %   Generator and a load
B1inv   = inv(RG*C/delta + 1);  % BC inverse -left
B1dir   =    (RG*C/delta - 1);  % BC direct  -left
B2inv   = inv(RL*C/delta + 1);  % BC inverse -right
B2dir   =    (RL*C/delta - 1);  % BC direct  -right

while n <= NT-1
    VG          = [Vg(n) + Vg(n+1)];                                %   BCs
    VL          = [0];                                              %   BCs
    Vnext(1)    = B1inv*( B1dir*Vpast(1)    - 2*RG*Ipast(1)  + VG); %   BCs
    Vnext(Nz+1) = B2inv*( B2dir*Vpast(Nz+1) + 2*RL*Ipast(Nz) + VL); %   BCs    

    Vnext(indv) = Vpast(indv) +  delta*(-1/C)*(Ipast(indv+0) - Ipast(indv-1));     
    Vnext(1)    = sin(2*pi*f*n*dt);
    Inext(indi) = Ipast(indi) +  delta*(-1/L)*(Vnext(indi+1) - Vnext(indi-0));      
    Vpast          = Vnext;
    Ipast          = Inext;
    n              = n + 1;  
    %------------------------------------------------------------------
    string = strcat(num2str(1e9*t(n)), ' ns'); 
    x = [min(zkv) max(zkv)]; y = [1 1];
    plot(zkv, Vnext(1, :), 'r', 'LineWidth', 2); grid on;    
    line(x, y, 'LineWidth', 2); line(x, -y, 'LineWidth', 2);
    axis([min(zkv), max(zkv) -3 3]); title (strcat('Line voltage at t=',string));     
    xlabel('TL length, m'); ylabel('Voltage instant. profile'); drawnow;
    %------------------------------------------------------------------  
end