% Eric J. Lundquist
% 1D FDTD Plots
% 2 MHz Low-pass Composite Filter Design

clear all
close all
clc

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

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

% Set Frequency Range
fmin=2e3;
fstep=10e3;
fmax=2.5e6;

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

maxT = 300 ; % maximum number of time steps (2*maxZ before wave hits the end of grid)
maxZ = 100 ; % 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 = 1.31e-6 + 3.582e-6 ;
        G0 = 0 ;
        C0 = L0/(75^2) ; % filter is designed for a 75 ohm line
    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 ;

% Filter Design

c1=636.5e-12;
l1=6.368e-6;
C(50) = c1/(1-l1*c1*w^2);

L(51) = 3.582e-6 + 5.97e-6 ;
C(51) = 2122e-12;

c2=465.8e-12;
l2=12.94e-6;
L(52) = 5.97e-6 + 1.31e-6 ;
C(52) = c2/(1-l2*c2*w^2);

L(53) = 1.31e-6 + 3.582e-6 ;
C(53) = C(50);

% C(50) = 1e-15;

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

% 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.
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;
%ttitle = sprintf(['Voltage Signal \n\n time_T_O_T_A_L = ' num2str(maxT*dt) ' seconds']) ;

% 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)/Z0 ;
    if movieon==1
        M(n)=getframe;
    end
    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(60:90));
end

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

end

if movieon==1
    movie2avi(M,'2 MHz FDTD Filter');
end

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

figure(2)
plot(fmin:fstep:fmax,VdB)
grid on
% axis([2e3 2.04e6 -20 5])
xlabel('frequency (Hz)')
ylabel('V_d_B')
title('2 MHz Low-pass Composite Filter Frequency Response')

figure(3)
Freq = [1 1.9 2   2.01 2.02 2.03 2.04 2.045 2.05 2.1 ]*1e6;  
S12  = [0 0   -1  -2   -3   -10  -20  -34   -35  -30 ];
plot(fmin:fstep:fmax,VdB,Freq,S12./2)
grid on
% axis([2e3 2.04e6 -20 5])
xlabel('frequency (Hz)')
ylabel('V_d_B')
title('2 MHz Low-pass Composite Filter Frequency Response')
legend('FDTD','Calculated')