function [x_est_PMU]=SE_WLS_PMU(m_VMA,m_IMA,m_IRI,flag_PMUSE,nob,line_pi)
%WLS Estimation (ONLY PMU)
%States are in rectangular form i.e. e+1i*f
x_est_PMU=[];

e_PMU=[];f_PMU=[];J_V=[];z_V=[];R_V=[];
m_V_polar=zeros(nob,5);%col in order V_abs V_ang R_V_abs R_V_ang
if isempty(m_VMA)==0
    for i=1:length(m_VMA(:,1))
        k=m_VMA(i,1);
        m=m_VMA(i,2);a=m_VMA(i,3);Rm=m_VMA(i,4);Ra=m_VMA(i,5);
        R_real=Rm*cos(a)*cos(a)+m*m*sin(a)*sin(a)*Ra;
        R_imag=Rm*sin(a)*sin(a)+m*m*cos(a)*cos(a)*Ra;
        
        m_V_polar(k,:)=[m a Rm Ra 1];
        if flag_PMUSE==2
            [e_PMU]=[e_PMU;k m*cos(a) 1/R_real];
            [f_PMU]=[f_PMU;k m*sin(a) 1/R_imag];
        else
            t_J_V=zeros(2,2*nob);t_J_V(1,k)=1;t_J_V(2,k+nob)=1;
            [J_V]=[J_V;t_J_V];
            [z_V]=[z_V;m*cos(a);m*sin(a)];
            [R_V]=[R_V;R_real;R_imag];
        end
    end
end

if flag_PMUSE==2
    %CURRENT PHASOR CONVERTED INTO VOLTAGE PHASOR(rectangular form)
    if isempty(m_IMA)==0
        for i=1:length(m_IMA(:,1))
            k=m_IMA(i,1);
            idx_vm=line_pi(k,1);
            if m_V_polar(idx_vm,5)==1
                VM=m_V_polar(idx_vm,1);VA=m_V_polar(idx_vm,2);
                RVM=m_V_polar(idx_vm,3);RVA=m_V_polar(idx_vm,4);
                
                IM=m_IMA(i,2);IA=m_IMA(i,3);RIM=m_IMA(i,4);RIA=m_IMA(i,5);
                
                t_gm=line_pi(k,3);t_bm=line_pi(k,4);
                t_g=line_pi(k,5);t_b=line_pi(k,6);
                
                temp=1/(t_g+1i*t_b);t_w=real(temp);t_v=imag(temp);
                
                temp=(t_g+t_gm+1i*(t_b+t_bm))/(t_g+1i*t_b);
                t_x=real(temp);t_y=imag(temp);
                
                t1=t_x*cos(VA)-t_y*sin(VA);t2=t_x*sin(VA)+t_y*cos(VA);
                t3=t_w*cos(IA)-t_v*sin(IA);t4=t_w*sin(IA)+t_v*cos(IA);
                
                t_e=VM*t1-IM*t3;t_f=VM*t2-IM*t4;
                t_re=RVM*t1^2+RVA*(t2*VM)^2+RIM*t3^2+RIA*(t4*IM)^2;
                t_rf=RVM*t2^2+RVA*(t1*VM)^2+RIM*t4^2+RIA*(t3*IM)^2;
                
                [e_PMU]=[e_PMU;line_pi(k,2) t_e 1/t_re];
                [f_PMU]=[f_PMU;line_pi(k,2) t_f 1/t_rf];
            end
        end
    end
    
    if isempty(m_IRI)==0
        for i=1:length(m_IRI(:,1))
            k=m_IRI(i,1);
            idx_vm=line_pi(k,1);
            if m_V_polar(idx_vm,5)==1
                VM=m_V_polar(idx_vm,1);VA=m_V_polar(idx_vm,2);
                RVM=m_V_polar(idx_vm,3);RVA=m_V_polar(idx_vm,4);
                
                IR=m_IRI(i,2);II=m_IRI(i,3);RIR=m_IRI(i,4);RII=m_IRI(i,5);
                
                t_gm=line_pi(k,3);t_bm=line_pi(k,4);
                t_g=line_pi(k,5);t_b=line_pi(k,6);
                
                temp=1/(t_g+1i*t_b);t_w=real(temp);t_v=imag(temp);
                
                temp=(t_g+t_gm+1i*(t_b+t_bm))/(t_g+1i*t_b);
                t_x=real(temp);t_y=imag(temp);
                
                t1=t_x*cos(VA)-t_y*sin(VA);t2=t_x*sin(VA)+t_y*cos(VA);
                t_e=VM*t1-t_w*IR+t_v*II;t_f=VM*t2-t_v*IR-t*w*II;
                t_re=RVM*t1^2+RVA*(t2*VM)^2+RIR*t_w^2+RII*t_v^2;
                t_rf=RVM*t2^2+RVA*(t1*VM)^2+RIR*t_v^2+RII*t_w^2;
                
                [e_PMU]=[e_PMU;line_pi(k,2) t_e 1/t_re];
                [f_PMU]=[f_PMU;line_pi(k,2) t_f 1/t_rf];
            end
        end
    end
    
    temp_e=zeros(nob,2);%col 2 sum of all weight correseponding to each e
    for i=1:length(e_PMU(:,1))
        k=e_PMU(i,1);
        temp_e(k,1)=temp_e(k,1)+e_PMU(i,2)*e_PMU(i,3);
        temp_e(k,2)=temp_e(k,2)+e_PMU(i,3);
    end
    
    temp_f=zeros(nob,2);%col 2 sum of all weight correseponding to each f
    for i=1:length(f_PMU(:,1))
        k=f_PMU(i,1);
        temp_f(k,1)=temp_f(k,1)+f_PMU(i,2)*f_PMU(i,3);
        temp_f(k,2)=temp_f(k,2)+f_PMU(i,3);
    end
    
    for i=1:nob
        t_V=temp_e(i,1)/temp_e(i,2)+1i*temp_f(i,1)/temp_f(i,2);
        [x_est_PMU]=[x_est_PMU;abs(t_V) angle(t_V) real(t_V) imag(t_V)];
    end
    
    
