% Staggered Leapfrog Script Template
clear;close all;


% build a cell-centered grid with N=200 cells
% on the interval x=0 to x=L, with L=1
L = 1;
sigma = 30;
N = 2000;
xmin = 0; xmax = L; h = (xmax-xmin)/N;
x = xmin-h/2:h:xmax+h/2;
%x=x';
y=x; yold=x; ynew=x;
j = 2:N+1;

% Set the values for parameters
%c=2; % wave speed
%c=2; % wave speed
c(1:(N+2)/2)= 2;
c( ((N+2)/2 + 1):(N+2)) = 4;

% define the initial displacement and velocity vs. x
%y = exp(-(x-L/2).^2*160/L^2)-exp(-(0-L/2).^2*160/L^2);
%y = exp(-(x-L/2).^2*900/L^2); %-exp(-(0-L/2).^2/(L/sigma)^2);
y = exp(-(x-L/4).^2*900/L^2); 
%vy = zeros(1,length(x)); %vy = vy';

% Choose a time step (Suggest that it is no larger than taulim)
taulim=h/max(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) - vy(j)*tau + c^2*tau^2/(2*h^2) * (y(j+1) - 2*y(j) + y(j-1));
%yold = exp(-((x+c*tau)-L/2).^2/(L/sigma)^2)-exp(-((0+c*tau)-L/2).^2/(L/sigma)^2);
yold = exp(-((x+c(1)*tau)-L/4).^2*900/L^2); 

% Apply the boundary conditions for yold(1) and yold(N+2)
yold(1) = -yold(2) + 2*0;
yold(N+2) = -yold(N+1) + 2*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)
plot(x,yold)
%xlabel('x');ylabel('v_y(x,0)');title('Initial Velocity')
xlabel('x');ylabel('yold');title('Older Displacement')
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

figure  % open a new frame for the animation
for n=1:nsteps
   time=n*tau;  % compute the time

   % 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) = 2*y(j) - yold(j) + c(j).^2*tau^2/h^2  .*(y(j+1) - 2*y(j) + y(j-1));
   ynew(1) = -ynew(2) + 2*0;
   ynew(N+2) = -ynew(N+1) + 2*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)])
      axis([min(x) max(x) -1 1]);
      pause(.1)
   end
end
