function [D,Ft] = d_matrix(IX,ju,jun,M,bctype);
% interpolation part of EJIIM for 2D Poisson
% use quadratic polynomial to find coefficients 
% IX : interface structure
% ju : [u]
% jun: [u_n]
%      for Dirichlet b.c. : [jun] = empty
%      for Neumann b.c.   : [u] = empty
% bctype: 'Dirichlet' or 'Neumann'
% M: indicator function at the grid points
% 21.04.04 by Vita
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global GRID  
  
IN = length(IX.type);
x_inters = find(IX.type=='x'); y_inters = find(IX.type=='y');
INx      = length(x_inters)  ; INy      = length(y_inters);

% initialise the memory
D  = sparse(IN*6,GRID.nx*GRID.ny); 
Ft = zeros(IN*6,1);

H = GRID.hx*GRID.hy;

S = spdiags([1 1 1 2 2 1]',0,6,6);

% x-intersections
for p = 1:IN
  i  = IX.anch(p,1) ; j  = IX.anch(p,2);
  ix = IX.coord(p,1); iy = IX.coord(p,2);
    
  if isequal(IX.type(p),'x'); % x-intersection
    hminus = GRID.x(i)-IX.coord(p,1);
  else
    hminus = GRID.y(j)-IX.coord(p,2);
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % STENCIL SELECTION
  % find the sphere of influence. If singular, enlarge k
  k_min = 5.00;
  k_max = 10;
  k = k_min;

  % to save time, search the influence sphere only inside the window
  % window in tripple indexation
  WIN_I = max(1,floor(IX.anch(p,1)-k_max)):min(GRID.nx,ceil(IX.anch(p,1)+k_max));
  WIN_J = max(1,floor(IX.anch(p,2)-k_max)):min(GRID.ny,ceil(IX.anch(p,2)+k_max));
 
  % window in single indexation for sub2ind function
  win_size_x = length(WIN_I);
  win_size_y = length(WIN_J);
  
  temp = ones(win_size_x,win_size_y);
  col1_matr = diag(WIN_I')*temp; 
  col1 = reshape(col1_matr',win_size_x*win_size_y,1);
  
  col2 = repmat(WIN_J',win_size_x,1);
  
  WIN = sub2ind([GRID.nx,GRID.ny],col1,col2);
  
  % x,y,z coordinates of the window
  X_WIN = GRID.X(WIN);
  Y_WIN = GRID.Y(WIN);
  M_WIN = M(WIN);
  
  LEFT = GRID.hy*GRID.hy*(X_WIN-ix).^2 + GRID.hx*GRID.hx*(Y_WIN-iy).^2;
       
  dx_max = k*GRID.hx; 
  dy_max = k*GRID.hy;
 
  % multiply with WIN to keep information about the global index
  STENCIL = ((LEFT<H^2*k^2)&(M_WIN==1)).*WIN; 

  clear X1 Y1 Z1 ST1
  
  [i,j,STENCIL] = find(STENCIL);
  
  hp = GRID.X(STENCIL)-ix;
  kp = GRID.Y(STENCIL)-iy;

  % set the ls-matrix
  d = sqrt(hp.^2 + kp.^2);
  ww = 1./(1+d/GRID.hx);
  W = diag(ww);  

  MM = [ones(size(hp)), hp, kp, hp.^2, kp.^2, hp.*kp];
  
  l = length(STENCIL);
  % multiply weights with the restriction matrix

  W2R = sparse([1:l]',STENCIL,ww.^2,l,GRID.nx*GRID.ny);
  B = S*(inv(MM'*W*W*MM))*MM'*W2R;

  if isequal(bctype,'Dirichlet')
    Ft(6*(p-1)+1)        = ju(p)       ; % [u]
    D(6*(p-1)+2,STENCIL) = B(2,STENCIL); % [ux]
    D(6*(p-1)+3,STENCIL) = B(3,STENCIL); % [uy]
  elseif isequal(bctype,'Neumann')
    D(6*(p-1)+1,STENCIL) = B(1,STENCIL); % [u]
    
    D(6*(p-1)+2,STENCIL) =  IX.t(p,1)^2*B(2,STENCIL) + ...
	                    IX.t(p,1)*IX.t(p,2)*B(3,STENCIL); % [ux]
    Ft(6*(p-1)+2) = IX.n(p,1)*jun(p)                        ; % [ux]
							      
    D(6*(p-1)+3,STENCIL) =  IX.t(p,1)*IX.t(p,2)*B(2,STENCIL) + ...
	                    IX.t(p,2)^2*B(3,STENCIL)        ; % [uy]
    Ft(6*(p-1)+3) = IX.n(p,2)*jun(p)                        ; % [uy]
  else
    error('only Dirichlet and Neumann b.c. implemented!')
  end
    
  D(6*(p-1)+4,STENCIL) = B(4,STENCIL); % [uxx]
  D(6*(p-1)+5,STENCIL) = B(5,STENCIL); % [uyy]
  D(6*(p-1)+6,STENCIL) = B(6,STENCIL); % [uxy]
  
end