function ABCD
% ABCD and Reflection Coefficient Evaluation

clc
clear all
close all

% User Input
wire_type={39.74,39.77,40.62,46.61,83.51,46.61,40.62,39.77,39.74};
% Input wire type such as 'RG58','RG59BU','RG62' or a number such as 50 to
% specify the desired Zo value which will be used to calculate RLGC values
% No damage to shield, Zo = 39.7433
% 5% shield damage, Zo = 39.7659
% 25% shield damage, Zo = 40.6166
% 45% shield damage, Zo = 46.6145
% 50% shield damage, Zo = 83.5104
wire_length=[10000,1,1,1,1,1,1,1,10000]*1e-3; % Length of each section, in meters 
Load_R=10^6; % Resistive Load (ohms)
Load_L=0; % Inductive Load (H)
Load_C=0; % Capacitive Load (F)
fixed_Zs=0; % 0: Zs=0.9*Z1 , 1: Zs=50 ohms

% Parameters
n=length(wire_length);  % number of sections
c=3e8;                  % speed of light
ns=1e-9;                % define nanosecond
b=500;             % scaling factor for freq steps
N=sum(wire_length)*b;   % total freq steps, scaled by b to account for VoP/c
                        % Explanation: The total length of the wires needs
                        % to be < 1/2 of simulated length, which ends up
                        % being scaled by VoP/c. Thus N, the # of frequency
                        % steps, needs to be > 20*2*sum(wire_length)/(VoP/c)
N=round(N/2)*2;         % make N a round and even number
samp_freq=1e9;          % sampling frequency
dt=1/samp_freq;         % sampling period
t=dt*(0:N-1);           % create t array for time-domain signal
df=1/(N*dt);            % freq steps
f=df*(0:N-1)+1e-6;      % freq distribution, start 10^-6 to avoid DC
w=2*pi*f;               % angular frequency
ZL=Load_R+i*w*Load_L+(1./(i*w*Load_C+10^-15)); % complex load, 10^-15 to prevent divide by zero

% Generate Pulse
pulsetype = 'Gaussian';
switch pulsetype
    case 'step'
        y=[zeros(1,0.5*N-1), .5, ones(1,0.5*N)]; % Step function with rise time of 2ns
    case 'SSTDR'
        dx = 40/N;
        X = 2;
        x = -X+dx:dx:X;
        s = sinc(x); % sinc function
        t = tripuls(x,2*X); % triangular pulse
        SSTDR = s.*t;
        z = round((N-length(SSTDR))/2);
        y=[zeros(1,z) , SSTDR , zeros(1,z)];
    case 'Gaussian'
        x = round(N/2);
        y = exp((-([-x+1:x]).^2)/(2*x));
    otherwise
        error('Pulse Type Unspecified!')
end
figure
plot(t,y)
title('Pulse Shape')
grid on
xlabel('time (s)')
ylabel('voltage (V)')

Y=fft(y); % Step function in freq domain

% Determine wire type, find Z & Gamma, Alpha & Beta
Zo=[];
for s=1:n
    if iscellstr(wire_type(s))==1 % determine whether wire type or Zo is specified
        wire=char(wire_type(s));
        wire=upper(wire);
    else
        wire='Zo_Specific'; % Zo has been specified
        Zo=cell2mat(wire_type(s));
    end
    switch wire
        % Obtain Zo and Gamma values for each wire section (s) as function of freq (:)
        case 'RG58'
            [Z(s,:),Gamma(s,:)]=RG58(f);
        case 'RG59BU'
            [Z(s,:),Gamma(s,:)]=RG59BU(f);
        case 'RG62'
            [Z(s,:),Gamma(s,:)]=RG62(f);
        case 'Zo_Specific'
            [Z(s,:),Gamma(s,:)]=Zo_Specific(f,Zo);
        otherwise
            error('Wire Type Specification Error!')
    end
    alpha(s,:)=real(Gamma(s,:)); % Attenuation constant
    beta(s,:)=imag(Gamma(s,:));  % Phase constant
