function [cbar] = getcbar(rhobar,distmat)
% get cbar from rhobar and distance matrix
% use bisection
c0 = 0.0000;
c1 = 1.0E+4;
vd = lvech(distmat);  

% Verify that c0 yields value gt rhobar
i1 = 0;
jj = 0;
while i1 == 0
  v = mean(exp(-c0*vd));
  i1 = v > rhobar;
  if i1 == 0
      c0 = c0/2;
      jj = jj+1;
  end
  if jj > 500
   error('rhobar too large');
  end
end

% Verify that c1 yields value lt rhobar
i1 = 0;
jj = 0;
while i1 == 0
  v = mean(exp(-c1*vd));
  i1 = v < rhobar;
  if i1 == 0
      c1 = 2*c1;
      jj = jj+1;
  end
  if jj > 500
   error('rhobar too small');
  end
end

% Bisection determination of cbar
while (c1-c0) > 0.001
   cm = (c0+c1)/2;
   v = mean(exp(-cm*vd));
   if v < rhobar
      c1 = cm;
   else
      c0 = cm;
   end
end

cbar = (c0+c1)/2;

end

