water1.m

声线图

clc;
clear all;
for theta=8:0.5:15
theta=theta*pi/180;
z1=50;
c1=1500-0.03*z1;     %声源深度的声速。
c2=1500-0.03*1000;   %1000m深度的声速。
 x1=c1/cos(theta)/0.03*abs(sin(theta)-sqrt(1-((1500-1000*0.03)/(1500-z1*0.03)*cos(theta))^2));%1000m处第一次的水平距离。
 theta2=acos((1500-1000*0.03)/(1500-z1*0.03)*cos(theta));      %第一次到1000m的射出角。
 x2=2*c1/cos(theta)*abs(sin(theta2))/0.018;
      x3=2*c1/cos(theta)*abs(sin(theta2))/0.03;
 if(-sqrt((c2/cos(theta2)/0.03)^2-(x1+x2+x3/2-tan(theta2)/(0.03/c2)-x1-x2)^2)+(1/(0.03/c2)+1000)<0)
  theta3=acos(1500/(1500-0.03*1000)*cos(theta2));
  x3=c1/cos(theta)*abs(sin(theta2)-sin(theta3))/0.03;
 end
  theta4=acos((1500-1000*0.03)/(1500)*cos(theta3));
  x4=2*c1/cos(theta)*abs(sin(theta2))/0.018;
     n=1;
     for x=1:10:x1+x2+x3+x3+x4+2*x3
    if x<=x1
        z(n)=-sqrt((c1/cos(theta)/0.03)^2-(x+tan(theta)/(0.03/c1))^2)+(1/(0.03/c1)+z1);  %第一个半圆。
    elseif x<=x1+x2
        z(n)=sqrt((c2/cos(theta2)/0.018)^2-(x-tan(theta2)/(0.018/c2)-x1)^2)+(-1/(0.018/c2)+1000);
    elseif x<=x1+x2+x3
            z(n)=-sqrt((c2/cos(theta2)/0.03)^2-(x-tan(theta2)/(0.03/c2)-x1-x2)^2)+(1/(0.03/c2)+1000);
    elseif x<=x1+x2+x3+x3
         z(n)=-sqrt((1500/cos(theta3)/0.03)^2-(x+tan(theta3)/(0.03/1500)-x1-x2-x3)^2)+(+1/(0.03/1500));
    elseif x<=x1+x2+x3+x3+x4
        z(n)=sqrt((c2/cos(theta2)/0.018)^2-(x-tan(theta2)/(0.018/c2)-x1-x2-x3-x3)^2)+(-1/(0.018/c2)+1000);
    elseif x<=x1+x2+x3+x3+x4+x3
            z(n)=-sqrt((c2/cos(theta2)/0.03)^2-(x-tan(theta2)/(0.03/c2)-x1-x2-2*x3-x4)^2)+(1/(0.03/c2)+1000);
    elseif x<=x1+x2+x3+x3+x4+2*x3
         z(n)=-sqrt((1500/cos(theta3)/0.03)^2-(x+tan(theta3)/(0.03/1500)-x1-x2-x3-2*x3-x4)^2)+(+1/(0.03/1500));
    end 
    n=n+1;
       end 
x=1:10:x1+x2+x3+x3+x4+2*x3;
plot(x*0.001,z,'r');
% axis([1 150000 0 3500]);
set(gca,'YDIR','Reverse');                                                                        %反转Y坐标
xlabel('水平距离范围(km)');
ylabel('水深(m)');
title('近表面声源汇聚区示意图(声源深度100m，角度为8--15度)');
hold on
end