end
L=wire_length; % Array for wire lengths

%Source Impedance
if fixed_Zs == 0
    Zs=0.9*Z(1,:); % Assume source impedance = RG58*0.9
else
    Zs=50*ones(1,N); % If fixed Zs, assume 50 ohms for now
end

% ABCD Parameters
for x=1:N
    M0=[1 Zs(x); 0 1];                      % Source ABCD parameter
    ML=[1 0; 1/ZL(x) 1];                    % Load ABCD parameter
    temp=[1 0; 0 1];                        % dummy variable - unit matrix
    temp2=1;                                % dummy variable
    for s=1:n
        M(s,:,:)=[cos(beta(s,x).*L(s)) j*Z(s,x)*sin(beta(s,x).*L(s)); j*sin(beta(s,x).*L(s))/Z(s,x) cos(beta(s,x).*L(s))]; % Wire_n ABCD parameter
        MM=reshape(M(s,:,:),2,2);           % Transform 3D matrix into 2D
        temp=temp*MM;                       % M1*M2...Mn
        temp2=temp2*exp(-alpha(s,x));       % Apply the attenuation constant
    end
    ABCD1(:,:,x)=temp*ML;                   % Reverse Transfer Function
    ABCD2(:,:,x)=M0*ABCD1(:,:,x);           % Forward Transfer Function (TDT)
    attn=temp2^2;                           % Attenuate twice for round trip
    H(x)=attn*ABCD1(1,1,x)/(ABCD2(1,1,x));  % TDR Transfer Funcion
end

YY=Y.*H;                            % TDR response in freq domain
yy=ifft(YY,'symmetric');            % TDR response in time domain
f1=yy(1:round(N/2));                % first half
f2=yy(round(N/2)+1:end);            % second half
TDR=f2-f1;                          % TDR Response
t_index=(0:length(TDR)-1)/10;       % time index     

p=1; % place holder
distance=ones(1,length(t_index)); % create distance vector
for s=1:n % scale each section
    Beta=imag(Gamma(s,:));
    VoP=w./Beta; % Velocity of Propagation
    VoP=mean(VoP);
    sf=VoP/(c*2/3); % scale by 2/3 the speed of light
    d=round(40*wire_length(s));
    if s==1 % first section of wire
        distance(p:p+d) = [1:d+1]*sf;
    elseif s==n % last section of wire
        distance(p:end) = [1:length(distance)-p+1]*sf+distance(p-1);
    else % middle wire sections
        distance(p:p+d) = [1:d+1]*sf+distance(p-1);
    end
    p=p+d+1;
end
distance=distance/10; % scale by 10 to correspond with time steps
if distance(end) < sum(wire_length)
    error('The total simulated distance is not long enough to account for the total length of wires. You can account for this by increasing Zo (which in turn scales VoP accordingly), increasing the number of frequency steps, or reducing the total length.')
end

figure
plot(distance,TDR,'b')
xlabel('Length (m)')
ylabel('Reflection Coefficient  (| \Gamma|)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Z,Gamma]=Zo_Specific(f,Zo)
 
cable = 1 ; % 1 is Lossless Line, 2 is Lossy Line, with user-specified Zo
switch cable
    case 1 % Lossless Line
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo^2) ;
    case 2 % Lossy Line with user-determined Zo
        R0 = 2e-2 ;
        L0 = 1e-6 ;
        G0 = 2e-4 ;
        C0 = L0/(Zo^2) ;
    otherwise
        warning('No cable selected! Lossless line selected by default.')
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo^2) ; % i.e. 2500 = (50 ohms)^2
end
w=2*pi*f; % angular frequency
Z=sqrt((R0+j.*w.*L0)./(G0+j.*w.*C0)); % characteristic impedance
Gamma=sqrt((R0+j.*w.*L0).*(G0+j.*w.*C0)); % propagation constant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RG58, RG59BU, & RG62 WIRE TYPE FUNCTIONS from Shang Chieh Wu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function evaluates the characteristic impedance (Z) and the
% propagation constant (Gamma) of RG58 coaxial cable.  The function has a
% format of [Z, Gamma]=RG58(f), where the frequency is in Hertz (Hz)
function [Z, Gamma]=RG58(f)

