%海浪建模
clear all;
clc;
dx=0.5;
dy=0.5;
Xm=-10:dx:10;
Yn=-10:dy:10;
M=size(Xm,2);
N=size(Yn,2);
Lx=20;
Ly=20;
r=3.3;
U=15; %10m处的风速为10m/s??
z=0; %风向
b=(9.8.*5000)./(U.^2);
a=0.076.*(b.^(-0.22));
w0=22.*(b.^(-0.33));
Q1=0;
Q2=0;
[xm,yn]=meshgrid(Xm,Yn);
for mk=(-M/2+1):1:(M/2)
Kmk=(2.*pi.*mk)./Lx;
for nk=(-N/2+1):1:(N/2)
Knk=(2.*pi.*nk)./Ly;
w=sqrt(9.8.*sqrt((Kmk).^2+(Knk).^2));
d=atan(Knk./Kmk);
if ((-pi/2)<(d-z)<(pi/2))
G=(2./pi).*((cos(d-z)).^2);
else
G=0;
end
if(w<=w0)
c=0.07;
else
c=0.09;
end
if (-M/2+1<=mk<0&&0<mk<(M/2)&&-N/2+1<=nk<0&&0<nk<(N/2))
A=(randn(M,N)+i.*randn(M,N))./sqrt(2);
else
A=randn(M,N);
end
S=a.*9.8.*(1./((sqrt(9.8.*sqrt(Kmk.^2+Knk.^2))).^5)).*exp(-1.25.*((w0./(sqrt(9.8.*sqrt(Kmk.^2+Knk.^2)))).^4)).*r.^(-((sqrt(9.8.*sqrt(Kmk.^2+Knk.^2))-w0).^2)./(2.*(c.^2).*(w0.^2))).*(9.8./(2.*sqrt(9.8.*sqrt(Kmk.^2+Knk.^2)))).*1/(Kmk*(Knk^2/Kmk^2 + 1)).*(Kmk/(Kmk^2 + Knk^2)^(1/2)).*G;
F=2*pi*(((Lx.*Ly).*S).^0.5).*exp(i.*w.*0).*A;
Q=F.*exp(i.*(Kmk.*xm+Knk.*yn));
Q1=Q1+Q;
end
Q2=Q2+Q1;
end
f=(1./(Lx*Ly)).*ifft(Q2);
f0=real(f);
surf(xm,yn,f0);
shading interp;
xlabel('x/m'),ylabel('y/m'),zlabel('z/m');
hold on
%光斑区域f数组计算
dx=0.1;
dy=0.1;
R=1; %光束半径
Xm=-R:dx:R;
Yn=-R:dy:R;
M=size(Xm,2);
N=size(Yn,2);
Lx=20;
Ly=20;
r=3.3;
U=15; %10m处的风速
z=0; %风向
b=(9.8.*5000)./(U.^2);
a=0.076.*(b.^(-0.22));
w0=22.*(b.^(-0.33));
Q1=0;
Q2=0;
[xm,yn]=meshgrid(Xm,Yn);
for mk=(-M/2+1):1:(M/2)
Kmk=(2.*pi.*mk)./Lx;
for nk=(-N/2+1):1:(N/2)
Knk=(2.*pi.*nk)./Ly;
w=sqrt(9.8.*sqrt((Kmk).^2+(Knk).^2));
d=atan(Knk./Kmk);
if ((-pi/2)<(d-z)<(pi/2))
G=(2./pi).*((cos(d-z)).^2);
else
G=0;
end
if(w<=w0)
c=0.07;
else
c=0.09;
end
if (-M/2+1<=mk<0&&0<mk<(M/2)&&-N/2+1<=nk<0&&0<nk<(N/2))
A=(randn(M,N)+i.*randn(M,N))./sqrt(2);
else
A=randn(M,N);
end
S=a.*9.8.*(1./((sqrt(9.8.*sqrt(Kmk.^2+Knk.^2))).^5)).*exp(-1.25.*((w0./(sqrt(9.8.*sqrt(Kmk.^2+Knk.^2)))).^4)).*r.^(-((sqrt(9.8.*sqrt(Kmk.^2+Knk.^2))-w0).^2)./(2.*(c.^2).*(w0.^2))).*(9.8./(2.*sqrt(9.8.*sqrt(Kmk.^2+Knk.^2)))).*1/(Kmk*(Knk^2/Kmk^2 + 1)).*(Kmk/(Kmk^2 + Knk^2)^(1/2)).*G;
F=2*pi*(((Lx.*Ly).*S).^0.5).*exp(i.*w.*0).*A;
Q=F.*exp(i.*(Kmk.*xm+Knk.*yn));
Q1=Q1+Q;
end
Q2=Q2+Q1;
end
f=(1./(Lx*Ly)).*ifft(Q2);
ff=real(f);
%光线绘制
a0=[1,0,0];
b0=[0,1,0];
c0=[0,0,1];
ct=pi./6; %入射角
P=M;
Q=N;
for n=1:Q-1
for m=1:P-1
x0=Xm(m);
y0=Yn(n);
f0=ff(m,n);
a1=[0.1,0,(ff(m+1,n)-ff(m,n))];
b1=[0,0.1,(ff(m,n+1)-ff(m,n))];
%a1=a1./norm(a1);
%b1=b1./norm(b1);
c1=cross(a1,b1);
c1=c1./norm(c1);
if c1(3)>0
n0=c1;
else
n0=-c1;
end
k=a0.*sin(ct).*cos(0)+b0.*sin(ct).*sin(ct)-c0.*cos(ct);
r0=k-2.*(k.*n0).*n0;
if r0(3)>0
r0=r0;
else
r0=-r0;
end
if ((y0^2)+(x0^2))<=R^2 %/((R/sin(ct))^2)
%绘制反射光线
t=0:0.5:10;
X=x0+r0(1).*t;
Y=y0+r0(2).*t;
Z=f0+(r0(3)).*t;
plot3(X,Y,Z,'r');
hold on;
%绘制入射光线
k=-k;
t=0:0.5:10;
X0=x0+k(1).*t;
Y0=y0+k(2).*t;
Z0=f0+(k(3)).*t;
plot3(X0,Y0,Z0,'g');
hold on
else
end
end
end
没有合适的资源?快使用搜索试试~ 我知道了~
资源详情
资源评论
资源推荐
收起资源包目录
roughsea.zip (2个子文件)
roughsea.m 1KB
OceanRefalction.m 4KB
共 2 条
- 1
小贝德罗
- 粉丝: 70
- 资源: 1万+
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功
评论4