• Newtonian mechanics
• Schwarzschild metric
• Kerr metric
• Embedding diagram
• Example 1
• Example 2
• Example 3

Relativistic Orbits and Black Holes

Frank Wang, fw@phys.columbia.edu

updated 5/20/2004

Using Maple to derive equations of motion from the Lagrangian, to solve differential equations numerically, and to form graphs based on numerical solutions, we present trajectories of a particle under gravity based on Newtonian theory and general relativity with the Schwarzschild metric and the Kerr metric.  

The Lagrangian for the Kerr solution in general relativity, which describes a spherical object of mass   and spin angular momentum  , is

where

  and  .

If  , namely the object does not rotate, the Kerr solution becomes the Schwarzschild solution.  

In the limit of   and  , the Schwarzschild solution further reduces to the Lagrangian for Newtonian mechanics, i.e., the difference between kinetic energy and potential energy.  

To find equations of motion, we simply employ the Euler-Lagrange equation

      -      =  0

where   is a coordinate and   denotes its derivative with respect to time:     =   (conventional notation for this quantity as   with a dot above this letter).

For additional information, see F. Wang, "Relativistic orbits with computer algebra," Am. J. Phys. 72, (2004 in press) and references therein; contact the author for preprint.   

By varying the initial conditions, one can discover many interesting orbits in relativity.  I thank Professor Roman Koniuk (koniuk@yorku.ca) for suggesting the last three sections.  

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

Newtonian mechanics

First define the kinetic energy, the potential energy, and the Lagrangian:

>	T := 1/2*(diff(r(t),t)^2 + r(t)^2*diff(phi(t),t)^2);

>	V := -G*M/r(t);

>	L := T - V;

>	L1 := subs({r(t)=var1, diff(r(t),t)=var2, phi(t)=var3, diff(phi(t),t)=var4}, L):

Employ the Euler-Lagrange equation for the Newtonian Lagrangian (Epr11 through Eq26)

>	Epr11 := diff(L1, var4):

>	Epr12 := diff(L1, var3):

>	Epr13 := subs({var1=r(t),var2=diff(r(t),t),var3=phi(t),var4=diff(phi(t),t)},Epr11):

>	Eq14 := Epr13 = l;

>	Epr21 := diff(L1, var2):

>	Epr22 := diff(L1, var1):

>	Epr23 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr21):

>	Epr24 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr22):

>	Epr25 := diff(Epr23,t):

>	Eq26 := Epr25 - Epr24 = 0;

Eq31 and Eq32 decouple the differential equations.

>	Eq31 := isolate(Eq14, diff(phi(t),t));

>	Eq32 := eval(Eq26, Eq31);

Provide the initial values, and use Maple to solve these differential equations numerically and to form plots based on numerical solutions.  

>	G:=1; M:=1;

>	Eq41 := r(0) = 26;

>	Eq42 := D(r)(0) = 0;

>	Eq43 := phi(0) = 0;

>	Eq44 := D(phi)(0) = 0.0071;

>	En := eval(T + V, {r(t)=rhs(Eq41), diff(r(t),t)=rhs(Eq42), diff(phi(t),t)=rhs(Eq44)});

>	l := eval(lhs(Eq14), {r(t)=rhs(Eq41), diff(phi(t),t)=rhs(Eq44)});

The eccentricity is not needed; it is listed as a reference.

>	epsilon := sqrt(1 + 2*En*l^2/((G*M)^2));

>	ini1 := Eq41, Eq42, Eq43;

>	Eq51:=dsolve({Eq31, Eq32, ini1}, {r(t), phi(t)}, numeric, output=listprocedure);

>	polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, legend="Newtonian");

>	pnewton := polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, color=black, legend="Newtonian", linestyle=3, thickness=2):

Schwarzschild metric

>	G := 'G': M := 'M':

The procedure for the Schwarzschild metric is almost identical to the Newtonian mechanics; we merely change the Lagrangian.

>	f := 1/2*(-(c^2 - 2*G*M/r(tau))*diff(t(tau),tau)^2      + 1/(1 - 2*G*M/(c^2*r(tau)))*diff(r(tau),tau)^2 + r(tau)^2*diff(phi(tau),tau)^2);

>	f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,            phi(tau)=var5, diff(phi(tau),tau)=var6}, f):

Employ the Euler-Lagrange equations (Epr11 through Eq34)

>	Epr11 := diff(f1, var2):

>	Epr12 := diff(f1, var1):

>	Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>	Eq14 := Epr13 = -a;

>	Epr21 := diff(f1, var4):

>	Epr22 := diff(f1, var3):

>	Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>	Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>	Epr25 := diff(Epr23, tau):

