clc
clear all
f=35;
c=300;
len=c/f;
d=5;
N=16;
k=2*pi/len;
%口径场
SEL=2;%自由输入1，泰勒分布2
PLL=2;%产生0~1的随机数1，正态分布2
SLR=30; % 泰勒分布副瓣电平要求
n=5; % n
Snx=zeros(1,n-1);
dh=zeros(1,N);
x=zeros(1,n-1);
znx=zeros(1,(n-1));%零点个数n-1
f2=zeros(1,(n-1));
f3=zeros(1,(n-1));
I=zeros(1,N);
% 单元布局
for t=1:N
dh(t)=(t-N/2-1/2)*d;
end
figure(1)
plot(dh,0,'.b');
if SEL==1;
I=input('\n请输入自定义电流值\n');
elseif SEL==2
%泰勒分布
%taylor分布电流计算
r0x=10^(SLR/20);%由副瓣电平求主副瓣电平比
Ax=acosh(r0x)/pi;
taox=n/sqrt(Ax^2+(n-1/2)^2);%波瓣展宽因子
for t=1:(n-1)
x(t)=1;
for t1=1:(n-1)
znx(t1)=taox*sqrt(Ax^2+(t1-.5)^2);%Xn（变量t1，维数为n-1）
x(t)=x(t)*(1-(t/znx(t1))^2);%叠乘
end
f1=factorial(n-1);
f2(t)=factorial(n+t-1);
f3(t)=factorial(n-1-t);
Snx(t)=f1^2*x(t)/f2(t)/f3(t);%阵列函数X（变量t，x维数为n-1）
end
for t=1:N
I(t)=1+2*sum(Snx.*cos(2*pi*(1:(n-1))*dh(t)/(N*d)));%激励系数（求和）
end
end
%误差
if PLL==1;
detx=rand(1,N);
phix=random('unif',0,0,1,N);
elseif PLL==2;
mu1=0;
sigma1=0.2;
detx=normrnd(mu1,sigma1,1,N);
mu2=0;
sigma2=6*pi/180;
phix=normrnd(mu2,sigma2,1,N);
a=zeros(1,N);
end
I=I+detx;
I=I/max(max(I));
for i=1:N
a(i)=i*20*pi/180;
end
%方向图
Theta=0;%扫描角
Theta=Theta*pi/180;
theta=[-90:0.2:90]*pi/180;
s1=zeros(1,numel(theta));
for t=1:numel(theta);
for u=1:N
s1(t)=abs(sum(I.*exp(j*k*dh*(sin(theta(t))-sin(Theta)))*exp(j*phix(u))))*exp(j*a(u));
end
end
s1=20*log10(s1/max(s1));
figure(2)
plot(theta*180/pi,s1,'b'),axis([-90,90,-90,0]);grid on
title(['线源方向图（' num2str(f) 'GHz,扫描' num2str(Theta*180/pi) '度）']);
%副瓣电平
a=findpeaks(s1,'sortstr','descend')
[min_value min_position]=min(a);%先得到最大值的数值和位置
a(min_position)=max(a);%将最大值的数值用向量最小值替代，这样第二大的值就变成了最大值，且所在位置不变
[min_value_2 min_position_2]=min(a)%这时取出的最大值就是我们所需要的第二大值了
a(min_position)=min_value;%记得将刚才的最大值复原，保持向量的完整性
times=100;
minvalue=zeros(1,times);
for m=1:times
if SEL==1;
I=input('\n请输入自定义电流值\n');
elseif SEL==2
%泰勒分布
%taylor分布电流计算
r0x=10^(SLR/20);%由副瓣电平求主副瓣电平比
Ax=acosh(r0x)/pi;
taox=n/sqrt(Ax^2+(n-1/2)^2);%波瓣展宽因子
for t=1:(n-1)
x(t)=1;
for t1=1:(n-1)
znx(t1)=taox*sqrt(Ax^2+(t1-.5)^2);%Xn（变量t1，维数为n-1）
x(t)=x(t)*(1-(t/znx(t1))^2);%叠乘
end
f1=factorial(n-1);
f2(t)=factorial(n+t-1);
f3(t)=factorial(n-1-t);
Snx(t)=f1^2*x(t)/f2(t)/f3(t);%阵列函数X（变量t，x维数为n-1）
end
for t=1:N
I(t)=1+2*sum(Snx.*cos(2*pi*(1:(n-1))*dh(t)/(N*d)));%激励系数（求和）
end
end
%误差
if PLL==1;
detx=rand(1,N);
phix=random('unif',0,0,1,N);
elseif PLL==2;
mu1=0;
sigma1=0.2;
detx=normrnd(mu1,sigma1,1,N);
mu2=0;
sigma2=6*pi/180;
phix=normrnd(mu2,sigma2,1,N);
a=zeros(1,N);
end
I=I+detx;
I=I/max(max(I));
for i=1:N
a(i)=i*20*pi/180;
end
%方向图
Theta=0;%扫描角
Theta=Theta*pi/180;
theta=[-90:0.2:90]*pi/180;
s1=zeros(1,numel(theta));
for t=1:numel(theta);
for u=1:N
s1(t)=abs(sum(I.*exp(j*k*dh*(sin(theta(t))-sin(Theta)))*exp(j*phix(u))))*exp(j*a(u));
end
end
s1=20*log10(s1/max(s1));
figure(2)
plot(theta*180/pi,s1,'b'),axis([-90,90,-90,0]);grid on
title(['线源方向图（' num2str(f) 'GHz,扫描' num2str(Theta*180/pi) '度）']);
%副瓣电平
a=findpeaks(s1,'sortstr','descend');
[min_value min_position]=min(a);%先得到最大值的数值和位置
a(min_position)=max(a);%将最大值的数值用向量最小值替代，这样第二大的值就变成了最大值，且所在位置不变
[min_value_2 min_position_2]=min(a);%这时取出的最大值就是我们所需要的第二大值了
a(min_position)=min_value;%记得将刚才的最大值复原，保持向量的完整性
minvalue(m)=min_value_2;
end
mean(minvalue)
% %立体方向图
% phi=(0:1:360)*pi/180;
% the=(-90:1:90)*pi/180;
% tt3=zeros(numel(phi),numel(the));
% tt4=zeros(numel(phi),numel(the));
% for t=1:numel(phi)
% for t1=1:numel(the)
% tt3(t,t1)=the(t1)*cos(phi(t));
% tt4(t,t1)=the(t1)*sin(phi(t));
% f1(t,t1)=abs(sum(I.*exp(j*k*dh*(sin(the(t1))-sin(Theta)))*exp(j*phix(u))));
% end
% end
% figure(3)
% mesh(tt3,tt4,10*log10(f1/max(max(f1))));
