clear all;close all;clc;
%%%--------------------------------------------------------------------------
%硬质均匀浅海声场的波束传播 （海面绝对软，海底绝对硬）
%参考资料：水声学原理
f=5;
j=(-1).^(1/2);
w=2*pi*f; %发射频率
z0=15; %声源深度
c0=1500; %海水中的声速
k=w/c0; %波数
H=1000; %海深度 N=round(H*w/(pi*c0)+1/2); %最大简正波阶次
f_min=(N-1/2)*c0/(2*H); %最小简正频率
r=5000;z=35; %测试点距离声原的长度以及距水面的深度
fs=300;
T=1;t=1/fs:1/fs:T;
%------------------固定r=500;z=35接收的各号简正波----------------------------
for n=1:N An=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
Kzn=(k.^2-An.^2).^(1/2); %本征值，波数k0的竖直分量
xishu=(-j*2/H*(2*pi/(An*r)).^(1/2))*sin(Kzn*z0)*exp(-j*(An*r-pi/4));%中间系数设置
p(n,:)=xishu*sin(Kzn*z)*exp(j*w*t); %各号简正波信号传播
plot(t,p(n,:));hold on; %简正波的时域波形
end
legend('各号简正波的时域波形');
%-----------------群速度设置------------------------------------------------
%各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小，越晚到达同一点
%--------------------------------------------------------------------------
A1=(k.^2-((1-1/2)*pi/H).^2).^(1/2);
theta_1=asin(A1/k);
Cg(1)=c0*A1/k; %阶数1的群速度
tao(1)=0;
A2=(k.^2-((2-1/2)*pi/H).^2).^(1/2);
theta_2=asin(A2/k);
Cg(2)=c0*A2/k; %阶数2的群速度
tao(2)=round(fs*abs(r/Cg(1)-r/Cg(2))); %阶数1,2的时间延迟点数
p2=[p(2,:)，zeros(1,tao(2))]; %阶数1的信号波形，因为考虑到矩阵维数相同才可相加
%--------------------------------------------------------------------------
%我其实对接下来这个循环抱有疑问，首先，不同阶数的简正波，阶数越大，时延点数越多,tao=[0,23,76,172,353,767,3498]
%p(n,:)原本是一个7*300的矩阵，代表7个阶数的简正波的传播。后面补零是合理的，因为理论上是n阶简正波，但是我们取了7阶，所以后面是零，补零是合适的
%自己这几天觉得自己好辣鸡，难过。jpg
for n=2:N
An=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
theta_n=asin(An/k); %各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小
Cg(n)=c0*An/k; %n阶群速度
tao(n)=round(fs*abs(r/Cg(1)-r/Cg(n))); %n阶波与n-1阶波的时延差
pn=[zeros(1,tao(n)),p(n,:)];
p2=[p2,zeros(1,length(pn)-length(p2))]+pn; %矩阵相加维数须相同
end
figure;plot(p2);title('各号简正波的时域波形到达同一点延时相加频率相位图');
figure;plot(1:3647,p2);title('各号简正波的时域波形到达同一点延时相加频率相位图');%把上图的x轴设为1：3647，得到不一样的图
figure;plot(abs(p2));title('各号简正波的时域波形到达同一点延时相加频率相位图');
自己按照这个程序的思路，编写了亮点模型的程序，将亮点看做两个点源，将两个点源的简正波时延后，相同阶数相加，得到结果，不知道对不对。
%水声传播简正波计算程序 %设置两个亮点，一个在（0，15），一个在（10，15），代表一个目标的两个亮点
clear all;close all;clc;
%%%--------------------------------------------------------------------------
%硬质均匀浅海声场的波束传播（海面绝对软海底绝对硬）
f=5;
w=2*pi*f; %发射频率
z0=15; %深度
c0=1500; %海水中的声速
k=w/c0; %波数
H=1000; %海深度
N=round(H*w/(pi*c0)+1/2); %最大简正波阶次
f_min=(N-1/2)*c0/(2*H); %最小简正频率
fs=300;T=1;
t=1/fs:1/fs:T;
%------------------1号亮点（0，15），固定r=5000;z=35接收的各号简正波----------------------------
figure(1); z01=15; %1号亮点在（r1,z01）=(0,15)
r1=5000;z1=35; %信号接收点距离声原的长度以及距水面的深度
for n=1:N
An=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
Kzn=(k.^2-An.^2).^(1/2); %本征值，波数k0的竖直分量
xishu1=(-j*2/H*(2*pi/(An*r1)).^(1/2))*sin(Kzn*z01)*exp(-j*(An*r1-pi/4));%中间系数设置
p1(n,:)=xishu1*sin(Kzn*z1)*exp(j*w*t); %各号简正波信号传播
plot(t,p1(n,:));hold on; %简正波的时域波形
end
legend('1号亮点（0，15），简正波阶数N=7，各号简正波的时域波形');
%------------------2号亮点（10，15），固定r=5000;z=35接收的各号简正波----------------------------
figure(2);
z02=15; %2号亮点在（r1,z01）=(0,15)
r2=4990;z2=35; %信号接收点距离声原的长度以及距水面的深度
for n=1:N
An=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
Kzn=(k.^2-An.^2).^(1/2); %本征值，波数k0的竖直分量
xishu2=(-j*2/H*(2*pi/(An*r2)).^(1/2))*sin(Kzn*z02)*exp(-j*(An*r2-pi/4));%中间系数设置
p2(n,:)=xishu2*sin(Kzn*z2)*exp(j*w*t); %各号简正波信号传播
plot(t,p2(n,:));hold on; %简正波的时域波形
end
legend('2号亮点（10，15），简正波阶数N=7，各号简正波的时域波形');
%%--求同一时刻的在信号接收点接受到的信号---------------------------------------
%对1号亮点进行处理----------------------------
%群速度设置：各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小，越晚到达同一点 %--------------------------------------------------------------------------
for n=1:N
A1(n)=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
theta1_(n)=asin(A1(n)/k); %各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小
Cg1(n)=c0*A1(n)/k; %n阶群速度
tao1(n)=round(fs*abs(r1/Cg1(1)-r1/Cg1(n))); %n阶波与n-1阶波的时延差
end
for n=1:N
pp1(n,:)=[zeros(1,tao1(n)) p1(n,:)
zeros(1,tao1(7)-tao1(n))];
end
%对2号亮点进行处理----------------------------
%群速度设置：各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小，越晚到达同一点 %--------------------------------------------------------------------------
for n=1:N
A2(n)=(k.^2-((n-1/2)*pi/H).^2).^(1/2); %波数k0的水平分量
theta2_(n)=asin(A2(n)/k); %各阶简正波的方向，波的阶数越高，则偏离z轴的角度越小，则群速度越小
Cg2(n)=c0*A2(n)/k; %n阶群速度
tao2(n)=round(fs*abs(r2/Cg2(1)-r2/Cg2(n))); %n阶波与n-1阶波的时延差
end
for n=1:N
pp2(n,:)=[zeros(1,tao2(n)) p2(n,:) zeros(1,tao2(7)-tao2(n))];
end
%-----------------相同阶数的简正波加在一起。
figure;
for n=1:N
pp_2(n,:)=[pp2(n,:)
zeros(1,length(pp1)-length(pp2))];
ppp(n,:)=pp1(n,:)+pp_2(n,:);
plot(1:3798,ppp(n,:));hold on; %简正波的时域波形
end
legend('两个亮点，简正波阶数N=7，各号简正波的时域波形'）；
虽然不知道对不对，可喜可贺，忙活了两天，终于懂了一些简正波的知识，去年学水声学原理的时候编写海洋分层的程序，简直要吐血，现在觉得还是很简单的，加油吧。
