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