clear all;close all;
tic
% z0=680;%入射深度
dz=1;%深度分辨率
w=800;%初始化角频率
z1=[0.0;150.0;305.0;533.0;610.0;680.0;762.0;1372.0;1829.0;3048.0;4000.0];%已知点
z=0:dz:4000;%按深度分辨率生成深度矩阵
c1=[1507.2;1498.1;1491.7;1480.7;1478.9;1478.0;1478.6;1483.2;1488.6;1507.5;1523.0];
%已知点速度
c= interp1(z1,c1,z,'linear');
k=w./c;%波束
% figure(1)
% plot(c,z);title('声速垂直分布');set(gca,'ydir','reverse');set(gca,'XAxisLocation','top');
[~,low]=min(c);%找声速最小值 low为声速最小值位置
z0=z(low);
kmax=k(low);
kmin=k(1);%若比k(1)还小 就不会在海面以下翻转
c0=c(z0/dz+1);%初始声速
ksi0=kmin;%恰好在上表面翻转的 是表面cosa=1 求得极限情况下的水平波束;
h1=1;
h2=find(k(low:end)<=ksi0,1,'first')+low-1;
kz=sqrt(k(h1:h2).^2-ksi0^2);
fl=2*sum(dz.*kz);%直接全范围一阶数值积分
N=1;%至少一阶 不考虑0阶
while fl>(-pi+2*(N+1)*pi)
N=N+1;
end%N不行的时候退出 即N+1 正好对应N+1个
ksi=zeros(1,N);
for n=1:N
ma=kmax;mi=kmin;e=(ma-mi)/2;
while e>0.000001 %二分法
ksi(n)=(ma+mi)/2;
h1=find(k(1:low)<=ksi(n),1,'last');
h2=find(k(low:end)<=ksi(n),1,'first')+low-1;%最靠近发射位置的两个
kz=sqrt(k(h1:h2).^2-ksi(n)^2);
fl=2*sum(dz.*kz);
if fl<(2*n-1)*pi
ma=ksi(n);
e=(ma-mi)/2;
elseif fl>(2*n-1)*pi
mi=ksi(n);
e=(ma-mi)/2;
else%考虑等于0 几率极小 可去掉可简化
e=0;
break
end
end
% sprintf('ξ(%d)=%f',n,ksi(n))
end
sprintf("新耗时:%f",toc)