《电磁场数值分析》（期末作业）--- 2019学年---学院：电子工程学院学号：姓名：联系方式：任课教师：2019年5月作业1模拟真空中二维TM 电磁波的传播，边界设置为一阶Mur 吸收边界，观察电磁波的传播过程。

波源为正弦函数：sin()sin(2)25z t cE t n t ωπ==∆➢ 代码：clc clearclose allxmesh =150; ymesh =150;mu0=4*pi*1.0E-7; eps0=8.85E-12;C= 3.0E8; dx=1.0; dt=0.7*dx/C; timestep=150; ez( 1:xmesh+1,1:ymesh+1 ) = 0.0; hx( 1:xmesh+1,1:ymesh ) = 0.0; hy( 1:xmesh,1:ymesh+1 ) = 0.0;coef1 = dt/( mu0 * dx ); coef2 =dt/( eps0 * dx );coef3=(C*dt-dx)/(C*dt+dx); ez1=ez;for now = 1 : timestep hx = hx - coef1 * ( ez( :, 2 : ymesh+1 ) - ez( :, 1 : ymesh ) ); hy = hy + coef1 * ( ez(2 : xmesh+1, : ) - ez(1 : xmesh, : )); ez( 2 : xmesh , 2 : ymesh ) = ez( 2 : xmesh , 2 : ymesh ) - ...coef2 * ( hx( 2 : xmesh, 2 : ymesh ) - hx( 2 : xmesh , 1 : ymesh - 1) ) + ...coef2 * ( hy( 2 : xmesh ,2 : ymesh ) - hy( 1 : xmesh - 1,2 : ymesh) );ez(1,:)=ez1(2,:)+coef3*(ez(2,:)-ez1(1,:));ez(xmesh+1,:)=ez1(xmesh,:)+coef3*(ez(xmesh,:)-ez1(xme sh+1,:));ez(:,1)=ez1(:,2)+coef3*(ez(:,2)-ez1(:,1));ez(:,ymesh+1)=ez1(:,ymesh)+coef3*(ez(:,ymesh)-ez1(:,y mesh+1));ez( xmesh/2+1, ymesh/2+1) = sin( now * dt * 2 * pi * C / 25.0 );mesh(ez);pause(0.05)ez1=ez;end➢结果与分析：第10时间步第100时间步第150时间步作业2基于Pocklington方程用MoM分析半波对称振子天线：观察天线线径和分段数目分别取不同值对天线阻抗和辐射特性的影响（半径分别取0.001λ，0.0001λ，0.00001λ，分段数取11,21,31，可列表说明）➢代码：clear all; close all; clc;% 初始化参数c=3e8; % 光速r=1 % 波长f=c/r; % 频率w=2*pi*f; % 角频率e0=8.85e-12; % 介电常数u0=4*pi*1e-7; % 磁导率a=0.00001*r; % 半径L=0.5*r; % 振子长度k=2*pi/r; % 波数N=11; % 分段数（奇数段）dl=L/(N+1); % 每段长度（分母中+1 为两头半段之和）l=L/2-dl/2; % 两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N); % 每一小段的起点坐标lzm=lz(1:N)+dl/2; % 每一小段的中点坐标lze=lz(2:N+1); % 每一小段的终点坐标%阻抗矩阵元素求解fi=log(dl/a)/(2*pi*dl)-k/(4*pi)*1i;fi_1=exp(-k*dl*1i)/(4*pi*dl);fi_2=exp(-k*2*dl*1i)/(8*pi*dl);z=ones(N,N);for m=1:Nfor n=1:Nif m==nfi1=fi;fi2=fi_1;fi3=fi_1;z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);elseif abs(m-n)==1fi1=fi_1;fi2=fi;fi3=fi_2;z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);elsefi1=exp(-k*abs(m-n)*dl*1i)/(4*pi*abs(m-n)*dl);fi2=exp(-k*abs(m+1-n)*dl*1i)/(4*pi*abs(m+1-n)*dl); fi3=exp(-k*abs(n+1-m)*dl*1i)/(4*pi*abs(n+1-m)*dl); z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);endendend%电压矩阵求解V=zeros(N,1);V((N+1)/2)=-1*(1i*w*e0);% 计算电流系数矩阵I=z\V;% 计算输入阻抗Z_in=1/I((N+1)/2);disp(['输入阻抗 = ',num2str(Z_in)]);% 计算振子上归一化电流分布I_amp=abs(I); Max=max(I_amp);Iunit2=[0;I_amp/Max(1);0]; % 两端零电流figure(1);h=0:dl/r:L/r;Ithe=sin(pi*h*r/L); % 半波振子电流解析值plot(h,Iunit2,'b',h,Ithe,'r','linewidth',2);legend('pocklinton','解析值');grid on;xlabel('电长度L/\lambda');ylabel('归一化电流');% 画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta));for n=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);theta_2=0:0.1:2*pi;Far_theory=abs((cos(k*(L/2)*cos(theta_2))-cos(k*L/2)) ./sin(theta_2));figure(2);polar(theta,Far_patten2,'-b');hold on;polar(theta_2,Far_theory,'or');hold off; legend('pocklinton','解析值');title('半波振子天线E面方向图');figure(3);polar(theta,ones(1,length(theta)),'-b');title('半波振子天线H面方向图');% 半波振子增益I_in=I((N+1)/2);A=(w*u0)^2/(4*pi*sqrt(u0/e0)*real(Z_in)*(abs(I_in))^2 );G_theta=A*abs_f.^2;Max_gain=max(G_theta)Max_gain_dB=10*log10(Max_gain);disp(['半波振子增益 = ',sprintf('%.4fdBi',Max_gain_dB)]);➢结果与分析：作业3基于电场积分方程用MoM分析对称振子天线：计算振子总长度分别为0.25λ ，0.5λ，λ，1.5λ时，振子的输入阻抗和E面方向图。

