% Eric J. Lundquist
% 1D FDTD
% Frequency Domain Testing

clear all
close all
clc

set(0,'DefaultAxesFontSize',18)
set(0,'DefaultLineLineWidth',2)
set(0,'DefaultTextFontSize',18)

% Characteristic Impedances of Line and Connector
Zo1 = 50 ; % Zo of line (ohms)
Zo2 = 49 ; % Zo of connector (ohms)
cc = 0.01 ; % connector length (m)

movieon = 0 ; % enter 1 to create movie (turning on movie does slow down simulation)
i=1;

% Set Frequency Range
fmin=1e9;
fstep=1e8;
fmax=3e10;

for f=fmin:fstep:fmax

% General Parameters
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

c = 3e8 ; % speed of light
maxZ = round(fmax/c*120*cc+1) ; % cells in the transmission line
maxT = 2*maxZ ; % maximum number of time steps (2*maxZ before wave hits the end of grid)
MR = round(20:30) ; % measurement region

% 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/(Zo1^2) ; % 50 ohm line
    otherwise
        warning('No cable selected! Lossless line selected by default.')
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo1^2) ;
end

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

% Compute Constants for End Terminal
loadZ = round(maxZ-1);
G(loadZ)=Zo1;

% Constants computed in previous section
wire_cell = 10 ; % cell number 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) ;

% Connector Parameters
ccc = round((cc/dz-1)/2) ;
if ccc==0
    warning('Connector size is smaller than 1 cell! Increase minimum frequency or connector size in order to correct this problem.')
end
cZ = round(maxZ/2)-ccc:round(maxZ/2)+ccc ; % location of connector
LC = 1e-6 ;
CC = LC/(Zo2^2) ;
R(cZ)=R0;
L(cZ)=LC;
G(cZ)=G0;
C(cZ)=CC;

% 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
V = zeros(1,maxZ) ;
I = zeros(1,maxZ) ;

ramp=1; % initialize parameter for ramped sine

% 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
    I(loadZ-1) = V(loadZ-1)/Zo1 ;
    ploton = 0 ; % enter 1 to play movie
    if ploton==1
        figure(1)
        plot(V)                             % plot the voltage all along the line at time N
        ttitle = ['Voltage for f = ' num2str(f/1e6) ' MHz' ] ;
        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
maxV(n) = max(V(MR));
end

maxVV(i)=max(maxV);
i=i+1;

end

VdB = 10.*log10(maxVV);
PdB = 20.*log10(maxVV);

figure(2)
plot(fmin:fstep:fmax,VdB)
grid on
xlabel('frequency (Hz)')
ylabel('V_d_B')
title('Frequency Response')

error=0.22142; % error calculated from lossless line simulations
figure(3)
plot(fmin:fstep:fmax,VdB-error)
grid on
xlabel('frequency (Hz)')
ylabel('V_d_B - error')
title('Frequency Response')