% To run, open matlab, so to this directory and enter cg

clear
f=figure;
slope=0.1;              %bed slope 
T=20;                   %period
l=slope*9.8*T^2/pi^2;   %characteristice length (lo in C&G58)
A=0.5;                  %amplitude
dd=0.1;                 %resolution of grid in the sigma-alpha plane

%define the potential phi according to the C&G solution (eqn 2.21 C&G58)
for sigma=[1:200];
for lambda=[1:200];
Phi(sigma,lambda)=A*besselj(0,(sigma-1)*dd)*cos((lambda-1)*dd);
end
end

%solve for x, t, u and eta (equations 2.16-2.19 C&G58)
for sigma=[2:199];
for lambda=[2:199];
u(sigma,lambda)=-(Phi(sigma+1,lambda)-Phi(sigma-1,lambda))/dd/2./(sigma*dd);
eta(sigma,lambda)=1/4.*(Phi(sigma,lambda+1)-Phi(sigma,lambda-1))/dd/2-1/2.*(u(sigma,lambda)).^2; % stage
t(sigma,lambda)=1/2.*(lambda*dd)+u(sigma,lambda);
x(sigma,lambda)=-1/4*(Phi(sigma,lambda+1)-Phi(sigma,lambda-1))/dd/2+1/16*(sigma*dd).^2+1/2*u(sigma, lambda).^2;
end
end

%plot using Matlab interpolation routines
figure
ti=-50:250;
tj=1:99;
[XI,YI] = meshgrid(ti/10,tj/10);
Z = griddata(x(:,:),t(:,:),eta(:,:),XI,YI);
mesh(XI,YI,Z)

%show animation using dimensional variables according to C&G58 scaling
figure(f)
a=0;
for tp=1:99
    plot(l/slope-l*XI(1,:),l-l*slope*XI(1,:),'k')
    hold on
    plot(l/slope-l*XI(1,:),l+l*slope*Z(tp,:),'b')
    text(l,l+l*slope/8,num2str(tp/10*T/pi))
    axis([0 1.1*l/slope l-l*slope/4 l+l*slope/4])
    pause(0.1)
    hold off
    a=a+1;
%     fname=[num2str(a) '.jpg'];
%     saveas(f,fname)
    F(a) = getframe;
end
movie(F,5)
% figure
% plot3(x, t, eta,'.')