%%load the relevant circuit parameters
clear all
close all
hr2
m = zeros(N-1,N-1,Nvals);
Current_Source = [1,zeros(1,N-2)]';

mm = sum(y,1);

for ii = 1:N-1

      m(ii,ii,:) = mm(1,ii,:);

      for jj = 1:N-1
            if  ii ~= jj 
                 m(ii,jj,:) = -y(jj,ii,:);
            end
      end
end

Volts = zeros(N-1,Nvals);
for ii = 1:Nvals
       Volts(:,ii) = inv(m(:,:,ii))*Current_Source;
end

%Plot out the voltage at node 1. Since this is the node where 1 Ampere is injected, we are in effect calculating the input impedance...

figure(1)
%subplot(3,1,1)
%Convert to Thevenin equivalent...factor of 2. as we have 6 dB into a 50 Ohm load
factor = 2./Zo;

plot(frequency,20*log10(abs(factor.*Volts(5,:))),'r')
axis([min(frequency) max(frequency)  -40 10]);
title('transfer function')
%remove the 50 Ohm in paralalle with the current source
break
Zc = Volts(1,:);
Zin = Zo*Zc./(Zo - Zc);

subplot(3,1,2)
Q = omega*L1/rloss;
plot(frequency,real(Zin/Zo),frequency,imag(Zin/Zo),'r',frequency,Q/1000,'g')
axis([min(frequency) max(frequency)  -2 2]);
title('Input Impedance/Zo')
subplot(3,1,3)

ref = (Zin - Zo)./(Zin + Zo);

VSWR = (1 + abs(ref).^2)./(1 - abs(ref).^2);
plot(frequency,VSWR,frequency, abs(ref),'ro')
axis([min(frequency) max(frequency)  0 3]);

shg