Er1=2.3;                        %dielectric constant, vop=0.6594
E0=8.854*10^-12;                %vacuum dielectric constant
mu=4*pi*10^-7;                  %permeability
a=0.032*2.54/200;               %conductor radius (m)
b=0.116*2.54/200;               %insulation radius (m)
sigma_c=5.813*10^7;             %conductor's conductivity
Rs=1.988*sqrt((f/10^6)/sigma_c);    %intrinsic resistance
R=Rs*(1/a+1/b);                %per unit reistance
L=mu*log(b/a)/(2*pi);          %per unit inductance
G=10^-12;                      %per unit conductance, very small number, but not zero
C=2*pi*Er1*E0/log(b/a);        %per unit capacitance
w=2*pi*f;
Z=sqrt((R+j*w*L)./(G+j*w*C));    %characteristic impedance of RG58
Gamma=sqrt((R+j*w*L).*(G+j*w*C)); %propagation constant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function evaluates the characteristic impedance (Z) and the
% propagation constant (Gamma) of RG59B/U coaxial cable.  The function has a
% format of [Z, Gamma]=RG59(f), where the frequency is in Hertz (Hz)
function [Z, Gamma]=RG59BU(f)
Er1=2.3;                            %dielectric constant, vop=0.6594
E0=8.854*10^-12;                    %vacuum dielectric constant
mu=4*pi*10^-7;                      %permeability
a=0.57404/2000;                     %conductor radius (m)
b=3.71/2000;                        %insulation radius (m)
sigma_c=5.813*10^7;                 %conductor's conductivity
Rs=1.988*sqrt((f/10^6)/sigma_c);    %intrinsic resistance
R=Rs*(1/a+1/b);                     %per unit reistance
L=mu*log(b/a)/(2*pi);               %per unit inductance
G=10^-12;                           %per unit conductance, very small number, but not zero
C=2*pi*Er1*E0/log(b/a);             %per unit capacitance
w=2*pi*f;                           %angular frequency
Z=sqrt((R+j*w*L)./(G+j*w*C));       %characteristic impedance of RG59
Gamma=sqrt((R+j*w*L).*(G+j*w*C));   %propagation constant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function evaluates the characteristic impedance (Z) and the
% propagation constant (Gamma) of RG62 coaxial cable.  The function has a
% format of [Z, Gamma]=RG62(f), where the frequency is in Hertz (Hz)
function [Z, Gamma]=RG62(f)
Er1=1.4173;                         %dielectric constant, vop=0.6594
E0=8.854*10^-12;                    %vacuum dielectric constant
mu=4*pi*10^-7;                      %permeability
a=0.64516/2000;                     %conductor radius (m)
b=3.66/2000;                        %insulation radius (m)
sigma_c=5.813*10^7;                 %conductor's conductivity
Rs=1.988*sqrt((f/10^6)/sigma_c);    %intrinsic resistance
R=Rs*(1/a+1/b);                     %per unit reistance
L=mu*log(b/a)/(2*pi);               %per unit inductance
G=10^-12;                           %per unit conductance, very small number, but not zero
C=2*pi*Er1*E0/log(b/a);             %per unit capacitance
w=2*pi*f;                           %angular frequency
Z=sqrt((R+j*w*L)./(G+j*w*C));       %characteristic impedance of RG58
Gamma=sqrt((R+j*w*L).*(G+j*w*C));   %propagation constant