% Matlab code to evaluate the RC lowpass filter response of the first (left) 2-port problem of
% HWb, page 2, Fig 1.
% Feb12, 2019: Fixed spelling error (Tram -> Tran

clear all;close all;

NFT=1024; NF=1+NFT/2;%number of non-negative frequencies
Fmax=1e4; Fs=2*Fmax; Fmin=Fs/NFT;
f=0:Fmin:Fmax; %sweep frequency from 0 to Fmax [Hz] in steps of Fmin [Hz]

R1=9.8e3; %in [Ohms]
R2=9.8e3;
C=47e-9; %[Farads]
T=zeros(NF,2,2);

for k=1:NF; 
s=2*pi*f(k)*1j; %Define complex frequency

%Define three matrices for the three elements in the network
T1=[1 R1;0 1];
T2=[1 0; s*C 1];
T3=[1 R2;0 1];
T(k,:,:)=T1*T2*T3; %avoid all the messy algebra

%disp(sprintf('Freq=%g',f(k)))

end %end for freq loop

%plot results:
kf=2:(NF-1);

%Transmission matrix elements
A=T(kf,1,1); B=T(kf,1,2);
C=T(kf,2,1); D=T(kf,2,2);

%compute Z matrix elements
Z1=D./C; Z12=1./C;
Z2=A./C; Z21=Z12; % Assumes reciprocity

figure(1)
%fig(1); %This may be replaced with figure(1) (line above), or simply removed.

%Plot |A|(f)
subplot(4,2,1)
loglog(f(kf)/1e3,abs(A));hold on;
%xlabel('Freq [Hz]')
ylabel('|A|(f)')
title('Transmission matrix |A|(f)');

subplot(4,2,3)
semilogx(f(kf)/1e3,1+angle(A)/pi,'g');
%xlabel('Freq [Hz]')
ylabel('\angle A(f)/pi')
%title('Transmission matrix |A|(f)');

%Plot |B|(f)
subplot(4,2,2)
%Plot |B|(f)
loglog(f(kf)/1e3,abs(B));
%xlabel('Freq [Hz]')
ylabel('B')
title('Transmission matrix |B|(f)');

subplot(4,2,4)
semilogx(f(kf)/1e3,1+angle(B)/pi,'g');
%xlabel('Freq [Hz]')
ylabel('\angle A(f)/pi')
%title('Transmission matrix |A|(f)');

%Plot |C|(f)
subplot(4,2,5)
%Plot |C|(f)
loglog(f(kf)/1e3,abs(C));
%xlabel('Freq [Hz]')
ylabel('C')
title('Transmission matrix |C|(f)');

subplot(4,2,7)
semilogx(f(kf)/1e3,1+angle(C)/pi,'g');
xlabel('Freq [Hz]')
ylabel('\angle A(f)/pi')
%title('Transmission matrix |A|(f)');

%Plot |C|(f)
subplot(4,2,6)
%Plot |D|(f)
loglog(f(kf)/1e3,abs(D))
%xlabel('Freq [Hz]')
ylabel('D')
title('Transmission matrix |D|(f)');

subplot(4,2,8)
semilogx(f(kf)/1e3,1+angle(D)/pi,'g');
xlabel('Freq [Hz]')
ylabel('\angle A(f)/pi')
%title('Transmission matrix |A|(f)');

%Next step, compute Z matrix Z1, Z2, Z12
%Z1=D/C, Z2=A/C, Z12=Z21=1/C (if Delta_T=1)
