function [k,c] = bintc(r,z,rx,zx)
% Bintc - consolidation for 4-point bilinear interpolation
% Example:
%   [k,c] = bintc(r,z,rx,zx)
%   Fi = bint(Fx,k,c)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
 
nrx = numel(rx); drx=rx(2)-rx(1);
nzx = numel(zx); dzx=zx(2)-zx(1);
  
% Interpolate flux
% 3 2   F = F0 + a dr + b dz + e dr dz =
% 0 1   F0(1-dr-dz+dr*dz) + F1(dr-dr*dz) + F2 dr*dz + F3(dz-dr*dz)
kr = max(min(floor((r-rx(1))/drx)+1,nrx-1),1);
kz = max(min(floor((z-zx(1))/dzx)+1,nzx-1),1);
k  = sub2ind([nzx,nrx],kz,kr);
dr = (r - rx(kr))/drx; dz = (z - zx(kz))/dzx;
c  = [1-dr-dz+dr.*dz , dr-dr.*dz , dr.*dz , dz-dr.*dz].';

end
