Seite 1 von 1

### Kerr Metric

Verfasst: Do 12. Apr 2018, 22:24

This is the english version.   Deutschsprachige Version auf kerr.yukterez.net und Yukipedia.

Here we use natural units of G=M=c=1, so lengths are in GM/c² and times in GM/c³. The metric signature is time-positive (+,-,-,-). a is the spin parameter (for black holes 0≤a≤M), M the mass equivalent of the total energy of the black hole, and Mirr its irreducible mass:

Shorthand terms:

Covariant metric coeffizients:

Contravariant components (superscripted letters are not powers, but indices):

The dimensionless spin parameter is a=Jc/G/M². Transformation into cartesian coordinates:

Line element in Boyer Lindquist coordinates:

Metric tensor (t,r,θ,Ф):

With the transformation:

where T is a finkelsteinlike time coordinate (radially infalling photons move with dr/dt=1) and ψ the flattened azimuthal angle:

the metric in Kerr Schild coordinates (T,r,θ,ψ) is:

Equations of motion in Boyer Lindquist coordinates:

Coordinate time by proper time (dt/dτ):

First proper time derivative of the radial coordinate (dr/dτ):

Derivative of the poloidial component of motion (dθ/dτ):

Derivative of the poloidial angular momentum (dpθ/dτ):

Axial angular momentum:

Derivative of the axial component of motion (dФ/dτ):

Axial angular momentum derivative (pФ/dτ):

Longitudinal component of the angular momentum:

Constant of motion, Carter's constant:

Constant of motion, Carter k:

Constant of motion, total energy:

Constant of motion, axial angular momentum:

Local 3-velocity component along the r-axis:

Local 3-velocity component along the θ-axis:

Local 3-velocity component along the Ф-axis:

Local 3-velocity, total:

For massive testparticles μ=-1 and for photons μ=-0. δ is the inclination angle. With α as the vertical launch anglel the components of the local velocity (relative to a ZAMO) are

Shapirodelayed and frame dragged velocity as observed at infinity:

Frame-Dragging angular velocity oberserved at infinity (dФ/dt):

Delayed Frame-Dragging transverse velocity at the equator of the outer horizon:

with the horizons and ergospheres (solution for r at Δ=0 and gtt=0):

r and θ dependend delayed Frame-Dragging transverse velocities:

at the equatorialen plane at θ=π/2:

r und θ dependend local Frame-Dragging transverse velocities (greater than c inside of the ergosphere):

at the equatorialen plane at θ=π/2:

Cartesian projection of the Frame-Dragging transverse velocity:

at the equatorialen plane at θ=π/2:

Gravitational time dilation component relative to a ZAMO (dt/dτ):

Axial and coaxial radius of gyration:

Axial and coaxial circumference:

For images and animations see the german version of this site.

### Kerr Metric

Verfasst: Mi 18. Apr 2018, 01:02

Simulator code, Boyer-Lindquist & Kerr-Schildcoordinates, particles & photons:

Code: Alles auswählen

