function wave2d_demo

clear -globals
Globals2D

N = 15; % cubic
FinalTime = 4;
CFL = .75;

% define mesh
K1D = 8;
[Nv, VX, VY, K, EToV] = unif_tri_mesh(K1D);
% [Nv, VX, VY, K, EToV] = MeshReaderGambit2D('Grid/Other/lshape.neu'); 

% computes coordinates, matrices, node mappings, and normals
StartUp2D;

% PlotMesh2D;axis on;return

% plotting nodes
[rp sp] = EquiNodes2D(15); 
Vp = Vandermonde2D(N,rp,sp)/V;
xp = Vp*x; yp = Vp*y;

% % modify node maps to be periodic
% dX = max(VX)-min(VX);
% dY = max(VY)-min(VY);
% BuildPeriodicMaps2D(dX,dY);

%% set wave parameters and initial condition

x0 = -.3; y0 = .3;
p = exp(-25^2*((x-x0).^2 + (y-y0).^2));
u = zeros(Np, K); 
v = zeros(Np, K);

global C2
C2 = ones(Np,K);
for e = 1:K    
    ya = mean(VY(EToV(e,:)));
    if ya > 0 
        C2(:,e) = 2;
    end
end

% % k = 1; % frequency of solution
% % W = (2*k-1)/2*pi;
% % p = cos(W*x).*cos(W*y);

%% run time-stepper

% Runge-Kutta residual storage
resp = zeros(Np,K);
resu = zeros(Np,K);
resv = zeros(Np,K);

% compute time step size
h = 2/K1D;
CN = (N+1)*(N+2)/2;
dt = CFL*h/CN;

Nsteps = ceil(FinalTime/dt);
dt = FinalTime/Nsteps;

% outer time step loop
for i = 1:Nsteps
    
    for INTRK = 1:5
                
        % increment Runge-Kutta residuals and solution
        [rhsp, rhsu, rhsv] = acousticsRHS2D(p,u,v);
        
        % initiate and increment Runge-Kutta residuals
        resp = rk4a(INTRK)*resp + dt*rhsp;
        resu = rk4a(INTRK)*resu + dt*rhsu;
        resv = rk4a(INTRK)*resv + dt*rhsv;
        
        % update fields
        u = u+rk4b(INTRK)*resu;
        v = v+rk4b(INTRK)*resv;
        p = p+rk4b(INTRK)*resp;
        
    end
    
    if mod(i,10)==0
        clf;
        up = Vp*p;
        color_line3(xp,yp,up,up,'.');
        axis([-1,1,-1,1])
        axis equal
        axis off        
%         caxis([0 .25])
%         view(3)
        drawnow
    end
    
end

return

function [rhsp, rhsu, rhsv] = acousticsRHS2D(p,u,v)

Globals2D;

% Define field differences at faces
pP = p(vmapP);
uP = u(vmapP);
vP = v(vmapP);

% way to enforce Dirichlet BCs: reflection BCs
pP(mapB) = -p(vmapB);
uP(mapB) = u(vmapB);
vP(mapB) = v(vmapB);

% % Neumann boundary conditions: dp/dn = 0, u \dot n = 0
% pP(mapB) = p(vmapB);
% uP(mapB) = -u(vmapB);
% vP(mapB) = -v(vmapB);
% 
% % Absorbing boundary conditions
% pP(mapB) = 0;
% uP(mapB) = 0;
% vP(mapB) = 0;

dp = zeros(Nfp*Nfaces,K); dp(:) = pP - p(vmapM);
du = zeros(Nfp*Nfaces,K); du(:) = uP - u(vmapM);
dv = zeros(Nfp*Nfaces,K); dv(:) = vP - v(vmapM);

% evaluate penalty fluxes
tau = 1;
ndotdU = nx.*du + ny.*dv;
fluxp =  tau*dp - ndotdU;
fluxu =  (tau*ndotdU - dp).*nx;
fluxv =  (tau*ndotdU - dp).*ny;

pr = Dr*p; ps = Ds*p;
dpdx = rx.*pr + sx.*ps;
dpdy = ry.*pr + sy.*ps;
divU = Dr*(u.*rx + v.*ry) + Ds*(u.*sx + v.*sy);

% compute right hand sides of the PDE's
rhsp =  -divU + LIFT*(Fscale.*fluxp)/2.0;
rhsu =  -dpdx + LIFT*(Fscale.*fluxu)/2.0;
rhsv =  -dpdy + LIFT*(Fscale.*fluxv)/2.0;

global C2
rhsp = C2.*rhsp;

return;

