function [mc,nc,lic]=CellData(N)
% [mc,nc,lic]=CellData(N)
% row ic of mc and nc define to and from nodes, respectively, for each 
% resistor bounding cell ic

% Colin Fox, 28 Jan 2003

nnps = N+1;
nn = nnps*nnps;     % number of nodes

mc = zeros(N*N,4);
nc = zeros(N*N,4);
lic = zeros(N*N,4); % linear index into DR

for ic = 1:N*N % loop ofer cells
    [i2,i1] = ind2sub([N N],ic); % indices of cell ic
    mc(ic,1) = sub2ind([nnps nnps],i2,i1);
    nc(ic,1) = sub2ind([nnps nnps],i2+1,i1);
    lic(ic,1) = sub2ind([nn nn],mc(ic,1),nc(ic,1));
    
    mc(ic,2) = sub2ind([nnps nnps],i2,i1);
    nc(ic,2) = sub2ind([nnps nnps],i2,i1+1);
    lic(ic,2) = sub2ind([nn nn],mc(ic,2),nc(ic,2));
    
    mc(ic,3) = sub2ind([nnps nnps],i2+1,i1);
    nc(ic,3) = sub2ind([nnps nnps],i2+1,i1+1);
    lic(ic,3) = sub2ind([nn nn],mc(ic,3),nc(ic,3));
    
    mc(ic,4) = sub2ind([nnps nnps],i2,i1+1);
    nc(ic,4) = sub2ind([nnps nnps],i2+1,i1+1);
    lic(ic,4) = sub2ind([nn nn],mc(ic,4),nc(ic,4));
end

mc=mc'; % resistors round cell ic now given by column ic
nc=nc';
lic=lic';