function y = Chi(x,L,delta,n)
% y = Chi(x,L)
% y = Chi(x,L,delta)
% y = Chi(x,L,delta,n)
%
% Chi(x,L) is 1 for |x|<=L/2, zero otherwise
% Chi(x,L,delta) is convolution with phi_delta
% Chi(x,L,delta,n) is convolution with phi_delta, n times (n=0,1,2,...)
%
% need 0<delta<=L/n
% Chi(x,L,delta,n) is zero outside of [(-L-n*delta)/2,(L+n*delta)]
%                  is 1            in [(-L+n*delta)/2,(L-n*delta)]
% The sum of Chi(x-j*L,L,delta,n) over all integers j is 1.
%
% The Fourier transform of Chi(x,L,delta,n) is L*sinc(pi*L*x)*sinc(pi*delta*x)^n

if nargin==2
  delta = 0;
  n = 0;
elseif nargin==3
  n = 1;
end
if n==0
  y = abs(x)<=L/2;
else
  x = (L/2-abs(x))/delta + n/2;
  X = repmat(x(:),1,n) - (0:n-1);
  y = zeros(size(X));
  y(X>=1) = 1;
  ind = (X>0 & X<1);
  y(ind) = X(ind);
  for k=2:n
    X = X(:,1:end-1);
    y = (X.*y(:,1:end-1) + (k-X).*y(:,2:end))/(k-1);
  end
  y = y/n;
end

