1 | % To run, open matlab, so to this directory and enter cg |
---|
2 | |
---|
3 | clear |
---|
4 | f=figure; |
---|
5 | slope=0.1; %bed slope |
---|
6 | T=20; %period |
---|
7 | l=slope*9.8*T^2/pi^2; %characteristice length (lo in C&G58) |
---|
8 | A=0.5; %amplitude |
---|
9 | dd=0.1; %resolution of grid in the sigma-alpha plane |
---|
10 | |
---|
11 | %define the potential phi according to the C&G solution (eqn 2.21 C&G58) |
---|
12 | for sigma=[1:200]; |
---|
13 | for lambda=[1:200]; |
---|
14 | Phi(sigma,lambda)=A*besselj(0,(sigma-1)*dd)*cos((lambda-1)*dd); |
---|
15 | end |
---|
16 | end |
---|
17 | |
---|
18 | %solve for x, t, u and eta (equations 2.16-2.19 C&G58) |
---|
19 | for sigma=[2:199]; |
---|
20 | for lambda=[2:199]; |
---|
21 | u(sigma,lambda)=-(Phi(sigma+1,lambda)-Phi(sigma-1,lambda))/dd/2./(sigma*dd); |
---|
22 | eta(sigma,lambda)=1/4.*(Phi(sigma,lambda+1)-Phi(sigma,lambda-1))/dd/2-1/2.*(u(sigma,lambda)).^2; % stage |
---|
23 | t(sigma,lambda)=1/2.*(lambda*dd)+u(sigma,lambda); |
---|
24 | 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; |
---|
25 | end |
---|
26 | end |
---|
27 | |
---|
28 | %plot using Matlab interpolation routines |
---|
29 | figure |
---|
30 | ti=-50:250; |
---|
31 | tj=1:99; |
---|
32 | [XI,YI] = meshgrid(ti/10,tj/10); |
---|
33 | Z = griddata(x(:,:),t(:,:),eta(:,:),XI,YI); |
---|
34 | mesh(XI,YI,Z) |
---|
35 | |
---|
36 | %show animation using dimensional variables according to C&G58 scaling |
---|
37 | figure(f) |
---|
38 | a=0; |
---|
39 | for tp=1:99 |
---|
40 | plot(l/slope-l*XI(1,:),l-l*slope*XI(1,:),'k') |
---|
41 | hold on |
---|
42 | plot(l/slope-l*XI(1,:),l+l*slope*Z(tp,:),'b') |
---|
43 | text(l,l+l*slope/8,num2str(tp/10*T/pi)) |
---|
44 | axis([0 1.1*l/slope l-l*slope/4 l+l*slope/4]) |
---|
45 | pause(0.1) |
---|
46 | hold off |
---|
47 | a=a+1; |
---|
48 | % fname=[num2str(a) '.jpg']; |
---|
49 | % saveas(f,fname) |
---|
50 | F(a) = getframe; |
---|
51 | end |
---|
52 | movie(F,5) |
---|
53 | % figure |
---|
54 | % plot3(x, t, eta,'.') |
---|
55 | |
---|