function Psi = corrections(IX)
% function Psi = corrections(IX,x,y) 
% input: interface structure array IX
%           the intersection points have to be ordered, at first 
%           come the x-type intersections, then y-type
%        x,y: grid points
% output: structure array Psi
% with fields Psi.xx and Psi.yy
% 
% 14.08.03 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);

Psi.xx = sparse(GRID.nx*GRID.ny,6*IN);
Psi.yy = sparse(GRID.nx*GRID.ny,6*IN);

SIZ = [GRID.nx,GRID.ny];

% numbering of jumps:
% u, ux, uy, uxx, uyy, uxy
% 1. 2.  3.  4.   5.   6.

% x-intersections
for p = 1:INx
  % anchor indices
  i = IX.anch(p,1);
  j = IX.anch(p,2);
  hminus = GRID.x(i)-(IX.coord(p,1));

  if hminus<=0 % intersection to the right from anchor
    hplus = GRID.x(i+1)-IX.coord(p,1);
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+1) = -1/GRID.hx/GRID.hx; %u
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+2) = -hplus/GRID.hx/GRID.hx; %ux   
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+4) = -hplus^2/2/GRID.hx/GRID.hx; %uxx

    if i+1<GRID.nx % boundary points need not to be corrected
      Psi.xx(sub2ind(SIZ,i+1,j),6*(p-1)+1) = 1/GRID.hx/GRID.hx;
      Psi.xx(sub2ind(SIZ,i+1,j),6*(p-1)+2) = hminus/GRID.hx/GRID.hx;
      Psi.xx(sub2ind(SIZ,i+1,j),6*(p-1)+4) = hminus^2/2/GRID.hx/GRID.hx;      
    end    
  else
    hplus = GRID.x(i-1)-IX.coord(p,1);
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+1) = -1/GRID.hx/GRID.hx; %u
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+2) = -hplus/GRID.hx/GRID.hx;   
    Psi.xx(sub2ind(SIZ,i,j),6*(p-1)+4) = -hplus^2/2/GRID.hx/GRID.hx;    
    
    if i-1>1
      Psi.xx(sub2ind(SIZ,i-1,j),6*(p-1)+1) = 1/GRID.hx/GRID.hx;
      Psi.xx(sub2ind(SIZ,i-1,j),6*(p-1)+2) = hminus/GRID.hx/GRID.hx;
      Psi.xx(sub2ind(SIZ,i-1,j),6*(p-1)+4) = hminus^2/2/GRID.hx/GRID.hx;        
    end
  end  
end

% y-intersections
for p = (INx+1):IN
  % anchor indices
  i = IX.anch(p,1);
  j = IX.anch(p,2);
  hminus = GRID.y(j)-(IX.coord(p,2));

  if hminus<=0 % intersection to the right from anchor
    hplus = GRID.y(j+1)-IX.coord(p,2);
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+1) = -1/GRID.hy/GRID.hy; %u
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+3) = -hplus/GRID.hy/GRID.hy; %uy   
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+5) = -hplus^2/2/GRID.hy/GRID.hy; %uyy   

    if j+1<GRID.ny
      Psi.yy(sub2ind(SIZ,i,j+1),6*(p-1)+1) = 1/GRID.hy/GRID.hy;
      Psi.yy(sub2ind(SIZ,i,j+1),6*(p-1)+3) = hminus/GRID.hy/GRID.hy;
      Psi.yy(sub2ind(SIZ,i,j+1),6*(p-1)+5) = hminus^2/2/GRID.hy/GRID.hy;      
    end
  else
    hplus = GRID.y(j-1)-IX.coord(p,2);
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+1) = -1/GRID.hy/GRID.hy; %u
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+3) = -hplus/GRID.hy/GRID.hy;    
    Psi.yy(sub2ind(SIZ,i,j),6*(p-1)+5) = -hplus^2/2/GRID.hy/GRID.hy;    
    
    if j-1>1
      Psi.yy(sub2ind(SIZ,i,j-1),6*(p-1)+1) = 1/GRID.hy/GRID.hy;
      Psi.yy(sub2ind(SIZ,i,j-1),6*(p-1)+3) = hminus/GRID.hy/GRID.hy;
      Psi.yy(sub2ind(SIZ,i,j-1),6*(p-1)+5) = hminus^2/2/GRID.hy/GRID.hy;      
    end    
  end
end