else
    %CONSIDER VOLTAGE & CURRENT PHASOR DIRECTLY in rectangular (BY DEFAULT)
    t_IRI=m_IRI;
    if isempty(m_IMA)==0
        for i=1:length(m_IMA(:,1))
            k=m_IMA(i,1);
            m=m_IMA(i,2);a=m_IMA(i,3);Rm=m_IMA(i,4);Ra=m_IMA(i,5);
            R_real=Rm*cos(a)*cos(a)+m*m*sin(a)*sin(a)*Ra;
            R_imag=Rm*sin(a)*sin(a)+m*m*cos(a)*cos(a)*Ra;
            [t_IRI]=[t_IRI;k m*cos(a) m*sin(a) R_real R_imag];
        end
    end
    
    z_I=[];J_I=[];R_I=[];
    if isempty(t_IRI)==0
        for i=1:length(t_IRI(:,1))
            k=t_IRI(i,1);idx_m=line_pi(k,1);idx_n=line_pi(k,2);
            [z_I]=[z_I;t_IRI(i,2);t_IRI(i,3)];
            [R_I]=[R_I;t_IRI(i,4);t_IRI(i,5)];
            
            t_gm=line_pi(k,3);t_bm=line_pi(k,4);
            t_g=line_pi(k,5);t_b=line_pi(k,6);
            
            temp=zeros(2,2*nob);
            %Jacobian element for real current measurement w.r.t. e & f
            temp(1,idx_m)=t_g+t_gm;temp(1,idx_n)=-t_g;
            temp(1,idx_m+nob)=-(t_b+t_bm);temp(1,idx_n+nob)=t_b;
            %Jacobian element for reactive current measurement w.r.t. e & f
            temp(2,idx_m)=t_b+t_bm;temp(2,idx_n)=-t_b;
            temp(2,idx_m+nob)=t_g+t_gm;temp(2,idx_n+nob)=-t_g;
            
            [J_I]=[J_I;temp];
        end
    end
    z_PMU=[z_V;z_I];
    R_PMU=diag([R_V;R_I]);
    J_PMU=[J_V;J_I];
    
    W_PMU=inv(R_PMU);
    t_x=inv(J_PMU'*W_PMU*J_PMU)*J_PMU'*W_PMU*z_PMU;
    for i=1:nob
        t_V=t_x(i,1)+1i*t_x(i+nob,1);
        [x_est_PMU]=[x_est_PMU;abs(t_V) angle(t_V) real(t_V) imag(t_V)];
    end
    
end
    

