非回転Blackhole の周りの軌道

1.有効ポテンシャル ((Subscript[r, 0,])Mを決めた後で，Manipulateして下さい)

ClearAll["Global`*"];
r0=3;
M=1;

(Subscript[r, max])=30M;(*maximum plot range of r*)
f[r_,l_]:=(1-((2M)/r))*(1+((l^2)/(r^2)));
f[r_]:=(1-((2M)/r))*(1+((l^2)/(r^2)));
(*f(r) is effective potential defined as (((dr/dτ))^2)=(e^2)-f(r).*)

Manipulate[
{cr1,cr2,fmax,fmin}=If[l/M>=Sqrt[12],{((l^2)/(2M))(1-(Sqrt[1-12(((M/l))^2)])),((l^2)/(2M))(1+(Sqrt[1-12(((M/l))^2)])),f[cr1],f[cr2]},{0,0,1,0}];
Show[
Plot[{f[r,l],1},{r,0,(Subscript[r, max])},
PlotStyle->{{Blue,Thickness[0.003]},{Gray,Dashing[0.01]}},PlotRange->{{0,(Subscript[r, max])},{0,1.5}}],(**)
Graphics[{
Dashed,Line[{{r0,0},{r0,2}}],
If[l/M>=Sqrt[12],{(*この部分は Mathematica6ではエラーが出ますが，問題ありません.*)
Text[Style[StringForm["(Subscript[r, knife])=``",N[cr1]],Italic,16],{0.5(Subscript[r, max]),1.5}],
Text[Style[StringForm["(Subscript[r, stable])=``",N[cr2]],Italic,16],{0.5(Subscript[r, max]),1.3}]}]
}]
]
,{{l,0,Style["L/m(l)",Italic,15]},0,10}
]

Manipulate[
l=r0 v Sqrt[ 1 /(1-(v^2))] Sin[θ];
e=(Sqrt[1-((2M)/r0)])/(Sqrt[1-(v^2)]);
Show[
Plot[{f[r,l],(e^2),1},{r,0,(Subscript[r, max])},
PlotStyle->{{Blue,Thickness[0.003]},Green,{Gray,Dashing[0.01]}},PlotRange->{{0,(Subscript[r, max])},{0,1.5}}],
Graphics[{Dashed,Line[{{r0,0},{r0,2}}],
Text[Style[StringForm["l=``",N[l]],Italic,15],{0.5(Subscript[r, max]),0.7}],
Text[Style[StringForm["(Subscript[θ, 0])=``°",N[θ*180/Pi] ],Italic,15],{0.5(Subscript[r, max]),0.5}],
Text[Style[StringForm["(e^2)=``",N[(e^2)]],Italic,15],{0.5(Subscript[r, max]),0.3}],
If[l/M>=Sqrt[12],(*この部分は Mathematica6ではエラーが出ますが，問題ありません.*)
{Text[Style[StringForm["(Subscript[r, knife])=``,",cr1],Italic,15],{0.45(Subscript[r, max]),1.8}],
Text[Style[StringForm["(Subscript[r, stable])=``,",cr2],Italic,15],{0.45(Subscript[r, max]),1.6}],
Text[Style[StringForm["f((Subscript[r, knife]))=``",fmax],Italic,15],{0.8(Subscript[r, max]),1.8}],
Text[Style[StringForm["f((Subscript[r, stable]))=``",fmin],Italic,15],{0.8(Subscript[r, max]),1.6}]}]
}]
]
,{{θ,(π/2),"(Subscript[θ, 0])"},(π/2),Pi},{{v,0,"(Subscript[v, 0])"},0,1}
]

2. 初期設定(上の グラフなどを参考にしながら,(Subscript[θ, 0]),(Subscript[v, 0]),M, (Subscript[r, 0])を決めてください)

M=1;
r0=15;
θ0=(π/2);
v0=0.265;

(*以下の4行は，変えないで下さい*)

(Subscript[γ, 0])=(1/(Sqrt[1-(v0^2)]));
e= (Subscript[γ, 0])(Sqrt[1-((2M)/r0)]);(*energy per 1 static mass*)
l=r0 v0 Sin[θ0](Subscript[γ, 0]) ;(*angular momentum per 1 static mass*)

b=(l/e);

3.物体の運動(宇宙船時間で測る)  [τmaxを適時変更してください]