`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* |||||||| Mathematica Syntax | http://kerr.yukterez.net | Version: 24.03.2018  |||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};mt2 = {"EventLocator", "Event"-> (r[τ]-1001/1000 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 = 1;                                    (* Boyer-Lindquist: crd=1, Kerr-Schild: crd=2 *)dsp = 1;                                                                 (* Display Modus *) tmax = 300;                                                                  (* Eigenzeit *)Tmax = 300;                                                            (* Koordinatenzeit *)TMax = Min[Tmax, т[plunge-1*^-3]]; tMax = Min[tmax, plunge];          (* Integrationsende *) r0 = 7;                                                                    (* Startradius *)θ0 = π/2;                                                                  (* Breitengrad *)φ0 = 0;                                                                     (* Längengrad *)a  = 9/10;                                                               (* Spinparameter *) v0 = 4/10;                                                      (* Anfangsgeschwindigkeit *)α0 = 0;                                                      (* vertikaler Abschusswinkel *)ψ0 = ArcTan[5/6];                                               (* 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   = Simplify[Limit[pθ0^2+(Lz^2 Csc[ϑ]^2-a^2 (ε^2+μ)) Cos[ϑ]^2, ϑ->θ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 *)rG = 1-Sqrt[1-a^2 Cos[θ]^2];                                         (* innere Ergosphäre *)rA = 1+Sqrt[1-a^2];                                                   (* äußerer Horizont *)rI = 1-Sqrt[1-a^2];                                                   (* innerer Horizont *) 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[A_, w1_, w2_] := Xyz[xyZ[{Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];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[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];ц       = 2r[τ]/Σ; ц0=2r0/Ξ; Σ       = 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 δ; т[τ_] := 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 *)Ω[τ_] := ω[τ] 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[Simplify[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*^-4, 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];   (* kartesisches R *) 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 *)   Mrec    = 100;                                           (* Parametric Plot Subdivisionen *)mrec    = 10; imgsize = 380;                                                               (* Bildgröße *)w1l     = 0;                                                 (* Startperspektiven, Winkel *)w2l     = 0;w1r     = 0;w2r     = 0; s[text_] := Style[text, FontSize->font]; font=11;                          (* Anzeigestil *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* |||||||| 11) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) display[T_] := Grid[{{s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},{s[" т coord"], " = ", s[n0[tk+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]},{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},{s[" γ total"], " = ", s[n0[If[dsp==1, γ[T]-2 If[crd==1, 0, R'[T] R[T]/(R[T]^2-2 R[T]+a^2)], γ[T]]]], s[If[dsp==1, "dt/dτ", "dт/dτ"]], s[dp]},{s[" ς gravt"], " = ", s[n0[ς[T]]], s["dt/dτ"], s[dp]},{s[" r coord"], " = ", s[n0[R[T]]], s["GM/c²"], s[dp]},{s[" φ longd"], " = ", s[n0[Φ[T] 180/π]], s["deg"], s[dp]},{s[" θ lattd"], " = ", s[n0[Θ[T] 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[T]]], s["GM/c²"], s[dp]},{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], s[dp]},{s[" E kinet"], " = ", s[n0[ekin[T]]], s["mc²"], s[dp]},{s[" E poten"], " = ", s[n0[epot[T]]], 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''[T]/.sol][[1]]]], s["c⁴/G/M"], s[dp]},{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},{s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c²"], s[dp]},{s[" ω fdrag"], " = ", s[n0[ω[T]]], s["c³/G/M"], s[dp]},{s[" v fdrag"], " = ", s[n0[й[T]]], s["c"], s[dp]},{s[" Ω fdrag"], " = ", s[n0[Ω[T]]], s["c"], s[dp]},{s[" v propr"], " = ", s[n0[v[T]/Sqrt[1-v[T]^2]]], s["c"], s[dp]},{s[" v obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},{s[" v escpe"], " = ", s[n0[ж[T]]], s["c"], s[dp]},{s[" v delay"], " = ", s[n0[sh[T]]], s["c"], s[dp]},{s[" v local"], " = ", s[n0[v[T]]], s["c"], s[dp]},{s[" "], s[" "], s["                   "], s["         "]}},Alignment-> Left, Spacings-> {0, 0}]; plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)Show[Graphics3D[{{PointSize[0.009], Red, Point[Xyz[xyZ[{x[T], y[T], z[T]}, 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 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],Cos[-φ0-ω0 tk+π/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 Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],Cos[-φ0-ц0 a Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],z0[A]}, w1], w2]]}}]]],If[crd==1, If[tk==0, {}, If[a==0, {},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, tk-199/100 π/ω0], tk},PlotStyle -> {Thickness[0.001], Dashed, Purple},PlotPoints-> Automatic,MaxRecursion-> mrec]]],If[tk==0, {}, If[a==0, {},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, tk-199/100 π/ω0], tk},PlotStyle -> {Thickness[0.001], Dashed, Purple},PlotPoints-> Automatic,MaxRecursion-> mrec]]]],If[tk==0, {},Block[{\$RecursionLimit = Mrec},ParametricPlot3D[Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, If[Tmax<0, Min[-1*^-16, T+d1/3], Max[1*^-16, T-d1/3]]},PlotStyle-> {Thickness[0.003], Gray},PlotPoints-> Automatic,MaxRecursion-> mrec]]],Block[{\$RecursionLimit = Mrec},If[tk==0, {},ParametricPlot3D[Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[TMax<0, Min[0, T+d1], Max[0, T-d1]], T},PlotStyle-> {Thickness[0.004]},ColorFunction-> Function[{x, y, z, t},Hue[0, 1, 0.5, If[TMax<0, Max[Min[(+T+(-t+d1))/d1, 1], 0], Max[Min[(-T+(t+d1))/d1, 1], 0]]]],ColorFunctionScaling-> False,PlotPoints-> Automatic,MaxRecursion-> mrec]]],ViewPoint-> {xx, yy, zz}]; Quiet[Do[Print[Rasterize[Grid[{{plot1a[{0, -Infinity, 0, tk, w1l, w2l}],plot1a[{0, 0, Infinity, tk, w1r, w2r}],display[Quiet[д[tk]]]}}, Alignment->Left]]],{tk, TMax, TMax, TMax}]] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* |||||||| 12) PLOT NACH EIGENZEIT ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) display[T_] := Grid[{{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},{s[" t coord"], " = ", s[n0[т[tp]]], s["GM/c³"], s[dp]},{s[" т coord"], " = ", 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[" γ total"], " = ", 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[If[dsp==1, "dt/dτ", "dт/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.009], 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, {},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-> Automatic,MaxRecursion-> 12]]],If[tk==0, {}, If[a==0, {},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-> Automatic,MaxRecursion-> 12]]]],If[tk==0, {},Block[{\$RecursionLimit = Mrec},ParametricPlot3D[Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},PlotStyle-> {Thickness[0.003], Gray},PlotPoints-> Automatic,MaxRecursion-> mrec]]],If[tk==0, {},Block[{\$RecursionLimit = Mrec},ParametricPlot3D[Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},PlotStyle-> {Thickness[0.004]},ColorFunction-> Function[{x, y, z, t},Hue[0, 1, 0.5, If[tp<0, Max[Min[(+tp+(-t+d1))/d1, 1], 0], 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, tMax, tMax, 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 ||||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)`

