      program wave

      parameter( mx=2000 )

      dimension u(mx), u_exact(Mx), r(mx)

      f(x) = sin( wave_no*(x-x_le) )
c      f(x) = 0.5*( 1 - cos( wave_no*(x-x_le) ) )

      pi = acos(-1.0)

      open(unit=1,file='frames.dat',form='unformatted')
      open(unit=3,file='error.dat')

      print *, 'enter Nw, cfl, and i_fd_type'
      read(5,*) Nw, cfl, i_fd_type

      if( i_fd_type .eq. 1 ) then
         a =  1.0
         b = -1.0
         c =  0.0
      else if( i_fd_type .eq.  0 ) then
         a =  0.5
         b =  0.0
         c = -0.5
      else if( i_fd_type .eq. -1 ) then
         a =  0.0
         b =  1.0
         c = -1.0
      else
         print *, 'bad input for i_fd_type'
         stop
      end if

      xL     = 1.0
      period = 0.2
      x_le   = 0.1
      x_te   = x_le + period
      wave_no = 2*pi/period

      t_final  = 0.5
      dt_print = 0.01

      dx = period/float(Nw)
      Nx = nint(xL/dx) + 1

      dt = CFL*dx
      nt = t_final/dt
      i_skip = nint(max(dt_print/dt, 1.0))
      n_frames = nt/i_skip + 1

      print *, n_frames

      do i=1, Nx
         x = (i-1)*dx
         if( x .lt. x_le .or. x .gt. x_te ) then
            u(i) = 0.0
         else
            u(i) = f(x)
         end if
         u_exact(i) = u(i)
      end do

      write(1) Nw, CFL, i_fd_type, 
     &         Nx, xL, dx, x_le, period,
     &         nt, t_final, dt, dt_print, i_skip, n_frames

      n = 0
      t = 0.0
      write(1) i_frame, n, t
      write(1) (u(i),i=1,Nx)
      write(1) (u(i),i=1,Nx)
      write(3,10) t, 0.0

      i_frame = 1

      do n=1, nt

         do i=2, Nx-1
            r(i) = -( a*u(i+1) + b*u(i) + c*u(i-1) )
         end do

         t = n*dt

         do i=2, Nx-1

            u(i) = u(i) + cfl*r(i)

            x = (i-1)*dx
            xw = x - t

            if( xw .lt. x_le .or. xw .gt. x_te ) then
               u_exact(i) = 0.0
            else
               u_exact(i) = f(xw)
            end if

         end do

         u(1 ) = 0.0
         u(Nx) = 0.0
         u_exact(1 ) = 0.0
         u_exact(Nx) = 0.0

         sum = 0.0
         do i=1, Nx
            sum = sum + ( u(i) - u_exact(i) )**2
         end do
         error = sqrt( sum / Nx )

         write(3,10) t, error
10       format(1p2e12.4)

         if( mod(n,i_skip) .eq. 0 ) then
            i_frame = i_frame + 1
            write(1) i_frame, n, t
            write(1) (u(      i),i=1,Nx)
            write(1) (u_exact(i),i=1,Nx)
         end if

      end do

      stop
      end
