电磁数值计算作业

1、作业13：Pocklington法求偶极子天线电流分布及方向图1011020804 章恺瑜问题描述：已知一偶极子天线长度为L，天线半径为a，试用Pocklington法求出电流分布图以及输入阻抗，辐射方向图以及功率增益，以及讨论当电长度改变时，上述的变化。图1 偶极子天线的结构模型计算结果：（1）当m，且天线为半波偶极子天线，即时，电流的分布图：图2 半波偶极子天线电流分布图（2）半波偶极子天线E面方向图：图3 半波偶极子天线E面方向图（3）半波偶极子天线H面方向图：图4 半波偶极子天线H面方向图（4）半波偶极子三维方向图：图5 半波偶极子天线三维方向图（5）半波偶极子天线功率增益图：图6 半

2、波偶极子天线功率增益图结果分析：由图2可以看出，用Pocklington法计算出来的电流分布图与之前作业中用解析法以及双位法计算得到的相近，故计算结果正确。由图3可以看出，用双位法计算出来的E面方向图与之前作业中用解析法计算得到的相近，故计算结果正确。由图4可以看出，用双位法计算出来的H面方向图仍然是一个圆，这与实际相符，故结果正确。图5所示的半波偶极子天线三维方向图与实际相符，故结果正确。计算得到的天线输入阻抗为80.23 + 42.395i，与理论值73.1+45.5i有一定误差，但可以看到误差不大，在可接受的范围内，故结果正确。计算得到的天线增益G=1.649，与理论值1.64接近，故结

3、果正确。图6所示的半波偶极子天线功率增益图与实际相符，故结果正确。程序代码：（1）当m，且天线为半波偶极子天线，即时，电流的分布图：clear;clc;miu=pi*4e-7;epsilon=8.854e-12;c=3*108;L=1;a=0.0003;lambda=2;f=c/lambda;omega=2*pi*f;k=2*pi/lambda;N=51;dl=L/(N+1);zpos=dl:dl:L-dl;zpos=0 zpos L;theta=linspace(0,2*pi,100);eta=120*pi; % psi(m,n)for m=1:N+2 for n=1:N+2 if(m=n)

