function x1=programengine(choice,nbus,nlines,nshunts,tolerance);
nbus=nbus;
nlines=nlines;
nshunts=nshunts;
tolerance=tolerance;

%reading data related to buses
fid=fopen('busdata.txt','r');
a=textread('busdata.txt');
fclose(fid);

%reading data pertaining to lines
fid=fopen('linedata.txt','r');
b=textread('linedata.txt');
fclose(fid);

%reading data pertaining to shunts
fid=fopen('shuntdata.txt','r');
c=textread('shuntdata.txt');
fclose(fid);

%reading data pertaining to weights
fid=fopen('busweightsdata.txt','r');
buswts=textread('busweightsdata.txt');
fclose(fid);

fid=fopen('lineweightsdata.txt','r');
linewts=textread('lineweightsdata.txt');
fclose(fid);

v_wts=buswts(:,1:nbus);
p_wts=buswts(:,nbus+1:2*nbus);
q_wts=buswts(:,2*nbus+1:3*nbus);
pij_wts=linewts(:,1:nlines);
qij_wts=linewts(:,nlines+1:2*nlines);
%pji_wts=linewts(:,2*nlines+1:3*nlines);
%qji_wts=linewts(:,3*nlines+1:4*nlines);

%STATE ESTIMATION PROGRAM

%Formation of Ybus
Ybus=zeros(nbus);
tap=b(:,6);
x=b(:,1);
y=b(:,2);
zreal=b(:,3);
zimag=b(:,4);
y1=b(:,5);

if size(c)~=0
    shun_bus=c(:,1);
    shun_val=c(:,2);
end

%effect of line impedances and shunt admittances
for i=1:length(x)
    tap_effect11=(1/(zreal(i)+zimag(i)*j))/(tap(i)^2);
    tap_effect12=(-1/(zreal(i)+zimag(i)*j))/(tap(i));
    tap_effect21=(-1/(zreal(i)+zimag(i)*j))/(tap(i));
    tap_effect22=(1/(zreal(i)+zimag(i)*j));
    Ybus(x(i),x(i))=Ybus(x(i),x(i))+tap_effect11+y1(i);
    Ybus(y(i),y(i))=Ybus(y(i),y(i))+tap_effect22+y1(i);
    Ybus(x(i),y(i))=Ybus(x(i),y(i))+tap_effect12;
    Ybus(y(i),x(i))=Ybus(y(i),x(i))+tap_effect21;
end
Ybus;

%inclusion of shunts
if size(c)~=0
    for i=1:length(shun_bus)
        Ybus(shun_bus(i),shun_bus(i))=Ybus(shun_bus(i),shun_bus(i))+shun_val(i);
    end
end

%separation into G and B
G=real(Ybus);
B=imag(Ybus);

%computation of measurement mismatch
%the measurement mismatch vector is of the form [V Pi Qi Pij Qij Pji Qji]' where each
%term represents a subvector

%state mismatch vector
%state mismatch vector is of the form [theta V] where term represents a
%subvector and the theta subvector excludes the slack bus angle


v_bus=a(:,2);
ang_bus=a(:,3);
v_meas=a(:,4);
p_meas=a(:,5);
q_meas=a(:,6);

pij_flow=b(:,7);
qij_flow=b(:,8);
%pji_flow=b(:,9);
%qji_flow=b(:,10);

count=0;
flag=1;
while flag>0
%for count=0:5
 %voltage calculation 
 v_calc=v_bus;

%Real power injection calculation
p_calc=[];
for i=1:nbus
    p_calc(i)=0;
    for j=1:nbus
        p_calc(i)=p_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)-ang_bus(j)));
    end
end

%reactive power injection calculations
q_calc=[];
for i=1:nbus
    q_calc(i)=0;
    for j=1:nbus
        q_calc(i)=q_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-ang_bus(j)));
    end
end

start_bus=x;
end_bus=y;

%i to j real flows calculation
pij_calc=[];
length(x)
for i=1:length(x)
    q=start_bus(i);
    r=end_bus(i);
    pij_calc(i)=-((v_bus(q))^2)*G(q,r)+v_bus(q)*v_bus(r)*(cos(ang_bus(q)-ang_bus(r))*G(q,r)+sin(ang_bus(q)-ang_bus(r))*B(q,r));
end

%i to j reactive flows calculation
qij_calc=[];
for i=1:length(x)
    q=start_bus(i);
    r=end_bus(i);
    qij_calc(i)=((v_bus(q))^2)*B(q,r)+v_bus(q)*v_bus(r)*(sin(ang_bus(q)-ang_bus(r))*G(q,r)-cos(ang_bus(q)-ang_bus(r))*B(q,r));
end


%the data for flows from bus j to bus i is not computed in this algorithm
%because the transformer taps are taken wrt bus i. hence since no separate
%calculations are performed for the tapping j to i flows are not considered
%j to i real flows calculation
%pji_calc=[];
%for i=1:length(x)
 %   q=start_bus(i);
  %  r=end_bus(i);
  %  pji_calc(i)=-((v_bus(r))^2)*G(q,r)+v_bus(q)*v_bus(r)*(cos(ang_bus(r)-ang_bus(q))*G(q,r)+sin(ang_bus(r)-ang_bus(q))*B(q,r));
  %end

