%%%%矩量法解雙位方程，求直線天線表面電流分布
clear all
clc
ft=3e8;%頻率
c=3e8;%光速
ep=(10^(-9))/(36*pi);
mu=4*pi*(10^(-7));
w=2*pi*ft;
lamda=c/ft;%波長
a0=0.0025*lamda;%天線半徑
L=1.77*lamda;%天線長度
k0=2*pi/lamda;%波數
z0=120*pi;%波阻抗
N=51;
xs0=1/2;
v0=L/N;
for m=1:(N-2)
    zs1=-xs0*L+m*L/N;%m-
    zs2=-xs0*L+(m+1)*L/N;%m+
    zs12=(zs1+zs2)/2;%m
    for n=1:(N-2)
        zs3=-xs0*L+n*L/N;%n-
        zs4=-xs0*L+(n+1)*L/N;%n+
        zs34=(zs3+zs4)/2;%n
        rn=sqrt((zs34-zs12)^2+a0^2);%m,n
        rn11=sqrt((zs2-zs4)^2+a0^2);%m+,n+
        rn12=sqrt((zs1-zs4)^2+a0^2);%m-,n+
        rn21=sqrt((zs2-zs3)^2+a0^2);%m+,n-
        rn22=sqrt((zs1-zs3)^2+a0^2);%m-,n-
        %aa=m-n;
        switch m-n
            case 0 
            phi=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi);
            phi11=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi);
            phi12=exp(-j*k0*rn12)/(4*pi*rn12);
            phi21=exp(-j*k0*rn21)/(4*pi*rn21);
            phi22=phi11;
            case 1
            phi=exp(-j*k0*rn)/(4*pi*rn);
            phi11=exp(-j*k0*rn11)/(4*pi*rn11);
            phi12=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi);
            phi21=exp(-j*k0*rn21)/(4*pi*rn21);
            phi22=exp(-j*k0*rn22)/(4*pi*rn22);
            case -1
            phi=exp(-j*k0*rn)/(4*pi*rn);
            phi11=exp(-j*k0*rn11)/(4*pi*rn11);
            phi12=exp(-j*k0*rn12)/(4*pi*rn12);
            phi21=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi);
            phi22=exp(-j*k0*rn22)/(4*pi*rn22);
            otherwise
            phi=exp(-j*k0*rn)/(4*pi*rn);
            phi11=exp(-j*k0*rn11)/(4*pi*rn11);
            phi12=exp(-j*k0*rn12)/(4*pi*rn12);
            phi21=exp(-j*k0*rn21)/(4*pi*rn21);
            phi22=exp(-j*k0*rn22)/(4*pi*rn22);
         end
         z(m,n)=j*w*mu*phi*(L/N)^2+(phi11-phi12-phi21+phi22)/(j*w*ep);
     end
end
%det_z=det(z)
v=zeros(N-2,1);
I1=zeros(1,N);
v(round((N-2)*xs0),1)=v0;
I2=inv(z)*v;
for n=1:(N-2)
    I1(1,n+1)=I2(n,1);
end
%I'
figure(1)
x=L/(2*N*lamda):L/(N*lamda):(2*N-1)*L/(2*N*lamda);
subplot(2,2,1);
plot(x,abs(I1),'k');
xlabel('L/lamda');
ylabel('Magnitude');
title('電流幅度');
axis tight

subplot(2,2,2);
plot(x,180*angle(I1)/pi,'k');
xlabel('L/lamda');
ylabel('Phase');
title('電流相位');
axis tight
%%%%%%%%%%%%%%%%%%%%%%畫E面方向圖%%%%%%%%%%%%%%%%%%%%%%%%
z=(1-N)*L/(2*N):L/N:(N-1)*L/(2*N);
thita=0:0.01:2*pi;
for n=1:length(thita)
    u=0;
    for m=1:N
        u=u+I1(1,m).*sin(thita(1,n)).*exp(j*k0*z(1,m).*cos(thita(1,n))).*L/(N*lamda);
    end
    f(1,n)=u;
end
subplot(2,2,3);
polar(thita,abs(f)./max(abs(f)),'k');
title('E面方向圖')

%s(1,round(N/2))
%max(s)
Zr=v0./I1(1,round(N*xs0))