%  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 density and velocity to the solution to the 2X2
% nonlinear Euler system (E3) with the condition p = c^2 rho using leap frog.
% rho is the density, c is the sound speed, m is the momentum (density times 
% velocity). 

      disp(' This code computes the pressure and velocity for the one dimensional Euler (E3) equation. ')
      disp(' m_t + (m^2/rho + c^2 rho)_x = 0 ; -L < x < L, 0 < t < T '  )
      disp(' rho_t + m_x = 0 ; -L < x < L, 0 < t < T '  )
      disp(' u(x,0) = q_u arctan(s_u*x)+r_u ')
      disp(' rho(x,0) = q_rho arctan(s_rho*x)+r_rho ')
      disp(' where m = rho u ')
      disp('  ')
      disp(' using leap frog. You will now be asked for the input data. ')

      disp(' Give the (constant) sound speed, c ')
      c = input(' ');
      disp(' Give the limit values for u(x,0) : u_left and u_right.  ')
      u_left = input(' u_left ');
      u_right = input(' u_right ');
      disp(' Give the value s for u(x,0).  ')
      s_u = input(' s_u ');
      disp(' Give the limit values for rho(x,0) : rho_left and rho_right.  ')
      rho_left = input(' rho_left ');
      rho_right = input(' rho_right ');
      disp(' Give the value s for rho(x,0).  ')
      s_rho = input(' s_rho ');
      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

      output_type =input(' Type a 0 for density frames and 1 for velocity frames. ');
      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;
        lambda = dt/dx;

%
% Set up the values for the initial conditions

        r_u = (u_left+u_right)/2;
        q_u = (u_right-u_left)/pi;

        r_rho = (rho_left+rho_right)/2;
        q_rho = (rho_right-rho_left)/pi;

        for j=1:jmax
          x(j) = -L+j*dx;
          rho0(j) = q_rho*atan(s_rho*x(j))+r_rho;
          rho1(j) = rho0(j);
          rho2(j) = rho0(j);
          m0(j) = rho0(j)*(q_u*atan(s_u*x(j))+r_u);
          m1(j) = m0(j);
          m2(j) = m0(j);
          a(j) = c^2;
        end 

        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. We use leap frog to get the first step.
%
         for j=2:jmax-1

           rho1(j) = rho0(j) - (dt/dx)*(m0(j+1)-m0(j-1))/2;

           m1(j) = m0(j) - (dt/dx)*(m0(j+1)^2/rho0(j+1)-m0(j-1)^2/rho0(j-1)+a(j+1)*rho0(j+1)-a(j-1)*rho0(j-1))/2;

         end     

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

         if (output_type)
        
           for j=2:jmax-1   
               
             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(:,1) = getframe;
           plot(x,m2);
           Mp(:,1) = getframe;  

         else

           for j = 1:jmax
           for k = 1:2
             z(k,j) = rho1(j);
           end
           end
           pcolor(z);
           Mf(:,1) = getframe;
           plot(x,rho1)
           axis([x(1) x(jmax) -0.5 ymax]);grid;
           Mp(:,1) = getframe;

         end

% tsteps is the number of time steps

         count = 0;

         for i=2:tsteps

           count=count+1;
           fprintf(' step = %d \n',i);
           
           for j=2:jmax-1 
 
             rho2(j) = rho0(j) - (dt/dx)*(m1(j+1)-m1(j-1));
             m2(j) = m0(j) - (dt/dx)*(m1(j+1)^2/rho1(j+1)-m1(j-1)^2/rho1(j-1)+a(j+1)*rho1(j+1)-a(j-1)*rho1(j-1));

           end

          if (count==interval) 
             
             count = 0;
           
             if (output_type) 
        
               for j=2:jmax-1
              
               
                 m0(j)=m2(j)/rho2(j);

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

             else

               for j = 1:jmax
               for k = 1:2
                 z(k,j) = rho2(j);
               end
               end
               pcolor(z);
               Mf(:,i) = getframe;  
               plot(x,rho2)
               axis([x(1) x(jmax) -0.5 ymax]);grid;
               Mp(:,i) = getframe;

             end

           end

%
%  Update the calculations in time
%
          for j=2:jmax-1
              
            rho0(j)=rho1(j);
            rho1(j)=rho2(j);
            m0(j)=m1(j);
            m1(j)=m2(j);

          end
         
      end

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