clc
clear all
f=3e9;   
N1=4;N2=8;N3=12;
a=pi/2; %饋電相位差
i=1;  %天線電流值
lambda=(3e8)/f; %lambda=c/f     波長
d=lambda/2;
beta=2.*pi/lambda;
W=-2*pi:0.001:2*pi;
 y1=sin((N1.*W./2))./(N1.*(sin(W./2)));  %歸一化陣因子
 y1=abs(y1);
 r1=max(y1);
 y2=sin((N2.*W./2))./(N2.*(sin(W./2)));  %歸一化陣因子
 y2=abs(y2);
 r2=max(y2);
 y3=sin((N3.*W./2))./(N3.*(sin(W./2)));  %歸一化陣因子
 y3=abs(y3);
 r3=max(y3);
 %歸一化陣因子繪圖程序,
 figure(1)
 subplot(311);plot(W,y1) ;  grid on;     %繪出N=4等幅等矩陣列的歸一化陣因子
  xlabel('f=3GHz,N=4,d=1/2波長,a=π/2')
 subplot(312);plot(W,y2) ;  grid on;     %繪出N=8等幅等矩陣列的歸一化陣因子
  xlabel('f=3GHz,N=8,d=1/2波長,a=π/2')
 subplot(313);plot(W,y3)  ;  grid on;    %繪出N=12等幅等矩陣列的歸一化陣因子
  xlabel('f=3GHz,N=12,d=1/2波長,a=π/2')
%---------------------
%只有參數N改變的天線方向圖
t=0:0.01:2*pi; 
W=a+(beta.*d.*cos(t));
 z1=(N1/2).*(W);
z2=(1/2).*(W);
W1=sin(z1)./(N1.*sin(z2));  %非歸一化的陣因子K1
K1=abs(W1);
%----------------------
W=a+(beta.*d.*cos(t));
 z3=(N2/2).*(W);
z4=(1/2).*(W);
W2=sin(z3)./(N2.*sin(z4));  %非歸一化的陣因子K2
K2=abs(W2);
%-------------------------
W=a+(beta.*d.*cos(t));
 z5=(N3/2).*(W);
z6=(1/2).*(W);
W3=sin(z5)./(N3.*sin(z6));  %非歸一化的陣因子K3
K3=abs(W3);
%--------------------繪圖函數
figure(2)
subplot(131);polar(t,K1);xlabel('f=3GHz,N=4,d=1/2波長,a=π/2');
subplot(132);polar(t,K2);xlabel('f=3GHz,N=8,d=1/2波長,a=π/2');
subplot(133);polar(t,K3);xlabel('f=3GHz,N=12,d=1/2波長,a=π/2');
%----------------------
%只有陣列單元方向角a改變的天線方向圖
N4=10;d1=lambda/4;a1=0;a2=pi/2;a3=pi/2+pi/10;
A=a1+(beta.*d1.*cos(t));
 x1=(N4/2).*(A);
x2=(1/2).*(A);
A1=sin(x1)./(N4.*sin(x2));  %非歸一化的陣因子K4
K4=abs(A1);
%---------------------------
B=a2+(beta.*d1.*cos(t));
 y_1=(N4/2).*(B);
y_2=(1/2).*(B);
B1=sin(y_1)./(N4.*sin(y_2));  %非歸一化的陣因子K5
K5=abs(B1);
%----------------------------
C=a3+(beta.*d1.*cos(t));
 v1=(N4/2).*(C);
