電磁數(shù)值計算作業(yè)

1、作業(yè)13：Pocklington法求偶極子天線電流分布及方向圖1011020804 章愷瑜問題描述：已知一偶極子天線長度為L，天線半徑為a，試用Pocklington法求出電流分布圖以及輸入阻抗，輻射方向圖以及功率增益，以及討論當(dāng)電長度改變時，上述的變化。圖1 偶極子天線的結(jié)構(gòu)模型計算結(jié)果：（1）當(dāng)m，且天線為半波偶極子天線，即時，電流的分布圖：圖2 半波偶極子天線電流分布圖（2）半波偶極子天線E面方向圖：圖3 半波偶極子天線E面方向圖（3）半波偶極子天線H面方向圖：圖4 半波偶極子天線H面方向圖（4）半波偶極子三維方向圖：圖5 半波偶極子天線三維方向圖（5）半波偶極子天線功率增益圖：圖6 半

2、波偶極子天線功率增益圖結(jié)果分析：由圖2可以看出，用Pocklington法計算出來的電流分布圖與之前作業(yè)中用解析法以及雙位法計算得到的相近，故計算結(jié)果正確。由圖3可以看出，用雙位法計算出來的E面方向圖與之前作業(yè)中用解析法計算得到的相近，故計算結(jié)果正確。由圖4可以看出，用雙位法計算出來的H面方向圖仍然是一個圓，這與實(shí)際相符，故結(jié)果正確。圖5所示的半波偶極子天線三維方向圖與實(shí)際相符，故結(jié)果正確。計算得到的天線輸入阻抗為80.23 + 42.395i，與理論值73.1+45.5i有一定誤差，但可以看到誤差不大，在可接受的范圍內(nèi)，故結(jié)果正確。計算得到的天線增益G=1.649，與理論值1.64接近，故結(jié)

3、果正確。圖6所示的半波偶極子天線功率增益圖與實(shí)際相符，故結(jié)果正確。程序代碼：（1）當(dāng)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(