% Model problem for the Finite Difference in Time Domain solver
% Transmission Line Problem, Telegraphic equations
% Sinusoidal generator

% Division in space and time
Nz = 100;
Nt = 50;

% Line parameters
L  = 5/3/1000/1000; % L=5/3 mH/km
C  = 1/((3e8)^2)/L; % C, light velocity, ideal line
Z0 = sqrt(L/C);     % Z0
l  = 100;           % Length

% Source data
Ug = 1;             % Thevenin generator
Rg = 500;
f  = 10e6;          % Frequency
w  = 2*pi*f;        % Angular frequency

% Load 
R2 = 1000;

% Reflexion, gamma and U2+
r  = (R2-Z0) / (R2+Z0);
b  = j * w/(3e8);
U2p= Ug / ( Rg/Z0 * (exp(b*l) - r*exp(-b*l)) + (exp(b*l) + r*exp(-b*l)) );

% Analytical formula
z  = [0:l/Nz:l];
Ua = U2p * ( exp(b*(l-z)) + r * exp(-b*(l-z)) );
Ia = U2p * ( exp(b*(l-z)) - r * exp(-b*(l-z)) ) / Z0;

u  = zeros(1,Nz+1);
i  = zeros(1,Nz);
dt = 1/f/Nt;
dz = l/Nz;

disp(['Criteria: ' num2str(3e8*dt/dz) ' <=1']);

% Simulation
for n = 1:10*Nt,
    t       = n*dt;
    ug      = 1 * cos(w*t);
    
    u(1)    = (1 - 2*dt/(Rg*C*dz)) * u(1)    - 2*dt/(C*dz) * i(1)    + 2*dt/(Rg*C*dz) * ug;
    u(Nz+1) = (1 - 2*dt/(R2*C*dz)) * u(Nz+1) + 2*dt/(C*dz) * i(Nz);    
    u(2:Nz) =                        u(2:Nz) - dt/(C*dz)   * (i(2:Nz) - i(1:Nz-1));
    
    i       = i - dt/(L*dz) * (u(2:Nz+1)-u(1:Nz));
    
    figure(1), plot(z,abs(Ua).*cos(w*t+phase(Ua)),'b',z,u,'r');
    xlabel('z [m]'); ylabel('u(z,t) [V]');
    axis([0 l -1 1]); axis square;
    legend('Analytical','FDTD');
    saveas(gcf,['pl' num2str(n)],'bmp');
    pause(0.01);
end;