Solver for the initial conditions;converseion between local velocity, proper time derivatives and constants of motion:

Code: Alles auswählen

`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 1: Lokale Geschwindigkeit nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *) ε  = (72 Sqrt[3/2136769])/7+5 Sqrt[3581/105087];Lz = (2 Sqrt[105087/61])/35;Q  = 700/183; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 1"Reduce[ε0==ε && L0==Lz && Q0==Q && vr0^2==v0^2-vφ0^2-vθ0^2 && v0>0, {v0,vr0,vφ0,vθ0}]N[%](* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* || *)    (* || *)        (* || *)            (* || *)                 (* || *)                      (* || *)                          (* || *)                              (* || *)                                  (* ||*)                                      (* || *)                                          (* || *)                                              (* || *)                                                   (* || *)                                                        (* || *)                                                            (* || *)                                                                (* || *)                                                                    (* || *)                                                                        (* ||*)                                                                            (* || *)                                                                                (* || *)                                                                                    (* || *)                                                                                   (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 2: Lokale Geschwindigkeit nach ersten Eigenzeitableitungen  |||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *) dt = (5 Sqrt[35029/10743])/7;dr = 0;dθ = 10/(7 Sqrt[1281]);dφ = (300 Sqrt[3/125438849])/7+40 Sqrt[3/2136769]; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Ξ=r0^2+a^2 Cos[θ0]^2;δ=r0^2-2r0+a^2;j[v_]:=Sqrt[1+μ v^2];pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];pθ0=vθ0 Sqrt[Ξ]/j[v0]; dT=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);dR=(pr0 δ)/Ξ;dΘ=pθ0/Ξ;dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 2"Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ && v0>0, {v0,vr0,vφ0,vθ0}]N[%] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || *)    (* || *)        (* || *)            (* || *)                 (* || *)                      (* || *)                          (* || *)                              (* || *)                                  (* ||*)                                      (* || *)                                          (* || *)                                              (* || *)                                                   (* || *)                                                        (* || *)                                                            (* || *)                                                                (* || *)                                                                    (* || *)                                                                        (* ||*)                                                                            (* || *)                                                                                (* || *)                                                                                    (* || *)                                                                                   (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 3: Erste Eigenzeitableitungen nach Erhaltungsgrößen ε, Lz, Q ||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Erhaltungsgrößen Gesamtenergie, axialer Drehimpuls & Carter Konstante eingeben |||| *) ε  = (72 Sqrt[3/2136769])/7+5 Sqrt[3581/105087];Lz = (2 Sqrt[105087/61])/35;Q  = 700/183; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Ξ=r0^2+a^2 Cos[θ0]^2;δ=r0^2-2r0+a^2;j[v_]:=Sqrt[1+μ v^2];pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];pθ0=vθ0 Sqrt[Ξ]/j[v0]; dT=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);dR=(pr0 δ)/Ξ;dΘ=pθ0/Ξ;dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 3"Solve[ε==ε0 && Lz==L0 && Q==Q0 && dt==dT && dr==dR && dθ==dΘ && dφ==dΦ && v0^2==vr0^2+vθ0^2+vφ0^2, {dt,dr,dθ,dφ}]N[%] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || *)    (* || *)        (* || *)            (* || *)                 (* || *)                      (* || *)                          (* || *)                              (* || *)                                  (* ||*)                                      (* || *)                                          (* || *)                                              (* || *)                                                   (* || *)                                                        (* || *)                                                            (* || *)                                                                (* || *)                                                                    (* || *)                                                                        (* ||*)                                                                            (* || *)                                                                                (* || *)                                                                                    (* || *)                                                                                   (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 4: Erste Eigenzeitableitungen nach lokalen Geschwindigkeiten ||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *) vr0 = 0;vθ0 = 2/Sqrt[61];vφ0 = 12/(5 Sqrt[61]);v0  = Sqrt[vr0^2+vθ0^2+vφ0^2]; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Ξ=r0^2+a^2 Cos[θ0]^2;δ=r0^2-2r0+a^2;j[v_]:=Sqrt[1+μ v^2];pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];pθ0=vθ0 Sqrt[Ξ]/j[v0]; dT=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);dR=(pr0 δ)/Ξ;dΘ=pθ0/Ξ;dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 4"Reduce[dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {dt,dr,dθ,dφ}]N[%] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || *)    (* || *)        (* || *)            (* || *)                 (* || *)                      (* || *)                          (* || *)                              (* || *)                                  (* ||*)                                      (* || *)                                          (* || *)                                              (* || *)                                                   (* || *)                                                        (* || *)                                                            (* || *)                                                                (* || *)                                                                    (* || *)                                                                        (* ||*)                                                                            (* || *)                                                                                (* || *)                                                                                    (* || *)                                                                                   (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 5: Erhaltungsgrößen ε, Lz, Q nach lokalen Geschwindigkeiten |||||||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *) vr0 = 0;vθ0 = 2/Sqrt[61];vφ0 = 12/(5 Sqrt[61]);v0  = Sqrt[vr0^2+vθ0^2+vφ0^2]; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Ξ=r0^2+a^2 Cos[θ0]^2;δ=r0^2-2r0+a^2;j[v_]:=Sqrt[1+μ v^2];pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];pθ0=vθ0 Sqrt[Ξ]/j[v0]; dT=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);dR=(pr0 δ)/Ξ;dΘ=pθ0/Ξ;dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 5"Reduce[ε==ε0 && Lz==L0 && Q==Q0, {ε,Lz,Q}]N[%] (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || *)    (* || *)        (* || *)            (* || *)                 (* || *)                      (* || *)                          (* || *)                              (* || *)                                  (* ||*)                                      (* || *)                                          (* || *)                                              (* || *)                                                   (* || *)                                                        (* || *)                                                            (* || *)                                                                (* || *)                                                                    (* || *)                                                                        (* ||*)                                                                            (* || *)                                                                                (* || *)                                                                                    (* || *)                                                                                   (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || CODE 6: Erhaltungsgrößen ε, Lz, Q  nach ersten Eigenzeitableitungen  |||||||||||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) ClearAll["Global`*"] (* || Startposition etc. eingeben  |||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 7;θ0 = π/2;φ0 = 0;a  = 9/10;μ  =-1; (* || Startwerte für die ersten Eigenzeitableitungen eingeben ||||||||||||||||||||||||||| *) dt = (5 Sqrt[35029/10743])/7;dr = 0;dθ = 10/(7 Sqrt[1281]);dφ = (300 Sqrt[3/125438849])/7+40 Sqrt[3/2136769]; (* || Gleichungen für Gesamtenergie, axialer Drehimpuls & Carter Konstante  ||||||||||||| *) ε0 = (Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))/Sqrt[1+μ v0^2];L0 = (vφ0 Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/Sqrt[1+μ v0^2];Q0 = (vθ0^2 (r0^2+a^2 Cos[θ0]^2)+(vφ0^2 Cos[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/(r0^2+a^2 Cos[θ0]^2)-a^2 Cos[θ0]^2 (-1+v0^2+(Sqrt[((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)]+(2 a r0 vφ0 Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]))^2))/(1+μ v0^2); (* || Benötigte Gleichungen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Ξ=r0^2+a^2 Cos[θ0]^2;δ=r0^2-2r0+a^2;j[v_]:=Sqrt[1+μ v^2];pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];pθ0=vθ0 Sqrt[Ξ]/j[v0]; dT=ε0+(2r0(r0^2+a^2)ε0-2 a r0 L0)/(Ξ δ);dR=(pr0 δ)/Ξ;dΘ=pθ0/Ξ;dΦ=(2 a r0 ε0+(Ξ-2 r0)L0 Csc[θ0]^2)/(Ξ δ); (* || Output: lokale Geschwindigkeitskomponenten  ||||||||||||||||||||||||||||||||||||||| *) "Code 6"Solve[ε==ε0 && Lz==L0 && Q==Q0 && dT==dt && dR==dr && dΘ==dθ && dΦ==dφ, {ε,Lz,Q}]N[%](* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* ||||| Mathematica Syntax |||| http://kerr.yukterez.net |||| Simon Tyran, Vienna  ||||| *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)`

