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)