(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| Mathematica Syntax | http://kerr.yukterez.net | Version: 24.02.2018 |||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] mt1={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}}; mt2={"EventLocator", "Event"-> (r[t]-1000001/1000000 rA)}; mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20}; mt4={"EquationSimplification"-> "Residual"}; mt0=Automatic; mta=mt0; wp=MachinePrecision; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) A=a; (* pseudosphärisch: A=0, kartesisch: A=a *) crd=2; (* Boyer-Lindquist: crd=1, Kerr-Schild: crd=2 *) dsp=1; (* Display Modus *) tmax=50; (* Eigenzeit *) Tmax=50; (* Koordinatenzeit *) TMax=Min[Tmax, т[plunge-1*^-4]]; tMax=Min[tmax, plunge-1*^-4]; (* Integrationsende *) r0=32/10; (* Startradius *) θ0=80 Pi/180; (* Breitengrad *) φ0=0; (* Längengrad *) a=998/1000; (* Spinparameter *) v0=ж0; (* Anfangsgeschwindigkeit *) α0=-Pi/2; (* vertikaler Abschusswinkel *) ψ0=0; (* Bahninklinationswinkel *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) vr0=v0 Sin[α0]; (* radiale Geschwindigkeitskomponente *) vθ0=v0 Cos[α0] Sin[ψ0]; (* longitudinale Geschwindigkeitskomponente *) vφ0=v0 Cos[α0] Cos[ψ0]; (* latitudinale Geschwindigkeitskomponente *) x0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0]; (* Anfangskoordinaten *) y0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0]; z0[A_]:=r0 Cos[θ0]; x0KS[A_]:=((r0 Cos[φ0]+A Sin[φ0]) Sin[θ0]); y0KS[A_]:=((r0 Sin[φ0]-A Cos[φ0]) Sin[θ0]); x0[A_]:=If[crd==1, x0BL[A], x0KS[A]]; y0[A_]:=If[crd==1, y0BL[A], y0KS[A]]; ε=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0; (* Energie und Drehimpulskomponenten *) Lz=vφ0 Ы/j[v0]; pθ0=vθ0 Sqrt[Ξ]/j[v0]; pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2]; Q=FullSimplify[Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0]]; (* Carter Q *) k=Q+Lz^2+a^2 (ε^2+μ); (* Carter k *) μ=If[v0==1, 0, -1]; (* Baryon: μ=-1, Photon: μ=0 *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 3) RADIUS NACH GESCHWINDIGKEIT UND VICE VERSA ||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) rPro=2 (1+Cos[2/3 ArcCos[-a]]); (* prograder Photonenorbitradius *) rRet=2 (1+Cos[2/3 ArcCos[+a]]); (* retrograder Photonenorbitradius *) rTeo=1+2 Sqrt[1-a^3/3] Cos[ArcCos[(1-a^2)/(1-a^2/3)^(3/2)]/3]; δp[r_,a_]:=Quiet[δi/.NSolve[(a^4(-1+r)+2(-3+r)r^4+a^2r(6+r(-5+3 r))+ (* Eq. Ink. Winkel *) 4a Sqrt[a^2+(-2+r)r](a^2+3 r^2)Cos[δi]-a^2(3+r)(a^2+(-2+r)r)Cos[2δi])/(2r(a^2+ (-2+r)r)(r^3+a^2(2+r)))==0&&δi<=π&&δi>=0,δi][[1]]]; vPro=(a^2-2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a+r0^(3/2))); (* Kreisgeschwindigkeit + *) vRet=(a^2+2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a-r0^(3/2))); (* Kreisgeschwindigkeit - *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) rE=1+Sqrt[1-a^2 Cos[θ]^2]; (* äußere Ergosphäre *) RE[A_, w1_, w2_]:=Xyz[xyZ[ {Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2]; rG=1-Sqrt[1-a^2 Cos[θ]^2]; (* innere Ergosphäre *) RG[A_, w1_, w2_]:=Xyz[xyZ[ {Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2]; rA=1+Sqrt[1-a^2]; (* äußerer Horizont *) RA[A_, w1_, w2_]:=Xyz[xyZ[ {Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2]; rI=1-Sqrt[1-a^2]; (* innerer Horizont *) RI[A_, w1_, w2_]:=Xyz[xyZ[ {Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2]; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) horizons[A_, mesh_, w1_, w2_]:=Show[ ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]], ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]], ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]], ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π}, Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]]; BLKS:=Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}]; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) j[v_]:=Sqrt[1+μ v^2]; (* Lorentzfaktor *) mirr=Sqrt[(Sqrt[1-a^2]+1)/2]; (* irreduzible Masse *) я=Sqrt[Χ/Σ]Sin[θ[τ]]; (* axialer Umfangsradius *) яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]]; Ы=Sqrt[χ/Ξ]Sin[θ0]; Σ=r[τ]^2+a^2 Cos[θ[τ]]^2; (* poloidialer Umfangsradius *) Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2; Ξ=r0^2+a^2 Cos[θ0]^2; Δ=r[τ]^2-2r[τ]+a^2; Δi[τ_]:=R[τ]^2-2R[τ]+a^2; δ=r0^2-2r0+a^2; Χ=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ; Χi[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ]; χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ; ц=2r[τ]/Σ; ц0=2r0/Ξ; т[τ_]:=Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit nach Eigenzeit *) д[ξ_]:=Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 0}]]; (* Eigenzeit nach Koordinatenzeit *) T:=Quiet[д[tk]]; ю[τ_]:=Evaluate[t'[τ]/.sol][[1]]; γ[τ_]:=If[μ==0, "Infinity", ю[τ]]; (* totale ZD *) R[τ_]:=Evaluate[r[τ]/.sol][[1]]; (* Boyer-Lindquist Radius *) Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]]; Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]]; ß[τ_]:=Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ]]; gs[τ_]:=(2 (R[τ]^2-a^2 Cos[Θ[τ]]^2) Sqrt[((a^2+R[τ]^2)^2-a^2 (a^2+(R[τ]-2) R[τ]) Sin[Θ[τ]]^2)/(a^2+2 R[τ]^2+a^2 Cos[2 Θ[τ]])^2])/(R[τ]^2+a^2 Cos[Θ[τ]]^2)^2; (* Gravitation *) ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0=Sqrt[χ/δ/Ξ]; (* gravitative ZD *) ω[τ_]:=2R[τ] a/Χi[τ]; ω0=2r0 a/χ; ωd=2r[τ] a/Χ; (* Frame Dragging Winkelgeschwindigkeit *) Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2]; (* Frame Dragging beobachtete Geschwindigkeit *) й[τ_]:=ω[τ] яi[τ] ς[τ]; й0=ω0 Ы ς0; (* Frame Dragging lokale Geschwindigkeit *) ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ]; ж0=Sqrt[ς0^2-1]/ς0; (* Fluchtgeschwindigkeit *) vd[τ_]:=Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2- r[τ]^4(ε-Lz ωd)^2+Δ(Σ+a^2 Sin[θ[τ]]^2 (ε- Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+ a^2 Sin[θ[τ]]^2 Δ](ε - Lz ωd)))]; v[τ_]:=If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]]; (* lokale Dreiergeschwindigkeit *) dst[τ_]:=Evaluate[str[τ]/.sol][[1]]; (* Strecke *) pΘ[τ_]:=Evaluate[pθ[τ] /. sol][[1]]; pΘks[τ_]:=Σi[τ] Θ'[τ]; (* Impuls *) pR[τ_]:=Evaluate[pr[τ] /. sol][[1]]; pRks[τ_]:=Σi[τ]/Δi[τ] R'[τ]; sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]]; epot[τ_]:=ε+μ-ekin[τ]; (* potentielle Energie *) ekin[τ_]:=If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1]; (* kinetische Energie *) (* beobachtete Inklination *) ink0:=б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]]; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_]:=Chop[N[z]]; (* Boyer-Lindquist-Koordinaten *) pr2τ[τ_]:=1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ; pθ2τ[τ_]:=(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ)); DG1={ t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ), t[0]==0, r'[τ]==(pr[τ] Δ)/Σ, r[0]==r0, θ'[τ]==pθ[τ]/Σ, θ[0]==θ0, φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ), φ[0]==φ0, pr'[τ]==pr2τ[τ], pr[0]==pr0, pθ'[τ]==pθ2τ[τ], pθ[0]==pθ0, str'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]], str[0]==0, vlt'[τ]==vd[τ], vlt[0]==0 }; (* Kerr-Schild-Koordinaten *) dr=(pr0 δ)/Ξ; dθ=pθ0/Ξ; dφ=(2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ; dΦ=If[θ0==0, 0, If[θ0==π, 0, dφ]]; φk=φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}]; DG2={ t''[τ]==(2 ((a^4 Cos[θ[τ]]^4+a^2 Cos[θ[τ]]^2 r[τ]-r[τ]^3-r[τ]^4) r'[τ]^2+r[τ] ((a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-2 a^3 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^3 θ'[τ] φ'[τ]+Sin[θ[τ]]^2 (r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]^2+a t'[τ] (a (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+2 (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]))+r'[τ] ((a^4 Cos[θ[τ]]^4+2 a^2 Cos[θ[τ]]^2 r[τ]-2 r[τ]^3-r[τ]^4) t'[τ]+a (a r[τ] (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+(-a^4 Cos[θ[τ]]^4-2 a^2 Cos[θ[τ]]^2 r[τ]+2 r[τ]^3+r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^3, t'[0]==Limit[(ц0 (-dr+a Sin[θ1]^2 dΦ))/(-1+ц0)+\[Sqrt]((1/((-1+ц0)^2 Ξ))(Ξ (dr^2+(-1+ц0) (μ-Ξ dθ^2))+Sin[θ1]^2 dΦ (-2a Ξ dr-(-1+ц0) χ dΦ+ц0^2 a^2 Ξ Sin[θ1]^2 dΦ))), θ1->θ0], t[0]==0, r''[τ]==(-8 (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2 Cos[2 θ[τ]]+r[τ] (2+r[τ])) r'[τ]^2+4 r'[τ] (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) (-2 r[τ]+a^2 Sin[θ[τ]]^2) t'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]-a Sin[θ[τ]]^2 (2 r[τ] (a^2 Cos[θ[τ]]^2 (-4+a^2+a^2 Cos[2 θ[τ]])+2 r[τ] ((2+a^2+a^2 Cos[2 θ[τ]]) r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+a^4 Sin[2 θ[τ]]^2) φ'[τ])+2 (a^2+(-2+r[τ]) r[τ]) (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+4 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+8 a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2))/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3), r'[0]==dr, r[0]==r0, θ''[τ]==(4 a^2 r[τ] Sin[2 θ[τ]] (r'[τ]+t'[τ])^2-8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]^2-8 a ((Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]]) r'[τ]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ]) φ'[τ]+(2 a^6 Cos[θ[τ]]^4+r[τ] (a^4 Cos[θ[τ]]^2 (5+Cos[2 θ[τ]]) r[τ]+2 a^2 (2+Cos[2 θ[τ]]) r[τ]^3+2 r[τ]^5+2 a^2 (a^2 (3+Cos[2 θ[τ]])+4 r[τ]^2) Sin[θ[τ]]^2)) Sin[2 θ[τ]] φ'[τ]^2)/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3), θ'[0]==dθ, θ[0]==θ0, φ''[τ]==If[θ[τ]==0, 0, (4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) r'[τ]^2+4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]^2+4 a r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-8 (a^2 Cos[θ[τ]]^2+r[τ]^2) (Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 r[τ] Sin[2 θ[τ]]) θ'[τ] φ'[τ]+a Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2+8 a t'[τ] (2 Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ]+a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ])+8 r'[τ] ((a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]+a Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ] (2+r[τ])) θ'[τ]-(r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]))/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3)], φ'[0]==dΦ, φ[0]==φk, str'[τ]==vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]], str[0]==0, vlt'[τ]==vd[τ], vlt[0]==0 }; DGL=If[crd==1, DG1, DG2]; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) sol= NDSolve[DGL, {t, r, θ, φ, str, vlt, pr, pθ}, {τ, 0, tmax}, WorkingPrecision-> wp, MaxSteps-> Infinity, Method-> mta, InterpolationOrder-> All, StepMonitor :> (laststep=plunge; plunge=τ; stepsize=plunge-laststep;), Method->{"EventLocator", "Event" :> (If[stepsize<1*^-5, 0, 1])}]; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) XBL[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]]; YBL[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]]; Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]]; XKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+a Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]]; YKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-a Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]]; X[τ_]:=If[crd==1, XBL[τ], XKS[τ]]; Y[τ_]:=If[crd==1, YBL[τ], YKS[τ]]; xBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]]; yBL[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]]; z[τ_]:=Z[τ]; xKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]]; yKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]]; x[τ_]:=If[crd==1, xBL[τ], xKS[τ]]; y[τ_]:=If[crd==1, yBL[τ], yKS[τ]]; XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2]; (* kartesischer Radius *) Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *) xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]}; xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]}; (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) PR=1.2r0; (* Plot Range *) d1=10; (* Schweiflänge *) plp=50; (* Flächenplot Details *) w1l=0; w2l=0; w1r=0; w2r=0; (* Startperspektiven *) Mrec=1000; mrec=15; (* Parametric Plot Subdivisionen *) imgsize=380; (* Bildgröße *) s[text_]:=Style[text, FontSize->font]; font=11; (* Anzeigestil *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 11) RADIUS NACH ZEIT |||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Plot[{R[τ], rI, rA, 0}, {τ,0,tMax}, PlotLabel->"y=r, x=τ (Eigenzeit)", PlotLabels->{None,"r-","r+","ring"}, Frame->True, ImageSize->500] Plot[{R[д[t]], rI, rA, 0}, {t,0,TMax}, PlotLabel->"y=r, x=т (Finkelsteinzeit)", PlotLabels->{None,"r-","r+","ring"}, Frame->True, ImageSize->500] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 12) PLOT NACH EIGENZEIT ||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) display[T_]:=Grid[{ {If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]}, {s[" т coord"], " = ", s[n0[т[tp]]], s["GM/c³"], s[dp]}, {s[" t statc"], " = ", s[n0[т[tp]+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]}, {s[" γ coord"], " = ", s[n0[If[dsp==2, γ[tp]-2 If[crd==1, 0, R'[tp] R[tp]/(R[tp]^2-2 R[tp]+a^2)], γ[T]]]], s["dt/dτ"], s[dp]}, {s[" γ statc"], " = ", s[n0[If[dsp==1, γ[tp]-2 If[crd==1, 0, R'[tp] R[tp]/(R[tp]^2-2 R[tp]+a^2)], γ[T]]]], s["dt/dτ"], s[dp]}, {s[" ς gravt"], " = ", s[n0[ς[tp]]], s["dt/dτ"], s[dp]}, {s[" r coord"], " = ", s[n0[R[tp]]], s["GM/c²"], s[dp]}, {s[" φ longd"], " = ", s[n0[Φ[tp] 180/π]], s["deg"], s[dp]}, {s[" θ lattd"], " = ", s[n0[Θ[tp] 180/π]], s["deg"], s[dp]}, {s[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]}, {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]}, {s[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], s["GMm/c"], s[dp]}, {s[" Ř crc.φ"], " = ", s[n0[яi[tp]]], s["GM/c²"], s[dp]}, {s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[tp]]]], s["GM/c²"], s[dp]}, {s[" E kinet"], " = ", s[n0[ekin[tp]]], s["mc²"], s[dp]}, {s[" E poten"], " = ", s[n0[epot[tp]]], s["mc²"], s[dp]}, {s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]}, {s[" CarterQ"], " = ", s[n0[Q]], s["GMm/c"], s[dp]}, {s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]}, {s[" d\.b2r/dτ\.b2"], " = ", s[n0[Evaluate[r''[T] /. sol][[1]]]], s["c⁴/G/M"], s[dp]}, {s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Evaluate[φ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]}, {s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Evaluate[θ''[T] /. sol][[1]]]], s["c⁶/G²/M²"], s[dp]}, {s[" α dv/dτ"], " = ", s[n0[Evaluate[vlt''[tp]/.sol][[1]]]], s["c⁴/G/M"], s[dp]}, {s[" R carts"], " = ", s[n0[XYZ[tp]]], s["GM/c²"], s[dp]}, {s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]}, {s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]}, {s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]}, {s[" s dstnc"], " = ", s[n0[dst[tp]]], s["GM/c²"], s[dp]}, {s[" ω fdrag"], " = ", s[n0[ω[tp]]], s["c³/G/M"], s[dp]}, {s[" v fdrag"], " = ", s[n0[й[tp]]], s["c"], s[dp]}, {s[" Ω fdrag"], " = ", s[n0[Ω[tp]]], s["c"], s[dp]}, {s[" v propr"], " = ", s[n0[v[tp]/Sqrt[1-v[tp]^2]]], s["c"], s[dp]}, {s[" v obsvd"], " = ", s[n0[ß[tp]]], s["c"], s[dp]}, {s[" v escpe"], " = ", s[n0[ж[tp]]], s["c"], s[dp]}, {s[" v delay"], " = ", s[n0[sh[tp]]], s["c"], s[dp]}, {s[" v local"], " = ", s[n0[v[tp]]], s["c"], s[dp]}, {s[" "], s[" "], s[" "], s[" "]}}, Alignment-> Left, Spacings-> {0, 0}]; plot1b[{xx_, yy_, zz_, tk_, w1_, w2_}]:= (* Animation *) Show[Graphics3D[{ {PointSize[0.012], Red, Point[ Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}}, ImageSize-> imgsize, PlotRange-> PR, SphericalRegion->False, ImagePadding-> 1], horizons[A, None, w1, w2], If[crd==1, If[a==0, {}, Graphics3D[{{PointSize[0.009], Purple, Point[ Xyz[xyZ[{ Sin[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2], Cos[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2], z0[A]}, w1], w2]]}}]], If[a==0, {}, Graphics3D[{{PointSize[0.009], Purple, Point[ Xyz[xyZ[{ Sin[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2], Cos[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2], z0[A]}, w1], w2]]}}]]], If[crd==1, If[tk==0, {}, If[a==0, {}, Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ Xyz[xyZ[{ Sin[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2], Cos[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2], z0[A]}, w1], w2], {tt, Max[0, д[т[tp]-199/100 π/ω0]], tp}, PlotStyle -> {Thickness[0.001], Dashed, Purple}, PlotPoints-> 5+Round[2т[tp]], MaxRecursion-> 12]]]], If[tk==0, {}, If[a==0, {}, Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ Xyz[xyZ[{ Sin[-φ0-ц0 a Ξ/χ т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2], Cos[-φ0-ц0 a Ξ/χ т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2], z0[A]}, w1], w2], {tt, Max[0, д[т[tp]-199/100 π/ω0]], tp}, PlotStyle -> {Thickness[0.001], Dashed, Purple}, PlotPoints-> 5+Round[2т[tp]], MaxRecursion-> 12]]]]], If[tk==0, {}, Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, Max[0, tp-d1], tp}, PlotStyle-> {Thickness[0.004]}, ColorFunction-> Function[{x, y, z, t}, Hue[0, 1, 0.5, Max[Min[(-tp+(t+d1))/d1, 1], 0]]], ColorFunctionScaling-> False, PlotPoints-> Automatic, MaxRecursion-> mrec]]], ViewPoint-> {xx, yy, zz}]; Quiet[Do[ Print[Rasterize[Grid[{{ plot1b[{0, -Infinity, 0, tp, w1l, w2l}], plot1b[{0, 0, +Infinity, tp, w1r, w2r}], display[tp] }}, Alignment->Left]]], {tp, {0, tau/.FindRoot[R[tau]-rA, {tau, 1}], tau/.FindRoot[R[tau]-rI, {tau, 2}], tau/.FindRoot[R[tau], {tau, 2.27}], tMax}}]] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||| 13) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* Export als HTML Dokument *) (* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *) (* Export direkt als Bildsequenz *) (* Do[Export["dateiname" <> ToString[tk] <> ".png", Rasterize[...] ], {tk, 0, 10, 5}] *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||| http://kerr.yukerez.net ||||| Simon Tyran, Vienna ||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)