>	Eq26 := Epr25 - Epr24 = 0;

>	Epr31 := diff(f1, var6):

>	Epr32 := diff(f1, var5):

>	Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>	Eq34 := Epr33 = h;

Decouple the differential equations in the following three lines.

>	Eq53 := isolate(Eq14, diff(t(tau),tau));

>	Eq54 := isolate(Eq34, diff(phi(tau),tau));

>	Eq55 := subs({Eq53, Eq54}, Eq26);

Provide the initial conditions, and plot the numerical solutions.

>	M := 1; G := 1; c := 1;

>	Eq77 := r(0) = 26;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = 0.0071;

>	h := eval(lhs(Eq34), {r(tau)=rhs(Eq77), diff(phi(tau),tau)=rhs(Eq80)});

>	a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2));

The effective potential illuminates the turning points.  

>	Vsq := sqrt((1 - 2*M/r)*(1 + h^2/r^2)):

>	plot([Vsq, a], r=2..40, 0.973..1.12, legend=["effective potential", "energy"]);

>	fsolve(Vsq=a, r, 10..20);

>	fsolve(Vsq=a, r, 20..30);

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,               output=listprocedure);

>	psch1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2300],          scaling=constrained, axesfont=[TIMES, ROMAN, 12], legend="Schwarzschild"):

>	psch2 := disk([15.47*cos(3.7797), 15.47*sin(3.7797)], 0.5, color=black):

>	psch3 := disk([15.47*cos(11.3399), 15.47*sin(11.3399)], 0.5, color=black):

>	psch4 := disk([15.47*cos(18.9008), 15.47*sin(18.9008)], 0.5, color=black):

>	psch5 := disk([15.47*cos(26.4221), 15.47*sin(26.4221)], 0.5, color=black):

>	pns := disk([0,0],5.8, color=orange):

>	display([psch1, psch2, psch3, psch4, psch5, pns, pnewton]);

>	pschwarz := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],          scaling=constrained, color=blue, legend=Schwarzschild):

>	pschwarz3d := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]),                 sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..2500],                 coords=cylindrical, color=black, numpoints=400):

Kerr metric

>	M := 'M': l := 'l':

The same procedure also applies to the Kerr metric; we only change the Lagrangian:

>	a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

>	f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2      + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));

>	f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,            phi(tau)=var5, diff(phi(tau),tau)=var6}, f):

Employ the Euler-Lagrange equations (Epr11 through Eq34)

>	Epr11 := diff(f1, var2):

>	Epr12 := diff(f1, var1):

>	Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>	Eq14 := Epr13 = -E;

>	Epr21 := diff(f1, var4):

>	Epr22 := diff(f1, var3):

>	Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>	Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>	Epr25 := diff(Epr23, tau):

>	Eq26 := Epr25 - Epr24 = 0;

>	Epr31 := diff(f1, var6):

>	Epr32 := diff(f1, var5):

>	Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>	Eq34 := Epr33 = l;

A similar decoupling manipulation (Eq53 through Eq55)

>	Eq53 := isolate(Eq14, diff(t(tau),tau));

>	Eq53p := subs(Eq53, Eq34);

>	Eq54 := isolate(Eq53p, diff(phi(tau),tau));

>	Eq54p := subs(Eq54, Eq53);

>	Eq55 := subs({Eq53, Eq54}, Eq26);

Solve the differential equations with initial conditions provided:

>	M := 1; J := 0.37;

>	l := 4.7996; E := 0.9772;

>	Eq77 := r(0) = 26;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,               output=listprocedure);

>	polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],          scaling=constrained, color=red, legend="Kerr direct");

>	pkerrd := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],          scaling=constrained, color=red, legend="Kerr direct"):

>	M := 1; J := -0.37;

>	l := 4.7996; E := 0.976797;

>	Eq77 := r(0) = 26;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,               output=listprocedure);

>	polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],          scaling=constrained, legend="Kerr retro");

>	pkerro := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],          scaling=constrained, color=green, legend="Kerr retro", axesfont=[TIMES, ROMAN, 12]):

>	pkerr1 := disk([16.31*cos(3.63467), 16.31*sin(3.63467)], 0.5, color=black):

>	pkerr2 := disk([14.53*cos(3.923), 14.53*sin(3.923)], 0.5, color=black):

>	pkerr3 := disk([16.31*cos(10.9765), 16.31*sin(10.9765)], 0.5, color=black):

>	pkerr4 := disk([14.53*cos(11.7778), 14.53*sin(11.778)], 0.5, color=black):

>	display([pnewton,pschwarz,pkerrd,pkerro,pkerr1,pkerr2,psch2,pkerr3,pkerr4,psch3,pns]);