%j to i reactive flows calculation
%qji_calc=[];
%for i=1:length(x)
 %   q=start_bus(i);
  %  r=end_bus(i);
   % qji_calc(i)=((v_bus(r))^2)*B(q,r)+v_bus(q)*v_bus(r)*(sin(ang_bus(r)-ang_bus(q))*G(q,r)-cos(ang_bus(r)-ang_bus(q))*B(q,r));
   %end

%DETERMINING MISMATCH VECTORS
%voltage mismatch
%v_meas
%v_calc
mmv=v_meas-v_calc;
%mmv
%real injections 
%p_meas
%p_calc
mmp=p_meas'-p_calc;
%mmp
%reactive injections
%q_meas
%q_calc
mmq=q_meas'-q_calc;
%mmq
%real ij flows
%pij_flow
%pij_calc
mmpij=pij_flow'-pij_calc;
%mmpij
%reactive ij flows
%qij_flow
%qij_calc
mmqij=qij_flow'-qij_calc;
%mmqij
%real ji flows
%pji_flow
%pji_calc
%mmpji=pji_flow'-pji_calc;
%mmpji
%reactive ji flows
%qji_flow
%qji_calc
%mmqji=qji_flow'-qji_calc;
%mmqji

%FINAL MISMATCH VECTOR
%mm=[mmv' mmp mmq mmpij mmqij mmpji mmqji]'
mm=[mmv' mmp mmq mmpij mmqij]';
nmeas=length(mm);

matrixstate=zeros(2*nbus);
state=(matrixstate(1,:))';
nstate=length(state);

%formation of the Jacobian
Hv=zeros(nbus,nstate);
Hp=zeros(nbus,nstate);
Hq=zeros(nbus,nstate);
Hpij=zeros(nlines,nstate);
Hqij=zeros(nlines,nstate);
%Hpji=zeros(nlines,nstate);
%Hqji=zeros(nlines,nstate);

bus_number=[];
for i=1:nbus
    bus_number=[bus_number i];
end

%The complete Jacobian H will be an agglomeration of these submatrices
line_tracker=[x';y'];
bus_number;
%formation of Hv submatrix corresponding to voltage mismatches
for i=1:length(mmv)
    for j=1:nstate
        if j~=(nbus+i)
            Hv(i,j)=0;
        end
        if j==(nbus+i)
            Hv(i,j)=1;
        end
    end
end
%Hv

%formation of Hp submatrix corresponding to real power injection mismatches
for i=1:length(mmp)
    for j=1:nstate
        if i==j
            Hp(i,j)=-q_calc(i)-B(i,i)*v_bus(i)*v_bus(i);
        end
        if j~=i & j<=nbus
            Hp(i,j)=v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-ang_bus(j)));
        end
        if j==i+nbus
            Hp(i,j)=p_calc(i)/v_bus(i)+G(i,i)*v_bus(i);
        end
        if j~=i+nbus & j>nbus
            Hp(i,j)=v_bus(i)*v_bus(j-nbus)*(G(i,j-nbus)*cos(ang_bus(i)-ang_bus(j-nbus))+B(i,j-nbus)*sin(ang_bus(i)-ang_bus(j-nbus)));
        end
    end
end
%Hp

%formation of Hq submatrix corresponding to reactive power injection
%mismatches
for i=1:length(mmq)
    for j=1:nstate
        if i==j
            Hq(i,j)=p_calc(i)-G(i,i)*v_bus(i)*v_bus(i);
        end
        if j~=i & j<=nbus
            Hq(i,j)=-v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)-ang_bus(j)));
        end
        if j==i+nbus
            Hq(i,j)=q_calc(i)/v_bus(i)-B(i,i)*v_bus(i);
        end
        if j~=i+nbus & j>nbus
            Hq(i,j)=v_bus(i)*v_bus(j-nbus)*(G(i,j-nbus)*sin(ang_bus(i)-ang_bus(j-nbus))-B(i,j-nbus)*cos(ang_bus(i)-ang_bus(j-nbus)));
        end
    end
end
%Hq

%formation of Hpij matrix corresponding to real power flow between bus i and j
%mismatches
for i=1:length(mmpij)
    for j=1:nstate
        temp=line_tracker(:,i);
        Hpij(i,j)=0;
        if j==temp(1)
            Hpij(i,j)=v_bus(temp(1))*v_bus(temp(2))*(-G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(2)
            Hpij(i,j)=v_bus(temp(1))*v_bus(temp(2))*(G(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2)))-B(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(1)+nbus
            Hpij(i,j)=-2*v_bus(temp(1))*G(temp(1),temp(2))+v_bus(temp(2))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
        end
        if j==temp(2)+nbus
            Hpij(i,j)=v_bus(temp(1))*(G(temp(1),temp(2))*cos(ang_bus(temp(1))-ang_bus(temp(2)))+B(temp(1),temp(2))*sin(ang_bus(temp(1))-ang_bus(temp(2))));
        end
    end
end
%Hpij

%formation of Hqij matrix corresponding to reactive power flow between bus i and j
%mismatches
for i=1:length(mmqij)
    for j=1:nstate
        temp=line_tracker(:,i);