水声工程声线轨迹图绘制

lhaozong

有没有哪位大神教一教水声中声线轨迹图的绘制。
帮小弟列一组数据参数和声线轨迹图的matlab绘制程序，求指导！

久病成欢

1. clc;
   
2. clear all;
   
3. for theta=8:0.5:15
   
4. theta=theta*pi/180;
   
5. z1=50;
   
6. c1=1500-0.03*z1;     %声源深度的声速。
   
7. c2=1500-0.03*1000;   %1000m深度的声速。
   
8. x1=c1/cos(theta)/0.03*abs(sin(theta)-sqrt(1-((1500-1000*0.03)/(1500-z1*0.03)*cos(theta))^2));%1000m处第一次的水平距离。
   
9. theta2=acos((1500-1000*0.03)/(1500-z1*0.03)*cos(theta));      %第一次到1000m的射出角。
   
10. x2=2*c1/cos(theta)*abs(sin(theta2))/0.018;
    
11.       x3=2*c1/cos(theta)*abs(sin(theta2))/0.03;
    
12. 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)
    
13.   theta3=acos(1500/(1500-0.03*1000)*cos(theta2));
    
14.   x3=c1/cos(theta)*abs(sin(theta2)-sin(theta3))/0.03;
    
15. end
    
16.   theta4=acos((1500-1000*0.03)/(1500)*cos(theta3));
    
17.   x4=2*c1/cos(theta)*abs(sin(theta2))/0.018;
    
18.      n=1;
    
19.      for x=1:10:x1+x2+x3+x3+x4+2*x3
    
20.     if x<=x1
    
21.         z(n)=-sqrt((c1/cos(theta)/0.03)^2-(x+tan(theta)/(0.03/c1))^2)+(1/(0.03/c1)+z1);  %第一个半圆。
    
22.     elseif x<=x1+x2
    
23.         z(n)=sqrt((c2/cos(theta2)/0.018)^2-(x-tan(theta2)/(0.018/c2)-x1)^2)+(-1/(0.018/c2)+1000);
    
24.     elseif x<=x1+x2+x3
    
25.             z(n)=-sqrt((c2/cos(theta2)/0.03)^2-(x-tan(theta2)/(0.03/c2)-x1-x2)^2)+(1/(0.03/c2)+1000);
    
26.     elseif x<=x1+x2+x3+x3
    
27.          z(n)=-sqrt((1500/cos(theta3)/0.03)^2-(x+tan(theta3)/(0.03/1500)-x1-x2-x3)^2)+(+1/(0.03/1500));
    
28.     elseif x<=x1+x2+x3+x3+x4
    
29.         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);
    
30.     elseif x<=x1+x2+x3+x3+x4+x3
    
31.             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);
    
32.     elseif x<=x1+x2+x3+x3+x4+2*x3
    
33.          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));
    
34.     end  
    
35.     n=n+1;
    
36.        end  
    
37. x=1:10:x1+x2+x3+x3+x4+2*x3;
    
38. plot(x*0.001,z,'r');
    
39. % axis([1 150000 0 3500]);
    
40. set(gca,'YDIR','Reverse');                                                                        %反转Y坐标
    
41. xlabel('水平距离范围(km)');
    
42. ylabel('水深(m)');
    
43. title('近表面声源汇聚区示意图(声源深度100m，角度为8--15度)');
    
44. hold on
    
45. end

复制代码

lhaozong

久病成欢 发表于 2014-5-23 09:24

先学习学习，十分感谢