!------------------------------------------------------------------------------!
! Main program to simulate Burgers equation
program main
implicit none

! grid and related constants
real*8, parameter :: xl = -5.; real*8, parameter :: xr = +5.;	! domain limits
integer, parameter :: N = 20; real*8, parameter :: dt = .1;
real*8, parameter :: zeta = 0.001;	 	! Note: this is "alpha" in the Burgers eqn
real*8, parameter :: T = 100;

! RK3 scheme constants
real*8, parameter :: a(1:3) = (/ 0., -5./9., -153./128. /);
real*8, parameter :: alpha(1:3) = (/ 0., 1./3., 3./4. /);
real*8, parameter :: b(1:3) = (/ 1./3., 15./16., 8./15. /);

! solution vectors and other variables
real*8 x, dx, u(1:N), h(1:N), f(1:N);
integer time_step, num_steps, i, rk;

dx = ( xr - xl ) / real(N) ;
num_steps = dnint(T / dt);
u = 0.; h = 0.; f = 0.;

! initial conditions
do i = 1,N
	x = xl + real(i-1) * dx;
    if (x>0 .and. x<=1) then; u(i) = x; end if
    if (x>1 .and. x<=2) then; u(i) = 2. - x; end if
end do

! open M-file for output
open (unit=10, file="output.m"); write(10,*) "u = [];";
print *, ""; print *, ""; 

! main time stepping loop
do time_step = 1, num_steps

	! for each RK step, compute solution
	do rk = 1,3
	
		! compute current du_dt in "f"
		call compute_du_dt(f, u, N, dx, zeta);

		! compute "h" 
		h = a(rk) * h + dt * f;
		
		! Then, update current solution
		u = u + b(rk) * h;
	
	end do
	
	! write current solution to M-file
	write (10,*) "u(end+1,:) = [ ", u," ];";

	if (mod(time_step, 100) == 0) then
		print *, time_step, " / ", num_steps;
	end if
	
end do

! write stuff to plot solution to M-file 
write (10,*) ""; write (10,*) ""; write (10,*) ""; 
write (10,*) "close all; figure; num_steps = ", num_steps, ";";
write (10,*) "for i = 1:size(u,1)"
write (10,*) "    clf; plot(u(i,:), 'o-'); axis([1 size(u,2) -0.4 +1]);"
write (10,*) "    title(sprintf('time step = %d / %d', i, num_steps));"
write (10,*) "    pause(.005); drawnow"
write (10,*) "end"
write (10,*) ""; write (10,*) ""; write (10,*) ""; 

! close output file
close (10);

print *, ""; print *, ""; 
print *, "Type 'output' in Matlab window to see plots"; 
print *, ""; print *, ""; 

end



!------------------------------------------------------------------------------!
! Subroutine to compute du_dt and store it in f
subroutine compute_du_dt(f, u, N, dx, zeta)
integer N; real*8 f(1:N), u(1:N), dx;
f = 0.;

! left periodic boundary condition
i = 1;
f(i) = - u(i) * ( u(i+1) - u(N) ) / ( 2 * dx ) &
	+ zeta * ( u(i+1) - 2*u(i) + u(N) ) / ( dx ** 2. );

! interior
do i = 2,N-1
	f(i) = - u(i) * ( u(i+1) - u(i-1) ) / ( 2 * dx ) &
		+ zeta * ( u(i+1) - 2*u(i) + u(i-1) ) / ( dx ** 2. );
end do

! right periodic boundary condition 
i = N;
f(i) = - u(i) * ( u(1) - u(i-1) ) / ( 2 * dx ) &
	+ zeta * ( u(1) - 2*u(i) + u(i-1) ) / ( dx ** 2. );

end