function [idxs, delta] = evalDelta( XIB, X, Y, N, h);
%
% [idxs delta] = evalDelta( XIB, X, Y, N, h );
%
%  Calculates the delta function at grid locations, idxs
%
%  Returns:
%     idxs  = NIB by 16 matrix, where NIB is the number of IB Pts
%             gives the indices in the Euclidean mesh of where
%             the weight should be applied
%     delta = the weights to apply at each index above
%     
%  Input:
%     XIB   = NIB by 2 array of the location of the IB Pts
%     X     = N*N length vector of all X spatial locations
%     Y     = " but in Y direction
%     N     = number of mesh points in each direction
%     h     = mesh width
%
%  NOTES:
%     MESH POINTS ARE ASSUME TO BE AT LOCATIONS (i,j)*h, i=0..N-1,
%     j=0..N-1
%


persistent isInit xL xR yD yU idxVecX idxVecY;

if( size( isInit ) == 0 )
  isInit = 1;
  
  xL = -2:1;
  xR = -1:2;
  yD = [-2:1];
  yU = [-1:2];
  
  idxVecX = zeros(1, 16);
  idxVecY = zeros(1, 16);
end


XIBLen = length(XIB(:,1));
idxs   = zeros(XIBLen, 16);
delta  = zeros(XIBLen, 16);

for( l = 1:XIBLen )
  
  % cell XIB lies in, assume Xmesh = j * h:
  fIdx    = XIB(l,:) / h;
  cellIdx = round( fIdx );

  % find quadrant of this cell, and hence mesh indices to use, note
  % indices start at zero not one!:
  if( cellIdx(1) > fIdx(1) )
    xIdx = xL + cellIdx(1);
  else
    xIdx = xR + cellIdx(1);    
  end
  
  if( cellIdx(2) > fIdx(2) )
    yIdx = yD + cellIdx(2);    
  else
    yIdx = yU + cellIdx(2);
  end
  
  idxVecX = [ xIdx xIdx xIdx xIdx ]';
  idxVecY = [yIdx(1); yIdx(1); yIdx(1); yIdx(1); yIdx(2); yIdx(2); yIdx(2); yIdx(2); yIdx(3);  yIdx(3);  yIdx(3); yIdx(3); yIdx(4); yIdx(4); yIdx(4); yIdx(4);];
  %idxVecY = reshape( repmat( yIdx, 4, 1 ), 16 , 1);
  
  
  delta(l,:) = (evalPhi( fIdx(1) - idxVecX ) .* evalPhi( fIdx(2) - idxVecY ) )' /(h*h);
  
  xIdx = mod(xIdx, N) + 1;
  yIdx = mod(yIdx, N);

  idxs(l,:) = [ (xIdx + N * yIdx(1)) (xIdx + N * yIdx(2)) (xIdx + N * yIdx(3)) (xIdx + N * yIdx(4))];  
  
end


function [phi] = evalPhi( r );
%
% [phi] = evalPhi( r );
%
%  Calculates phi( r ), phi the discrete delta function (without scaling)
% 
%  Returns:
%     phi = value of unscaled discrete delta function at each point, r
%
%  Input:
%     r   = the points to evaluate the discrete delta function at
%



%idx1 = find ( r < -1 );
%idx2 = find(  (-1 <= r) & (r <= 0) );
%idx3 = find( (1 > r) & (r > 0) );
%idx4 = find( r >= 1 );
%idx2 = setdiff( idx2, idx1 );
%idx3 = setdiff( idx3, idx4 );

phi  = zeros(length(r), 1);

%rT        = r(idx1);
%phi(idx1) = 5 + 2*rT - sqrt( -7 - 12*rT - 4*rT.*rT );

%rT        = r(idx2);
%phi(idx2) = 3 + 2*rT + sqrt( 1 - 4*rT - 4*rT.*rT ); 

%rT        = r(idx3);
%phi(idx3) = 3 - 2*rT + sqrt( 1 + 4*rT - 4*rT.*rT );

%rT        = r(idx4);
%phi(idx4) = 5 - 2*rT - sqrt( -7 + 12*rT - 4*rT.*rT );


aR   = abs( r );
idx3 = find( aR < 1 );
idx4 = find( aR >= 1 );

rT        = aR(idx3);
phi(idx3) = 3 - 2*rT + sqrt( 1 + 4*rT - 4*rT.*rT );

rT        = aR(idx4);
phi(idx4) = 5 - 2*rT - sqrt( -7 + 12*rT - 4*rT.*rT );


phi = phi / 8;