function [WB SumOfAll fd FirstOrder]=AllOrder(f0,Vp,Q,QN,Q0,Reso)
%海浪散射系数仿真程序 个参数下面有说明 可以计算单基地雷达和双基地雷达的海浪散射系数
DopVal=zeros(1,401);
% f0=15; Vp=15; Q=0*pi/180;QN=90*pi/180.0;Q0=0*pi/180.0;Reso=2000;
g=9.8;
A=4/(3*pi);
s=4;
K0=2*pi*f0/300; %//无线电波长
WB=sqrt(2*g*K0*cos(Q0)); % //Bragg 基频
%/* 参数声明
% Vp---输入风速19.5(m/s)
% f0---雷达操作频率
% K0---雷达波数
% K1---一阶海洋波数
% K2---二阶海洋波数
% Wd---多普勒频率
% WB---Bragg多普勒频率
% Q ---风向（角度)
% QN---散射椭圆法线矢量方向
% Q0---半双基地角
% Q1---一阶海洋波方向
% Q2---二阶海洋波方向
% s ---扩展因子 */
% //计算二阶散射截面积
% //-3WB<Wd<WB部分，步长0.025WB
Singular1=0;
Singular2=0;
Singular3=0;
Singular4=0;
Singular5=0;
Singular6=0;
Singular7=0;
Singular8=0;
Singular9=0;
Singular10=0;
for N=1:160
Wd=(-5.0+0.025*(N-1))*WB;
if(Wd+sqrt(2.0)*WB>0) %// Wd>=-sqrt(2)*WB部分
if(N==145)
Wd=-sqrt(2)*WB; %// 奇异值
Singular1=N;
end
if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular2=N;
%//电场耦合系数产生的奇异值
end
DopVal(N)=dif(QN-pi,QN,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
else %// Wd<-sqrt(2)WB
if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular3=N;
%//电场耦合系数产生的奇异值
end
if(abs(Wd+(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=-(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB;
Singular4=N;
end
DopVal(N)=dif(QN-pi,QN-acos(2*(WB/Wd)*(WB/Wd)),1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN+acos(2*(WB/Wd)*(WB/Wd)),QN+pi,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
end
end
% //计算-WB<Wd<0部分
for N=162:200
Wd=(-5+0.025*(N-1))*WB;
if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular5=N;
end
DopVal(N)=dif(QN-pi,QN,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
end
N=N+1;
if(N==201) % //Wd=0
Wd=(-5+0.025*(N-1))*WB;
DopVal(N)=dif(QN-0.5*pi,QN+0.5*pi,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
end
for N=202:240
Wd=(-5+0.025*(N-1))*WB;
if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular6=N;
%//电场耦合系数产生的奇异值
end
DopVal(N)=dif(QN-pi,QN,1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
%// printf("%d,%e\n",N,DopVal[N]);
end
%//计算WB<Wd<3WB部分
for N=242:401
Wd=(-5+0.025*(N-1))*WB;
if(Wd-sqrt(2.0)*WB<0) %// Wd>=sqrt(2)*WB部分
if(N==257)
Wd=sqrt(2)*WB; %// 奇异值
Singular7=N;
end
if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular8=N;
end
%//电场耦合系数产生的奇异值
DopVal(1,N)=dif(QN-pi,QN,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
else % //Wd>sqrt(2)*WB
if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB;
Singular9=N;
end
if(abs(Wd-(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB)
Wd=(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB;
Singular10=N;
end
DopVal(N)=dif(QN-pi,QN-acos(2*(WB/Wd)*(WB/Wd)),-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN+acos(2*(WB/Wd)*(WB/Wd)),QN+pi,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0);
end
end
DopVal(201)=(DopVal(200)+DopVal(202))/2;
SecondDopVal=DopVal;
SecondDopVal1=DopVal;
SecondDopVal=SmoothScecondOrderData(SecondDopVal,Singular1,Singular2,Singular3,Singular4,Singular5,Singular6,Singular7,Singular8,Singular9,Singular10);
%计算一阶海杂波
FirstOrder=zeros(1,401);
for N=1:401
if(N-1>200)
m=-1;
else
m=1;
end
K=((-5+0.025*(N-1))*WB)^2/g;
if(N-1==160|N-1==240)
Val=(2^4)*pi*(K0^2)*0.00405*A*(K^-4)*exp(-0.74*g*g*(K^-2)*(Vp^-4))*(cos(((1-m)/2*pi+QN-Q)/2.0)^s)*(K^2.5)*cos(Q0)/sqrt(g)*Reso*(cos(Reso/2*(K/cos(Q0)-2*K0))^2);
else
Val=(2^4)*pi*(K0^2)*0.00405*A*(K^-4)*exp(-0.74*g*g*(K^-2)*(Vp^-4))*(cos(((1-m)/2*pi+QN-Q)/2.0)^s)*(K^2.5)*cos(Q0)/sqrt(g)*Reso*(sin(Reso/2*(K/cos(Q0)-2*K0))/(Reso/2*(K/cos(Q0)-2*K0))^2);
end
if(Val>=FirstOrder(N))
FirstOrder(N)=Val;
else
end
end
FirstOrder(201)=(FirstOrder(200)+FirstOrder(202))/2;
FirstOrder=dBConverter(FirstOrder);
FirstOrder=SmoothData(FirstOrder);
FirstOrder=10.^(FirstOrder/10);
%给出一阶和二阶海杂波RCS的和
SumOfAll=SecondDopVal+FirstOrder;
%给出多普勒坐标fd
N=1:401;
Wd=(-5.0+0.025*(N-1))*WB;
fd=Wd;
%plot(fd,10*log10(SumOfAll));
%%定义求数组转化成dB表示的函数
function dBVal=dBConverter(AA)
for L=1:length(AA)
if(AA(L)<=0.000000)
AA(L)=-90.0;
else
AA(L)=10*log10(AA(L));
end
if(AA(L)<-90.0)
AA(L)=-90.0;
end
end
dBVal=AA;
%定义平滑曲线的函数
function SmoothVal=SmoothData(B)
%求极大值点，比把极大值点坐标保存在jj数组里面，之后把极大值点数据保存在flav里面
jj=[];
for j=2:length(B)-1
if((B(j)>(B(j-1)))&(B(j)>B(j+1)))
jj(end+1)=j;
end
end
jj=[1 jj length(B)];
jjj=length(jj);
flatv=[];
for i=1:jjj
flatv(i)=B(jj(i));
end
%插值进行曲线平滑
i=1:length(B);
SmoothVal=interp1(jj,flatv,i,'spline');
%定义功能函数平滑二阶数据
function OutSmoothResult=SmoothScecondOrderData(InputSecondOderData,Singular1,Singular2,Singular3,Singular4,Singular5,Singular6,Singular7,Singular8,Singular9,Singular10)
NN=length(InputSecondOderData);
XX=InputSecondOderData;
for i=1:NN
if(XX(i)<=0)
XX(i)=-90;
else
XX(i)=10*log10(XX(i));
end
if(XX(i)<-90)
XX(i)=-90;
end
end
for i=2:NN-1
if(i==Singular1|i==Singular2|i==Singular3|i==Singular4|i==Singular5|i==Singular6|i==Singular7|i==Singular8|i==Singular9|i==Singular10)
if(XX(i)>XX(i-1)+15)
XX(i)=XX(i-1)+15;
end
else
if((XX(i)>XX(i-1)&XX(i)>XX(i+1))&(XX(i)>(max(XX(i-1),XX(i+1))+2)))
XX(i)=(XX(i-1)+XX(i+1))/2;
end
end
end
XX=10.^(XX./10);
OutSmoothResult=XX;
%定义求二阶海杂波积分的函数
function val=dif(LowLimit,UpLimit,m1,m2,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)
C=3.0e8; %光速
g=9.8;
A=4/(3*pi);
s=4;
val=0;
Sigma2=0.0; %初始化积分值
n3=m1*m2;
if(abs(Wd)>sqrt(2)*WB)
Beta=acos(2*(WB/Wd)*(WB/Wd));
else
Beta=0;
end
h=0.5*pi/180; %Theta=0.5的积分步长
Q1=LowLimit; %积分下限
while(Q1<=UpLimit)
Sigma=0.0; %初始积分值
% 在给定Wd,Q1时，用牛顿法计算海洋波数K1*/
K_sqrt=0.00000001; %猜测的K1的初始均方根值
fK_sqrt=10;
while((abs(fK_sqrt))>=0.000001)
if(K_sqrt<0)
K_sqrt=-K_sqrt;
else
K_sqrt=K_sqrt;
end
fK_sqrt=Wd+sqrt(g)*(m1*K_sqrt+m2*((((K_sqrt^4)+((2*K0*cos(Q0))^2)-4*K_sqrt*K_sqrt*K0*cos(Q0)*cos(Q1-QN)))^0.25));
% fK_sqrt=Wd+m1*sqrt(g*K1)+m2*sqrt(g*K2)
fK_sqrtP=sqrt(g)*(m1+m2*K_sqrt*(K_sqrt*K_sqrt-2*K0*cos(Q0)*cos(Q1-QN))*((((K_sqrt^4)+((2*K0*cos(Q0))^2)-4*K_sqrt*K_sqrt*K0*cos(Q0)*cos(Q1-QN)))^-0.75));
% fK_sqrtP 是fK_sqrt的微分
K_sqrt=K_sqrt-fK_sqrt/fK_sqrtP; %K_sqrt的新值
end
K1=K_sqrt*K_sqrt; %一阶海洋波
