% Staggered Leapfrog  with damping
% Problem 5.4B

clear;close all;

% Set the values for parameters
c=2; % wave speed
gamma = 0.2; %damping

% build a grid on the interval x=0 to x=L, with L=1
L = 1;
N = 500;
xmin = 0; xmax = L; h = (xmax-xmin)/N;  %cell-centered
x = xmin-h/2:h:xmax+h/2;    %cell-centered



j = 2:N+1;  % trick from manual, index to cover interior points

% define the initial displacement and velocity vs. x
vy = exp(-(x-L/2).^2*160/L^2)-exp(-(0-L/2).^2*160/L^2);
y = x*0;  % zero vector of proper length
yold=y; ynew=y;  % just to ensure same dimension of all arrays; they'll be set later

% Choose a time step (Suggest that it is no larger than taulim)
taulim=h/c;
fprintf(' Courant time step limit %g \n',taulim)
tau=input(' Enter the time step - ')

% Get the initial value of yold from the initial y and vy
yold(j) = y(j) - (2 + gamma*tau)* vy(j)*tau/2 +    c^2*tau^2/(2*h^2) * (y(j+1) - 2*y(j) + y(j-1));

% Apply the boundary conditions for yold(1) and yold(N+2)
yold(1) = -yold(2) + 2*0;     %Dirichlet, bndy value = 0
yold(N+2) = -yold(N+1) +2*0;  %Dirichlet, bndy value = 0

% plot the initial conditions and pause to look at them
subplot(2,1,1)
plot(x,y)
xlabel('x');ylabel('y(x,0)');title('Initial Displacement')
subplot(2,1,2)
plot(x,vy)
xlabel('x');ylabel('v_y(x,0)');title('Initial Velocity')
pause;

% Choose how long to run and when to plot
tfinal=input(' Enter tfinal - ')
skip=input(' Enter # of steps to skip between plots (faster) - ')
nsteps=tfinal/tau;

% here is the loop that steps the solution along

time = zeros(1,nsteps);   %pre-allocate memory for time array
maxamp = zeros(1,nsteps); %pre-allocate memory for maximum amplitude graph
figure  % open a new frame for the animation
for n=1:nsteps
   time(n)=n*tau;  % compute the time. Index added so we can use time array after loop ends

   % Use leapfrog and the boundary conditions to load ynew with y
   % at the next time step using y and yold, i.e., ynew(2:N+1)=...
   % Be sure to use colon commands so it will run fast.
   ynew(j) = 1/(2 + gamma*tau) * (4*y(j) - 2*yold(j) + gamma*tau*yold(j)   + 2*c^2*tau^2/h^2 * (y(j+1) - 2*y(j) + y(j-1))); %Eqn 5.17
   
   ynew(1) = -ynew(2) + 2*0;     %Dirichlet, bndy value = 0
   ynew(N+2) = -ynew(N+1) + 2*0; %Dirichlet, bndy value = 0

   %update yold and y
   yold=y;y=ynew;

% make plots every skip time steps
   if mod(n,skip)==0
      plot(x,y,'b-')
      xlabel('x');ylabel('y');
      title(['Staggered Leapfrog Wave: time=' num2str(time(n))])
      axis([min(x) max(x) -.05 .05]);
      pause(.1)
   end
   
   maxamp(n) = max(abs(y)); %to make max amplitude plot
end

figure;
plot(time,maxamp)   %plotting maximum amplitude vs. time
hold on
plot(time,max(maxamp)*exp(-gamma*time/2),'r-') %plotting expected decay in red, on same graph
title('Max Amplitude vs. Time')
