function d0=minque(B,L,T,d0)
% variance-ccovariance estimation based on MINQUE
% E(L)=B*X
% e=d1*T1+d2*T2+...+d3*T3;
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
for i=1:m
   T(:,:,i)=T(:,:,i)*d(i);
      TT=TT+T(:,:,i);
end
C=TT^-1-TT^-1*B*(B'*TT^-1*B)^-1*B'*TT^-1;
for i=1:m
    w(i)=L'*C*T(:,:,i)*C*L;
    for j=1:m
      s(i,j)=trace(C*T(:,:,i)*C*T(:,:,j)); s(j,i)=s(i,j); 
    end
end
d=s^-1*w';
d0=d0.*d;
if norm(d-1)<0.0001, break;  end;
end