%%?開始
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