Embedding diagram

We derive the embedding formula, see Misner, Thorne, and Wheeler, Gravitation, p. 613 for details.  A star is assumed to have uniform density:

>	M := 'M': m := rp^3/R^3*M;

>	zin := int(sqrt(1/(1 - 2*m/rp) -1), rp=2*M..r) + C;

>	zout := int(sqrt(1/(1 - 2*M/rp) -1), rp=2*M..r);

>	Eq201 := eval(zin, r=R) = eval(zout, r=R):

>	C := solve(Eq201, C):

>	pin := plot3d([r*cos(phi), r*sin(phi), eval(zin, {M=1, R=5.8})], phi=0..2*Pi, r=0..5.8, grid=[25,6], orientation=[88,135], shading=zhue):

>	pout := plot3d([r*cos(phi), r*sin(phi), eval(zout, {M=1, R=5.8})], phi=0..2*Pi, r=5.8..27, grid=[25,12]):

We replot the trajectory in the section "Schwarzschild metric" in this embedding diagram:

>	display([pin, pout, pschwarz3d], style=hidden);

>

Example 1

This example demonstrates an infalling particle with zero angular momentum encountering a Kerr black hole with maximal spin angular momentum ( ).

>

restart:

>	a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

>	f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2      + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));

>	f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,            phi(tau)=var5, diff(phi(tau),tau)=var6}, f):

>	Epr11 := diff(f1, var2):

>	Epr12 := diff(f1, var1):

>	Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>	Eq14 := Epr13 = -E;

>	Epr21 := diff(f1, var4):

>	Epr22 := diff(f1, var3):

>	Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>	Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>	Epr25 := diff(Epr23, tau):

>	Eq26 := Epr25 - Epr24 = 0;

>	Epr31 := diff(f1, var6):

>	Epr32 := diff(f1, var5):

>	Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>	Eq34 := Epr33 = l;

>	Eq53 := isolate(Eq14, diff(t(tau),tau));

>	Eq53p := subs(Eq53, Eq34);

>	Eq54 := isolate(Eq53p, diff(phi(tau),tau));

>	Eq54p := subs(Eq54, Eq53);

>	Eq55 := subs({Eq53, Eq54}, Eq26);

>	M := 1; J := 1;

>	l := 0; E := 1;

>	Eq77 := r(0) = 5;

>	Eq78 := D(r)(0) = -0.64;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,               output=listprocedure);

>	with(plots):

Warning, the name changecoords has been redefined

>	polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..4.49],          scaling=constrained, numpoints=2400);

>

Example 2

This example demonstrates the behavior of a particle near the unstable circular orbit.  The particle orbits several times, then flies out.  An animation of the motion is included.    

>

restart:

>	f := 1/2*(-(c^2 - 2*G*M/r(s))*diff(t(s),s)^2 + 1/(1 - 2*G*M/(c^2*r(s)))*diff(r(s),s)^2 + r(s)^2*diff(phi(s),s)^2);

>	f1 := subs({t(s)=var1, diff(t(s),s)=var2, r(s)=var3, diff(r(s),s)=var4, phi(s)=var5, diff(phi(s),s)=var6}, f):

>	Epr11 := diff(f1, var2):

>	Epr12 := diff(f1, var1):

>	Epr13 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr11):

>	Eq14 := Epr13 = -a/c;

>	Epr21 := diff(f1, var4):

>	Epr22 := diff(f1, var3):

>	Epr23 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr21):

>	Epr24 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr22):

>	Epr25 := diff(Epr23, s):

>	Eq26 := Epr25 - Epr24 = 0;

>	Epr31 := diff(f1, var6):

>	Epr32 := diff(f1, var5):

>	Epr33 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr31):

>	Eq34 := Epr33 = h/c;

>	Eq53 := isolate(Eq14, diff(t(s),s));

>	Eq54 := isolate(Eq34, diff(phi(s),s));

>	Eq55 := eval(Eq26, {Eq53, Eq54});

>	M := 1; G := 1; c := 1;

>	Eq77 := r(0) = 3.76777;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = 0.30290054;

>	h := eval(lhs(Eq34), {r(s)=rhs(Eq77), diff(phi(s),s)=rhs(Eq80)});

>	a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2));

>	Vsq := sqrt((1 - 2*M/x)*(1 + h^2/x^2)):

>	plot([Vsq, a], x=2..25, a-0.1*a..a+0.1*a, legend=["effective potential", "energy"]);

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(s), phi(s)}, numeric, output=listprocedure);

>	with(plots):

Warning, the name changecoords has been redefined