➢代码：clear all; close all; clc;% 初始化参数c=3e8; % 光速r=1; % 波长f=c/r; % 频率w=2*pi*f; % 角频率e0=8.85e-12; % 介电常数u0=4*pi*1e-7; % 磁导率a=0.001*r; % 半径L=1.5*r; % 振子长度k=2*pi/r; % 波数N=11; % 分段数（奇数段）dl=L/(N+1); % 每段长度（分母中+1 为两头半段之和）l=L/2-dl/2; % 两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N); % 每一小段的起点坐标lzm=lz(1:N)+dl/2; % 每一小段的中点坐标lze=lz(2:N+1); % 每一小段的终点坐标%阻抗矩阵元素求解fi=2*log(dl/a)/dl-k*1i;fi_1=exp(-k*dl*1i)/dl;fi_2=exp(-k*2*dl*1i)/(2*dl);for m=1:Nfor n=1:Nif m==nfi11=fi;fi12=fi;fi13=fi;fi2=fi_1;fi3=fi_1;elseif abs(m-n)==1fi11=fi_1;fi12=fi_1;fi13=fi_1;if n>mfi2=fi_2;fi3=fi;elsefi3=fi_2; fi2=fi;endelsefi11=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);fi12=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);fi13=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);if n>mfi2=exp(-k*(abs(n-m+1)*dl*1i))/((abs(n-m+1)*dl));fi3=exp(-k*(abs(n-m-1)*dl*1i))/((abs(n-m-1)*dl));elsefi2=exp(-k*(abs(m-n-1)*dl*1i))/((abs(m-n-1)*dl));fi3=exp(-k*(abs(m-n+1)*dl*1i))/((abs(m-n+1)*dl)); endendz(m,n)=1i*w*u0/(4*pi)*dl*dl*fi11+(1/(1i*4*pi*w*e0))*( fi12-fi3-fi2+fi13);endend%电压矩阵求解V=zeros(N,1);V((N+1)/2)=1;% 计算电流系数矩阵I=z\V;% 计算输入阻抗Z_in=1/I((N+1)/2);disp(['输入阻抗 = ',num2str(Z_in)]);% 画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta));for n=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i); endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);polar(theta,Far_patten2,'-b');title('振子天线E面方向图');➢结果与分析：0.25λ输入阻抗= = 12.76051-440.6278iE面方向图0.5λ输入阻抗= 81.1696+36.6291iE面方向图1λ输入阻抗=1076.2926-1067.7693i E面方向图1.5λ输入阻抗=121.0432+101.4787i E面方向图作业4请你谈谈对课程的建议、自己的收获等等。