Derivation from the line element, Boyer-Lindquist coordinates:

Code: Alles auswählen

`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* | Mathematica Syntax | GEODESIC SOLVER | geodesics.yukterez.net | Version 25.09.2017 | *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)ClearAll["Global`*"]℧=0; q=0;                                               (* Input: kovariante metrische Komponenten *)gtt=1-(2r-℧^2)/Σ;grr=-Σ/Δ;gθθ=-Σ;gφφ=-Χ/Σ Sin[θ]^2;gtφ=a (2r-℧^2) Sin[θ]^2/Σ;                                                                           (* Abkürzungen *)Σ=r^2+a^2 Cos[θ]^2;Δ=r^2-2 r+a^2+℧^2;Χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;                    (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung *)x={t, r, θ, φ}; n=4; P=0; Ω=℧;                                                                         "Metrischer Tensor"metrik={{gtt, 0, 0, gtφ}, {0, grr, 0, 0}, {0, 0, gθθ, 0}, {gtφ, 0, 0, gφφ}};Subscript["g", μσ] -> MatrixForm[metrik]inversemetrik=Simplify[Inverse[metrik]];"g"^μσ -> MatrixForm[inversemetrik]                                                                            "Maxwell Tensor"A={(Ω r)/Σ+P/Σ Cos[θ] a, 0, 0, -(Ω r/Σ) Sin[θ]^2 a-P/Σ Cos[θ](r^2+a^2)};F=Simplify[Table[((D[A[[j]], x[[k]]]-D[A[[k]], x[[j]]])), {j, 1, 4}, {k, 1, 4}], Reals];Subscript["F", μσ] -> MatrixForm[F]f=FullSimplify[Table[Sum[inversemetrik[[i, k]] inversemetrik[[j, l]] F[[k, l]], {k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];"F"^μσ -> MatrixForm[f]                                                                           "Einstein Tensor"G=FullSimplify[Table[Sum[-2(inversemetrik[[k,l]] F[[i, k]] F[[j, l]] - metrik[[i, j]] F[[k, l]] f[[k, l]]), {k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];Subscript["G", μσ] -> MatrixForm[G]H=FullSimplify[Table[Sum[inversemetrik[[i, k]] inversemetrik[[j, l]] G[[k, l]], {k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals];"G"^μσ -> MatrixForm[H]                                                                       "Christoffel Symbole"christoffel:=Simplify[Table[(1/2)Sum[(inversemetrik[[i, s]])(D[metrik[[s, j]], x[[k]]]+D[metrik[[s, k]], x[[j]]] -D[metrik[[j, k]], x[[s]]]), {s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}]];christoffelsymbole=Table[If[UnsameQ[christoffel[[i, j, k]], 0], {ToString[Γ[i, j, k]], christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}];rplc[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]-> t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]->θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ];rple[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->  t'[τ])/. r\.b4->r'[τ])/. θ\.b4->θ'[τ])/. φ\.b4->φ'[τ];list[y_]:=y[[1]]==y[[2]];list[christoffelsymbole[[1]][[2]][[1]]]list[christoffelsymbole[[1]][[3]][[1]]]list[christoffelsymbole[[1]][[4]][[2]]]list[christoffelsymbole[[1]][[4]][[3]]]list[christoffelsymbole[[2]][[1]][[1]]]list[christoffelsymbole[[2]][[2]][[2]]]list[christoffelsymbole[[2]][[3]][[2]]]list[christoffelsymbole[[2]][[3]][[3]]]list[christoffelsymbole[[2]][[4]][[1]]]list[christoffelsymbole[[2]][[4]][[4]]]list[christoffelsymbole[[3]][[1]][[1]]]list[christoffelsymbole[[3]][[2]][[2]]]list[christoffelsymbole[[3]][[3]][[2]]]list[christoffelsymbole[[3]][[3]][[3]]]list[christoffelsymbole[[3]][[4]][[1]]]list[christoffelsymbole[[3]][[4]][[4]]]list[christoffelsymbole[[4]][[2]][[1]]]list[christoffelsymbole[[4]][[3]][[1]]]list[christoffelsymbole[[4]][[4]][[2]]]list[christoffelsymbole[[4]][[4]][[3]]]                                                                      "Bewegungsgleichungen"geodäsie=Simplify[Table[-Sum[christoffel[[i, j, k]] x[[j]]' x[[k]]'+q f[[i, k]] x[[j]]' metrik[[j, k]], {j, 1, n}, {k, 1, n}], {i, 1, n}]];bewegungsgleichung=Table[{x[[i]]''[τ]==FullSimplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];bewegungsgleichung[[1]][[1]]bewegungsgleichung[[2]][[1]]bewegungsgleichung[[3]][[1]]bewegungsgleichung[[4]][[1]]ClearAll[Σ, Δ, Χ];                                                                            "Zeitdilatation"t'[τ]==rple[Simplify[Normal[Solve[-μ==gtt Т^2+grr r\.b4^2+gθθ θ\.b4^2+gφφ φ\.b4^2 + 2 gtφ Т φ\.b4, Т, Reals]]][[2]][[1]][[2]]]`