>	polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..100], scaling=constrained);

>	p1 := plot3d([x*cos(t), x*sin(t), sqrt(8*(x-2))], x=2..14, t=0..2*Pi, style=hidden):

>	p2 := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..100], coords=cylindrical, color=black, numpoints=400):

>	display([p1, p2]);

>	N := 50; nstep := 2;

>

for i from 0 by nstep to N*nstep do  pl[i] := polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=i..i+nstep], numpoints=5):  pl3d[i] := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=i..i+nstep], coords=cylindrical, numpoints=5, color=black):  plsum[i] := display(seq(pl[j*nstep], j=0..i/nstep)):  pl3dsum[i] := display({p1, seq(pl3d[j*nstep], j=0..i/nstep)}): end do:

>	display([seq(plsum[i*nstep], i=0..N)], insequence=true, scaling=constrained);

>	display([seq(pl3dsum[i*nstep], i=0..N)], insequence=true);

>

Example 3

Literature on numerical integrations for relativistic orbits can be found in the early 1970's; see M. Johnston and R. Ruffini, "Generalized Wilkins effect and selected orbits in a Kerr-Newman geometry," Phys. Rev. D 10, 2324-2329 (1974).  This example uses the same initial conditions for Fig. 6 of that paper (a particle initially orbits against a maximally rotating black hole).  Our calculation shows interesting details of the last orbits that was omitted in that paper perhaps due to limited computational power of that era.      

>

restart:

>	a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

>	f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2      + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));

>	f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,            phi(tau)=var5, diff(phi(tau),tau)=var6}, f):

>	Epr11 := diff(f1, var2):

>	Epr12 := diff(f1, var1):

>	Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>	Eq14 := Epr13 = -E;

>	Epr21 := diff(f1, var4):

>	Epr22 := diff(f1, var3):

>	Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>	Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>	Epr25 := diff(Epr23, tau):

>	Eq26 := Epr25 - Epr24 = 0;

>	Epr31 := diff(f1, var6):

>	Epr32 := diff(f1, var5):

>	Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>	Eq34 := Epr33 = l;

>	Eq53 := isolate(Eq14, diff(t(tau),tau));

>	Eq53p := subs(Eq53, Eq34);

>	Eq54 := isolate(Eq53p, diff(phi(tau),tau));

>	Eq54p := subs(Eq54, Eq53);

>	Eq55 := subs({Eq53, Eq54}, Eq26);

>	M := 1; J := -1;

>	l := 4; E := 0.97;

>	A := (r^3 + a^2*r + 2*M*a^2);

>	B := -4*a*M*l;

>	C := -(r - 2*M)*l^2 - r*(r^2 - 2*M*r + a^2);

>	Veff := (-B + sqrt(B^2 - 4*A*C))/(2*A);

>	plot({Veff, E}, r=1..30, 0..1.2);

>	fsolve(Veff=E, r=20..30);

>	Eq77 := r(0) = %;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = -5*Pi/4;

>	Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,               output=listprocedure);

>	with(plots): with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

>	pbh := disk([0,0], 1, color=black):

>	pergo := polarplot(2, phi=0..2*Pi, color=black, thickness=2):

>	porb1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=-180..207.4], scaling=constrained, axes=boxed):

This plot resembles Fig. 6 of Johnston and Ruffini: the particle initially orbits counterclockwise, and the black hole spins clockwise.  There is an indication of reversal of direction when the particle approaches the black hole, but we can see more details in the following.  

>	display([pbh, porb1]);

>	pergo3d := sphereplot(M + sqrt(M^2 - a^2*cos(theta)^2), phi=0..2*Pi, theta=0..Pi, scaling=constrained, style=wireframe, grid=[16,16]):

>	pkerrbh := sphere([0,0,0], 1, color=black):

This plot shows the ergosphere of a Kerr black hole: the region between the static limit described by and the event horizon described by .  Within the ergosphere, a particle (even a photon) must corotate with the black hole.  

>	display([pergo3d, pkerrbh]);

>	porb3d1 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=200..207.4], coords=cylindrical, scaling=constrained, color=red, thickness=2):

>	porb3d2 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=207.4..207.648], coords=cylindrical, scaling=constrained, numpoints=800, color=red, thickness=2):

>	porb2 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=200..207.4], scaling=constrained, axes=boxed):

>	porb3 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=207.4..207.648], scaling=constrained, numpoints=400, axes=boxed):

The numerical solution shows that the particle reverses its direction and corotates with the black hole as it spirals toward it.  

>	display([pkerrbh, pergo3d, porb3d1, porb3d2]);

Viewing from the top:

>	display([pbh, pergo, porb2, porb3]);

>