不同深度的水下声源声线绘制

• 第一次写这个，主要是发现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老师您布置的作业是真滴难！！！