矩量法分析八木天线

如何用基于hallen方程的MoM分析八木天线？clearlamda=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.499而不是0.5L2=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;forkm=1:Mifkm=N1A(km).y=-d;A(km).z=lzm1(km);elseif(km=N1+N2)&(kmN1)A(km).y=0;A(km).z=lzm2(km-N1);elseA(km).y=d;A(km).z=lzm3(km-N1-N2);endendform=1:Mifm=N1flagm=1;Vim=0;hm.y=-d;hm.z=L1/2;elseifm=N1+N2flagm=2;Vim=1;hm.y=0;hm.z=L2/2;elseflagm=3;Vim=0;hm.y=d;hm.z=L3/2;endforn=1:Mifn=N1flagn=1;dL=dL1;elseifn=N1+N2flagn=2;dL=dL2;elseflagn=3;dL=dL3;endifm==nRhmn=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;elseifflagm==flagnRhmn=sqrt((hm.z-A(n).z)^2+a^2);Rmn=sqrt((A(m).z-A(n).z)^2+a^2);elseRhmn=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);endFhmn=exp(-j*k*Rhmn)/(4*pi*Rhmn)*dL;Fmn=exp(-j*k*Rmn)/(4*pi*Rmn)*dL;endZ(m,n)=cos(k*hm.z)*Fmn-cos(k*A(m).z)*Fhmn;endV(m,1)=j*Vim*sin(k*(hm.z-abs(A(m).z)))/(2*120*pi);endI=linsolve(Z,V);figure(1)plot(abs(I),'m:','linewidth',2)fedp=N1+(N2+1)/2;Zin=1/I(fedp)gridontheta=0:pi/100:2*pi;%画E面方向图pha=pi/2;form=1:length(theta)F(m)=0;forn=1:Mifn=N1dL=dL1;elseifnN1&n=N1+N2dL=dL2;elsedL=dL3;endF(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)));endendF=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;form=1:length(pha)F11(m)=0;forn=1:Mifn=N1dL=dL1;elseifnN1&n=N1+N2dL=dL2;elsedL=dL3;endF11(m)=F11(m)+I(n)*exp(j*k*(A(n).y*sin(theta)*sin(pha(m))+A(n).z*cos(theta)))*dL*(-sin(theta));endendF21=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));谢谢哦~~我已经知道了。那个如果反射线取0.500，会产生奇点的。矩量法分析八木天线clearlamda=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;forkm=1:Mifkm=N1A(km).y=-d;A(km).z=lzm1(km);elseifkm=N1+N2A(km).y=0;A(km).z=lzm2(km-N1);elseA(km).y=d;A(km).z=lzm3(km-N1-N2);endendform=1:Mforn=1:Mifn==mFmnmm=(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);elseifabs(n-m)==1Rmm=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);ifmnRmnse=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);elseFmnes=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);endelseifmnRmnes=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);elseFmnse=(1/(2*pi*dL))*log(dL/a)-j*k/(4*pi);endendelseRmm=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);endz(m,n)=j*arg*mu*dL*dL*Fmnmm+(1/(j*arg*me))*(Fmnee-Fmnes-Fmnse+Fmnss);endendV=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);gridon