function [G, H] = engnorm(n)

% function [G, H] = engnorm(n)
% 
% This function constructs matrices that define the energy norm
% for the wave equation with Dirichlet boundary conditions on [0,1].
%
% n: discretization parameter (default: n=32).
% G, H: matrices of dimension n-1 such that 
%    energy = u'*G*u + v'*H*v approx int_0^1 |u'(x)|^2 + |v(x)|^2 dx,
% where u stands for displacement, and v stands for velocity.
% This corresponds to the linearization [O inv(M); K D] = [0 I; K D]
% using the mass, stiffness, and damping matrices from ccdamp.m
% 
% Uses Trefethen's "cheb.m" (http://www.comlab.ox.ac.uk/nick.trefethen/cheb.m )
% and "clencurt.m" (http://www.comlab.ox.ac.uk/nick.trefethen/clencurt.m ).

 if nargin<1, n=32; end

 [Dx,x] = cheb(n);                       % Trefethen's cheb.m for [-1,1]
 [x,w]  = clencurt(n);                   % Trefethen's clencurt.m for [-1,1]

 Dx = Dx*2; x = (x+1)/2; w=w/2;          % scale domain to [0,1]

 Dhat = Dx(:,2:n);
 G = Dhat'*diag(w)*Dhat;                 % (d/dx)*(quadrature weights)*(d/dx)
 H = diag(w(2:n));                       % (quadrature weights)