v2=(1/2).*(C);
C1=sin(v1)./(N4.*sin(v2));  %非歸一化的陣因子K6
K6=abs(C1);
%--------------------------繪圖函數
figure(3)
subplot(131);polar(t,K4);xlabel('f=3GHz,N=10,d=1/4波長,a=0');
subplot(132);polar(t,K5);xlabel('f=3GHz,N=10,d=1/4波長,a=π/2');
subplot(133);polar(t,K6);xlabel('f=3GHz,N=10,d=1/4波長,a=π/2+π/10');
%----------------------------------------------------------------------
%只有陣列單元間隔d改變的天線方向圖
N5=20;d2=lambda/4;d3=lambda/2;d4=0.7*lambda;a4=pi/2;
D=a4+(beta.*d2.*cos(t));
p1=(N5/2).*(D);
p2=(1/2).*(D);
D1=sin(p1)./(N5.*sin(p2));  %非歸一化的陣因子K7
K7=abs(D1);
%------------------------------
E=a4+(beta.*d3.*cos(t));
q1=(N5/2).*(E);
q2=(1/2).*(E);
E1=sin(q1)./(N5.*sin(q2));  %非歸一化的陣因子K8
K8=abs(E1);
%-------------------------------
F=a4+(beta.*d4.*cos(t));
r_1=(N5/2).*(F);
r_2=(1/2).*(F);
F1=sin(r_1)./(N5.*sin(r_2));  %非歸一化的陣因子K9
K9=abs(F1);
%-----------------------繪圖函數
figure(4)
subplot(131);polar(t,K7);xlabel('f=3GHz,N=20,d=1/4波長,a=π/2');
subplot(132);polar(t,K8);xlabel('f=3GHz,N=20,d=1/2波長,a=π/2');
subplot(133);polar(t,K9);xlabel('f=3GHz,N=20,d=0.7波長,a=π/2');

%--------------------------------------------------------------------------
%---------------------------------------3D-天線方向圖
n_tehta = 130; %-------------------- 采樣視角點的仰角
n_phi = 130;   %--------------------采樣點的方向角 
[tehta,phi]=meshgrid(eps:pi./(n_tehta-1):pi,...    %meshgrid函數為矩形區域的設定范圍是epf<tehta<π   0<phi<2π
            0:2*pi./(n_phi-1):2*pi) ; 
 t3=tehta;
%-------------只有參數N改變的天線方向3D圖
M=a+(beta.*d.*cos(t3));    %----N1=4;N2=8;N3=12;
 z_1=(N1/2).*(M);
z_2=(1/2).*(M);
M1=sin(z_1)./(N1.*sin(z_2));  %非歸一化的陣因子K1
K_1=abs(M1);
radio_1 =K_1; 
X1=radio_1.*sin(tehta).*cos(phi); 
Y1=radio_1.*sin(tehta).*sin(phi); 
Z1=radio_1.*cos(tehta);
%-----------------------------------
M=a+(beta.*d.*cos(t3));
 z_3=(N2/2).*(M);
z_4=(1/2).*(M);
M2=sin(z_3)./(N2.*sin(z_4));  %非歸一化的陣因子K2
K_2=abs(M2);
radio_2 =K_2; 
X2=radio_2.*sin(tehta).*cos(phi); 
Y2=radio_2.*sin(tehta).*sin(phi); 
Z2=radio_2.*cos(tehta);
%-------------------------------------
M=a+(beta.*d.*cos(t3));
 z_5=(N3/2).*(M);
z_6=(1/2).*(M);
M3=sin(z_5)./(N3.*sin(z_6));  %非歸一化的陣因子K3
K_3=abs(M3);
radio_3 =K_3; 
X3=radio_3.*sin(tehta).*cos(phi); 
Y3=radio_3.*sin(tehta).*sin(phi); 
Z3=radio_3.*cos(tehta);
%------------------------------------3D繪圖函數
figure(5)
surf(X1,Y1,Z1);   %三維繪圖函數surf，采用偽彩色表示曲面的高度
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=4,d=1/2波長,a=π/2');
figure(6)
surf(X2,Y2,Z2);   
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=8,d=1/2波長,a=π/2');
figure(7)
surf(X3,Y3,Z3)  
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=12,d=1/2波長,a=π/2');
%----------------------------------------------
%--------------------只有陣列單元方向角a改變的天線方向圖
N4=10;d1=lambda/4;a1=0;a2=pi/2;a3=pi/2+pi/10;
A_3d=a1+(beta.*d1.*cos(t3));
 x_1=(N4/2).*(A_3d);
