function FDTD_Connector3

% Eric J. Lundquist
% FDTD Connector Characteristic Impedance / Reflection Simulation
% 3 Regions

clear all
close all
clc

medium = 'air' ; % enter air, teflon, or sea water
cable = 3 ; % 1 is Measured 3300 Cable (predetermined Zo), 2 is RG58 Coaxial Cable (predetermined Zo), 3 is Lossless Line, 4 is Lossy Line with user-controlled Zo
src = 'gaussian' ; % enter desired voltage source type, choose from pulse, sine, rsine, spulse, gaussian
pulsetype = 1 ; % if Gaussian pulse selected, 1 is raised only, 2 is monopulse, and 3 is Gaussian-modulated sinusoidal pulse

% Zo (Characteristic Impedance) Changes
Zo1 = 50 ; % ohms
Zo2 = 55 ; % ohms
Zo3 = 50 ;
Zo = [Zo1,Zo2,Zo3];

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

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

% Statistical Randomness
srA = randn(1,maxZ)/100+1;
srB = randn(1,maxZ)/100+1;
srC = randn(1,maxZ)/100+1;
srD = randn(1,maxZ)/100+1;

% Compute Constants
lossyon = 0 ; % lossyon=0 means the line is lossless, otherwise it is lossy  
switch cable
    case 1 % Measured 3300 Cable
        R0 = 3.6996 ;
        L0 = 2.7630e-007 ;
        G0 = 4.5602e-004 ;
        C0 = 9.0846e-011 ;
        lossyon=1;
    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) ;
        lossyon=1;
    case 3 % Lossless Line
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo1^2) ;
        C0_max = L0/(max(Zo)^2);
    case 4 % Lossy Line with user-determined Zo
        R0 = 4 ;
        L0 = 1e-6 ;
        G0 = 4e-4 ;
        C0 = L0/(Zo1^2) ;
        C0_max = L0/(max(Zo)^2);
    otherwise
        warning('No cable selected! Lossless line selected by default.')
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo1^2) ; % 2500 = (50 ohms)^2
        C0_max = L0/(max(Zo)^2);
end
R(1:maxZ) = R0 ;
L(1:maxZ) = L0 ;
G(1:maxZ) = G0 ;
C(1:maxZ) = C0 ;

% Constants computed in previous section
gamma = sqrt((R+1i.*w.*L).*(G+1i.*w.*C)) ;
alpha = real(gamma) ;
beta = imag(gamma) ;
vp = w./beta ;
lambda = vp./f ;

gamma_max = sqrt((R0+1i*w*L0)*(G0+1i*w*C0_max)) ;
beta_max = imag(gamma_max) ;
vp_max = w/beta_max ;
lambda_max = vp_max/f ;
dz = lambda_max/20 ;
dt = 0.5*dz/vp_max ;
length_total = dz*maxZ

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

% Pulse Parameters
ramp = 1 ; % initialize parameter for ramped sine case
t = -dt*maxT/2:dt:dt*maxT/2;
switch pulsetype
    case 1 % raised only
        width = length(gmonopuls(t,f));
        gp = exp((-([-width:width]).^2)/(2*width));
    case 2 % monopulse
        gp = gmonopuls(t,f);
    case 3 % Gaussian-modulated sinusoidal pulse
        bw = .5 ; % bandwidth (i.e., 0.6=60%)
        gp = gauspuls(t,f,bw);
    otherwise
        warning('Gaussian Pulse Type Unspecified! Monopulse selected by default.')
        gp = gmonopuls(t,f);
end
x = length(gp);
s = round(x*0.4); % shorten Guassian pulse by 40%
g = gp(s:end);
gp = zeros(1,x);
gp(1:x-s+1)=g;

% Connector Parameters
cc = 2 ; % number of cells from center point of Region 2
connector_length = dz*(2*cc+1)
ccc = cc*5 ; % number of cells from center where reflection measurements take place

% Initialize voltage and current values of transmission line
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;
Vo = zeros(1,maxT); % original signal
Vr = zeros(1,maxT); % signal with reflections

% FDTD Loop 1 (without reflection)
for n=1:maxT
    switch src
        case 'pulse'
            if n<=25 % length of pulse
                V(1)=1;
            else
                V(1)=-1;
            end
        case 'sine'
            V(1) = sin(2*pi*f*n*dt); % sine wave source
        case 'rsine'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*f*n*dt);
            else
                V(1) = sin(2*pi*f*n*dt)*ramp;
                if ramp>0
                    ramp = ramp-.01;
                end
            end
        case 'spulse'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*f*n*dt);
            else
                V(1)=0;
            end
        case 'gaussian'
            V(1) = gp(n);
        otherwise
            error('Voltage Source Type Unspecified!')
    end
    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(round(maxZ/2)-ccc);
    I(maxZ-1) = V(maxZ-1)/Zo1 ; % absorbing boundary (not needed if simulation stops before hitting end; maxT=2*maxZ)
    movieon = 0 ;
    if movieon==1
        x = 50 ; % movie rate, ie. x=2 causes movie to play every other frame 
        if logical(x*floor(n/x)==n) == 1 % test for xth values of n (shorter movie times, plays only every xth calculation)
            figure(1)
            plot(V)
            title('Original Signal')
            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
    end
