function [M, K, C] = ielem_mkc(dim, nb, Map, mat, P, G, prop)
% IELEM_MKC Assemble infinite element mass, stiffness, and damping matrices
% 
% Inputs:
%   dim - dimensionality of the problem
%   nb - number of base points
%   Map - coordinates of mapping nodes
%   prop - element properties
%       if prop(1) is set to 1, zero mass matrix is forced 
%       (This is useful for time domain formulations.)

% Number of gaussian quadrature points
ngauss = length(G.w);

% The average distance from the mapping center
a = sqrt(dot(Map(1:nb,:)-Map(nb+1:2*nb,:),Map(1:nb,:)-Map(nb+1:2*nb,:),2));

% Take the distance into account
dMu = reshape(reshape(G.dMu,dim*ngauss*nb*P,nb)*a,dim*ngauss,nb*P); 

% Secondary variables for coordinate transform
J = G.dL*Map;                   % Jacobian of the transformation
dNx = zeros(size(G.dN));        % Gradient of shape functions in xyz space
dDx = zeros(size(G.dD));        % Gradient of weighting functions in xyz space
dMux = zeros(size(dMu));        % Gradient of phase functions in xyz space
dMuxN = zeros(size(dMu));       % Gradient of phase functions x N
dDxN = zeros(size(G.dD));       % Gradient of weight functions x N
dNxD = zeros(size(G.dD));       % Gradient of shape functions x D
dMu2N = zeros(size(G.N));       % Grad Mu dot Grad Mu x N
dNdMu = zeros(size(G.N));

% Coordinbte transform
j = zeros(ngauss,1);
for n = 1 : ngauss
    ind = dim*(n-1)+(1:dim);                % Current indices
    Jac = (J(ind,:)).';                     % Current jacobian matrix
    dNx(ind,:)  = (Jac.') \ G.dN(ind,:);    
    dDx(ind,:)  = (Jac.') \ G.dD(ind,:);    
    dMux(ind,:) = (Jac.') \ dMu(ind,:);
    diagN = diag(G.N(n,:)); 
    dMu2N(n,:)   = dot(dMux(ind,:),dMux(ind,:),1).*G.N(n,:);
    dNdMu(n,:)   = dot(dNx(ind,:),dMux(ind,:),1);
    dDxN(ind,:)  = dDx(ind,:).*repmat(G.N(n,:),dim,1);
    dMuxN(ind,:) = dMux(ind,:)*diagN;
    dNxD(ind,:)  = dNx(ind,:).*repmat(G.D(n,:),dim,1);
    j(n) = abs(det(Jac));
end

% Assemble element matrices

% Numerical integration with Gaussian quadrature
q = G.w.*j;
qd = reshape(repmat(q.',dim, 1),dim*ngauss,1);


K = ((dDxN.'*diag(qd)*dNx + dNx.'*diag(qd)*dNxD))*(mat(1)*mat(2)^2);

C = (((G.D.*G.N).'*diag(q)*dNdMu)-(dMuxN.'*diag(qd)*dDxN)-(dNxD.'*diag(qd)*dMuxN))*mat(1)*mat(2);

% If the property is set to one, then force zero mass matrix
switch prop(1)
    case 1
        M = zeros(size(K));
    otherwise
        M = (G.N.'*diag(q)*(G.D.*G.N) - dMu2N.'*diag(q)*(G.D.*G.N))*mat(1);
end
