% 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]')