ClearAll[r,ϕ,τ,c1];
τmax= 1000.52;
c1=NDSolve[{ r''[τ]==-(M/((r[τ])^2))+((l^2)/((r[τ])^3))-3((M (l^2))/((r[τ])^4)),ϕ'[τ]==(l/((r[τ])^2)),
ϕ[0]==0,r[0]==r0,r'[0]==-v0 Cos[θ0](Subscript[γ, 0])(Sqrt[1-2M/(Subscript[r, 0])])},{r,ϕ},{τ,0,τmax} ];
r1[τ_]=Evaluate[r[τ]/.c1][[1]];ϕ1[τ_]=Evaluate[ϕ[τ]/.c1][[1]];
ParametricPlot[{r1[τ]Cos[ϕ1[τ]],r1[τ]Sin[ϕ1[τ]]},{τ,0,τmax}]
xmax=NMaximize[{r1[τ]Cos[ϕ1[τ]],0≤τ<=τmax},τ][[1]];(*以下は Manipulate用の計算*)
ymax=NMaximize[{r1[τ]Sin[ϕ1[τ]],0≤τ<=τmax},τ][[1]];
xmin=NMinimize[{r1[τ]Cos[ϕ1[τ]],0≤τ<=τmax},τ][[1]];
ymin=NMinimize[{r1[τ]Sin[ϕ1[τ]],0≤τ<=τmax},τ][[1]];
xmax=Max[xmax,1.5r0];
ymax=Max[ymax,1.5r0];
xmin=Min[xmin,-1.5r0];ymin=Min[ymin,-1.5r0];

Manipulate[
knifeedge:={Dashed,Circle[{0,0},cr1]};
stable:={Dashed,Circle[{0,0},cr2]};
spaceship1:=Point[{r1[τ] Cos[ϕ1[τ]],r1[τ] Sin[ϕ1[τ]]}];
blackhole=Disk[{0,0},2];
Show[
Graphics[{Gray,blackhole,Red,knifeedge,Green,stable,Blue,spaceship1}
,PlotRange->{{xmin,xmax},{ymin,ymax}},Axes->True
],ParametricPlot[{r1[u]Cos[ϕ1[u]],r1[u]Sin[ϕ1[u]]},{u,0,τ}]],{τ,0.001,τmax}
]

[有効ポテンシャルのグラフとの重ね合わせ]

potential=ParametricPlot3D[{r Cos[u],r Sin[u],f[r]},{r,2,100M},{u,0,2Pi},Mesh->{100,12},PlotRange->{{xmin,xmax},{ymin,ymax},{0,fmax}},MeshStyle->Gray,PlotStyle->LightBlue];
Manipulate[
spaceship3=Point[{r1[τ] Cos[ϕ1[τ]],r1[τ] Sin[ϕ1[τ]],f[r1[τ]]}];
Show[potential,
Graphics3D[{Blue,PointSize[Large],spaceship3}],Boxed->False],{τ,0,τmax}
]

4.物体の運動 (schwarzschild 時間で測る) [tmaxを適時変更してください]

ClearAll[r,ϕ,t];
tmax=2000;
c2=If[(e^2)>fmax,NDSolve[{ r[t]ϕ'[t]==b/r[t](1-2M/r[t]),r'[t]==-(1-2 M/r[t])Sqrt[1-(1-2 M/r[t])( 1/(e^2)+(b^2)/((r[t])^2))],ϕ[0]==0,r[0]==r0},{r[t],ϕ[t]},{t,0,tmax} ,PrecisionGoal->∞],NDSolve[{ r''[t]==(1/((e^2) ((r[t])^6)))(-2 M+r[t]) (M ((r[t])^2) (6 M+(-3+2 (e^2)) r[t])+(l^2) (10 (M^2)-7 M r[t]+((r[t])^2))),ϕ'[t]==(((1-2M/r[t])l)/(((r[t])^2)e)),
ϕ[0]==0,r[0]==r0,r'[0]==-v0 Cos[θ0](1-2M/r0)},{r,ϕ},{t,0,tmax},PrecisionGoal->∞ ]];
r2[t_]=Evaluate[r[t]/.c2][[1]];ϕ2[t_]=Evaluate[ϕ[t]/.c2][[1]];
ParametricPlot[{r2[t]Cos[ϕ2[t]],r2[t]Sin[ϕ2[t]]},{t,0,tmax}]

Manipulate[
knifeedge:={Dashed,Circle[{0,0},cr1]};
stable:={Dashed,Circle[{0,0},cr2]};
spaceship2=Point[{r2[t] Cos[ϕ2[t]],r2[t] Sin[ϕ2[t]]}];
blackhole=Disk[{0,0},2];
Show[
Graphics[{Gray,blackhole,Red,knifeedge,Green,stable,Blue,spaceship2}
,PlotRange->{{xmin,xmax},{ymin,ymax}},Axes->True
],ParametricPlot[{r2[u]Cos[ϕ2[u]],r2[u]Sin[ϕ2[u]]},{u,0,t}]
],{t,0.001,tmax}
]