% SIMPLE script to do illustrate sound channel % These parameters can be changed w=47.4; % acoustic frequency (radian/sec) a w=[47:.1:47.8]; lambda=200; % horizontal wavelength (m) zini=-1200; % depth at which to place the sound source dt=.01; % time step (s) T=100; % total integration time (s) % Grid dz=5; z=[0:dz:4000]'; %set up sound speed profile, this particular one % is the Munk profile but any will work. B=1200; z0=1200.0; c0=1492.0; ep=0.006; eta=2*(z-z0)/B; c=c0*(1+ep*(eta+exp(-eta) -1)); clear B z0 c0 ep eta z=-z; z=flipud(z); c=flipud(c); dcdz=gradient(c,dz); iz0=find(z==zini); nt=T/dt; k=2*pi./lambda; % horizontal wavenumber (radian/m) nw=length(w); xw=repmat(0,[nt nw]); zw=repmat(0,[nt nw]); kw=repmat(0,[nt nw]); mw=repmat(0,[nt nw]); zw(1,:)=z(iz0); kw(1,:)=k; cz=interp1(z,c,z(iz0)); % local sound speed m=sqrt((w.^2-(cz^2)*(k.^2))/(cz^2)); % local vertical wave number if any(~isreal(m)) error('m is not real: w must be >= c*k') end mw(1,:)=m; for i=1:nt-1 cz=interp1(z,c,zw(i,:)); % local speed dc=interp1(z,dcdz,zw(i,:)); % local gradient of c kappa=sqrt(k.^2+m.^2); cgx=cz.*k./kappa; cgz=cz.*m./kappa; DkDt=0; DmDt=-dc.*kappa; % integrate ray equations with a simple Euler step k=k+DkDt*dt; m=m+DmDt*dt; xw(i+1,:)=xw(i,:)+cgx*dt; zw(i+1,:)=zw(i,:)+cgz*dt; kw(i+1,:)=k; mw(i+1,:)=m; end figure subplot(1,2,1) plot(c,z) hold on %plot(w/kw(1,1),z,'r') grid xlabel('Sound speed [m/s]') ylabel('z [m]') subplot(1,2,2) plot(xw/1e3,zw) set(gca,'ylim',[min(z) max(z)]) grid title('Ray path') xlabel('horizontal distance [km]')