function [L_unob,flag_obs]=check_obs_CSE(line,connect_m,Pindex,Pflowindex,nob)

%"A Direct Numerical Method for Observability Analysis", Bei Gou and Ali Abur
%IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 15, NO. 2, MAY 2000
%If flag_obs==1, observability passed, L_unob=[]
%If flag_obs==0, observability failed, calculate L_unob
nol=length(line(:,1));

Hto=zeros(length(Pindex)+length(Pflowindex),nob);
for i=1:length(Pindex)
    t=Pindex(i);
    Hto(i,:)=Hto(i,:)-connect_m(t,:);
    Hto(i,t)=sum(connect_m(t,:));
end
for i=1:length(Pflowindex)
    t=Pflowindex(i);
    if t>nol
        t_from=line(t-nol,2);t_to=line(t-nol,1);
    else
        t_from=line(t,1);t_to=line(t,2);
    end
    Hto(i+length(Pindex),t_from)=1;
    Hto(i+length(Pindex),t_to)=-1;
end
clear i t t_from t_to;
% nob=4;nol=4;
% Hto=[2 -1 0 -1;-1 2 -1 0;1 -1 0 0];
% line=[1 2;2 3;1 3;3 4];
G=Hto'*Hto;%Formation of gain matrix

%Evaluate L and D such that G=L*D*L'
D=eye(nob);L=eye(nob);D_z=[];
for j=1:nob
    if j==1
        D(j,j)=G(j,j);
    else
        temp=0;
        for k=1:j-1
            temp=temp+D(k,k)*L(j,k)*L(j,k);
        end
        D(j,j)=G(j,j)-temp;
    end
    if abs(D(j,j))<1e-12
        [D_z]=[D_z;j];
    end
    for i=j+1:nob
        temp=0;
        for k=1:j-1
            temp=temp+L(i,k)*L(j,k)*D(k,k);
        end
        L(i,j)=(G(i,j)-temp)/D(j,j);
    end
end

%Formation of branch to node matrix
B2N_mat=zeros(nol,nob);
for i=1:nol
    B2N_mat(i,line(i,1))=1;
    B2N_mat(i,line(i,2))=-1;
end

%Form matrix W consisting of the rows of inverse of L; 
%rows decided by zero value of the diagonal element of matrix D
W_temp=inv(L);
W=[];

L_unob=[];
if isempty(D_z)==0
    for i=1:length(D_z)
        [W]=[W;W_temp(D_z(i),:)];
    end
    if length(D_z)==1
        flag_obs=1;
    else
        flag_obs=0;
        C=B2N_mat*W';
        
        for i=1:nol
            temp=0;
            for j=1:length(D_z)
                if C(i,j)<1e-5
                    C(i,j)=0;
                end
                if C(i,j)~=0
                    temp=temp+1;
                end
            end
            if temp~=0
                [L_unob]=[L_unob;i];
            end
        end
    end
else
    flag_obs=1;
end



