% shallow_water_rays.m
%
% traces the ray traveled by a shallow water wave packet in a field of 
% varying depth, H(y)
%
%   omega^2 = c^2 (k^2+l^2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

g = 9.8;                % gravity

ymax=300;               % y dimension of domain (m)
nt = 100;               % integration length (s)
dt = 1;                 % timestep (s)

k = 2*pi/10;            % 10m x wavelength
l = 2*pi/10;            % 10m initial y wavelength

x = nan(nt,1);
y = nan(nt,1);
x(1) = 0;              % initial x
y(1) = 0;              % initial y

H0 = 2;                 % depth scale
H = H0*(1-y(1)/ymax);   % depth at initial location of wave packet
c = sqrt(g*H);
omega = c*sqrt(k^2+l^2);    % omega is constant

for i = 1:nt
    
    H = H0*(1-y(i)/ymax);   % calc local depth
    c = sqrt(g*H);          % calc local c
    
    l = sqrt( omega^2 / c^2 - k^2 );
    
    cg_x = c*k / sqrt(k^2+l^2);
    cg_y = c*l / sqrt(k^2+l^2);
    
    x(i+1) = x(i) + cg_x*dt;    
    y(i+1) = y(i) + cg_y*dt;
    
end

yvec = linspace(0,300,100);
Hvec = H0*(1-yvec/ymax);

figure(1); 
subplot(1,3,1:2)
    plot(x,y)
    ylim([0,ymax])
    grid on;
    title('Ray path','FontSize',14)
    xlabel('x','FontSize',14); ylabel('y','FontSize',14)
subplot(1,3,3)
    plot(Hvec,yvec)
    ylim([0,ymax])
    grid on;
    title('Depth','FontSize',14)
    xlabel('H(y)','FontSize',14)
    