!
!  This fortran 90 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 
!  Compuatational Science Program between JMU and North Carolina Central
!  University. His info is
!  
!  Jim Sochacki
!  Department of Mathematics
!  James Madison University
!  Harrisonburg, VA   22807
!  jim@math.jmu.edu


      program density_velocity_lw

! This code calculates the Lax Wendroff approximations for the
! density and velocity to the one dimensional Euler equations (E3) of
! gas dynamics. The equations are those given in

      implicit none

      real :: rho_left,rho_right,u_left,u_right
      real :: s_u,s_rho
      real :: c,dx,dt,L,x,time
      
      real,dimension(:),allocatable :: m0,m1,a
      real,dimension(:),allocatable :: rho0,rho1
      
      integer :: count,interval,tstep,jmax
      integer :: output_type,i,j
      
      character*6 type_d,type_v
      character*11 filename

      print*,' This code computes the density and velocity for '
      print*,' the one dimensional nonlinear Euler equations  '
      print*,' m_t + (m**2/rho + c**2 rho)_x = 0 ; -L < x < L, 0 < t < T '  
      print*,' rho_t + m_x = 0 ; -L < x < L, 0 < t < T '  
      print*,' u(x,0) = q_u arctan(s_u*x)+r_u '
      print*,' rho(x,0) = q_rho arctan(s_rho*x)+r_rho '
      print*,' where m = rho u '
      print*,' using Lax Wendroff. '

       print*,' Give c '
       read*,c
       print*,' Give u_left the value of u(x,0) at -infinity '
       read*,u_left
       print*,' Give u_right the value of u(x,0) at infinity '
       read*,u_right
       print*,' Give s_u the "slope" of u(0,0) '
       read*,s_u
       print*,' Give rho_left the value of rho(x,0) at -infinity '
       read*,rho_left
       print*,' Give rho_right the value of rho(x,0) at infinity '
       read*,rho_right
       print*,' Give s_rho the "slope" of rho(0,0) '
       read*,s_rho
       print*,' Give the time to run the model '
       read*,time
       print*,' Give the time step '
       read*, dt
       print*,' Give L '
       read*,L
       print*,' Give the grid size '
       read*, dx


      print*,' Type a 0 for density or a 1  for velocity output '
      print*,' or a 2 or greater for both to be outputted. '
      read*,output_type
      if (output_type==0 .or. output_type>=2) then
        type_d = 'densty'
      end if
      if (output_type>=1) then
        type_v = 'velcty'
      end if
      print*,' Give an integer indicating how often a snapshot '
      print*,' of the profile should be taken. '
      read*,interval

! jmax is the number of grid points in the model
      
      jmax=(2*L)/dx+1
      tstep=time/dt

! Set up the memory

        allocate(m0(0:jmax+1))
        allocate(m1(0:jmax+1))
        allocate(rho0(0:jmax+1))
        allocate(rho1(0:jmax+1))
        allocate(a(0:jmax+1))

        x = -L

!
! Set up the initial values for the leap frog scheme.
!

        do j = 0,jmax+1

          rho0(j) = p(rho_left,rho_right,s_rho,x)
          m0(j) = p(rho_left,rho_right,s_rho,x)*p(u_left,u_right,s_u,x)
          a(j) = c**2
          x = x+dx

        end do
 
        rho1(0) = rho0(0)
        rho1(jmax+1) = rho0(jmax+1)
        m1(0) = m0(0)
        m1(jmax+1) = m0(jmax+1)

! output the initial condition

        i = 0

        if (output_type==0 .or. output_type>=2) then
        
          write (filename,'(a6,i5.5)') type_d,i

          open (6,file=filename,status='new')

          do j=1,jmax

            write (6,'(f14.8)') rho0(j)

          end do

          close (6)

        end if

        if (output_type==1 .or. output_type>=2) then

          write (filename,'(a6,i5.5)') type_v,i

          open (6,file=filename,status='new')

          do j=1,jmax

            write (6,'(f14.8)') m0(j)/rho0(j)

          end do

          close (6)

        end if

!
! Start the calculations.
! tstep is the number of time steps

        count = 0

        do i=1,tstep

          count=count+1
          if (count==interval) print*,' step = ',i

          do j = 1 , jmax           

            rho1(j) = rho0(j) - 1./2.*dt*(-m0(j-1)+m0(j+1))/dx &
+1./4.*dt**2*(-4*m0(j)**2/rho0(j) - &
4*a(j)*rho0(j)+2*m0(j+1)**2/rho0(j+1)+2*a(j+1)*rho0(j+1) + &
2*m0(j-1)**2/rho0(j-1)+2*a(j-1)*rho0(j-1))/(dx**2)

             m1(j) = m0(j) - &
1./2.*dt*(-m0(j-1)**2/rho0(j-1)-a(j-1)*rho0(j-1) + &
m0(j+1)**2/rho0(j+1)+a(j-1)*rho0(j+1))/dx + &
1./4.*dt**2*((2*m0(j+1)/rho0(j+1) + &
2*m0(j)/rho0(j))*(-m0(j)**2/rho0(j)-a(j)*rho0(j)+m0(j+1)**2/rho0(j+1) + &
a(j)*rho0(j+1))+(-m0(j+1)**2/(rho0(j+1)**2)+2*a(j) - &
m0(j)**2/(rho0(j)**2))*(-m0(j)+m0(j+1))-(2*m0(j)/rho0(j) + &
2*m0(j-1)/rho0(j-1))*(-m0(j-1)**2/rho0(j-1)-a(j-1)*rho0(j-1) + &
m0(j)**2/rho0(j)+a(j)*rho0(j))-(-m0(j)**2/(rho0(j)**2) + &
2*a(j-1)-m0(j-1)**2/(rho0(j-1)**2))*(-m0(j-1)+m0(j)))/(dx**2)

          end do

! Create the snapshots.

          if (count==interval)then

            count = 0

            if (output_type==0 .or. output_type>=2) then
        
              write (filename,'(a6,i5.5)') type_d,i

              open (6,file=filename,status='new')

              do j=1,jmax

                write (6,'(f14.8)') rho1(j)

              end do

              close (6)

            end if

            if (output_type>=1) then

              write (filename,'(a6,i5.5)') type_v,i

              open (6,file=filename,status='new')

              do j=1,jmax

                write (6,'(f14.8)') m1(j)/rho1(j)

              end do

              close (6)

            end if

          end if
!
!  Update the calculations in time
!

        do j = 1 , jmax 

          m0(j)=m1(j)
          rho0(j)=rho1(j)
          
        end do

      end do

     if (output_type==0 .or. output_type>=2) then

        print*,' the snapshots are in file: ',type_d

     end if

     if (output_type>=1) then

        print*,' the snapshots are in file: ',type_v

     end if

  contains

  function p(a_l,a_r,s,x)

    real :: x,a_l,a_r,s,p,q,r

    r = (a_l+a_r)/2;
    q = (a_r-a_l)/(4*atan(1.0));

    p = q*atan(s*x)+r

    return

  end function


end  program density_velocity_lw