function Q = GetInterpMatrix(rec,mesh,indelem2)
%% Field interpolation function

% rec  = [0 0 0.20];
% rec  = randn(1000,3);

t = tic;
mu0 = 4*pi*1e-7  ;  % V·s/ (A·m)
nrec = size(rec,1);
if nargin==2
    TR = triangulation(mesh.elem2node,mesh.node2coord);
    ID = pointLocation(TR,rec);
elseif nargin==3
    ID = indelem2;
end

a           = mesh.abcd(ID,[1   2  3  4]);
b           = mesh.abcd(ID,[5   6  7  8]);
c           = mesh.abcd(ID,[9  10 11 12]);
d           = mesh.abcd(ID,[13 14 15 16]);
ve          = mesh.vol(ID,1);
Sign        = mesh.sign(ID,:);
elem2edge   = mesh.elem2edge(ID,:);
elemEdgeLength  = mesh.elemEdgeLength(ID,:);
local_i    = [1 2;1 3;1 4;2 3;4 2;3 4];

[Nie, CurlNie] = getEit1(rec(:,1),rec(:,2),rec(:,3),a,b,c,d,ve,local_i,Sign,elemEdgeLength);


Q.Ex=sparse(nrec,mesh.nedge);
Q.Ey=sparse(nrec,mesh.nedge);
Q.Ez=sparse(nrec,mesh.nedge);
Q.Hx=sparse(nrec,mesh.nedge);
Q.Hy=sparse(nrec,mesh.nedge);
Q.Hz=sparse(nrec,mesh.nedge);

ii=repmat((1:nrec)',1,6);

Q.Ex=sparse(ii,elem2edge,squeeze(Nie(:,1,:)),nrec,mesh.nedge);
Q.Ey=sparse(ii,elem2edge,squeeze(Nie(:,2,:)),nrec,mesh.nedge);
Q.Ez=sparse(ii,elem2edge,squeeze(Nie(:,3,:)),nrec,mesh.nedge);
Q.Hx=sparse(ii,elem2edge,squeeze(CurlNie(:,1,:)),nrec,mesh.nedge);
Q.Hy=sparse(ii,elem2edge,squeeze(CurlNie(:,2,:)),nrec,mesh.nedge);
Q.Hz=sparse(ii,elem2edge,squeeze(CurlNie(:,3,:)),nrec,mesh.nedge);

Q.Hx = -Q.Hx;
Q.Hy = -Q.Hy;
Q.Hz = -Q.Hz;

Q.Hx=Q.Hx/mu0;
Q.Hy=Q.Hy/mu0;
Q.Hz=Q.Hz/mu0;

GetInterpMatrixTime = toc(t)
end



function [Nie, CurlNie]=getEit1(xc,yc,zc,a,b,c,d,ve,local_i,TE_Sign,TE_length)
Ne          = size(xc,1);
Nie         = zeros(Ne,3,6);
CurlNie     = zeros(Ne,3,6);

for kt1 = 1:6
    i1              = local_i(kt1,1);
    i2              = local_i(kt1,2);
    Sign_i          = TE_Sign(:,kt1);
    ai1             = a(:,i1);
    ai2             = a(:,i2);
    bi1             = b(:,i1);
    bi2             = b(:,i2);
    ci1             = c(:,i1);
    ci2             = c(:,i2);
    di1             = d(:,i1);
    di2             = d(:,i2);
    l               = TE_length(:,kt1);
    Li1             = (ai1 + bi1.*xc + ci1.*yc + di1.*zc)./(6*ve);
    Li2             = (ai2 + bi2.*xc + ci2.*yc + di2.*zc)./(6*ve);
    GradLi1         = [bi1,ci1,di1]./repmat(6*ve,1,3);
    GradLi2         = [bi2,ci2,di2]./repmat(6*ve,1,3);
    Nie(:,:,kt1)    = repmat(Sign_i,1,3).*(repmat(Li1,1,3).*GradLi2 - repmat(Li2,1,3).*GradLi1).*repmat(l,1,3);
    CurlNie(:,:,kt1) =  repmat(Sign_i,1,3).*[ci1.*di2-di1.*ci2, di1.*bi2-bi1.*di2, bi1.*ci2-ci1.*bi2].*[2*l./((6*ve).^2),2*l./((6*ve).^2),2*l./((6*ve).^2)];
end

end