%function rbisect(n0)
function rbisect(T0)
%n = 6; 
%n = n0;
global m n T k
T = T0;
n = size(T,1);
m=2;
k = 2;
%T = zeros(n);
%T(2:n,1) = ones(n-1,1);
%T=makecycle(n);
if 1
  rr = [];
  aa = [];
  for r=0.001:.01:n/2
    aa(end+1) = evalb(r);
    rr(end+1) = r;
  plot(rr,aa);
  drawnow
  end
  return
end
  
LB = 1;
RB = n/2;
while (RB - LB > 1e-9)
  b = (LB + RB)/2;
  b1 = (LB + b)/2;
  b2 = (RB + b)/2;
  a1 = evalb(b1);
  a2 = evalb(b2);
  if (a1 > a2)
    RB = b;
    a = a1;
  else
    LB = b;
    a = a2;
  end
  fprintf('b = %1.9f   a = %1.9f \n',b,a);
end
%b = (LB + RB)/2
%a = evalb(b)
return

function a = evalb(b)
global m n T
psi = ones(2);
%psi(1,2) = b;
%psi(2,1) = b;
psi(1,1) = 1.7;
psi(2,2) = b;
R = max(psi(:))/min(psi(:));
P = treefillprob(m,n,T,psi);
[hh,Ht,H] = gethhn(P,m,n);
%a = log(1-(H-1)/(n-1))/(1-b);
neigh = max(sum(T));
a = H/R/neigh;
return

function T = makecycle(n)
T = zeros(n);
for i=2:n
  T(i,i-1) = 1;
end
T(n,1) = 1;
T = T + T';
return