[5649] | 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 | |
---|