水声物理中波数积分方法的MATLAB实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% 实验三 波数积分方法的编程实现 %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc; clear all; close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 设置参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H=100; H0=200; zs=36; f=20; zr=46; %设置接收点！ c1=1500; ro1=1000; c2=1800; ro2=1800; cmin=c1; k1=1.05*pi*f/c1; k2=2*pi*f/c2; krmax=0.9*2*pi*f/cmin; %对波数的最大值作了一定的调整！ krmin=10^(-8); Rmax=5000; %对距离的最大值作了一定的调整！ R=2*Rmax; M=2^round((log2(Rmax*krmax*2/pi)+0.5)); %kr和r的采样点数，亦即FFT的采样点数 phai=zeros(M); del_kr=(krmax-krmin)/(M-1); j=0:1:M-1; kr=krmin+del_kr.*j; del_r=2*pi/(M*del_kr); rmin=del_r; r=rmin+j.*del_r; u=3*(krmax-krmin)/(2*pi*(M-1)*log10(exp(1))); %计算复积分围线时的偏移单位 beita1=sqrt(k1^2-kr.^2); beita2=sqrt(k2^2-kr.^2); b=ro1/ro2; A = (2./beita1).*((beita1.*cos(beita1*(H-zs))-1i*b.*beita2.*sin(beita1*(H-zs)))./(beita1.*cos(beita1*H)-1i*b.*beita2.*sin(beita1*H))); B = 2*sin(beita1*zs)./beita1.*(beita1.*sin(beita1*H)+1i*b*beita2.*cos(beita1*H))./(beita1.*cos(beita1*H)-1i*b.*beita2.*sin(beita1*H)); C = 2.*sin(beita1*zs)./beita1; D = (2*b*sin(beita1*zs).*exp(-1i*beita2*H))./(beita1.*cos(beita1*H)-1i*b*beita2.*sin(beita1*H)); %%%%% 画出接收点zr深度处积分核函数随水平波数的变化曲线、声传播损失曲线 %%%% if (zr<=zs) phai=A.*sin(beita1.*zr); %积分核函数随水平波数变化 p=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*pi*0.25-i*krmin*r)*M.*ifft(phai.*sqrt(kr).*exp(-1i*rmin.*(0:M-1)*del_kr)); end if ((zr<=H)&&(zr>zs)) phai=B.*sin(beita1*zr)+C.*cos(beita1*zr); p=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*0.25-i*krmin*r)*M.*ifft(phai.*sqrt(kr).*exp(-1i*rmin.*(0:M-1)*del_kr)); end if (zr>H) phai=D.*exp(1i*beita2*zr); %积分核函数随水平波数变化 p=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*0.25-i*krmin*r)*M.*ifft(phai.*sqrt(kr).*exp(-1i*rmin.*(0:M-1)*del_kr)); end figure(1); plot(kr,abs(phai)); title('积分核函数随水平波数的变化曲线') xlabel('kr'); ylabel('Phai'); figure(2) plot(r,-20*log10(abs(p))); axis ([r(1) Rmax 0 100]); axis ij; title('声传播损失曲线') xlabel('Range(m)'); ylabel('TL(dB)'); %%%%%%%%%%%%%%%%%% 积分核函数伪彩图、绘制声压损失伪彩图 %%%%%%%%%%%%%%%%%%% for j=1:1:zs z1(j)=0+j-1; Phai1(j,:)=A.*sin(beita1*z1(j)); end for n=1:1:zs p(n,:)=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*pi*0.25-i*krmin*r)*M.*ifft(Phai1(n,:).*sqrt(kr).*exp(1i*rmin.*(0:M-1)*del_kr)); end for j=1:1:(H-zs) z2(j)=zs+j-1; Phai2(j,:)=B.*sin(beita1*z2(j))+C.*cos(beita1*z2(j)); end for n=1:1:(H-zs) p(n+zs,:)=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*pi*0.25-i*krmin*r)*M.*ifft(Phai2(n,:).*sqrt(kr).*exp(1i*rmin.*(0:M-1)*del_kr)); end for j=1:1:H0-H z3(j)=H+j-1; Phai3(j,:)=D.*exp(1i*beita2*z3(j)); end for n=1:1:H0-H p(n+H,:)=del_kr*sqrt(1./(2*pi.*r)).*exp(u*r+1i*pi*0.25-i*krmin*r)*M.*ifft(Phai3(n,:).*sqrt(kr).*exp(1i*rmin.*(0:M-1)*del_kr)); end z=[z1 z2 z3]; phai=[Phai1; Phai2; Phai3]; figure(3) he=pcolor(kr,z,20*log10(abs(phai))); set(he,'Linestyle','None'); %绘制积分核函数伪彩图 title('积分核函数'); xlabel('kr'); ylabel('Phai'); axis ([kr(1) kr(M) z(1) z(H0)]); axis ij; colormap(jet(250)); colorbar; figure(4) sheng=pcolor(r,z,-20*log10(abs(p))); set(sheng,'Linestyle','None'); %绘制传播损失伪彩图 title('传播损失分布'); xlabel('Range(m)'); ylabel('TL(dB)'); axis ([r(1) Rmax z(1) z(H0)]); axis ij; caxis([0 150]); colormap(jet(250)); colormap(flipud(colormap)); %反色 colorbar; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%