% Eric J. Lundquist
% 1D FDTD Plots

clear all
close all
clc

% General Parameters
f = 1e9; 
w = 2*pi*f; 
eo = 8.854e-12;
u0 = 4e-7*pi;

% Copper Conductor
sigma_c =5.8e7;
u_c = u0;

% Parameters for Surrounding Medium
medium = 'air' ; % enter air, teflon, or sea water
switch medium
    case 'air'
        sigma = 0 ;
        E_r = 1 ;
        u_r = u0 ;
    case 'teflon'
        sigma = 0 ;
        E_r = 2.1 ;
        u_r = u0 ;
    case 'sea water'
        sigma = 5 ;
        E_r = 81 ;
        u_r = u0 ;
    otherwise
        warning('Surrounding medium has not been specified! Air selected by default.')
        sigma = 0 ;
        E_r = 1 ;
        u_r = u0 ;
end

maxT = 400 ; % maximum number of time steps (2*maxZ before wave hits the end of grid)
maxZ = 400 ; % cells in the transmission line

% Compute Constants

cable = 3 ;
switch cable
    case 1 % Measured 3300 Cable
        R0 = 3.6996;
        L0 = 2.7630e-007 ;
        G0 = 4.5602e-004 ;
        C0 = 9.0846e-011 ;
    case 2 % RG58 Coaxial Cable
        a = 0.445e-3 ;
        b = 1.765e-3 ;
        rs = sqrt(pi*f*u_c/sigma_c) ;
        R0 = (rs/(2*pi))*(1/a + 1/b) ;
        L0 = (u_r/(2*pi))*log(b/a) ;
        G0 = (2*pi*sigma)/log(b/a) ;
        C0 = (2*pi*E_r*eo)/log(b/a) ;
    case 3 % Lossless Line
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/2500 ; % 2500 = (50 ohms)^2
    otherwise
        warning('No cable selected! Lossless line selected by default.')
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/2500 ; % 2500 = (50 ohms)^2
end

Z0 = sqrt(L0/C0)
R(1:maxZ) = R0 ;
L(1:maxZ) = L0 ;
G(1:maxZ) = G0 ;
C(1:maxZ) = C0 ;

% Constants computed in previous section
wire_cell = 10 ; % number of cell where velocity of propagation and lambda will be calculated (should be on wire segment)
gamma = sqrt((R+j.*w.*L).*(G+j.*w.*C)) ;
alpha = real(gamma) ;
beta = imag(gamma) ;
vp = w./beta ;
lambda = vp./f ;
dz = lambda(wire_cell)/20 ;
dt = 0.5*dz/vp(wire_cell) ;

% Additional Constants
A = -1./(dz.*(R./2 + L./dt));
B = -1.*((R./2 - L./dt)./(R./2 + L./dt));
D = -1./(dz.*(G./2 + C./dt));
E = -1.*(G./2 - C./dt)./(G./2+C./dt);

% Initialize voltage and current values of transmission line
maxZ=200;
maxT=400;
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;

% FDTD loops
for n=1:maxT;
	V(1) = sin(2*pi*f*n*dt);                % sine wave source
	for k=2:maxZ;                           % find voltage everywhere on the line
		V(k)= D(k)*(I(k)-I(k-1))+E(k)*V(k);
    end
    for k=1:maxZ-1;                         % find current everywhere on the line
		I(k)= A(k)*(V(k+1)-V(k))+B(k)*I(k);
    end
    Vo(n)=V(90);
end

% Initialize voltage and current values of transmission line
maxZ=100;
maxT=400;
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;

% Compute Constants for End Terminal
ZL = Z0
loadZ = maxZ-1;
R(loadZ)=0;
L(loadZ)=0;
G(loadZ)=ZL;
C(loadZ)=0;

% FDTD loops
for n=1:maxT;
	V(1) = sin(2*pi*f*n*dt);                % sine wave source
	for k=2:maxZ;                           % find voltage everywhere on the line
		V(k)= D(k)*(I(k)-I(k-1))+E(k)*V(k);
    end
    V(loadZ+1) = V(loadZ);
    for k=1:maxZ-1;                         % find current everywhere on the line
		I(k)= A(k)*(V(k+1)-V(k))+B(k)*I(k);
    end
    I(loadZ) = V(loadZ)/ZL ;    
    Vr(n)=V(90);
% plot(V)
% ttitle = ['Voltage for n = ' int2str(maxT) ' Time Steps'] ;
% title(ttitle)
% xlabel('Distance')
% ylabel('Voltage')
% axis([0 maxZ -2 2])                         % control the axis for uniform pictures
% pause(.0001);                               % give the program time to plot to screen
end

Vrefl = Vr-Vo;
ReflCo = max(Vrefl); % does not need to be divided by max(Vo) because max(Vo)=1 for this source
tvalue = sprintf(['Gamma = '  num2str(ReflCo)]) ;

figure(2)
plot(Vrefl)
title('Reflected Voltage')
xlabel(tvalue)

if ZL>=Z0
    Simulated = ReflCo-.0108 % ~0.0125 Error
else
    Simulated = -(ReflCo-.0108) % ~0.0125 Error
end
Analytic = (ZL-Z0)/(ZL+Z0)