clc
clear;
f0=1000;   %信號頻率
fs=80*f0;   %采樣頻率 
Ns=2048;   %采樣點數 
N1=5;      %陣元數
c=1500;    %聲速
d=0.5*c/f0;      %陣元間距 
fai=0;     %信號入射方向 
w=2*pi*f0;
ph1=-180;                   
ph2=180;                   
phn=1;                  
ax=(ph2-ph1)/phn+1;
t0=d*sin(fai*pi/180)/c;
for n=1:N1;
    for k=1:Ns;
        p(n,k)=exp(-j*(2*pi*f0*(k-1)/fs-2*(n-1)*pi*f0*t0));
        vx(n,k)=p(n,k)*cos(fai*pi/180);
        vy(n,k)=p(n,k)*sin(fai*pi/180);  
    end
end

for i=1:ax;     
    num=ph1+(i-1)*phn;
    fain=num*pi/180;
    Xp1=0;
    Xv1=0;
    for n=1:N1;
        Xp(n,1,i)=p(n,1)*exp(-j*(n-1)*2*pi*d*f0*sin(fain)/c); 
        Xp1=Xp1+Xp(n,1,i);
         Xv(n,1,i)=Xp(n,1,i)*(cos(fai*pi/180)*cos(fain)+sin(fai*pi/180)*sin(fain));
         Xv1=Xv1+Xv(n,1,i);
    end
%        Rp(i)=abs((Xp1+Xv1)*Xv1);
       Rp(i)=abs(Xp1*Xv1);
end 
i=1:ax;
maxRp=max(Rp);
Rpp=20*log10(Rp(i)/maxRp);
figure
plot(ph1+(i-1)*phn,Rpp);
grid on
hold on
axis([-180 180 -100 0])
title('直線矢量陣指向性')
xlabel('角度')
ylabel('dB')