x_2=(1/2).*(A_3d);
A_1=sin(x_1)./(N4.*sin(x_2));  %非歸一化的陣因子K4
K_4=abs(A_1);
radio_4 =K_4; 
X4=radio_4.*sin(tehta).*cos(phi); 
Y4=radio_4.*sin(tehta).*sin(phi); 
Z4=radio_4.*cos(tehta);
%-----------------------------------
B_3d=a2+(beta.*d1.*cos(t3));
 y_1_3d=(N4/2).*(B_3d);
y_2_3d=(1/2).*(B_3d);
B_1=sin(y_1_3d)./(N4.*sin(y_2_3d));  %非歸一化的陣因子K5
K_5=abs(B_1);
radio_5 =K_5; 
X5=radio_5.*sin(tehta).*cos(phi); 
Y5=radio_5.*sin(tehta).*sin(phi); 
Z5=radio_5.*cos(tehta);
%------------------------------------
C_3d=a3+(beta.*d1.*cos(t3));
 v_1=(N4/2).*(C_3d);
v_2=(1/2).*(C_3d);
C_1=sin(v_1)./(N4.*sin(v_2));  %非歸一化的陣因子K6
K_6=abs(C_1);
radio_6 =K_6; 
X6=radio_6.*sin(tehta).*cos(phi); 
Y6=radio_6.*sin(tehta).*sin(phi); 
Z6=radio_6.*cos(tehta);
%-----------------------------------
figure(8)
surf(X4,Y4,Z4);   %三維繪圖函數surf，采用偽彩色表示曲面的高度
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=10,d=1/4波長,a=0');
figure(9)
surf(X5,Y5,Z5);   
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=10,d=1/4波長,a=π/2');
figure(10)
surf(X6,Y6,Z6)  
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=10,d=1/4波長,a=π/2+π/10');
%---------------------------------------------
%-------------------只有陣列單元間隔d改變的天線方向3D圖
N5=20;d2=lambda/4;d3=lambda/2;d4=0.7*lambda;a4=pi/2;
D_3d=a4+(beta.*d2.*cos(t3));
p_1=(N5/2).*(D_3d);
p_2=(1/2).*(D_3d);
D_1=sin(p_1)./(N5.*sin(p_2));  %非歸一化的陣因子K7
K_7=abs(D_1);
radio_7 =K_7; 
X7=radio_7.*sin(tehta).*cos(phi); 
Y7=radio_7.*sin(tehta).*sin(phi); 
Z7=radio_7.*cos(tehta);
%---------------------------------------
E_3d=a4+(beta.*d3.*cos(t3));
q_1=(N5/2).*(E_3d);
q_2=(1/2).*(E_3d);
E_1=sin(q_1)./(N5.*sin(q_2));  %非歸一化的陣因子K8
K_8=abs(E_1);
radio_8 =K_8; 
X8=radio_8.*sin(tehta).*cos(phi); 
Y8=radio_8.*sin(tehta).*sin(phi); 
Z8=radio_8.*cos(tehta);
%------------------------------------------
F_3d=a4+(beta.*d4.*cos(t3));
r_1_3d=(N5/2).*(F_3d);
r_2_3d=(1/2).*(F_3d);
F_1=sin(r_1_3d)./(N5.*sin(r_2_3d));  %非歸一化的陣因子K9
K_9=abs(F_1);
radio_9 =K_9; 
X9=radio_9.*sin(tehta).*cos(phi); 
Y9=radio_9.*sin(tehta).*sin(phi); 
Z9=radio_9.*cos(tehta);
%-----------------------------------
figure(11)
surf(X7,Y7,Z7);   %三維繪圖函數surf，采用偽彩色表示曲面的高度
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=20,d=1/4波長,a=π/2');
figure(12)
surf(X8,Y8,Z8);   
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=20,d=1/2波長,a=π/2');
figure(13)
surf(X9,Y9,Z9)  
camlight right
light
colorbar
axis image 
rotate3D on
title('f=3GHz,N=20,d=0.7波長,a=π/2');