【分享】温度梯度场中的声线路径

久病成欢

相关理论基础见附件，下面是matlab的实现代码

1. function ret=T1(x,y)
   
2. ret=2000/(0.05*(x^2+y^2)+1);

复制代码

1. clear all;clc;close all
   
2. 
   
3. xmin=-5;xmax=5;ymin=-5;ymax=5;
   
4. [X Y]=meshgrid(xmin:.1:xmax,ymin:.1:ymax);Z=2000./(0.05*(X.^2+Y.^2)+1);[C,h]=contour(X,Y,Z);text_handle = clabel(C,h);
   
5. view(0,90)
   
6. 
   
7. jj=0;
   
8. for ii=1:5:89
   
9. hold on
   
10. sita=-ii/180*pi;
    
11. jj=jj+1;
    
12. 
    
13. x_in=0;y_in=5;
    
14. 
    
15. xori=x_in;yori=y_in;
    
16. a_x=cos(sita);a_y=sin(sita);h= 0.1;t=0.1;flag=1;
    
17. 
    
18. D=0;MD=[];Time=0;MT=[];MP=[];
    
19. 
    
20. while (flag)
    
21. norm=sqrt(a_x^2+a_y^2);a_x=a_x/norm;a_y=a_y/norm;
    
22. 
    
23. xh=h*a_x+x_in;yh=h*a_y+y_in;
    
24. b1_x=a_y;b1_y=-a_x;b2_x=-a_y;b2_y=a_x;
    
25. x2=xh+t/2*b1_x;y2=yh+t/2*b1_y;x3=xh+t/2*b2_x;y3=yh+t/2*b2_y;
    
26. %%
    
27. 
    
28. z=20.08;c1=z*sqrt(T1(x_in,y_in));c2=z*sqrt(T1(x2,y2));c3=z*sqrt(T1(x3,y3));
    
29. % if(c2<c3)
    
30. %     tmp=x2;
    
31. %     x2=x3;
    
32. %     x3=tmp;
    
33. %     tmp=y2;
    
34. %     y2=y3;
    
35. %     y3=tmp;
    
36. % end
    
37. %     
    
38. %     c1=z*sqrt(T(x_in,y_in));c2=z*sqrt(T(x2,y2));c3=z*sqrt(T(x3,y3));
    
39. %     
    
40. %%
    
41. 
    
42. A=[x_in y_in 1;x2 y2 1; x3 y3 1];Co=A\[c1;c2;c3];
    
43. cx=Co(1);cy=Co(2);c0=Co(3);
    
44. 
    
45. R=(c1*t/(c2-c3));
    
46.     d=(R-sqrt(R^2-h^2));
    
47. % if(R<0)
    
48. %     R=-R;
    
49. %     d=(R-sqrt(R^2-h^2));
    
50. % end
    
51. 
    
52. l=-a_y;m=a_x;
    
53. xout=xh+d*l/(sqrt(l^2+m^2));
    
54. yout=yh+d*m/(sqrt(l^2+m^2));
    
55. Rm=1/R;
    
56. bx=h*Rm*l/(sqrt(l^2+m^2))+(l-d*Rm)*a_x;
    
57. by=h*Rm*m/(sqrt(l^2+m^2))+(l-d*Rm)*a_y;
    
58. a_x=bx;a_y=by;
    
59. 
    
60. %%
    
61. dleD=sqrt((x_in-xout)^2+(y_in-yout)^2);
    
62. cs1=x_in/sqrt(x_in^2+y_in^2);
    
63. cs2=xout/sqrt(xout^2+yout^2);
    
64. delT=1/(2*sqrt(cx^2+cy^2))*abs(log((1+cs1)/(1-cs1)*(1-cs2)/(1+cs2)));
    
65. 
    
66. D=D+dleD;
    
67. MD=[MD; D];
    
68. Time=Time+delT;
    
69. %%
    
70. MT=[MT;Time];
    
71. 
    
72. hold on
    
73. figure(1)
    
74. x_in=xout;
    
75. y_in=yout;
    
76. 
    
77. 
    
78. if(x_in>xmax  || x_in<xmin  ||y_in>ymax  ||y_in<ymin)
    
79.     flag=0;
    
80. end
    
81. 
    
82. if(flag)
    
83.     MP=[MP;x_in y_in];
    
84. plot(MP(:,1),MP(:,2),'linewidth',2.3);
    
85. set(gca,'fontsize',14)
    
86. drawnow;
    
87. 
    
88. end
    
89. 
    
90. end
    
91. plot(-MP(:,1),MP(:,2),'linewidth',2.3);
    
92. 
    
93. MDD(ii)=sqrt((MP(end,1)-xori)^2+ (MP(end,2)-yori)^2);
    
94. MHD(ii)=MD(end);
    
95. MMT(ii)=MT(end);
    
96. disp(['第' num2str(jj) '次飞行的时间、弧线距离、直线距离分别为：'   num2str(MT(end))          '          '  num2str(MHD(end))        '       '  num2str(MDD(end))   ])
    
97. end
    

复制代码