4、 psi(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(m-n)*dl; psi(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end; end;end;for m=2:N+1 for n=2:N+1 z(m,n)=(2-k2*dl2)*psi(m,n)-psi(m+1,n)-psi(m-1,n)/(1j*omega*epsilon); endendz=z(2:N+1,2:N+1);one=zeros(N,1);one(N+1)/2)=1;I=zone;figure(1);plot(zpos-L/2,real(0;I;

5、0)/max(real(I),*-g);hold on;zz=linspace(-L/2,L/2);i=sin(k*(L/2-abs(zz);plot(zz,i,x-b);hold on;N=101;mid=ceil(N/2);dl=L/(N+1);zpos=dl:dl:L-dl;theta=linspace(0,2*pi,100);eta=120*pi; z=zeros(N);psi11=z;psi01=z;psi10=z; for m=1:N for n=1:N if(m=n) psi11(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs

6、(m-n)*dl; psi11(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end; end;end;for m=1:N for n=1:N if n-m=1 psi10(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(m-n+1)*dl; psi10(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end endendfor m=1:N for n=1:N if m-n=1 psi01(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(n-m+1)

7、*dl; psi01(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end endendfor m=1:N for n=1:N z(m,n)=1j*omega*miu*dl*dl*psi11(m,n)+(psi11(m,n)-psi10(m,n)-psi01(m,n)+psi11(m,n)/(1j*omega*epsilon); endendV=zeros(N,1);I=zeros(N,1);for m=1:N if m=mid V(m,1)=1; else V(m,1)=0; end endI=inv(z)*V;x=-L/2:dl:L/2;plot(x,real(0;I;

8、0)/max(real(I),-r);grid on;xlabel(电长度);ylabel(归一化电流);legend(Pocklington法,解析法,双位法)（2）半波偶极子天线E面方向图：clear;clc;miu=pi*4e-7;epsilon=8.854e-12;c=3e8;L=1;a=0.0003;lambda=2;f=c/lambda;omega=2*pi*f;k=2*pi/lambda;N=51;dl=L/(N+1);zpos=dl:dl:L-dl;zpos=0 zpos L;theta=linspace(0,2*pi,100);eta=120*pi; % psi(m,n)fo

9、r m=1:N+2 for n=1:N+2 if(m=n) psi(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(m-n)*dl; psi(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end; end;end;for m=2:N+1 for n=2:N+1 z(m,n)=(2-k2*dl2)*psi(m,n)-psi(m+1,n)-psi(m-1,n)/(1j*omega*epsilon); endendz=z(2:N+1,2:N+1);one=zeros(N,1);one(N+1)/2)=1;I=zone;t=ze

10、ros(N,1);E=zeros(N,1);t=linspace(0,2*pi,N);for q=1:Nfor p=1:N E(p)=dl*sin(t(q)*exp(1j*k*p*dl*cos(t(q);endEE(q)=I*E;endhold off;for p=1:NF(p)=abs(EE(p)/max(EE);endpolar(t,F,-*r);hold onlam=1;L1=lam/2;k1=2*pi/lam;for i=1:360 g1(1,i)=abs(cos(k1*L1/2)*cos(2*pi/360*(i-1)-cos(k1*L1/2)/sin(2*pi/360*(i-1);e

11、ndfor i=1:360 f1(1,i)=g1(1,i)/max(g1); j1(1,i)=2*pi/360*(i-1); polar(j1,f1,-b); end legend(Pocklington法,解析法) title(半波偶极子E面方向图); （3）半波偶极子三维方向图：clear;clc;miu=pi*4e-7;epsilon=8.854e-12;c=3*108;L=1;a=0.0003;lambda=2;f=c/lambda;omega=2*pi*f;k=2*pi/lambda;N=51;dl=L/(N+1);zpos=dl:dl:L-dl;zpos=0 zpos L;thet

12、a=linspace(0,2*pi,100);eta=120*pi; % psi(m,n)for m=1:N+2 for n=1:N+2 if(m=n) psi(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(m-n)*dl; psi(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end; end;end;%zfor m=2:N+1 for n=2:N+1 z(m,n)=(2-k2*dl2)*psi(m,n)-psi(m+1,n)-psi(m-1,n)/(1j*omega*epsilon); endendz=z(2:N+

13、1,2:N+1);%one=zeros(N,1);one(N+1)/2)=1;I=zone;%t=zeros(N,1);E=zeros(N,1);t=linspace(0,2*pi,N);for q=1:Nfor p=1:N E(p)=dl*sin(t(q)*exp(1j*k*p*dl*cos(t(q);endEE(q)=I*E;endhold off;for p=1:NF(p)=abs(EE(p)/max(EE);endtheta=linspace(pi/2,-pi/2,N/2);psi=linspace(0,2*pi,N/2);for i=1:N/2 for j=1:N/2 X1(i,j)

14、=real(F(i)*cos(theta(i)*cos(psi(j); Y1(i,j)=real(F(i)*cos(theta(i)*sin(psi(j); Z1(i,j)=real(F(i)*sin(theta(i); endendsurf(X1,Y1,Z1);axis(-1,1,-1,1,-0.5,0.5);title(Pocklington半波振子的三维方向图)（4）半波偶极子天线功率增益图：clear;clc;miu=pi*4e-7;epsilon=8.854*10-12;c=3*108;L=1;a=0.0003;lambda=2;f=c/lambda;omega=2*pi*f;k=2

15、*pi/lambda;N=51;dl=L/(N+1);zpos=dl:dl:L-dl;zpos=0 zpos L;theta=linspace(0,2*pi,100);eta=120*pi; % psi(m,n)for m=1:N+2 for n=1:N+2 if(m=n) psi(m,n)=1/(2*pi*dl)*log(dl/a)-1j*k/4/pi; else Rmn=abs(m-n)*dl; psi(m,n)=(exp(-1j*k*Rmn)/(4*pi*Rmn); end; end;end;for m=2:N+1 for n=2:N+1 z(m,n)=(2-k2*dl2)*psi(m,n)-psi(m+1,n)-psi(m-1,n)/(1j*omega*epsilon); endendz=z(2:N+1,2:N+1);one=zeros(N,1);one(N+1)/2)=1;I=zone;Zin=1/I(N+1)/2); Ftheta=zeros(1,100);Etheta=Ftheta; for ii=1:100 temp=0; for jj=1:N temp=temp+I(jj)*dl*(cos(k*cos(theta(