function [ao,mu,d]=dvaughan(beta,a1,a2,b,q11,q12,R,w1,w2,mceq1,mceq2);
[n12,nc]=size(b);
[row,ns3]=size(a2);
ns=n12+ns3;
nsc=ns+nc;
m1=[eye(n12) -a2 -b zeros(n12);
    zeros(nc,n12) w2' R b';
    zeros(n12) beta*q12 beta*w1 beta*a1';
    mceq1];
m2=[a1 zeros(n12,ns+nc);
    -w1' zeros(nc,ns+nc);
    -beta*q11 zeros(n12,ns3+nc) eye(n12);
    mceq2];

if rank(m1)>=ns+nc+n12;
  [V,D]=eig(m1\m2);
  d=diag(D);
else;
  [V,D]=eig(m2,m1);
  d=diag(D);
  rdd=find(real(d)>=inf);
  if ~isempty(rdd);
    d(rdd)=real(d(rdd));
  end;
end;
[d,index]=sort(abs(d));
r=1:nsc;
s=index(1:nsc);
t=nsc+1:nsc+n12;
mu=V(t,s)/V(r,s);
Y=inv(V);
ao=V(r,s)*D(s,s)*(Y(s,r)+Y(s,t)*mu);
