clc;
clear all;
close all;
    %     From      To      R       X           Link identification
    %      Bus      Bus     p.u     p.u         link=1 branch=0
linedata=  [0       1       0       0.2         0;
            2       0       0       0.4         0;
            1       2       0       0.8         1;
            1       3       0       0.4         0;            
            2       3       0       0.4         1];
nl=linedata(:,1);
nr=linedata(:,2);
R=linedata(:,3);
X=linedata(:,4);
ll=linedata(:,5);
nbr=length(nl);
nbus=max(max(nl),max(nr));
ZB=R+j*X;
Zbus=zeros(nbus,nbus);
% Rule - 1
for I=1:nbr
    if nl(I)==0 | nr(I)==0
        if nl(I)==0
            K=nr(I);
        elseif nr(I)==0
            K=nl(I);
        end
        Zbus(K,K)=ZB(I);
    end
end
% Rule -2
for J=1:nbr
    if nl(J)>0 & nr(J) >0
        if ll(J) == 0
          Zbus(:,nr(J))=Zbus(:,nl(J));
          Zbus(nr(J),:)=Zbus(nl(J),:);
          Zbus(nr(J),nr(J))= Zbus(nr(J),nr(J))+ZB(J);
        end
    end
end
%Rule-3
for J=1:nbr
    if nl(J)>0 & nr(J) >0
        if ll(J) == 1
          delta_Z=Zbus(:,nr(J))-Zbus(:,nl(J));
          Zll=ZB(J)+Zbus(nr(J),nr(J))+Zbus(nl(J),nl(J))-2*Zbus(nr(J),nl(J));
          P=(delta_Z*(-delta_Z)')./Zll;
          Zbus=Zbus-P;
        end
    end
end
Zbus



