- 第一次写这个,主要是发现CSDN上有一些理论知识,
但不多 ,然而并没有关于这方面的代码,为了方便以后的学弟学妹,现将代码陈列如下:
clear;clc;close all;
z1=[0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0 900.0 1000.0];
z2=[1100.0 2000.0 3000.0 4000.0 5000.0];
c1=[1520.0 1517.0 1514.0 1511.0 1508.0 1505.0 1502.0 1499.0 1496.0 1493.0 1490.0];
c2=[1490.5 1495.0 1500.0 1505.0 1510.0];
z=[z1 z2];
c=[c1 c2];
z0=4900;%声源深度
a0=[-20:5:20]*pi/180; %入射角度
for j=1:length(a0)
ct=interp1(z,c,z0,'linear')/cos(a0(j)); %翻转深度的声速
if ct>=1520
zt1=0;
zt2=5000;
elseif ct<1520&&ct>=1510
zt1=interp1(c1,z1,ct);
zt2=5000;
elseif ct<1510
zt1=interp1(c1,z1,ct);
zt2=interp1(c2,z2,ct); %获得上下翻转深度 zt1 zt2
end
if a0(j)>=0
s=1; %判断声线开始的出射方向是偏上还是偏下
else
s=-1;
end
fc1=zt1:5:z0;
ds1=5;
fc2=z0:5:zt2;
ds2=5; %出射深度向上向下分别分层
zk=z0;
a=a0(j);
x=0;
for i=1:6000
ds=5;
zk(i+1)=zk(i)-s*ds;
ck(i+1)=interp1(z,c,zk(i+1));
a(i+1)=acos(ck(i+1)/ct);
dx=abs(zk(i)-zk(i+1))/tan((abs(a(i))+abs(a(i+1)))/2);
x(i+1)=x(i)+dx;
if zk(i+1)<=zt1||zk(i+1)>=zt2
s=-s;
end
if(x(i+1)>=1e5)
break;
end
%求每点的坐标
end
% subplot 122
plot(x,zk);
set(gca,'ydir','reverse');
% set(gca,'position',[0.1,0.1,0.2,0.85]);
axis([0,1e5,0,5000]);
xlabel('水平距离(m)');
ylabel('depth(m)');
title('z0=4900m时声线轨迹图');
grid;
hold on;
end
hold off;
% subplot 121
% plot(c,z);grid;
% set(gca,'ydir','reverse');
% % set(gca,'position',[0.1,0.1,0.2,0.85]
% xlabel('sound speed(m/s)');
% ylabel('depth(m)');
% title('声速分布图');
% xlim([min(c)-10 max(c)+10]);
**这些代码认真看的话大家都是能够看得懂的,其中一些核心的计算公式在《水声学原理》这本书上是有的。顺便吐槽一下,xl老师您布置的作业是真滴难!!!