function u = backwave(data, uend, c, DX, DY, grid, time, u0max)
%BACKAVE: solves a backward wave equation on a bounded domain with prescribed 
% boundary conditions and final conditions


% load grid parameters
N = grid.N;
bsq = grid.bsq;
xxsq = grid.xxsq;
yysq = grid.yysq;

% load time parameters
T = time.T;
dtdisp = time.dtdisp;
dt = time.dt;
nt = time.nt;
ntdisp = time.ntdisp;

% set the initial conditions
u = uend;

vx = zeros( (N+1)^2, 1);
vy = zeros( (N+1)^2, 1);

% sound speed map
c2 = c.^2; % square it now to avoid repeating the task
c2 = reshape(c2, (N+1)^2,1);

figure(10); clf;

% iterations in time
for t=1:nt,
    
    % reshape the solution u(t) as a vector
    uvec = reshape(u, (N+1)^2, 1);
    
    vx = vx - dt*DX*uvec;
    vy = vy - dt*DY*uvec;
    
    uvec = uvec - dt*c2.*(DX*vx + DY*vy);
    
    % implement the boundary condition
    uvec(bsq) = data(end-t,:);
   
    % reshape u(t) back into a matrix for visualization
    u = reshape(uvec, N+1, N+1);

    % display every _ntdisp_ steps
    if ~mod(t,ntdisp), 
        figure(10);        
        nu = norm(u,'fro'); 
        tit = ['u, t= ', num2str(t*dt), ', L^2 norm = ', num2str(nu)]; 
        myImagesc(u,tit);
    end    
    
end

