function dgclaw_shallowwater
%DGCLAW_SHALLOWWATER
%   Shallow water equations, to be solved by the program
%   dgclaw_compute.
%
%   2012/06/18 by Benjamin Seibold

%================================================================
% Parameters
%================================================================
par = struct(...
'name','Shallow water equations',... % name of example
'len',1,...                % length of computational domain
'f',@f,...                 % flux function
'feigs',@feigs,...         % eigenvalues of Jacobian
'fL',@fL,'fR',@fR,...      % basis change matrices of Jacobian
'g',@g,...                 % source term
'ic',@ic,...               % initial conditions
'tfinal',1,...             % final time
'nx',500,...               % number of grid cells
'np',2,...                 % polynomial order
'M',10,...                 % upper bound on second derivative at extrema
'steps',200,...            % number of output steps
'flag_plot',1,...          % if yes, animation gets plotted in figure
'flag_savefig',0,...       % if yes, output figure gets saves as png file
'flag_save',0      );      % if yes, solution gets stored in data files

%================================================================
% Run Computation
%================================================================
dgclaw_compute(par)

%================================================================
% Functions
%================================================================
function y = ic(par,x) % Initial condition
y = x;
y(1,:,:,:) = .1+.4*(.2<x&x<.3)+.6*(.6<x&x<.7);
y(2,:,:,:) = x*0;

function z = f(par,u,x,t) % Flux function
z = [u(2,:,:);u(2,:,:).^2./u(1,:,:)+.5*u(1,:,:).^2];

function z = feigs(par,u,x,t) % Eigenvalues of flux function
v = u(2,:)./u(1,:);
w = sqrt(u(1,:));
z = cat(1,v-w,v+w);

function z = fL(par,u,w) % Transform quantity w to
% characteristic variables, based on state u.
la = feigs(par,u);
c = 1./(la(2,:)-la(1,:));
z = [(la(2,:).*w(1,:)-w(2,:)).*c;...
     (-la(1,:).*w(1,:)+w(2,:)).*c];

 function z = fR(par,u,w) % Transform quantity w back
% to physical variables based on state u.
la = feigs(par,u);
z = cat(1,w(1,:)+w(2,:),...
    la(1,:).*w(1,:)+la(2,:).*w(2,:));

function z = g(par,u,x,t) % Source term
z = u*0;
