function psiab(m0)
global m mrat
%m0 = 2;
m = m0;
mrat = 0;

xdim = 3*m;
R = 1;
%r = 1e-10;
r = .01;
LB = zeros(xdim,1)+r;
UB = zeros(xdim,1)+R;

%A = zeros(0,xdim);
A = assembleblock({ones(1,m),ones(1,m)});
B = ones(2,1);

%LB = LB+.1;
%UB = UB-.1;

while 1
  %for i=1
  if 0
    abc = rand(m,3);
    x = abc(:);
    x = x/min(x(:));
    [x,fval] = fmincon(@thabc,x,[],[],[],[],0*x+1,[]);
    abc = reshape(x,[m 3]);
    a = abc(:,1);
    b = abc(:,2);
    c = abc(:,3);
  else
    ab = rand(m,2);
    x = ab(:);
    x = x/min(x(:));
    [x,fval] = fminunc(@thab,x);
    ab = reshape(x,[m 2]);
    a = ab(:,1);
    b = ab(:,2);
  end
end

return


function y = thab(x)
global m
ab = reshape(x,[m 2]);
ab = abs(ab);
a = ab(:,1);
b = ab(:,2);

aa = a/sum(a);
bb = b/sum(b);
th = .5*sum(abs(aa-bb));

ab = [a b];
R = max(ab(:));
r = min(ab(:));

%B=(R-r)/(R+r*(m-1));
B=(R-r)/(R+r);

r = th/B;

y = -r;
global mrat mab
if r > mrat
  mrat = r;
  mab = ab;
end
fprintf('psiab m=%d   MRAT = %1.9f\n',m,mrat);
return



function y = thabc(x)
global m
abc = reshape(x,[m 3]);
abc = abs(abc);
a = abc(:,1);
b = abc(:,2);
c = abc(:,3);

aa = a/sum(a);
bb = b/sum(b);
th0 = .5*sum(abs(aa-bb));

ac = a.*c; ac = ac/sum(ac);
bc = b.*c; bc = bc/sum(bc);
thc = .5*sum(abs(ac-bc));

ab = [a b];
R = max(ab(:));
r = min(ab(:));

B=(R-r)/(R+r);

r = thc/B;
%r = abs(th0-thc);
%r = thc/th0;

y = -r;
global mrat mabc
if r > mrat
  mrat = r;
  mabc = abc;
end
fprintf('psiab m=%d   MRAT = %1.9f\n',m,mrat);
return