Derivation from the line element, Kerr-Schild coordinates:

Code: Alles auswählen

`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* | Mathematica Syntax | GEODESIC SOLVER | geodesics.yukterez.net | Version 25.09.2017 | *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)ClearAll["Global`*"]                                               (* Input: kovariante metrische Komponenten *)gtt=-(ц-1);gtr=-ц;gtφ=ц a Sin[θ]^2; gφφ=-χ/Σ Sin[θ]^2;grφ=a(1+ц)Sin[θ]^2; grr=-(1+ц);gθθ=-Σ;                    (* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung *)x={t, r, θ, φ}; n=4;                                                                         "Metrischer Tensor"metrik={{gtt, gtr, 0, gtφ}, {gtr, grr, 0, grφ}, {0, 0, gθθ, 0}, {gtφ, grφ, 0, gφφ}};Subscript["g", μσ] -> MatrixForm[metrik]inversemetrik=Simplify[Inverse[metrik]];"g"^μσ -> MatrixForm[inversemetrik]                                                                           (* Abkürzungen *)Σ=r^2+a^2 Cos[θ]^2;Δ=r^2-2 r+a^2;χ=(r^2+a^2)^2-a^2 Sin[θ]^2 Δ;ц=2r/Σ;                                                                       "Christoffel Symbole"christoffel:=Simplify[Table[(1/2)Sum[(inversemetrik[[i, s]])(D[metrik[[s, j]], x[[k]]]+D[metrik[[s, k]], x[[j]]] -D[metrik[[j, k]], x[[s]]]), {s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}]];christoffelsymbole=Table[If[UnsameQ[christoffel[[i, j, k]], 0], {ToString[Γ[i, j, k]], christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}]rplc[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]-> t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]->θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ];rple[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->  t'[τ])/. r\.b4->r'[τ])/. θ\.b4->θ'[τ])/. φ\.b4->φ'[τ];                                                                      "Bewegungsgleichungen"geodäsie=Simplify[Table[-Sum[christoffel[[i, j, k]] x[[j]]' x[[k]]', {j, 1, n}, {k, 1, n}], {i, 1, n}]];bewegungsgleichung=Table[{x[[i]]''[τ]==FullSimplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];bewegungsgleichung[[1]][[1]]bewegungsgleichung[[2]][[1]]bewegungsgleichung[[3]][[1]]bewegungsgleichung[[4]][[1]]ClearAll[Σ, Δ, χ, ц];                                                                            "Zeitdilatation"t'[τ]==rple[FullSimplify[Normal[Solve[-μ==gtt Т^2+grr r\.b4^2+gθθ θ\.b4^2+gφφ φ\.b4^2+2 gtr Т r\.b4+2 gtφ Т φ\.b4+2 grφ r\.b4 φ\.b4,Т,Reals]]][[2]][[1]][[2]]]`

