clear
lamda=1;
a=0.0025;
me=8.85e-12;
mu=4*pi*(1e-7);
arg=2*pi*(3e8)/lamda;
k=2*pi/lamda;
L1=0.499*lamda;%取0.5防止奇点出现
L2=0.460*lamda;
L3=0.419*lamda;
d=0.20*lamda;
N1=59;
N2=59;
N3=59;
dL1=L1/(N1+1);
dL2=L2/(N2+1);
dL3=L3/(N3+1);
lz1=-L1/2:dL1:L1/2;
lzm1=lz1(2:N1+1);
lz2=-L2/2:dL2:L2/2;
lzm2=lz2(2:N2+1);
lz3=-L3/2:dL3:L3/2;
lzm3=lz3(2:N3+1);
M=N1+N2+N3;
for km=1:M
if km<=N1
A(km).y=-d;
A(km).z=lzm1(km);
elseif (km<=N1+N2)&(km>N1)
A(km).y=0;
A(km).z=lzm2(km-N1);
else
A(km).y=d;
A(km).z=lzm3(km-N1-N2);
end
end
for m=1:M
if m<=N1
flagm=1;
Vim=0;
hm.y=-d;
hm.z=L1/2;
elseif m<=N1+N2
flagm=2;
Vim=1;
hm.y=0;
hm.z=L2/2;
else
flagm=3;
Vim=0;
hm.y=d;
hm.z=L3/2;
end
for n=1:M
if n<=N1
flagn=1;
dL=dL1;
elseif n<=N1+N2
flagn=2;
dL=dL2;
else
flagn=3;
dL=dL3;
end
if m==n
Rhmn=sqrt((hm.z-A(n).z)^2+a^2);
Fmn=(1/(2*pi))*log(dL/a)-j*k*dL/(4*pi);
Fhmn=exp(-j*k*Rhmn)/(4*pi*Rhmn)*dL;
else
if flagm==flagn
Rhmn=sqrt((hm.z-A(n).z)^2+a^2);
Rmn=sqrt((A(m).z-A(n).z)^2+a^2);
else
Rhmn=sqrt((hm.y-A(n).y)^2+(hm.z-A(n).z)^2);
Rmn=sqrt((A(m).y-A(n).y)^2+(A(m).z-A(n).z)^2);
end
Fhmn=exp(-j*k*Rhmn)/(4*pi*Rhmn)*dL;
Fmn=exp(-j*k*Rmn)/(4*pi*Rmn)*dL;
end
Z(m,n)=cos(k*hm.z)*Fmn-cos(k*A(m).z)*Fhmn;
end
V(m,1)=j*Vim*sin(k*(hm.z-abs(A(m).z)))/(2*120*pi);
end
I=linsolve(Z,V);
figure(1)
plot(abs(I),'m:','linewidth',2)
fedp=N1+(N2+1)/2;
Zin=1/I(fedp)
grid on
theta=0:pi/100:2*pi;%画E面方向图
pha=pi/2;
for m=1:length(theta)
F(m)=0;
for n=1:M
if n<=N1
dL=dL1;
elseif n>N1&n<=N1+N2
dL=dL2;
else
dL=dL3;
end
F(m)=F(m)+I(n)*exp(j*k*(A(n).y*sin(theta(m))*sin(pha)+A(n).z*cos(theta(m))))*dL*(-sin(theta(m)));
end
end
F=abs(F);
F1=F/max(F);
figure(2)
polar(theta,F1,'mo')
FBR=20*log10(F(51)/F(151));
theta=pi/2;%画H面方向图
pha=0:pi/100:2*pi;
for m=1:length(pha)
F11(m)=0;
for n=1:M
if n<=N1
dL=dL1;
elseif n>N1&n<=N1+N2
dL=dL2;
else
dL=dL3;
end
F11(m)=F11(m)+I(n)*exp(j*k*(A(n).y*sin(theta)*sin(pha(m))+A(n).z*cos(theta)))*dL*(-sin(theta));
end
end
F21=abs(F11);%方向图函数
F21=F21/max(F21);%归一化方向图函数
figure(3)
polar(pha,F21,'mo');
Rin=real(Zin);
G=(1/(480*pi^2))*arg^2*mu^2*max(F)^2/(Rin*abs(I(fedp))^2);
T=(Zin-50)/(Zin+50);
SWR=(1+abs(T))/(1-abs(T));
clear
lamda=1;%波长
a=0.0025;%振子的半径
me=8.85e-12;%介电常数
mu=4*pi*(1e-7);%磁导率
arg=2*pi*(3e8)/lamda;%角频率
k=2*pi/lamda;%波数
N2=51;
L1=0.5*lamda;
L2=0.460*lamda;
L3=0.419*lamda;
d=0.2*lamda;
dL=L2/(N2+1);
N1=fix(L1/dL);
dL1=rem(L1,dL);
N3=fix(L3/dL);
dL3=rem(L3,dL);
l1=L1/2-dL1/2;
l2=L2/2-dL/2;
l3=L3/2-dL3/2;
lzintal1=-l1:dL:l1;
lzm1=lzintal1(1:N1)+dL/2;
lzintal2=-l2:dL:l2;
lzm2=lzintal2(1:N2)+dL/2;
lzintal3=-l3:dL:l3;
lzm3=lzintal3(1:N3)+dL/2;
M=N1+N2+N3;
for km=1:M
if km<=N1
A(km).y=-d;
A(km).z=lzm1(km);
elseif km<=N1+N2
A(km).y=0;
A(km).z=lzm2(km-N1);
else
A(km).y=d;
A(km).z=lzm3(km-N1-N2);
end
end
for m=1:M
for n=1:M
if n==m
Fmnmm=(1/(2*pi*dL))*log(dL/a)- j*k/(4*pi);
Fmnee=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);
Fmnss=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);
Fmnse=exp(-j*k*dL)/(4*pi*dL);
Fmnes=exp(-j*k*dL)/(4*pi*dL);
elseif abs(n-m)==1
Rmm=sqrt((A(m).y-A(n).y)^2+(A(m).z-A(n).z)^2);
Fmnmm=exp(-j*k*Rmm)/(4*pi*Rmm);
Fmnee=exp(-j*k*Rmm)/(4*pi*Rmm);
Fmnss=exp(-j*k*Rmm)/(4*pi*Rmm);
if m<n
Rmnse=sqrt((A(m).y-A(n).y)^2+((A(m).z-dL/2)-(A(n).z+dL/2))^2);
Fmnse=exp(-j*k*Rmnse)/(4*pi*Rmnse);
if (m==N1)||(m==N1+N2)
Rmnes=sqrt((A(m).y-A(n).y)^2+((A(m).z+dL/2)-(A(n).z-dL/2))^2);
Fmnes=exp(-j*k*Rmnes)/(4*pi*Rmnes);
else
Fmnes=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);
end
elseif m>n
Rmnes=sqrt((A(m).y-A(n).y)^2+((A(m).z+dL/2)-(A(n).z-dL/2))^2);
Fmnes=exp(-j*k*Rmnes)/(4*pi*Rmnes);
if (n==N1)||(n==N1+N2)
Rmnse=sqrt((A(m).y-A(n).y)^2+((A(m).z-dL/2)-(A(n).z+dL/2))^2);
Fmnse=exp(-j*k*Rmnse)/(4*pi*Rmnse);
else
Fmnse=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);
end
end
else
Rmm=sqrt((A(m).y-A(n).y)^2+(A(m).z-A(n).z)^2);
Rmnse=sqrt((A(m).y-A(n).y)^2+((A(m).z-dL/2)-(A(n).z+dL/2))^2);
Rmnes=sqrt((A(m).y-A(n).y)^2+((A(m).z+dL/2)-(A(n).z-dL/2))^2);
Fmnmm=exp(-j*k*Rmm)/(4*pi*Rmm);
Fmnee=exp(-j*k*Rmm)/(4*pi*Rmm);
Fmnss=exp(-j*k*Rmm)/(4*pi*Rmm);
Fmnse=exp(-j*k*Rmnse)/(4*pi*Rmnse);
Fmnes=exp(-j*k*Rmnes)/(4*pi*Rmnes);
end
z(m,n)=j*arg*mu*dL*dL*Fmnmm+(1/(j*arg*me))*(Fmnee-Fmnes-Fmnse+Fmnss);
end
end
V=zeros(M,1);
fedp=N1+(N2+1)/2;%馈电点的位置
V(fedp)=1*dL;%馈电的电位
I=linsolve(z,V);
Zin=V(fedp)/I(fedp)%输入阻抗
I1=abs(I);
plot(I1);grid on