create_matrix

function create_matrix
global zeng2
zeng2.setup.n_var      =length(zeng2.dummyvar);
zeng2.setup.lamda_dt   =zeros(1,zeng2.setup.n_var);
zeng2.setup.nseg       =zeros(zeng2.setup.n_var,1);
zeng2.setup.nseg_length=zeros(zeng2.setup.n_var,2);

for i=1:zeng2.setup.n_var
   zeng2.setup.lamda_dt(i)= zeng2.var{i}.setting.ga/((zeng2.dummyvar{i}.cm*1e-3)*zeng2.var{i}.setting.nodalarea); %1/(mohm.F)
   zeng2.setup.nseg(i)    = zeng2.dummyvar{i}.nseg;
end

zeng2.setup.nseg_total=sum(zeng2.setup.nseg);
zeng2.setup.nseg_length(1,1)=1;
zeng2.setup.nseg_length(:,2)=cumsum(zeng2.setup.nseg);
zeng2.setup.nseg_length(2:zeng2.setup.n_var,1)=zeng2.setup.nseg_length(1:(zeng2.setup.n_var-1),2)+1;

%set L U
zeng2.setup.lamda_dt_matrix=zeros(zeng2.setup.nseg_total);
for i=1:(zeng2.setup.n_var)
   	zeng2.setup.lamda_dt_matrix = zeng2.setup.lamda_dt_matrix -   diag([ zeros(1,zeng2.setup.nseg_length(i,2)-zeng2.setup.nseg(i) )      zeng2.setup.lamda_dt(i)*ones(1,zeng2.setup.nseg(i)-1)        zeros(1,zeng2.setup.nseg_total-zeng2.setup.nseg_length(i,2)) ],-1) - diag([ zeros(1,zeng2.setup.nseg_length(i,2)-zeng2.setup.nseg(i) )       zeng2.setup.lamda_dt(i)*ones(1,zeng2.setup.nseg(i)-1)        zeros(1,zeng2.setup.nseg_total-zeng2.setup.nseg_length(i,2) )],1);
end

%Set connection
if ~isempty(zeng2.setup.connect)
    for i=1:length(zeng2.setup.connect(:,1))
        node(1)=round(zeng2.setup.connect(i,2)*(zeng2.setup.nseg(zeng2.setup.connect(i,1))-1)+1);
        node(2)=round(zeng2.setup.connect(i,4)*(zeng2.setup.nseg(zeng2.setup.connect(i,3))-1)+1);

        node_all(1)=zeng2.setup.nseg_length(zeng2.setup.connect(i,1),2)-zeng2.setup.nseg(zeng2.setup.connect(i,1))+node(1);
        node_all(2)=zeng2.setup.nseg_length(zeng2.setup.connect(i,3),2)-zeng2.setup.nseg(zeng2.setup.connect(i,3))+node(2);
    
        ra1=(zeng2.dummyvar{zeng2.setup.connect(i,1)}.ra*zeng2.var{zeng2.setup.connect(i,1)}.setting.lseg*1e-4)/zeng2.var{zeng2.setup.connect(i,1)}.setting.axialarea;
        ra2=(zeng2.dummyvar{zeng2.setup.connect(i,3)}.ra*zeng2.var{zeng2.setup.connect(i,3)}.setting.lseg*1e-4)/zeng2.var{zeng2.setup.connect(i,3)}.setting.axialarea;
        ga=2/(ra1+ra2); %1/ohm
        zeng2.setup.lamda_dt_matrix(node_all(1),node_all(2))=-ga/((zeng2.dummyvar{zeng2.setup.connect(i,1)}.cm*1e-3)*zeng2.var{zeng2.setup.connect(i,1)}.setting.nodalarea);      %1/(mohm.F)
        zeng2.setup.lamda_dt_matrix(node_all(2),node_all(1))=-ga/((zeng2.dummyvar{zeng2.setup.connect(i,3)}.cm*1e-3)*zeng2.var{zeng2.setup.connect(i,3)}.setting.nodalarea);      %1/(mohm.F)
    
    end
end

%set the diag
for i=1:zeng2.setup.nseg_total;
    zeng2.setup.lamda_dt_matrix(i,i)=-sum(zeng2.setup.lamda_dt_matrix(i,:));
end

%Gamma-Omega
zeng2.setup.Gamma_dt=zeros(zeng2.setup.nseg_total,1);
zeng2.setup.Omega_dt=zeng2.setup.Gamma_dt;
for p=1:zeng2.setup.n_var
   if isempty(zeng2.dummyvar{p}.model)
      zeng2.setup.Gamma_dt([zeng2.setup.nseg_length(p,1):zeng2.setup.nseg_length(p,2)])=1/zeng2.dummyvar{p}.rm/zeng2.dummyvar{p}.cm;
      zeng2.setup.Omega_dt([zeng2.setup.nseg_length(p,1):zeng2.setup.nseg_length(p,2)])=zeng2.setup.Gamma_dt(zeng2.setup.nseg_length(p,1))*zeng2.dummyvar{p}.vini;
   end
end