function [J]=Jacobian(Node,Element,Agrad,U,U0,rho,style);% computes the Jacobian for 2D EIT when elementwise basis is used.%% INPUT% Node = 剖分節點結構數組% Element =剖分單元結構數組% Agrad = \int_{Element(ii) \nabla\phi_i\cdot\nabla\phi_j% U = 電極處的電壓% U0 = 計算場域內的電壓% rho = 電阻率或電導率向量% style = 重構的電阻率或電導率為實數或復數 %% OUTPUT% J = JacobianNNode=max(size(Node));       %剖分的節點數NElement=max(size(Element)); %剖分的單元數if nargin<3 Agrad=sparse(NNode^2,NElement);%初始化單元的整體雅可比矩陣(梯度陣) %******************************************************* for ii=1:NElement          %對每個單元進行處理  ind=Element(ii).Topology; %提取該單元的拓樸結構（三個相關節點）  gg=reshape([Node(ind).Coordinate],2,max(size(ind)))';   % 三個相關節點的坐標陣3*2  if max(size(ind))==3      %判斷是否為一階單元   anis=grinprodgaus(gg,1); %調用該單元三節點的梯度計算  else   anis=grinprodgausquad(gg,1);% 二階單元處理  end  Aa=sparse(NNode,NNode);   %初始化雅可比矩陣的子矩陣    Aa(ind,ind)=anis;         %賦值雅可比矩陣的子矩陣  Agrad(:,ii)=Aa(:);        %由子矩陣封裝構成雅可比整體矩陣 end J=Agrad; %*******************************************************elseif style=='comp'            %復數計算則如下  J=zeros(size(U,2)*size(U0,2),2*size(Agrad,2)); for ii=1:size(Agrad,2);    Agrad1=reshape(Agrad(:,ii),NNode,NNode);    AgU1=Agrad1*U(1:NNode,:);    AgU2=Agrad1*U(NNode+1:2*NNode,:);    JJ=-U0.'*[AgU1;AgU2];    JJ=JJ(:);    J(:,ii)=JJ;    JJ=-U0.'*[-AgU2;AgU1];    J(:,ii+size(Agrad,2))=JJ(:); endelseif style=='real'       %實數計算則如下 J=zeros(size(U,2)*size(U0,2),size(Agrad,2)); for ii=1:size(Agrad,2);    JJ=U0.'*1/rho(ii)^2*reshape(Agrad(:,ii),NNode,NNode)*U;  JJ=JJ(:);  J(:,ii)=JJ; endendend