function [Es,tsd,sige,sigspm,h_rms]=mom_sin_HH_PM(theta_i,t,R,theta_s,nsa,L,N,g,U,phase1,phase2,direction) %入射角,时间,水波波数与布拉格波数比值,水波幅度,雷达距离,雷达方位,散射方位计算点数,积分区域长度,积分细分点数
% function mom_sin_HH2
% theta_i = 60;t=0;amp=0.02;
% N=500;
% L=10;
% g=L/4;
euler=1.78107;
w_wave=13000000*2*pi; %入射波频率
miu_0=4*pi*10^-7;
e0=8.85418781*10^-12;
k = sqrt(w_wave^2*miu_0*e0); %入射波数
tai=theta_i*pi/180; %入射角
% k1=n_water*k*sin(tai); %水波波数
% k2=4/3*k*sin(tai); %水波波数
% k3=2*k*sin(tai); %水波波数
ti=tan(tai);
ci=cos(tai);
denom=1+2*ti^2;
denom=denom/(2*(k*g*ci)^2);
denom=8*pi*k*g*sqrt(pi/2)*ci*(1-denom);
% x=(-L/2+0.5*L/N:L/N:L/2);
dx=L/N;
% f=2*amp*sin(k1*x-sqrt(9.8*k1)*t)+amp*cos(k3*x-sqrt(9.8*k3)*t); %水波轮廓
% df=2*k1*amp*cos(k1*x-sqrt(9.8*k1)*t)+k3*amp*cos(k3*x-sqrt(9.8*k3)*t); %水波轮廓导数
% M=3000;%傅里叶级数总数
% L=10000;%傅里叶级数展开长度
% m=1:M;
% w=sqrt(9.8*m*2*pi/L);
% % phase = 2*pi*rand(1,M);
% F_PM = 8.1e-3*9.8^2*exp(-0.74*(9.8/U./w).^4)./w.^5;%%Pierson-Moskwitz 无向谱
% f=zeros(1,N);
% for n=1:N
% f(n)=f(n)+sum(sqrt(F_PM*(2*pi/L)).*exp(1i*2*pi*x(n)*m/L + 1i*phase).*(exp(1i*sqrt(9.8*2*pi*m/L)*t)+exp(-1i*sqrt(9.8*2*pi*m/L)*t)) + sqrt(F_PM*(2*pi/L)).*exp(-1i*2*pi*x(n)*m/L - 1i*phase).*(exp(1i*sqrt(9.8*2*pi*m/L)*t)+exp(-1i*sqrt(9.8*2*pi*m/L)*t)));
% end
% df=diff(f);
% df(N)=df(N-1);
[f,df,ddf,x,h_rms]=rsgeno_t1(N,L,0.01,2,U,phase1,phase2,t,direction);
b=incid(k,x,f,tai,g); %产生矩阵b
[xm,xn]=meshgrid(x);
[fm,fn]=meshgrid(f);
arg=k*sqrt((xm-xn).^2+(fm-fn).^2);
for m=1:N
arg(m,m)=1;
end
A=dx*i*besselh(0,1,arg)/4; %产生矩阵A
for m=1:N
dl=sqrt(1+df(m)*df(m))*dx;
A(m,m)=(i/4)*dx*(1+(2*i/pi)*(log(euler*k*dl/4)-1));
end
b=b.';
u=A\b;
u=u.';
fac=((x+0*ti)/g).^2;
kg=(k*g*ci)^2;
w=(2*fac-1)/kg;
% u=u.^(1+w).*exp(-fac);
% u=u.*exp(-fac);
% u = u.*hanning(N,'periodic')';
dan=180/(nsa+1);
for m=1:nsa
tsd(m)=-90+m*dan;
tas=tsd(m)*pi/180;
ss=sin(tas);
cs=cos(tas);
csa(m)=cs^2;
kxq=k*(ss-sin(tai));
kx1(m)=kxq;
integ=exp(-1i*k*(ss*x+f*cs));
integ=integ.*u*dx;
psis=sum(integ);
sige(m)=abs(psis)^2/denom;
% sig(m)=(sige+(ir-1)*sig(m))/ir;
end
% sigspm=csa.*(amp^2)*4*k^3*cos(tai);
sigspm=csa.*wkos(kx1,U)*4*k^3*cos(tai);
tas=theta_s*pi/180;
ss=sin(tas);
cs=cos(tas);
integ=exp(-1i*k*(ss*x+f*cs));
integ=integ.*u*dx;
psis=sum(integ);
Es=1i*sqrt(2/pi/k/R)*exp(-1i*pi/4+1i*k*R)/4*psis;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% incid.m %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b=incid(k,x,z,tai,g)
% generates the spatial tapered incident wave
ti=tan(tai);
ci=cos(tai);
si=sin(tai);
fac=((x+z*ti)/g).^2;
kg=(k*g*ci)^2;
w=(2*fac-1)/kg;
% b=1*exp(i*k*(x*si-z*ci)); %使用平面波
b=1*exp(i*k*(x*si-z*ci).*(1+w)-fac); %使用锥形波
%%%%%%%%%%%%%%%%%%%%
% 1D ocean spectrum %
%%%%%%%%%%%%%%%%%%%%%
function y=wkos(kx,us)
kr=abs(kx);
%a0=0.008;
a0=0.004;
a=0.225;
b=1.25;
kj=2;
gs=9.81+7.25e-5*kx.^2;
z0=6.84e-5/us+4.28e-3*us^2-4.43e-4;
U195=us/0.4*log(19.5/z0);
kc=9.81/U195^2;
i1=find(kr>kj);
kr1=kr(i1);
y(i1)=(a0./kr1.^3).*(b*kr1*us^2./gs(i1)).^(a*log10(kr1/kj));
i2=find(kr<=kj);
kr2=kr(i2);
y(i2)=(a0./kr2.^3).*exp(-0.74*(kc./kr2).^2);
评论7