`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || Mathematica Syntax | BLACK HOLE SHADOWS | kerr.yukterez.net | Version 03.10.2017 || *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)rE=M+Sqrt[M^2-a^2 Cos[θ]^2]; re={Sqrt[rE^2+a^2] Sin[θ],rE Cos[θ]};rG=M-Sqrt[M^2-a^2 Cos[θ]^2]; rg={Sqrt[rG^2+a^2] Sin[θ],rG Cos[θ]};rA=M+Sqrt[M^2-a^2]; ra={Sqrt[rA^2+a^2] Sin[θ],rA Cos[θ]};rI=M-Sqrt[M^2-a^2]; ri={Sqrt[rI^2+a^2] Sin[θ],rI Cos[θ]};A=Interpolation[{{0,0},{0.1,0.21},{0.2,0.4},{0.3,0.61},{0.4,0.82},{0.5,1.01},{0.6,1.24},{0.7,1.48},{0.8,1.71},{0.866,1.9},{0.9,2.01},{0.95,2.17},{0.98,2.3},{0.99,2.38},{0.999,2.47},{1,2.5}},InterpolationOrder -> 1]; B=Interpolation[{{0,Sqrt[27]},{0.1,5.19},{0.2,5.18},{0.3,5.16},{0.4,5.14},{0.5,5.11},{0.6,5.07},{0.7,5.02},{0.8,4.93},{0.866,4.86},{0.9,4.82},{0.95,4.74},{0.98,4.68},{0.99,4.65},{0.999,4.62},{1,4.6}},InterpolationOrder -> 1]; pp:=ParametricPlot[{re,rg,ra,ri},{θ,0,2π},Frame->True,ImageSize->400,PlotRange->{{-8,8},{-6,6}},PlotStyle->{{Black,Thick}}]; cp:=ContourPlot[(x^2+y^2+A[a] x)^2-B[a]^2 (x^2+y^2)==0,{X,-8,8},{y,-6,6},ContourStyle->{Red,Thick}]; a=Sin[w]; M=1; x=-X; Do[Print[Rasterize[Grid[{{Show[pp,cp]},{" a"->N[a]}},Alignment->Left]]],{w,0,π/2,π/200}]`
`(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)(* || Mathematica Syntax | BLACK HOLE SHADOWS | kerr.yukterez.net | Version 03.10.2017 || *)(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)Manipulate[rE=1+Sqrt[1-a^2 Cos[θ]^2]; re={Sqrt[rE^2+a^2] Sin[θ] Cos[φ], Sqrt[rE^2+a^2] Sin[θ] Sin[φ], rE Cos[θ]};rG=1-Sqrt[1-a^2 Cos[θ]^2]; rg={Sqrt[rG^2+a^2] Sin[θ] Cos[φ], Sqrt[rG^2+a^2] Sin[θ] Sin[φ], rG Cos[θ]};rA=1+Sqrt[1-a^2]; ra={Sqrt[rA^2+a^2] Sin[θ] Cos[φ], Sqrt[rA^2+a^2] Sin[θ] Sin[φ], rA Cos[θ]};rI=1-Sqrt[1-a^2]; ri={Sqrt[rI^2+a^2] Sin[θ] Cos[φ], Sqrt[rI^2+a^2] Sin[θ] Sin[φ], rI Cos[θ]};rot[{x_, y_, z_}, w_]:={x, y Sin[w]-z Cos[w], y Cos[w]+z Sin[w]};pp:=Show[ParametricPlot3D[rot[re, ϑ], {φ, 0, 2 π}, {θ, 0, π},Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.10]],ImageSize->400, PlotRange->{{-8, 8}, {-8, 8}, {-6, 6}}, ViewPoint->{0, 20, 0}, ImagePadding->1],ParametricPlot3D[rot[ra, ϑ], {φ, 0, 2 π}, {θ, 0, π},Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.15]]],ParametricPlot3D[rot[ri, ϑ], {φ, 0, 2 π}, {θ, 0, π},Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.25]]],ParametricPlot3D[rot[rg, ϑ], {φ, 0, 2 π}, {θ, 0, π},Mesh -> None, PlotPoints -> 20, PlotStyle -> Directive[Gray, Opacity[0.35]]]];                       α=Interpolation[{{0,0},{0.1,0.21},{0.2,0.4},{0.3,0.61},{0.4,0.82},{0.5,1.01},{0.6,1.24},{0.7,1.48},{0.8,1.71},{0.866,1.9},{0.9,2.01},{0.95,2.17},{0.98,2.3},{0.99,2.38},{0.999,2.47},{1,2.5}},InterpolationOrder -> 1]; β=Interpolation[{{0,Sqrt[27]},{0.1,5.19},{0.2,5.18},{0.3,5.16},{0.4,5.14},{0.5,5.11},{0.6,5.07},{0.7,5.02},{0.8,4.93},{0.866,4.86},{0.9,4.82},{0.95,4.74},{0.98,4.68},{0.99,4.65},{0.999,4.62},{1,4.6}},InterpolationOrder -> 1]; A[a_]:=α[a] Sin[ϑ]+a Sin[ϑ]^3 Cos[ϑ]^2/5;B[a_]:=β[a]+0.23 Cos[ϑ]^4(1-Sqrt[1-a^4]);cp:=ContourPlot3D[(x^2+y^2+A[a] x)^2-B[a]^2 (x^2+y^2)==0, {x, -8, 8}, {z, -3, 3}, {y, -6, 6}, ContourStyle->{Black, Thick}, PlotPoints -> 20];Rasterize[Grid[{{Show[pp, cp]}, {" a"->N[a]}}, Alignment->Left]], {{ϑ, π/2}, 0, π/2}, {{a, 1}, 0, 1}]`