%                     Lab Assignment #1
%  
%   Program the following Matlab code. The best way is to make a 
%   file (filename: examp1.m) of the commands as they are shown below.  
%   Startup matlab on your system and at the matlab prompt 
%   (typically >) type: examp1
%   The program should start up and prompt you for input.
%
%   The system you are solving is describe below in the fprintf statements.  
%   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
%
%   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: (Optional) 
%          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


%% Start of program.........................................
      fprintf('u + u  = 0, Periodic BC\n');
      fprintf(' t   x  \n');
      fprintf('\n');
      fprintf('Solved by \n');
      fprintf('\n');
      fprintf(' n+1    n          n       n      n\n');
      fprintf('u    = u  - CFL(a u   + b u  + c u   )\n');
      fprintf(' j      j          j-1     j      j+1 \n');

% $$$       jmax = 51;       
% $$$       jmax = input('Enter the number points D= ');
      m5 = findobj(gcf,'Tag','jmax_in');
      jmax = str2num(get(m5,'String'));

      x = [1:jmax];

% $$$       ic = 0;       
% $$$       ic = input('Enter 0 for spike ic, 1 for jump ic D= ');
m0 = findobj(gcf,'Tag','ic_in');
ic = get(m0,'Value')-1;

m6 = findobj(gcf,'Tag','meth_in');
meth = get(m6,'Value')-1;

%  Initial Condition
      u  = zeros(jmax,1);   
      un = zeros(jmax,1);   

      js = (jmax-1)/2;
      jl = js - 5;
      jr = js + 5;

      if ic==0
        u(js) = 1.0;
      elseif ic == 1
	j = [jl:jr];
	u(j) = 1.0;
      elseif ic == 2
	j=[2:jmax-1];
	dpi = 2*pi/(jmax-2);
	u(j) = sin((j-2)*dpi);
	u(1) = u(jmax-1);
      end

      if ic==1
      end 

      v = u;
      un = u;
      unm1 = u;
%  Input
      %  Input
      m1 = findobj(gcf,'Tag','cfl_in');
      CFL = str2num(get(m1,'String'));
      m2 = findobj(gcf,'Tag','n_in');
      nmax = str2num(get(m2,'String'));
      m3 = findobj(gcf,'Tag','a_in');
      a = str2num(get(m3,'String'));
      m4 = findobj(gcf,'Tag','b_in');
      b = str2num(get(m4,'String'));
      m5 = findobj(gcf,'Tag','c_in');
      c = str2num(get(m5,'String'));

      for iter = 1:nmax 
	j=[2:jmax-1];
	if meth == 0 | iter == 1
	  un(j) = u(j) - CFL * (a*u(j-1) + b*u(j) + c*u(j+1));
	elseif meth == 1 & iter > 1
	  un(j) = unm1(j) - 2*CFL * (a*u(j-1) + b*u(j) + c*u(j+1));
	end
	un(1)    = un(jmax-1);  
	un(jmax) = un(2);
	unm1     = u;
	u        = un;

	if ic == 0 
	  % exact shift of spike
	  js = js + 1;  jsm1 = js-1; if js == jmax js = 2; jsm1 = jmax-1; end
	  v(jsm1) = 0.0;
	  v(js) = 1.0;
	elseif ic == 1
	  jl = jl + 1; if jl == jmax jl = 2; end
	  jr = jr + 1; if jr == jmax jr = 2; end
	  jj = [jl:jr];
	  v([1:jmax]) = 0.0;
	  v(jj) = 1.0;
	elseif ic == 2
	  j=[2:jmax-1];
	  dpi = 2*pi/(jmax-2);
	  v(j) = sin((j-2-iter)*dpi);
	  v(1) = v(jmax-1);
	  v(jmax) = v(2);
	end
%  If you have plotting Capability use the next statement 
	plot(x,u,'b');
	hold on;
	plot(x,v,'r-.')
	hold off;
%  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 (1).
         pause(1);

      end