end
if movieon==1 % option available for visual validation purposes in programming
    Vo(round(maxT/2):end)=0; % clear out slight reflections caused by imperfect absorbing boundaries
    figure(3)
    plot(Vo)
    title('Pure Signal')
end

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

% Compute Constants for Connector Impendances (Region 2)
loadC = round(maxZ/2)-cc:round(maxZ/2)+cc; % connector region
L(loadC)=1e-6;
C(loadC)=(1e-6)/(Zo2^2);
if lossyon==0
    R(loadC)=0;
    G(loadC)=0;
end

% Compute Constants for Cable 2 (Region 3)
r3=round(maxZ/2)+cc+1; % starting point after connector
L(r3:maxZ) = 1e-6 ;
C(r3:maxZ) = (1e-6)/(Zo3^2) ;
if lossyon==0
    R(r3:maxZ) = 0 ;
    G(r3:maxZ) = 0 ;
end

% Constants computed in previous section
gamma = sqrt((R+j.*w.*L).*(G+j.*w.*C)) ;
alpha = real(gamma) ;
beta = imag(gamma) ;
vp = w./beta ;
lambda = vp./f ;

% Compute Constants for End Terminal
loadZ = maxZ-1;
R(loadZ)=Zo3;
R(2)=Zo1;

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

% FDTD Loop 2 (with reflections)
for n=1:maxT
    switch src
        case 'pulse'
            if n<=25 % length of pulse
                V(1)=1;
            else
                V(1)=-1;
            end
        case 'sine'
            V(1) = sin(2*pi*f*n*dt); % sine wave source
        case 'rsine'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*f*n*dt);
            else
                V(1) = sin(2*pi*f*n*dt)*ramp;
                if ramp>0
                    ramp = ramp-.01;
                end
            end
        case 'spulse'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*f*n*dt);
            else
                V(1)=0;
            end
        case 'gaussian'
            V(1) = gp(n);
        otherwise
            error('Voltage Source Type Unspecified!')
    end
    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) = V(loadZ)/Zo3 ; % right absorbing boundary
    Vr(n)=V(round(maxZ/2)-ccc);
    movieon = 1 ;
    if movieon==1
        x = 20 ; % movie rate, ie. x=2 causes movie to play every other frame 
        if logical(x*floor(n/x)==n) == 1 % test for xth values of n (shorter movie times, plays only every xth calculation)
            figure(2)
            subplot(2,1,1)
            plot(V)
            title('Signal through Varying Cables and Connectors')
            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
            subplot(2,1,2)
            bar = zeros(1,maxT);
            bar(1,1:n)=1;
            status = sprintf(['\n\nSimulation is ',num2str(round(n/maxT*100)),'%% Complete']);
            time_status = sprintf(['\n\nTime Step ',num2str(n),' of ',num2str(maxT),'\n']);
            imagesc(bar)
            colormap(winter)
            daspect([8,100/maxT,1]) % x/y aspect ratio
            axis off
            title([status,time_status])
        end
    end
end

Vrefl = Vr-Vo;
ReflCo = max(Vrefl); % does not need to be divided by max(Vo) because max(Vo)=1 for this source

figure(4)
plot(Vrefl)
title('Signal Reflections')
xlabel('time steps')
ylabel('Voltage')

if Zo2>=Zo1
    Simulated = ReflCo;
else
    Simulated = -ReflCo;
end
% Analytic = (Zo2-Zo1)/(Zo2+Zo1) % applies for 2 regions

function [yc,ys,ye] = gauspuls(t,fc,bw,bwr,tpr)

% Check input parameters:
error(nargchk(1,5,nargin));
if nargin<5, tpr = []; end
if nargin<4, bwr = []; end
if nargin<3, bw = []; end
if nargin<2, fc = []; end
if isempty(tpr), tpr = -60; end
if isempty(bwr), bwr = -6; end
if isempty(bw), bw = 0.5; end
if isempty(fc), fc = 1E3; end

if bw<=0, error('Fractional bandwidth BW must be > 0.'); end
if fc<0, error('Frequency FC must be >= 0 Hz.'); end
if bwr>=0, error('Bandwidth reference level BWR must be < 0 dB.'); end

% Determine Gaussian mean and variance in the
% frequency domain to match specifications:
r = 10.^(bwr/20);             % Ref level (fraction of max peak)
fv = -bw*bw*fc*fc/(8*log(r)); % variance is fv, mean is fc
% Determine corresponding time-domain parameters:
tv = 1/(4*pi*pi*fv);  % variance is tv, mean is 0

compute_cutoff = strncmp(lower(t),'cutoff',length(t));

if compute_cutoff,
  % Compute pulse cutoff time:
  if nargout>1, error('Too many return parameters.'); end

  % Determine extent (pulse length) of time-domain envelope:
  delta = 10.^(tpr/20);        % Ref level (fraction of max peak)
  yc = sqrt(-2*tv*log(delta)); % Pulse cutoff time

else
  % Return RF pulses:
  if nargin>4, error('Too many input arguments.'); end
  if isempty(t),
    ye=[]; yc=[]; ys=[];
    return
  end

  % Compute time-domain pulse envelope, normalized by sqrt(2*pi*tv):
  ye = exp(-t.*t/(2*tv));

  % Modulate envelope to form in-phase and quadrature components:
  yc = ye .* cos(2*pi*fc*t);    % In-phase
  if nargout>1,
    ys = ye .* sin(2*pi*fc*t);  % Quadrature
  end
end