clear;clc;close all;
np = 129; %pulses
mu = 4*pi*(10^(-7));
eps = 8.854*(10^(-12));
eta = 120*pi;
acc = 99; %accuracy/number of frequency subdivisions, added acc's around
%to just easily change one number
length = .7; %antenna length
h = length/2; %height
lambda = linspace(length*pi/1,length*pi/3.5,acc); %LAMBDA DEPENDS ON L IN
%THIS EXPRESSION, AND ITS THE INVERSE OF WHAT WE HAD BEFORE
Y = zeros(acc,2);
for f = 1:acc
    Y(f,1) = pi*length/lambda(f); %to get x-axis right, do this instead
    r = 0.007022*lambda(f); %set value of changing radius
    delta3 = 2*h/(np-1);
    delta4 = 2*h/(2*(np-1));
    k = (2*pi)/lambda(f);
    Z = zeros(np,np);
    Z(1,1) = 1;
    Z(np,np) = 1;
    pos2 = delta3;
    for s = 2:(np-1)
        poss = (s-1)*delta3;
        Rr = sqrt(((pos2-poss)^2)+r^2);
        term1 = (1j*eta/k)*(exp(-1j*k*Rr)/(4*pi*((Rr)^5)));
        term2 = ((1+(1j*k*Rr))*((2*((Rr)^2))-(3*((r)^2))))+((k*r*Rr)^2);
        Z(2,s) = term1*term2*delta3;
    end
    for p = 3:(np-1)
        for g = 2:(np-1)
            if(p >= g)
                Z(p,g) = Z(2,(p-g+2));
            else
                Z(p,g) = Z(2,(g-p+2));
            end
        end
    end
    V = zeros(np,1);
    V((np+1)/2) = 1;
    I = (Z\V);%*10; %to use this, remove the ;% in between (Z\V) and *10
    Y(f,2) = I((np+1)/2);
end
x = Y(:,1);
y = Y(:,2);
plot(x,real(y));
hold on
plot(x,imag(y));
legend('Conductance','Susceptance');
ylabel('G and B in mmhos'); %added labels
xlabel('Beta*L/2'); %labels