function tm_bnd_g
% Schwarzchild null geodesic
% critical orbit
M = 2/14;

phi = 0:0.1:15;
phi_0 =  -2*pi/3;    %2*atanh(1/sqrt(3));
u1 = -1/6/M +tanh((phi-phi_0)/2).^2/2/M;
u2 = 1/3/M + 2*exp(phi)/M./(exp(phi-1)).^2;

polar(phi, 1./u1)
hold on
polar(phi, 1./u2)
title('Special case  {\it D} = 3{\it M}3^{1/2}', 'Fontsize'...
    , 15);