%                     Lab Assignment #1
%  
%   Program the following Matlab code. The best way is to make a 
%   file (filename: lab1.m) of the commands as they are shown below.  
%   Startup matlab on your system and at the matlab prompt 
%   (typically >) type: lab1
%   The program should start up and prompt you for input.
%
%   The system you are solving is 
%                n+1    n          n       n      n
%               u    = u  - CFL(a u   + b u  + c u   )
%                j      j          j-1     j      j+1
%   with periodic boundary conditions
%
%   It is the simple wave equation on a periodic domain.  We have chosen a 
%   three point differencing in space and explicit Euler time differencing.
%             
%   Problem #1:  For nmax = 10, CFL = 1.0, run the program with
%      a)  a = -0.5, b = 0.0, c = 0.5 and describe what happens.
%          Short descriptions please.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) for a = -1, b = 1, c = 0
%      c)  Repeat a) for CFL = 0.7
%      d)  Repeat a) for CFL = 1.4
%
%   Problem #2:  For nmax = 10, a = -1, b = 1, c = 0  run the program with
%      a)  CFL = 0.7 and describe what happens.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) CFL = 1.4
%        
%   Problem #3: 
%          For nmax = 10, CFL=1.0, run the program with
%      a)  a = -0.9, b = 0.8, c = 0.1 and describe what happens.
%          Short descriptions please.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) for a = -0.99, b = 0.98, c = 0.01

help lab1;

%% Start of program ::::::::::::::::::::::::::::::::::::::::::::

jmax = 51;       

% The grid (simply the points 1 to jmax)
x = [1:jmax];

%  Initial Conditions
u  = zeros(jmax,1);   
un = zeros(jmax,1);   
%  set spike at middle
u((jmax-1)/2) = 1.0;

%  Input
nmax = 10; 
nmax = input('Enter the number of steps to take :: ');

CFL = 1.0;  
CFL = input('Enter CFL :: ');

a = -0.5;  
b = 0.; 
c = 0.5;
fprintf('Enter a,b,c of B(a,b,c) \n');
a = input('Enter of  a :: ');
b = input('Enter of  b :: ');
c = input('Enter of  c :: ');

%  iterate nmax times 
for iter = 1:nmax 
  
  % interior scheme
  %  The points we are updating (integrating) are j = 2 to j = jmax-1
  %  The boundary points j=1 and j=jmax are overloaded 
  %  with j=jmax-1 and j=2 respectively to simulate periodicty, i.e. to
  %  elliminate end effects

  j=[2:jmax-1]; % the interior index array

  %  load new values into array un
  un(j) = u(j) - CFL * (a*u(j-1) + b*u(j) + c*u(j+1));

  % Note the above two lines are equivalent to
  %    for j=2:jmax-1
  %       un(j) = u(j) - CFL * (a*u(j-1) + b*u(j) + c*u(j+1));
  %    end
  % But, using the vector form is faster (much!)
  
  % boundary scheme enforces periodicity
  un(1)    = un(jmax-1);  
  un(jmax) = un(2);
  
  % update u
  u        = un;

  %  If you have plotting capability use the next statement 
  
  if iter==5 
      clf;
      subplot(1,2,1)
      plot(x,u);
      title(strcat('iter=',num2str(iter)));
      xlabel('x');ylabel('u');
      axis([1 jmax min(u) max(max(u),1)]);
      %axis_ob=get(gca);
      set(gca,'Position',[0.0875,0.11,0.4,0.7]);
  end
  if iter==10
      subplot(1,2,2)
      plot(x,u);
      title(strcat('iter=',num2str(iter)));
      xlabel('x');ylabel('u');
      axis([1 jmax min(u) max(max(u),1)]);
      text_ob=text(0,0,strcat('Solution to wave equation with Nmax=',num2str(nmax),...
          ' CFL=',num2str(CFL),' a=',num2str(a),' b=',num2str(b),' c=', num2str(c)))
      set(text_ob,'units','normal');
      set(text_ob,'position',[-.2,1.2,0]);
      set(text_ob,'horizontalalignment','center');
      set(gca,'Position',[0.5785,0.11,0.4,0.7]);
      fig_ob=get(gcf); 
      
  end
  %  Note!! the scales on the plot may be changing dynamically
  
  %  The next statement outputs data to the screen, the lack 
  %  of a trailing ";" is correct, causing the output.
  u
  %  The next statement slows down the output.  You can increase the
  %  pause time by increasing the (0.1).
  
end