function [V] = time_pot(l, r)
% Schwarzschild time-like geodesic
% Potential V for parametr l = L/M (L denotes angular momentum
% about an axis normal to the invariant plane, M is mass)
% r is in mass units.
V = sqrt((1-2./r) .*(1+l^2./r.^2));