%  This Matlab program is used in 
%  MATH 365 - Computational Fluid Dynamics at 
%  James Madison University. This program was
%  written by James Sochacki of the 
%  Department of Mathematics at James Madison 
%  University as part of the NSF Grant - 
%  A Collaborative Computational Science Program 
%  between JMU and North Carolina Central University. 
  
%  Jim Sochacki
%  Department of Mathematics
%  James Madison University
%  Harrisonburg, VA   22807
%  jim@math.jmu.edu



% This code calculates the solution to the nonlinear burger's equation in
% conservation form using Lax-Wendroff. The initial condition is a front based
% on arctan(x). 

     
      disp(' This code computes the solution for the burger equation ')
      disp(' u_t + (0.5*u^2)_x = 0 ; -L < x < L, 0 < t < T '  )
      disp(' u(x,0) = q arctan(s*x)+r ')
      disp(' using Lax-Wendroff. ')
      disp('  ')
      disp(' Give the value for the initial condition at ')
      disp(' a_left and a_right. ')
    
      a_left = input(' a_left ');
      a_right = input(' a_right ');
      cmax=max(abs(a_left),abs(a_right));

      s = input(' Give s ');

      L = input(' Give  L. ');

    T = input(' Give the time T to run the model. ');
 
    dx = input(' Give the grid size. ');
    dt = input(' Give the time step. ');

    tsteps = T/dt;  % number of time steps

      disp(' Give an integer indicating how often a snapshot ')
      disp(' of the profile should be taken. ')
      interval=input(' ');
      
%
% Set up the parameters for the Lax-Wendroff scheme.
%

        jmax = (2*L)/dx;
        lambda1 = dt/dx;
        lambda2 = dt/(2.*dx);
        r = (a_left+a_right)/2;
        q = (a_right-a_left)/pi;

        for j=1:jmax
          x(j) = -L+j*dx;
          u0(j) = q*atan(s*x(j))+r;
          u1(j) = u0(j);
          u2(j) = u1(j);
        end      

        Mf = moviein(tsteps);
        Mp = moviein(tsteps);

        for j = 1:jmax
        for k = 1:2
          z(k,j) = u0(j);
        end
        end
        pcolor(z);
        Mf(:,1) = getframe;
        plot(x,u0);
        Mp(:,1) = getframe;
        pause(0.25);

%
% Start the calculations.
%
         count = 0;

         for i=1:tsteps

           count=count+1;
           
           for j=2:jmax-1

             

              uxf = (u0(j+1)+u0(j))/2;
              uxb = (u0(j)+u0(j-1))/2;
              u1(j) = u0(j)-(u0(j+1)^2-u0(j-1)^2)*lambda1+(uxf*(u0(j+1)^2-u0(j)^2)-uxb*(u0(j)^2-u0(j-1)^2))*lambda2;

           end
%
% Create the snapshots.

          if (count == interval)

            for j = 1:jmax
            for k = 1:2
              z(k,j) = u2(j);
            end
            end
        pcolor(z);
        Mf(:,i) = getframe;
        plot(x,u2);
        Mp(:,i) = getframe;     
            count = 0;

            plot(x,u1);
            pause(0.25);
            count = 0;
 
          end
%
%  Update the calculations in time
%
          for j=1:jmax-1
            u0(j)=u1(j);
          end

      end

movie(Mf,0)
movie(Mp,0)

        
