矩量法计算Pocklington积分方程

%% 开始
function pocklington_rgtx
%% 初始化参数
clear all; close all; clc;
c=3e8;              % 光速
f=1e9;              % 频率
w=2*pi*f;           % 角频率
e0=8.85e-12;        % 介电常数
u0=4*pi*1e-7;       % 磁导率
r=c/f;              % 波长
a=0.0001*r;         % 半径
L=0.5*r;            % 振子长度
k=2*pi/r;           % 波数
N=31;               % 分段数
dl=L/（N+1）;         % 每段长度
l=L/2-dl/2;
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（NN）;
for m=1:N
    for n=1:N
        if m==n 
            fi1=fi;fi2=fi_1;fi3=fi_1;
            z（mn）=（（2-k^2*dl^2）*fi1-fi2-fi3）/（1i*w*e0）;
        elseif   abs（m-n）==1
            fi1=fi_1;fi2=fi;fi3=fi_2;
            z（mn）=（（2-k^2*dl^2）*fi1-fi2-fi3）/（1i*w*e0）;
        else
            fi1=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（mn）=（（2-k^2*dl^2）*fi1-fi2-fi3）/（1i*w*e0）;
        end  
    end
end

%% 计算电流
V=zeros（N1）;
V（（N+1）/2）=1;
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（hIunit2‘b‘hIthe‘r‘‘linewidth‘2）;
legend（‘pocklinton‘‘解析值‘）;grid on;
xlabel（‘电长度L/\lambda‘）;ylabel（‘归一化电流‘）;

%% 画方向图
theta=0:0.01:2*pi;
abs_f=zeros（1length（theta））;
for n=1:1:N
    abs_f=abs_f+I（n）*exp（k*（n-（N+1）/2）*dl*cos（theta）*1i）;
end
abs_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（（eps+cos（pi/2*cos（theta_2）））./（sin（theta_2）+eps））;
figure（2）;
polar（thetaFar_patten2‘-b‘）;
hold on;polar（theta_2Far_theory‘or‘）;hold off;
legend（‘pocklinton‘‘解析值‘）;
% title（‘半波振子天线E面方向图‘）;
figure（3）;
polar（thetaones（1length（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）/mean（G_theta）;
Max_gain_dB=10*log10（Max_gain）;
disp（[‘半波振子增益 = ‘sprintf（‘%.4fdB‘ Max_gain_dB）]）;

%% 绕杆天线方向图
theta=0:0.01:2*pi;
abs_fx=zeros（1length（theta））;
abs_fy=zeros（1length（theta））;
for n=1:1:N
    abs_fx=abs_fx+I（n）*exp（k*（n-（N+1）/2）*dl*cos（theta）*1i）;
    abs_fy=abs_fy+I（n）*exp（k*（n-（N+1）/2）*dl*cos（theta-pi/2）*1i）;
end
abs_fx=abs（sin（theta）*dl.*abs_fx）;
abs_fy=abs（sin（theta-pi/2）*dl.*abs_fy）;
F_x=abs_fx/Max_f（1）;
F_y=abs_fy/Max_f（1）;
figure（4）;
polar（thetaF_x‘-b‘）;hold on;
polar（thetaF_y‘-r‘）;
legend（‘沿x轴‘‘沿y轴‘）;

figure（5）;
FF_xy = zeros（length（0:pi/20:2*pi）1）;
nn = 0;
for iii=0:pi/20:2*pi
    F_xy = F_x*cos（iii）+F_y*sin（iii）;
    nn = nn + 1;
    FF_xy（nn） = max（abs（F_xy））;
    polar（thetaF_xy）;hold on;
end
ti