%**** ABCD_method.m simulates the TDR response of cascaded transmission lines
%using ABCD method.  It also compares the simulated result with the
%measured one. ****

clear all; close all;
%---------------  User Input Area  ---------------
wire_type={'RG58', 'RG59BU', 'RG58', 'RG62', 'RG58'}; %wire type
wire_length=[17.2, 2.77, 0.33, 2.2*(66/84), 3.68];   %length in meters (VOP normalize to RG58)
Load_R=99999;                       %Real Part (ohms)
Load_I=10^-15;                      %Imaginary part
%------------  End of User Input Area  -------------

n=length(wire_length);      %number of sections
c=3*10^8;                   %speed of light
ns=10^-9;                   %nano second;
N=5000;                     %total freq steps
samp_freq=10^9;             %sampling freq
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)+10^-6;         %freq distribution, start 10^-6 to avoid DC
y=[zeros(1,0.5*N-1), .5, ones(1,0.5*N)]; %a step function with rise time of 2ns

Y=fft(y);                   %Step function in freq domain
w=2*pi*f;                   %angular frequency
ZL=Load_R+i.*w*Load_I;      %Load impedance

for s=1:n                                           %determine wire type, find Z & Gamma, Alpha & Beta
    wire=char(wire_type(s));
    switch upper(wire)                      
        case 'RG58'
            [Z(s,:),Gamma(s,:)]=RG58(f);
        case 'RG59BU'
            [Z(s,:),Gamma(s,:)]=RG59BU(f);
        case 'RG62'
            [Z(s,:),Gamma(s,:)]=RG62(f);
    end
end
L=wire_length;                                       %Array for wire lengths

Zs=50*ones(1,N);                                 %if fixed Zs, assume 50 ohms for now

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,:,:)=[cosh(Gamma(s,x).*L(s)), Z(s,x)*sinh(Gamma(s,x).*L(s)); sinh(Gamma(s,x).*L(s))/Z(s,x) cosh(Gamma(s,x).*L(s))];
        MM=reshape(M(s,:,:),2,2);           %transform 3D matrix into 2D
        temp=temp*MM;                       %M1*M2...Mn
    end
    ABCD1(:,:,x)=temp*ML;                   %Reverse Transfer Function
    ABCD2(:,:,x)=M0*ABCD1(:,:,x);           %Forward Transfer Function (TDT)
    H(x)=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_response=f2-f1;                 %Raw TDR Response
TDR_offset=TDR_response(1);         %offset of the result
TDR=TDR_response-TDR_offset;        %final result
t_index=(0:length(TDR)-1)/10;       %x-axis index

%Load the measured TDR data
Tdata=load('TDR3.DAT');
avg=Tdata(1);
vop=Tdata(2);
points=Tdata(3);
start=Tdata(4);
test_length=Tdata(5);
probe_length=Tdata(6);
offset=Tdata(7)+0.36;                   %TDR100 has 0.36m of offset
mTDR=(Tdata(8:end))';
dx=(test_length+offset)/(points-1);
m_index=(-offset:dx:test_length);     

%plot the figure
plot(m_index,mTDR,'r',t_index,TDR,'b');
legend('measured','simulated', 'Location','Best');
axis([0 100 -1.5 1.5]);
xlabel('Length(m)');
ylabel('Magnitude');
title('L=[17.2, 2.77, 0.33, 2.2, 3.68], Z=[50,75,50,93,50], Open-ended');