搜索
声波传播路径程序说明一、程序流程图读入温度数据设定初始值设定前向三角形形状参数设置打靶入射角度的y方向分量k（j）J=3求取本次k（j）的值求解入射角度及三角形其他2个顶点坐标判断三角形是否越出边界求取三角形出射角坐标及方向求出本次打靶结束点与指定点的距离清除部分变量的值得到能打到指定点的入射角度依照之前方法得到声线创博路径上各点坐标计算相邻点之间声线传播的路程及时间累积出总路程和总时间，并在命令框中显示；画图K(1),k(2)是否打到指定点求取三角形底边中点坐标及3个顶点的温度、声速大小声线是否到达边界求取本次打靶结束时点的坐标求取三角形内声速分布函数及声线圆弧圆心坐标是否是是否否是否程序流程图二、变量含义说明1、求解结果必须用到的变量：[z]：读入温度数据时的温度向量；tt：存储温度数据的51*51的温度矩阵；T：将温度转化为声速大小的函数；s：某次打靶结束时点与指定点的距离；e：循环终止判定的阈值；i，j：记录循环次数的变量；x_int，y_int：某一声线的入射点坐标；h，t：前向展开三角形的高和底边长度；thelt1：前向展开三角形的顶角半角的大小；z：温度转换为声速的常数；x_over，y_over：指定点的坐标（声线结束点的坐标）；[xx_over]，[yy_over]：记录每次打靶结束时点坐标的向量；[x]，[y]：用于存储声波传播路径上离散点的坐标；[k]：每次打靶时初始入射向量的y方向的分量大小；[ax]，[ay]：声波传播路径上离散点对应的入（出）射方向向量；[thelt2]：记录声波传播路径上离散点的入（出）射方向向量所对应的角度；[x1]，[y1]，[x2]，[y2]：声波传播路径上离散点对应三角形的其他2个顶点的坐标；xh，yh：三角形中底边中点的坐标；[c]，[c1]，[c2]：三角形3个顶点对应的声速大小；[cx]，[cy]，[c0]：某段声线圆弧曲线对应声速分布函数的系数；[x0]，[y0]：某段声线圆弧曲线对应的圆心坐标；r：某段声线圆弧曲线对应的曲率半径大小；d：三角形底边中心与三角形出射点的距离；d1：声波传播结束前倒数第二个离散点与边界线的距离；（1：是数字）dl：某三角形内声波实际传播路径的距离；（l是字母）dc：声波传播路径上相邻两离散点的声速差值；[dtime]：某三角形内声波实际传播的时间；D0：两点之间直线距离；D1：声波实际传播路径的总距离；time：声波实际传播的总时间；2、求解结果时不需要用到的变量（可舍弃）：[thelt3]：某三角形对应的声速分布函数的某条直线与x轴正方向的夹角；[thelt4]，[thelt5]：某三角形的入射点、出射点和圆心的连线与x轴正方向的夹角；[dthelt1]、[dthelt2]：某三角形的入射点、出射点和圆心的连线与声速分布函数的某条直线的夹角；三、程序语句具体说明1、主程序-------------------------设置图片的界线clearfiguregridaxis([-3.7853.785-3.7853.785])holdon--------------------------读入温度数据，并将其存入51*51的tt矩阵中[z]=xlsread('Temperature','c2:c2602');globaltt;N=51;fori=1:Nforj=1:Ntt(i,j)=z((i-1)*N+j);endend---------------------------部分参数的初始化s=100;e=0.01;j=1;x_int=2.0482;y_int=3.785;-----------------------打靶算法开始whileseh=0.01;%以下为部分变量的初始化；t=0.01;thelt1=atand(t/2/h);z=19.08;x_over=3.785;y_over=1.5304;x(1)=x_int;y(1)=y_int;k(1)=-1;k(2)=-2;ifj=3%当j=3时，采用线性差值得出下次打靶时初始的入射方向矢量；k(j)=k(j-1)+(k(j-1)-k(j-2))*(y_over-yy_over(j-1))/(yy_over(j-1)-yy_over(j-2));endax(1)=1;ay(1)=k(j);i=1;whileabs(x(i))=3.785&&abs(y(i))=3.785%判断声线是否超出边界；thelt2(i)=atand(ay(i)/ax(i));%某个三角形的入射角度（在-90~90之间）；ifax(i)0&&ay(i)0%以下两个if语句用于修正入射角度（-180~180之间）；thelt2(i)=thelt2(i)+180;endifax(i)0&&ay(i)0thelt2(i)=thelt2(i)-180;end%以下4个式子得到三角形另外2顶点的坐标；x1(i)=x(i)+sqrt(h^2+t^2/4)*cosd(thelt1+thelt2(i));y1(i)=y(i)+sqrt(h^2+t^2/4)*sind(thelt1+thelt2(i));x2(i)=x(i)+sqrt(h^2+t^2/4)*cosd(thelt2(i)-thelt1);y2(i)=y(i)+sqrt(h^2+t^2/4)*sind(thelt2(i)-thelt1);ifabs(x1(i))3.785||abs(x2(i))3.785||abs(y1(i))3.785||abs(y2(i))3.785%判断三角形是否越界；并且求出越界时最后一个离散点的坐标（此时可知道x或y的坐标，先判断再利用出射方向得到另外的坐标；ifabs(x(i))=abs(y(i))ify(i)0d1=abs(3.785-y(i));y(i+1)=3.785;x(i+1)=x(i)+d1/ay(i)*ax(i);elsed1=abs(3.785+y(i));y(i+1)=-3.785;x(i+1)=x(i)+d1/abs(ay(i))*ax(i);endelseifx(i)0d1=abs(3.785-x(i));x(i+1)=3.785;y(i+1)=y(i)+d1/ax(i)*ay(i);elsed1=abs(3.785+x(i));x(i+1)=-3.785;y(i+1)=y(i)+d1/abs(ax(i))*ay(i);endendxx_over(j)=x(i+1);yy_over(j)=y(i+1);breakend%以下几个式子可以得到三角形底边中点坐标，顶点对应的声速；xh=(x1(i)+x2(i))/2;yh=(y1(i)+y2(i))/2;c(i)=z*sqrt(T(x(i),y(i)));%从温度得到声速的转换关系；c1(i)=z*sqrt(T(x1(i),y1(i)));c2(i)=z*sqrt(T(x2(i),y2(i)));%通过解方程得到三角形对应的声速场参数；A1=[x(i)y(i)1;x1(i)y1(i)1;x2(i)y2(i)1];B1=[c(i)c1(i)c2(i)]';X1=A1\B1;cx(i)=X1(1);cy(i)=X1(2);c0(i)=X1(3);%同样通过解方程得到三角形对应的圆弧曲线的圆心坐标；A2=[cx(i)cy(i);ax(i)ay(i)];B2=[-c0(i)ax(i)*x(i)+ay(i)*y(i)]';X2=A2\B2;x0(i)=X2(1);y0(i)=X2(2);%以下判断声线实际的弯曲方向，并利用几何知识求半径、出射点与底边中点距离、出射点坐标、出射方向矢量；ifabs(c1(i)-c2(i))==0r=0;d=0;x(i+1)=xh;y(i+1)=yh;ax(i+1)=ax(i);ay(i+1)=ay(i);elser=c(i)*t/abs((c1(i)-c2(i)));d=r-sqrt(r^2-h^2);end%以下两种情况需要注意符号，并且涉及到角度的正负和加减，还需用到和差公式：cos(𝜃1±𝜃2)或者sin(𝜃1±𝜃2)其中对于𝜃1，三角形入射的角度，cos(𝜃1)=𝑎𝑥√𝑎𝑥2+𝑎𝑦2；对于𝜃2，转过的角度，经推导：sin(𝜃2)=2𝑑ℎ(ℎ2+𝑑2)，cos(𝜃2)=ℎ2−𝑑2(ℎ2+𝑑2)；ifc1(i)c2(i)x(i+1)=xh+d*(-ay(i))/norm([-ay(i)ax(i)]);y(i+1)=yh+d*(ax(i))/norm([-ay(i)ax(i)]);ax(i+1)=(ax(i)*(h^2-d^2)-ay(i)*2*d*h)/(sqrt(ax(i)^2+ay(i)^2)*(h^2+d^2));ay(i+1)=(ay(i)*(h^2-d^2)+ax(i)*2*d*h)/(sqrt(ax(i)^2+ay(i)^2)*(h^2+d^2));endifc1(i)c2(i)x(i+1)=xh+d*(ay(i))/norm([ay(i)-ax(i)]);y(i+1)=yh+d*(-ax(i))/norm([ay(i)-ax(i)]);ax(i+1)=(ax(i)*(h^2-d^2)+ay(i)*2*d*h)/(sqrt(ax(i)^2+ay(i)^2)*(h^2+d^2));ay(i+1)=(ay(i)*(h^2-d^2)-ax(i)*2*d*h)/(sqrt(ax(i)^2+ay(i)^2)*(h^2+d^2));endi=i+1;ends=norm([xx_over(j)-x_overyy_over(j)-y_over]);%某段声线终点与指定点的距离；j=j+1;clearvars-exceptskejxx_overyy_overhtzthelt1x_inty_intttx_overy_over;end------------------------以下为打靶结束后，重复上述工作作出路径并求取传播路程和时间；ax(1)=1;ay(1)=k(j-1);x(1)=x_int;y(1)=y_int;D0=norm([x_int-x_overy_int-y_over]);%声线始终点直线距离；D1=0;time=0;i=1;%以下大部分程序与打靶法中的程序相同；whileabs(x(i))=3.785&&abs(y(i))=3.785thelt2(i)=atand(ay(i)/ax(i));ifax(i)0&&ay(i)0thelt2(i)=thelt2(i)+180;endifax(i)0&&ay(i)0thelt2(i)=thelt2(i)-180;endx1(i)=x(i)+sqrt(h^2+t^2/2)*cosd(thelt1+thelt2(i));y1(i)=y(i)+sqrt(h^2+t^2/2)*sind(thelt1+thelt2(i));x2(i)=x(i)+sqrt(h^2+t^2/2)*cosd(thelt2(i)-thelt1);y2(i)=y(i)+sqrt(h^2+t^2/2)*sind(thelt2(i)-thelt1);ifabs(x1(i))3.785||abs(x2(i))3.785||abs(y1(i))3.785||abs(y2(i))3.785ifabs(x(i))=abs(y(i))ify(i)0d1=abs(3.785-y(i));y(i+1)=3.785;x(i+1)=x(i)+d1/ay(i)*ax(i);elsed1=abs(3.785+y(i));y(i+1)=-3.785;x(i+1)=x(i)+d1/abs(ay(i))*ax(i);endelseifx(i)0d1=abs(3.785-x(i));x(i+1)=3.785;y(i+1)=y(i)+d1/ax(i)*ay(i);elsed1=abs(3.785+x(i));x(i+1)=-3.785;y(i