function [L_plus] = L_plus(r, M, a, e, delta_1)
% Kerr time-like geodesics
% angular momemntum in r-L plots
% M is mass of a black hole,
% a is an angular momentum per unit mass of a black hole,
% e is an escape energy
% delta_1 is an integral of motion
L_plus = 2*M*a*sqrt(e +delta_1.*r +sqrt(r.*(r.^2 -2*M*r ...
    +a^2).*(e*r.^3 +2*M*delta_1*r.^2 -Q*r +2*M*Q)));