clear;
clc;
Ybus=zeros(4,4);
Za=%i*.1;
Z12=%i*.25;
Z14=%i*.4;
Z43=%i*.1;
Z42=%i*.2;
Z32=%i*.1;
Z13=%i*.5;
Zg=%i*4;

V1=1;
V2=complex(.9787,-.0513);
V3=complex(1.019,.0445);
V4=complex(1.0098,-.0176);
Ybus(1,1)=1/Zg+1/Z12+1/Z13+1/Z14;
Ybus(1,2)=-1/Z12;
Ybus(2,1)=Ybus(1,2);
Ybus(1,3)=-1/Z13;
Ybus(3,1)=Ybus(1,3);
Ybus(1,4)=-1/Z14;
Ybus(4,1)=Ybus(1,4);
Ybus(2,2)=1/Z12+1/Z32+1/Z42;
Ybus(2,3)=-1/Z32;
Ybus(3,2)=Ybus(2,3);
Ybus(2,4)=-1/Z42;
Ybus(4,2)=Ybus(2,4);
Ybus(3,3)=1/Z13+1/Z43+1/Z32;
Ybus(3,4)=-1/Z43;
Ybus(4,3)=Ybus(3,4);
Ybus(4,4)=1/Z14+1/Z43+1/Z42;
Z=inv(Ybus);
Zt34=Z(3,3)+Z(4,4)-2*Z(3,4);
V=[V1;V2;V3;V4];
I12=(V1-V2)/Z12;
I13=(V1-V3)/Z13;
I14=(V1-V4)/Z14;
I23=(V2-V3)/Z32;
I24=(V2-V4)/Z42;
I34=(V3-V4)/Z43;
Zaa=-Za;
Zden=Zt34-Zaa;
L1234=Zaa/Z12*((Z(1,3)-Z(1,4))-(Z(2,3)-Z(2,4)))/Zden;
L1334=Zaa/Z13*((Z(1,3)-Z(1,4))-(Z(3,3)-Z(3,4)))/Zden;
L1434=Zaa/Z14*((Z(1,3)-Z(1,4))-(Z(4,3)-Z(4,4)))/Zden;
L2334=Zaa/Z32*((Z(2,3)-Z(2,4))-(Z(3,3)-Z(3,4)))/Zden;
L2434=Zaa/Z42*((Z(2,3)-Z(2,4))-(Z(4,3)-Z(4,4)))/Zden;
I112=I12+L1234*I34;
I113=I13+L1334*I34;
I114=I14+L1434*I34;
I223=I23+L2334*I34;
I224=I24+L2434*I34;
Iinj=Ybus*V;
S1=V(1)*conj(Iinj(1));
S2=V(2)*conj(Iinj(2));
S3=V(3)*conj(Iinj(3));
S4=V(4)*conj(Iinj(4));
L3234=Zaa/Z43*((Z(3,3)-Z(3,4))-(Z(2,3)-Z(2,4)))/Zden;
K323=(Z(3,3)-Z(2,3))/Z43;
K321=(Z(3,1)-Z(2,1))/Z43;
K341=(Z(3,1)-Z(4,1))/Z43;
K343=(Z(3,3)-Z(4,3))/Z43;
K323d=K323+L3234*K343;
K321d=K321+L3234*K341;
delP32=K323d*(-1.6760/2)+K321d*(1.670/2);
mprintf("Change in power flow is %.4f",delP32);