% Wave erosion of a beach with a rectangular extension
% of the planform of distance Y on (-L/2,L/2), zero 
% elsewhwere. The longshore diffusivity is G which is a rather
% complicated function of the tidal currents and other parameters.
% Ref. 'Beach Nourishment: Theory and Practice' by Robert G. Dean
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
xoverL=linspace(0,1.5,91);
sn=[];
nterm=10;
ftime=1;
ntime=101;
tau=linspace(0,ftime,ntime);%tau=G*t/L^2, scaled time
tau(1,1)=.00001;
one=ones(1,91);
%T1=(ones(1,101)-xoverL);
framedata=[];
%framedata=[ones(1,31) zeros(1,60)];
for k=1:ntime
    den=4*sqrt(tau(1,k));
    V=0.5*(erf((one-2*xoverL)/den)+erf((one+2*xoverL)/den));
    framedata=[framedata;V];
end
framedata2=[fliplr(framedata(:,2:91)) framedata];
xoverL2=linspace(-1.5,1.5,181);

figure(1);clf;% Plot y(x,t) vs x for various t
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Beach Extension,y(x,t)/Y')
for k=1:4:101
    plot(xoverL2,framedata2(k,:),'k')
    hold on
    axis([-1.5 1.5 -.1 1.2]);
    hold on
    %pause
end
%pause
text(-.45,.05,['Gt/L^2 = ' num2str(0)])
%text(.16,.2,num2str(0.02))
text(-.2 ,.23,['Gt/L^2 = ' num2str(1)])
%text(1.04,.8,'T = 0')
%text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(-1.25,1.1,'Press Enter to Continue')
pause

figure(2);clf;%Plot y(x,t)/Y vs t for various x
%plot(tau,ones(1,101))
hold on
for j=1:12:85
    plot(tau,framedata(:,j))
    hold on
    axis([0 1 0 1]);
    hold on
end
box on
text(.2,.62,['x/L = ' num2str(0)])
text(.05,.5,['x/L = ' num2str(0.4)])
text(.05,.3,['x/L = ' num2str(.8)])
text(0.15,0.14,['x/L = ' num2str(1.2)])
xlabel('Dimensionless Time, Gt/L^2')
ylabel('Beach Extension, y(x,t)/Y')
text(.2,.9,'Press Enter to Continue')
pause

figure(3);clf;%now do animation
axis([-1.5 1.5 -.1  1.2])
hold on
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Beach Extension, y(x,t)/Y')
hold on
L=plot(xoverL2, framedata2(1,:),'k','EraseMode','xor');
hold on
texthandl=text(-1.25, 1.1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:101
    set(L,'Ydata',framedata2(i,:))
    pause(.05)
end
hold on
plot(xoverL2,framedata2(1,:),'k')
hold on
set(texthandl,'String','Press Enter to Continue');
pause

figure(4);clf;%Animate in color
axis([-1.5 1.5 -.1  1.2])
hold on
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Beach Extension, y(x,t)/Y')
hold on
xpatch=[xoverL2 3 -3 -3];
ypatch=[framedata2(1,:) 1.2 1.2 framedata2(1,1)]; 
L1=patch(xpatch, ypatch,'b','EraseMode','background');
hold on
box on
texthandl=text(-1.25, 1.1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:101
    ypatch=[framedata2(i,:) 1.2 1.2 framedata2(i,1)];
    set(L1,'Ydata',ypatch)
    pause(.05)
end
hold on
%plot(xoverL2,framedata2(1,:),'k')
%hold on
set(texthandl,'String','Press Enter to Continue')
pause
    
figure(5);clf%3-D Plot of T(x,t) vs x and t
%tau1=linspace(0,.1,51);
[X,Y]=meshgrid(tau,(xoverL2));
mesh(X,Y,framedata2')
colormap(cool)
text(1,-.6,-.5,'Distance, x/L','Rotation',12)
text(.1, -1.5,-.5,'Time, Gt/L^2','Rotation',-32)
text(0.1,1.4,.8,'Press Enter to Continue','Rotation',-32)
zlabel('Beach Extension, y(x,t)/Y')
view(60,30)
axis([0 1 -1.5 1.5  -.2 1.2])