function c = orderstatvec(l,wvec,gbar)
a = 1e-6;
% a = 1-.01/max(gbar);
alpha = -log(2)*l./wvec;
xmax = max(-wvec.*log2(alpha./gbar/log(a))+wvec/log(2).*alpha./gbar/log(a)-wvec/log(2));
c = zeros(size(wvec));
for m = 1 : length(c)
  c = quad(@(x) myfun(x,wvec,gbar,l),0,xmax);
end
end

function y = myfun(x,wvec,gbar,l)
alpha = -log(2)*l./wvec;
% X = ones(length(wvec),1)*x;
y = x.*prod(1-exp(-((1./gbar.*alpha)'*ones(1,length(x)))./mylambertw(-exp(-(log(2)./wvec)'*x-1)))).*...
  sum((((1./gbar)'*ones(1,length(x))).*exp(-((1./gbar.*alpha)'*ones(1,length(x)))./...
  mylambertw(-exp(-(log(2)./wvec)'*x-1))).*(-((alpha.^2)'*ones(1,length(x)))/l./...
  (mylambertw(-exp(-(log(2)./wvec)'*x-1)).*(1+mylambertw(-exp(-(log(2)./wvec)'*x-1))))))./...
  (1-exp(-((1./gbar.*alpha)'*ones(1,length(x)))./mylambertw(-exp(-(log(2)./wvec)'*x-1)))));
y(isnan(y)) = 0;
if sum(isnan(y))>0  
  keyboard
  y(isnan(y)) = 0;
end;
if sum(isinf(y))>0
  keyboard
end
% y = zeros(size(x));
% for i  = 1 : length(x)
%   y(i) = x(i).*prod(1-exp(-1./gbar.*alpha./mylambertw(-exp(-x(i)*log(2)./wvec-1)))).*...
%     sum((1./gbar.*exp(-1./gbar.*alpha./mylambertw(-exp(-x(i)*log(2)./wvec-1))).*...
%     (-alpha.^2/l./(mylambertw(-exp(-x(i)*log(2)./wvec-1)).*(1+mylambertw(-exp(-x(i)*log(2)./wvec-1))))))./...
%     (1-exp(-1./gbar.*alpha./mylambertw(-exp(-x(i)*log(2)./wvec-1)))));
% end
end
%
% function y = gcdf(x,)
% y = 1-exp(-1./gbar.*alpha./mylambertw(-exp(-x*log(2)./wvec-1)));