wire_scattering.m

矩量法计算直天线电流及其方向图

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  05微波 许金根  学号：P200502040
%  作业：矩量法/滞后位积分方程计算直天线电流及其方向图/基函数：脉冲函数/点选配
%  完成时间：2006/5/14

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
clear
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 初始化
wavelength=0.1;   % 波长
k=2*pi/wavelength;   % 波数
co=3.0*10^8;   % 光速
omiga=k*co;   % 角频率
u0=4*pi*10^(-7);   % 真空磁导率   
e0=1/(36*pi*10^9);   % 真空介电常数
l=2*wavelength;   % 天线长度
a=0.5/74.2*l;   % 天线截面半径
N_point=63;   % 节点数   N_point+1为段数
N=N_point+2;   % 包含边界的总节点数
del_l=l/(N_point+1.0);   % 每段步长
% cc=1/(2*pi*del_l)*log(del_l/a)-i*k/(4*pi);   % 奇异点处fi的常数值
dd=u0/e0*k^2/(4*pi);   % 系数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 节点位置
for n=1:N_point
    point_pos(n)=-1/2*l+n*del_l;
    point_pos1(n)=point_pos(n)-1/2*del_l;
    point_pos2(n)=point_pos(n)+1/2*del_l;
end

% 源点到场点矢径长度
for m=1:N_point
    for n=1:N_point
%         R(m,n)=abs(point_pos(m)-point_pos(n));   % 第一种取法
%         R11(m,n)=abs(point_pos1(m)-point_pos1(n));
%         R12(m,n)=abs(point_pos1(m)-point_pos2(n));
%         R21(m,n)=abs(point_pos2(m)-point_pos1(n));
%         R22(m,n)=abs(point_pos2(m)-point_pos2(n));
        
        R(m,n)=sqrt(a^2+(point_pos(m)-point_pos(n))^2);   % 第二种取法
        R11(m,n)=sqrt(a^2+(point_pos1(m)-point_pos1(n))^2);
        R12(m,n)=sqrt(a^2+(point_pos1(m)-point_pos2(n))^2);
        R21(m,n)=sqrt(a^2+(point_pos2(m)-point_pos1(n))^2);
        R22(m,n)=sqrt(a^2+(point_pos2(m)-point_pos2(n))^2);
    end
end

% 阻抗矩阵计算
for m=1:N_point
    for n=1:N_point
%         if m==n
%            fi(m,n)=cc;
%            fi11(m,n)=cc;
%            fi22(m,n)=cc;
%         else
           fi(m,n)=exp(-i*k*R(m,n))/(4*pi*R(m,n));
           fi11(m,n)=exp(-i*k*R11(m,n))/(4*pi*R11(m,n));
           fi22(m,n)=exp(-i*k*R22(m,n))/(4*pi*R22(m,n));
%         end
%         if m==n+1
%            fi12(m,n)=cc;
%         else
           fi12(m,n)=exp(-i*k*R12(m,n))/(4*pi*R12(m,n));
%         end
%         if m==n-1
%            fi21(m,n)=cc;
%         else
           fi21(m,n)=exp(-i*k*R21(m,n))/(4*pi*R21(m,n));
%         end
        Z(m,n)=i*omiga*u0*del_l^2*fi(m,n)+1/(i*omiga*e0)*(fi22(m,n)-fi21(m,n)-fi12(m,n)+fi11(m,n));
    end
end

% 导纳矩阵
Y=inv(Z);

% 电压矩阵
delta=pi/2;
for n=1:N_point
        Vst(n)=del_l*exp(i*k*point_pos(n)*cos(pi-delta))*sin(delta);
end

% 电流
I=Y*Vst';

% 边界点加入
for n=N_point:(-1):1
    point_pos(n+1)=point_pos(n);
    I(n+1)=I(n);
end
point_pos(1)=-1/2*l;
point_pos(N)=1/2*l;
I(1)=0;
I(N)=0;

% E面辐射场方向图和增益方向图计算
M=361;
cta=linspace(0,2*pi,M);
for m=1:M
    for n=1:N
        Vrt(m,n)=del_l*exp(-i*k*point_pos(n)*cos(pi-cta(m)))*sin(cta(m));
    end
end
w=Vrt*I;
cgm=dd*(abs(w)).^2;   % 微分散射宽度

% for m=1:M
%     if cgm(m)<10^(-6)
%        cgm(m)=10^(-6);
%     end
% end
% rcs_double=10*log10(cgm)+60;   % 双站RCS
% RCS_double=10*log10(cgm/wavelength)+50;   % 双站RCS归一化

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 结果显示

figure(1) 
[ax,h1,h2]=plotyy(point_pos,abs(I),point_pos,180/pi*angle(I));
legend([h1,h2],'电流幅度','电流相位');
title('直导线电流幅度和相位图');

figure(2)
polar(cta,cgm');
title('双站RCS分布图');

figure(3)
polar(cta,(cgm/wavelength^2)');
title('归一化双站RCS分布图');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc