function zzz=main

%  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 leap frog approximations for the
% density, velocity and energy to the one dimensional Euler equations of
% gas dynamics. 

      disp(' This code computes the density, velocity and energy for ')
      disp(' the one dimensional nonlinear Euler equations using ')
      disp(' leap frog. ')

      disp(' Give the gas constant, gamma ')
      gamma = input(' ');
      
      disp(' Give the limit values for the initial velocity u(x,0) : u_left and u_right.  ')
      u_left = input(' u_left ');
      u_right = input(' u_right ');
      disp(' Give the slope value s for u(x,0) between u_left and u_right.  ')
      s_u = input(' s_u ');
      disp(' Give the limit values for the initial density rho(x,0) : rho_left and rho_right.  ')
      rho_left = input(' rho_left ');
      rho_right = input(' rho_right ');
      disp(' Give the slope value s for rho(x,0) between rho_left and rho_right. ')
      s_rho = input(' s_rho ');
      disp(' Give the limit values for the initial energy E(x,0) : E_left and E_right.  ')
      E_left = input(' E_left ');
      E_right = input(' E_right ');
      disp(' Give the slope value s for E(x,0) between E_left and E_right.  ')
      s_E = input(' s_E ');
      disp(' ')

      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(' Type a 0 for density output, a 1 for velocity output or ')
      disp(' a 2 or larger for energy output to be outputted.')
      output_type = input(' ');
      interval = input(' Give an integer indicating how often a frame should be made ( i.e., 2 means every other time step). ');
      
%
% Set up the parameters for the leap frog scheme.
%

        jmax = (2*L)/dx;

% jmax is the number of grid points in the model
%
% Set up the initial values for the leap frog scheme.
%

        for j=1:jmax

          x(j) = -L+(j-1)*dx;
          rho0(j) = p(rho_left,rho_right,s_rho,x(j));
          m0(j) = p(rho_left,rho_right,s_rho,x(j))*p(u_left,u_right,s_u,x(j));
          E0(j) = p(E_left,E_right,s_E,x(j));

        end
 
        E1(1) = E0(1);
        E1(jmax) = E0(jmax);
        rho1(1) = rho0(1);
        rho1(jmax) = rho0(jmax);
        m1(1) = m0(1);
        m1(jmax) = m0(jmax); 

        ymax = rho_left;
        if (rho_left < rho_right) 

          ymax = rho_right;

        end

        plot(x,rho0)
        axis([x(1) x(jmax) -0.5 ymax]);grid;

%
% Start the calculations.
% tsteps is the number of time steps
%

        count = 0;

        for i=1:tsteps

          count=count+1;
          if (count==interval) 

            disp(' step = ')
            i

          end

          for j=2:jmax-1
 
m1(j) = m0(j)-1/2.*dt*(-(gamma-1)*E0(j-1)-(3/2.-1/2.*gamma)*m0(j-1)^2./rho0(j-1)+(gamma-1)*E0(j+1)+(3/2.-1/2.*gamma)*m0(j+1)^2./rho0(j+1))/dx+1/4.*dt^2.*((2.*(3/2.-1/2.*gamma)*m0(j+1)/rho0(j+1)+2.*(3/2.-1/2.*gamma)*m0(j)/rho0(j))*(-(gamma-1)*E0(j)-(3/2.-1/2.*gamma)*m0(j)^2./rho0(j)+(gamma-1)*E0(j+1)+(3/2.-1/2.*gamma)*m0(j+1)^2./rho0(j+1))+(-(3/2.-1/2.*gamma)*m0(j+1)^2./(rho0(j+1)^2.)-(3/2.-1/2.*gamma)*m0(j)^2./(rho0(j)^2.))*(-m0(j)+m0(j+1))+(2.*gamma-2.)*(-m0(j)*gamma*E0(j)/rho0(j)-1/2.*(1-gamma)*m0(j)^3/(rho0(j)^2.)+m0(j+1)*gamma*E0(j+1)/rho0(j+1)+1/2.*(1-gamma)*m0(j+1)^3/(rho0(j+1)^2.))-(2.*(3/2.-1/2.*gamma)*m0(j)/rho0(j)+2.*(3/2.-1/2.*gamma)*m0(j-1)/rho0(j-1))*(-(gamma-1)*E0(j-1)-(3/2.-1/2.*gamma)*m0(j-1)^2./rho0(j-1)+(gamma-1)*E0(j)+(3/2.-1/2.*gamma)*m0(j)^2./rho0(j))-(-(3/2.-1/2.*gamma)*m0(j)^2./(rho0(j)^2.)-(3/2.-1/2.*gamma)*m0(j-1)^2./(rho0(j-1)^2.))*(-m0(j-1)+m0(j))-(2.*gamma-2.)*(-m0(j-1)*gamma*E0(j-1)/rho0(j-1)-1/2.*(1-gamma)*m0(j-1)^3/(rho0(j-1)^2.)+m0(j)*gamma*E0(j)/rho0(j)+1/2.*(1-gamma)*m0(j)^3/(rho0(j)^2.)))/(dx^2.);

rho1(j) = rho0(j)-1/2.*dt*(-m0(j-1)+m0(j+1))/dx+1/4.*dt^2.*(-4*(gamma-1)*E0(j)-4*(3/2.-1/2.*gamma)*m0(j)^2./rho0(j)+2.*(gamma-1)*E0(j+1)+2.*(3/2.-1/2.*gamma)*m0(j+1)^2./rho0(j+1)+2.*(gamma-1)*E0(j-1)+2.*(3/2.-1/2.*gamma)*m0(j-1)^2./rho0(j-1))/(dx^2.);
 
E1(j) = E0(j)-1/2.*dt*(-m0(j-1)*gamma*E0(j-1)/rho0(j-1)-1/2.*(1-gamma)*m0(j-1)^3/(rho0(j-1)^2.)+m0(j+1)*gamma*E0(j+1)/rho0(j+1)+1/2.*(1-gamma)*m0(j+1)^3/(rho0(j+1)^2.))/dx+1/4.*dt^2.*((gamma*E0(j+1)/rho0(j+1)+3/2.*(1-gamma)*m0(j+1)^2./(rho0(j+1)^2.)+gamma*E0(j)/rho0(j)+3/2.*(1-gamma)*m0(j)^2./(rho0(j)^2.))*(-(gamma-1)*E0(j)-(3/2.-1/2.*gamma)*m0(j)^2./rho0(j)+(gamma-1)*E0(j+1)+(3/2.-1/2.*gamma)*m0(j+1)^2./rho0(j+1))+(-m0(j+1)*gamma*E0(j+1)/(rho0(j+1)^2.)-(1-gamma)*m0(j+1)^3/(rho0(j+1)^3)-m0(j)*gamma*E0(j)/(rho0(j)^2.)-(1-gamma)*m0(j)^3/(rho0(j)^3))*(-m0(j)+m0(j+1))+(m0(j+1)*gamma/rho0(j+1)+m0(j)*gamma/rho0(j))*(-m0(j)*gamma*E0(j)/rho0(j)-1/2.*(1-gamma)*m0(j)^3/(rho0(j)^2.)+m0(j+1)*gamma*E0(j+1)/rho0(j+1)+1/2.*(1-gamma)*m0(j+1)^3/(rho0(j+1)^2.))-(gamma*E0(j)/rho0(j)+3/2.*(1-gamma)*m0(j)^2./(rho0(j)^2.)+gamma*E0(j-1)/rho0(j-1)+3/2.*(1-gamma)*m0(j-1)^2./(rho0(j-1)^2.))*(-(gamma-1)*E0(j-1)-(3/2.-1/2.*gamma)*m0(j-1)^2./rho0(j-1)+(gamma-1)*E0(j)+(3/2.-1/2.*gamma)*m0(j)^2./rho0(j))-(-m0(j)*gamma*E0(j)/(rho0(j)^2.)-(1-gamma)*m0(j)^3/(rho0(j)^3)-m0(j-1)*gamma*E0(j-1)/(rho0(j-1)^2.)-(1-gamma)*m0(j-1)^3/(rho0(j-1)^3))*(-m0(j-1)+m0(j))-(m0(j)*gamma/rho0(j)+m0(j-1)*gamma/rho0(j-1))*(-m0(j-1)*gamma*E0(j-1)/rho0(j-1)-1/2.*(1-gamma)*m0(j-1)^3/(rho0(j-1)^2.)+m0(j)*gamma*E0(j)/rho0(j)+1/2.*(1-gamma)*m0(j)^3/(rho0(j)^2.)))/(dx^2.);

          end

% Create the snapshots.

          if (count==interval) 
             
             count = 0;
           
             if (output_type==0) 

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

        ymax = rho_left;
        if (rho_left < rho_right) 

          ymax = rho_right;

        end

        plot(x,rho0)
        axis([x(1) x(jmax) -0.5 ymax]);grid;
               Mp(:,i) = getframe;  

             elseif (output_type==1) 
        
               for j=1:jmax
              
               
                 m2(j)=m1(j)/rho1(j);

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

             else

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

             end

           end


%
%  Update the calculations in time
%

        for j=2:jmax-1

          m0(j)=m1(j);
          rho0(j)=rho1(j);
          E0(j)=E1(j);

        end

      end

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


function f=p(left,right,s,x)

  
  a = (right-left)/pi;
  b = (left+right)/2;

  f = a*atan(s*x)+b;












