function solution(type)

%************************************************
% matlab solution file for assignment 2, Hyperbolic Systems, Spring 2001
% by Hans De Sterck
% (starting from Ben Bergen's solution file)
%************************************************
%       type: shock (type=1) or rarefaction (type=2)
%************************************************
%to run: e.g., solution(1)
%************************************************

%output format (for debugging only)
	format long;

%final time, number of spatial points, midpoint, CFL constant
	t_f = 0.3;
	n = 80;
	mid = n/2;
	CFL = 0.4;
        label = 'n=80, CFL=0.4, t=0.3';

%spatial interval and spatial range
	dx = 2/(n-1);
	x = [-1:dx:1];

%initialization
	t = 0;
	tmp_t = 0;
	frame = 2;

%four schemes are implemented, corresponding to the problems in assignment 2:
%FOU (1a), LW (1b), non-conservative (1c), SOU (2a)

%two problems (initial conditions), shock (type 1) or rarefaction (type 2)
	if(type==1)
		left = 2;
		right = 1;
	else
		left = 1;
		right = 2;
        end

% there are n 'physical' cells: numbered from 1 to n
% there are 2 'ghost cells' on every side, numbered -1:0 and n+1:n+2
% u(.,prob) contains the solution in the n cell centers and ghost cells for scheme prob
% x() contains the spatial coordinates of the n cell centers (from -1 to 1)

% matlab does not allow for negative matrix indices, so we have to shift
% all values in the solution vector u() by 2 to the right to accommodate for the two
% ghost cells on the left
% e.g.: u(2+[1:n],prob)=u([3:n+2],prob) gives the solution in the n physical cell centers
%       u(2+[-1:0],prob)=u([1:2],prob) gives the solution in the two ghost cells on the left

%set the boundary conditions in the ghost cells, for the 4 schemes
	u(2+[-1:0],1:4) = left;
	u(2+[n+1:n+2],1:4) = right;

%set the initial condition
	for i=1:n
		if(i<=mid)
			u(2+i,1:4) = left;
		else
			u(2+i,1:4) = right;
		end
	end

%plot the initial condition
	plot(x,u(2+[1:n],1));
	xlabel('x');
	ylabel('u');
	axis([-1 1 0 4]);
%get this picture as the first frame for a movie
	M(1) = getframe;

        %loop through time
	while(t<=t_f)
	
            %determine timestep as minimum over physical cells 
	    dt = min(dx./u(2+[1:n]))*CFL;

            %for most schemes, at every cell interface we need
            %a left and a right state, fluxes corresponding to this
            %state, and average wave speed a_tilde, and finally the
            %numerical flux at the interface
            %
            %there are n+1 interfaces surrounding the n cell centers
            %that need to be updated
            %the intefaces are numbered starting from 1 to the left
            %of cell (center) 1, until n+1 to the right of cell (center) n

            %******************
            %*** FOU scheme ***
            %******************

            %first order scheme: just take the values of u in the
            %cell centers to the left and right of the interface
            %for the left and right values at the interface
            u_left(1:n+1)=u(2+[1:n+1]-1,1);
            u_right(1:n+1)=u(2+[1:n+1],1);

            %the Burgers flux
            flux_left=u_left.^2/2;
            flux_right=u_right.^2/2;

            %a simple average for the wave speed at the interface
            %(wave speed=derivative of flux function)
            a_tilde=(u_left+u_right)/2;

            %general expression for the upwind numerical flux function
            %for a nonlinear equation
            numflux=(flux_left+flux_right)/2-1/2*abs(a_tilde).*(u_right-u_left);

            %the final update in the conservative scheme
	    u(2+[1:n],1) = u(2+[1:n],1) - dt/dx*reshape((numflux([1:n]+1)-numflux([1:n])),n,1);

            %*****************
            %*** LW scheme ***
            %*****************

            %we need the values of u in the
            %cell centers to the left and right of the interface
            u_left(1:n+1)=u(2+[1:n+1]-1,2);
            u_right(1:n+1)=u(2+[1:n+1],2);

            %the Burgers flux
            flux_left=u_left.^2/2;
            flux_right=u_right.^2/2;

            %a simple average for the wave speed at the interface
            %(wave speed=derivative of flux function)
            a_tilde=(u_left+u_right)/2;

            %general expression for the LW numerical flux function
            %for a nonlinear equation
            numflux=(flux_left+flux_right)/2-1/2*dt/dx*a_tilde.^2.*(u_right-u_left);

            %the final update in the conservative scheme
	    u(2+[1:n],2) = u(2+[1:n],2) - dt/dx*reshape((numflux([1:n]+1)-numflux([1:n])),n,1);

            %*******************************
            %*** non-conservative scheme ***
            %*******************************

            %a 'naive' non-conservative generalization of the
            %linear FOU scheme for positive u in the Burgers equation:
            %the final update
	    u(2+[1:n],3) = u(2+[1:n],3) - dt/dx*(u(2+[1:n],3).*(u(2+[1:n],3)-u(2+[1:n]-1,3)));

            %******************
            %*** SOU scheme *** (second order upwind)
            %******************

            %second order scheme: at the cell interfaces, take for left
            %and right states the values of u reconstructed using the slopes
            %and using a limiter to avoid oscillations

            %we need two values on the left and two on the right for every interface
            u_ll(1:n+1)=u(2+[1:n+1]-2,4);
            u_l(1:n+1)=u(2+[1:n+1]-1,4);
            u_r(1:n+1)=u(2+[1:n+1],4);
            u_rr(1:n+1)=u(2+[1:n+1]+1,4);

            %the slope over the interface
            d1=u_r-u_l;
            %the slope to the right of the interface
            d2=u_rr-u_r;
            %the slope to the left of the interface
            d0=u_l-u_ll;

            %determine the left state
            %lim is the limited slope
            lim=sign(d0).*max(0.,min(sign(d0).*d1,abs(d0)));
            u_left=u_l+.5*lim;
            %determine the right state
            lim=sign(d2).*max(0,min(sign(d2).*d1,abs(d2)));
            u_right=u_r-.5*lim;

            %the Burgers flux
            flux_left=u_left.^2/2;
            flux_right=u_right.^2/2;

            %a simple average for the wave speed at the interface
            %(wave speed=derivative of flux function)
            a_tilde=(u_left+u_right)/2;

            %general expression for the upwind numerical flux function
            %for a nonlinear equation
            numflux=(flux_left+flux_right)/2-1/2*abs(a_tilde).*(u_right-u_left);

            %the final update in the conservative scheme
	    u(2+[1:n],4) = u(2+[1:n],4) - dt/dx*reshape((numflux([1:n]+1)-numflux([1:n])),n,1);

            %make a plot of the solution in the present timestep, for the 4 schemes
	    plot(x,u(2+[1:n],1),'ro-');
            hold on;
	    plot(x,u(2+[1:n],2),'gs-');
	    plot(x,u(2+[1:n],3),'b^-');
	    plot(x,u(2+[1:n],4),'k*-');
	    
            %plot the analytically obtained location of a shock
	    tmp_u(1) = left;
	    tmp_u(2) = right;
            %1.5 is the shock speed obtained from the RH relations
	    tmp_t = tmp_t + 1.5*dt;
	    tmp_x(1) = tmp_t;
	    tmp_x(2) = tmp_t;
	    plot(tmp_x, tmp_u, 'r');
	    title(label);
            xlabel('x');
	    ylabel('u');
	    axis([-1 1 0.5 2.5]);
	    hold off;

            %put it in a movie frame
	    M(frame) = getframe;

            %update frame number and time
	    frame = frame+1;
	    t = t+dt;

        end

	movie(M);

