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 |
|
---|