function tlmgcirc
% Kerr timelike marginally-bound (E = 1) critical orbits
a = .8;
r_s = 2.154;      % r_s = 1.85 gives negative values for eta_s
M = 1;

Delta = r_s^2 -2*M*r_s +a^2;
eta_s = r_s^2/a^2/(r_s-M) ...
    *(r_s*(4*a^2*M-r_s*(r_s-3*M)^2)/(r_s-M) ...
    +r_s*(r_s-2*M)^2 -a^2*M ...
    -2*r_s*M/(r_s-M)*Delta*(1-sqrt(M/r_s)))
xi_s = (M*(r_s^2-a^2) -r_s*Delta*sqrt(M/r_s)) /a /(r_s-M);
gama = -a^2*eta_s/r_s^2 +2*M*r_s
beta = 2*M;
mu_max = sqrt(eta_s/(xi_s^2+eta_s));

r = 0:.1:15;
x = 1./(r - r_s);
F = beta*x +gama*x.^2;
mu = sqrt(eta_s/(xi_s^2+eta_s)) ...
    *cos(log(abs(beta+2*gama*x +2*sqrt(gama*F))));
plot(r, r.*mu, 'LineWidth',2)
hold on
ezpolar(num2str(r_s), [-pi/2, pi/2])
line([0 15], [0 15*mu_max], 'LineStyle', ':', 'Color', 'k')
line([0 15], [0 -15*mu_max], 'LineStyle', ':', 'Color', 'k')
axis equal
axis([0 15 -3 3])
title(...
'Timelike geodesic: projection on the (r, \theta) plane', ...
    'Fontsize', 15);