function[v,yb,s,Eg]=lf(data,line,load,xa,b,l)
% y bus formation
yb=0;
for x=1:l
yb(line(x,1),line(x,2))=-1/(line(x,3));
yb(line(x,2),line(x,1))=-1/(line(x,3));
yb(line(x,1),line(x,1))=yb(line(x,1),line(x,1))+1/(line(x,3))+(line(x,4));
yb(line(x,2),line(x,2))=yb(line(x,2),line(x,2))+1/(line(x,3))+(line(x,4));
end
ybz=yb;
for x=1:b
yy(x,:)=ybz(x,x);
ybz(x,x)=0;
end
v=data(:,3);
s=data(:,4)+data(:,5)*sqrt(-1);
q=data(:,5);
% initaliazation of voltages and reactive power
vx=[1.04;0;1];
qx=[0,4,1];
% load flow
while max(abs(v-vx))>.01 ||abs(max(qx-q))>.01
vx=v;
qx=q;
for x=1:b
if data(x,2)==2
v(x)=( ( (conj(s(x))/conj(v(x)) ) -(ybz(x,:)*v) )/yy(x));
end
if data(x,2)==1
vmag=abs(v(x));
q(x)=-sqrt(-1)*imag( conj(v(x))*v.'*yb(:,x));
s(x)=real(s(x))+q(x);
v(x)=( ( (conj(s(x))/conj(v(x)) ) -(ybz(x,:)*v) )/yy(x)) ;
v(x)=vmag*v(x)/abs(v(x));
end
end
end
Eg=v+yb*v.*xa;