function [L,W] = get_l(n,d)
%GET_L Compute discrete derivative operators.
%
% [L,W] = get_l(n,d)
%
% Computes the discrete approximation L to the derivative operator
% of order d on a regular grid with n points, i.e. L is (n-d)-by-n.
%
% L is stored as a sparse matrix.
%
% Also computes W, an orthonormal basis for the null space of L.

% Per Christian Hansen, UNI-C, 05/26/93.

% Initialization.
if (d<1), error ('Order d must be positive'), end
nd = n-d;

% Compute L.
c = [1,-1,zeros(1,d-1)];
for i=2:d, c = [c(1:d),0] - [0,c(1:d)]; end
L = sparse(nd,n);
for i=1:d+1
  L = L + sparse(1:nd,[1:nd]+i-1,c(i)*ones(1,nd),nd,n);
end

% If required, compute the null vectors W.
if (nargout==2)
  W = zeros(n,d);
  W(:,1) = ones(n,1);
  for i=2:d, W(:,i) = W(:,i-1).*[1:n]'; end
  W = mgs(W);
end
