function d0=minqe(B,L,T,d0)
%variance-ccovariance estimation based on MINQUE
b=size(T);
m=b(length(b));%subtract the number of variance components 
TT=zeros(length(L));
d=d0;
while 1
%for t=1:2000
d1=d;
for i=1:m
   T(:,:,i)=T(:,:,i)*d(i);
      TT=TT+T(:,:,i);
end
v=(eye(length(L))-B*(B'*TT^-1*B)^-1*B'*TT^-1)*L
C=TT^-1-TT^-1*B*(B'*TT^-1*B)^-1*B'*TT^-1;%
W=TT^-1;

for i=1:m
    %w(i)=L'*C*T(:,:,i)*C*L;
    w(i)=v'*TT^-1*T(:,:,i)*TT^-1*v;
    for j=1:m
      s(i,j)=trace(W*T(:,:,i)*W*T(:,:,j)); s(j,i)=s(i,j); 
    end
end
d=s^-1*w'
d0=d0.*d;
if norm(d-d1)<0.001, break;  end;
end