function [Agrad,Kb,M,S,C]=FemMatrix(Node,Element,z);

%FemMatrix Computes the blocks of the system matrix for 2D EIT with linear and quadratic basis
% Function [Agrad,Kb,M,S,C]=FemMatrix(Node,Element,z); 
% computes the matrices needed in the finite element approximation of the 2D EIT forward problem. 
%
% INPUT
% Node = nodal data structure
% Element = element data structure
% z = a vector of (complex) contact impedances
%
% OUTPUT
% Agrad = the gradient part of the system matrix 
% Kb,M and S = other blocks of the system matrix 
% C = voltage reference matrix

Nel=max(size(z));                           %電極數
NNode=max(size(Node));                      %節點數
NElement=max(size(Element));                %單元數
M=sparse(NNode,Nel);            %以稀疏陣初始化M
Kb=sparse(NNode,NNode);         %以稀疏陣初始化Kb
Agrad=sparse(NNode^2,NElement); %以稀疏陣初始化Agrad
s=zeros(Nel,1);                 %電極的電導率列向量
g=reshape([Node.Coordinate],2,NNode)';      %節點坐標矩陣

for ii=1:NElement               %對每個單元進行
  A=sparse(NNode,NNode);        %Agrad的預處理矩陣
  ind=(Element(ii).Topology);   %第ii單元的拓樸結構 節點號
  gg=g(ind,:);                  %第ii單元的拓樸節點坐標矩陣 3x2 or 6x2
  if max(size(gg))==3           % 一階單元
   grint=grinprodgaus(gg,1);    % 一階單元梯度計算
  else
   grint=grinprodgausquad(gg,1);% 二階單元梯度計算
  end

 if any([Element(ii).Face{:,3}]),
     %如果ii單元在電極下(Element.Face中的信息)則：
    [In,Jn,InE]=find([Element(ii).Face{:,3}]); 
    bind=Element(ii).Face{In,1};% 單元在邊界上的邊的節點
    ab=g(bind(:),:);            % 單元在邊界上的邊的節點坐標

    if max(size(bind))==2       % 一階單元邊？
     bb1=bound1([ab]);Bb1=zeros(max(size(ind)),1);%計算bound1
     bb2=bound2([ab]);Bb2=zeros(max(size(ind)));  %計算bound2

     s(InE)=s(InE)+1/z(InE)*2*bb1; % 2*bb1 = length of the electrode.
        
     eind=[find(bind(1)==ind),find(bind(2)==ind)];   
    else				    % Second order basis
      bb1=boundquad1([ab]);Bb1=zeros(max(size(ind)),1);
      bb2=boundquad2([ab]);Bb2=zeros(max(size(ind)));

      s(InE)=s(InE)+1/z(InE)*electrlen([ab]); 

      eind=[find(bind(1)==ind),find(bind(2)==ind),find(bind(3)==ind)];
    end 

    Bb1(eind)=bb1;   
    M(ind,InE)=M(ind,InE)-1/z(InE)*Bb1;%計算M矩陣
    
    Bb2(eind,eind)=bb2;
    A(ind,ind)=grint;
    Agrad(:,ii)=A(:);                  %封裝構成Agrad矩陣
    Kb(ind,ind)=Kb(ind,ind)+1/z(InE)*Bb2; %計算Kb矩陣
 else                                     %如果單元不在邊界上則：計算。。.
    A(ind,ind) = grint; 
    Agrad(:,ii)=A(:); 
 end
end  

S=sparse(diag(s));		 %計算S矩陣           

[II1,C]=Current(Nel,NNode,'adj'); %為計算電壓參考C矩陣
C=C(:,1:Nel-1);                   
C=sparse(C(:,1:Nel-1));                             

S=C'*S*C;
M=M*C;