function wmed=wmedian(l,P);

%WMEDIAN computes weighted median value of any observation datas.
%
%        wmed = wmedian(l,P) 
%
%        INPUTS : observation vector (l) and Diagonal vector of weight matrix (P).
%
%        OUTPUT: Weighted Median Value (wmed)       
%
%        Reference:
%        Wicki (1999), Robuste Schätzverfahren für die parameterschätzung in geodätischen Netzen,
%        Institut für Geodäsie und Photogrammetrie an der ETH, Zürich, Mitt. Nr. 67.
%        -----------------------------------------------------------------------------      
%        Prepared by R Cuneyt Erenoglu, 2005, Yildiz Technical University, Geodesy Division
%        e-mail:ceren@yildiz.edu.tr
%        -----------------------------------------------------------------------------      

[n,u]=size(l');

lP=[l' P'];
lP1=sortrows(lP);

fark=0;
fark1=0;
for i=1:n
    sum1=0;
    for j=1:i
        sum1=sum1+P(j);
    end
    
    sum2=0;
    for k=1:n
       sum2=sum2+P(k);
    end
    
    dif(i,1)=abs(2*sum1-sum2);
    dif(i,2)=i;
    
    sum3=0;
    for j=1:i
        sum3=sum3+P(j);
    end
    
    sum4=0;
    for k=i+1:n
       sum4=sum4+P(k);
    end
    
    dif1(i,1)=abs(sum3-sum4);
    dif1(i,2)=i;
    
    
    
end

if dif~=dif1
    break
end

dif;

siralidif=sortrows(dif);
smallest=siralidif(1,2);

if dif(smallest,1)==dif((smallest+1),1)
    wmed=l(dif(smallest+1,2));  % weighted median value
else
    sum5=0;
    for i=1:smallest
        sum5=sum5+P(i);
    end
    sum6=0;
    for i=(smallest+1):n
        sum6=sum6+P(i);
    end
    sum7=0;
    for i=1:n
        sum7=sum7+P(i);
    end
    wmed=(sum5*l(smallest)+sum6*l(smallest+1))/sum7;   % weighted median value
end