function idn = Idn(ndx,n)
%usage: Idn(ndx,n)
% description: Creates the necessary sinc identity matrices I^(i),
% with exponent i=-1,0,1,2.  Be warned that this proc is of a split
% personality.  For i=0,1,2, the matrices of Lund and Bowers [LB,
% p.106] are used, while for i=-1, Stenger [St,p.176] is used...-1 is
% in working order, but very slow.  One could improve this by doing
% table lookup for the necessary values of Si(k).

% Programmer: t. shores	    Last Rev: 4/6/96

% local variables:
% idn: temporary storage for returned matrix
% k,m: loop indices
% inputs:
% ndx: exponent of requested identity
% n: dimension of the requested identity
% outputs:
% return floating point value of the identity.

idn = eye(n);
c = zeros(n,1);
if (ndx == 1)
  for k = 2:n 
     c(k,1) = (-1)^k/(k-1);
  end
  idn = toeplitz(c,-c);
elseif (ndx == -1)
% use sinc integration on intervals [k,k+1] here
% assuming value 0 at endpoints, so do c(2,1) separately
  c(2,1) = 0.5894898722360839;
  M = 64;
  h = pi/sqrt(2*M);
  nodes = h*[-M:M];
  enodes = exp(nodes);
  for k = 3:n
    a = k-2;
    b = k-1;
    inodes = (b*enodes + a) ./ (enodes + 1);
    c(k,1) = c(k-1,1)+h*sum((sin(pi*inodes) ./ (pi*inodes)) .* ...
 (inodes - a) .* (b - inodes));
  end
  idn = 0.5*ones(n,n) + toeplitz(c,-c);
elseif (ndx == 2)
  c(1,1) = -pi^2/3;
  for k = 2:n
      c(k,1) = -2*(-1)^(k-1)/(k-1)^2;
  end
  idn = toeplitz(c,c);
end