Flammsches Paraboloid

Deutschsprachige Version
Benutzeravatar
Yukterez
Administrator
Beiträge: 194
Registriert: Mi 21. Okt 2015, 02:16

Flammsches Paraboloid

Beitragvon Yukterez » So 21. Jan 2018, 17:02

Bild
Bild

Radiale und zeitliche metrische Komponenten für die innere und äußere Schwarzschildlösung; ist die Ausdehnung der Masse in umfanggetreuen Schwarzschildkoordinaten:



Auf die Fläche als Kurvenlänge projizierte physikalische Radialdistanz :



Zeitdilatation (Eigenzeit des lokalen Schalenbeobachters durch Koordinatenzeit des feldfreien Buchhalters):


Bild

Plot A1: Isometrische Einbettung der Radialdistanz; die blaue Kurvenlänge im 1D Plot () stellt den tatsächlichen radialen Abstand im Verhältnis zum auf der -Achse dargestellten Umfangradius dar. Die -Achse ist eine räumliche Hilfsdimension.

Plot A2: Die rote Funktion () zeigt die Zeitdilation eines lokalen Schalenbeobachters; im inneren des schwarzen Lochs läuft die Eigenzeit eines stationären Schalenbeobachters relativ zur Koordinatenzeit des externen Koordinatenbuchalters wieder rückwärts - was ein mathematisches Artefakt ist, da es innerhalb des Horizonts keine stationären Beobachter mehr gibt, und es außerdem bereits eine unendliche Koordinatenzeit (wenn auch endliche Eigenzeit) dauert damit die Masse kleiner als ihr Ereignishorizont werden kann.

Plot A; Animationsparameter: Massenausdehnung :


Bild

Standbild A-Ⅰ: Ausdehnung der Masse bis zum zehnfachen Schwarzschildradius:

Bild

Standbild A-Ⅱ: Ausdehnung der Masse bis zum dreifachen Schwarzschildradius:

Bild

Standbild A-Ⅲ: Ausdehnung der Masse bis zur Photonensphäre bei :

Bild

Standbild A-Ⅳ: ; ab dieser Ausdehnung ist der Gravitationskollaps unvermeidlich, unendliche Zeitdilatation im Zentrum:

Bild

Standbild A-Ⅴ: ; der umfanggetreue Koordinatenradius entspricht dem metrischen Radius , unendliche ZD am Horizont:

Bild

Bild

Plot B zeigt das Paraboloid aus der umfanggetreuen 2D () Vogelperspektive; alle Schalen haben einen physikalischen Abstand von Animationsparameter: Massenausdehnung :


Bild

Standbild B-Ⅰ: euklidscher Raum, keine Masse:

Bild

Standbild B-Ⅱ: Ausdehnung bis :

Bild

Standbild B-Ⅲ: Ausdehnung bis zur Photonensphäre bei :

Bild

Standbild B-Ⅳ: Ausdehnung der Masse konvergiert auf den Ereignishorizont zu:

Bild

Standbild B-Ⅴ: Vakuumlösung (); die gesamte Masse befindet sich in der Singularität im Zentrum:

Bild

Bild

Diese Seite ist ein Unterkapitel von schwarzschild.yukterez.net; eine genauere Beschreibung und weitere Bilder folgen demnächst. Code (Mathematica Syntax):

Code: Alles auswählen

(* schwarzschild.yukterez.net || flamm.yukterez.net *)

set= {"GlobalAdaptive", "MaxErrorIncreases" -> 100, Method -> "GaussKronrodRule"};

grr[r_, R_]:=If[Abs[r]<R, (1-(2Abs[r]^2)/R^3)^-1, (1-2/Abs[r])^-1 ];
gtt[r_, R_]:=If[Abs[r]<R, 4(3Sqrt[1-2/R]-Sqrt[1-(2Abs[r]^2)/R^3])^-2, (1-2/Abs[r])^-1];

Ȓ[я_, R_]:=NIntegrate[Sqrt[Abs[grr[r, R]]], {r, 0, Abs[я]}, Method -> set, MaxRecursion -> 100];
Ř[я_, R_]:=NIntegrate[Sqrt[Abs[grr[r, R]-1]], {r, Abs[я], 24}, Method -> set, MaxRecursion -> 100];
ř[я_, R_]:=Ř[0, R]-Ř[я, R];
г[я_, R_]:=x/.FindRoot[Ȓ[x, R]-я, {x, 1}];

ť[я_, R_]:=1/Sqrt[gtt[я, R]];
ṫ[я_, R_]:=If[R>2.25, Limit[ť[r, R], r->я],
If[Abs[я]<f/.FindRoot[ť[f, R], {f, R}], -Limit[ť[r, R], r->я], Limit[ť[r, R], r->я]]];

s[text_]:=Style[text, FontSize->font];font=11;

plot[R_]:=Grid[{{
s["   \[Square]: flamm's paraboloid, \[Square]: time dilation dт/dt, r"]},
{Show[
Plot[ř[я, R]-ř[24, R], {я, -24, 24}, GridLines->{{-R, R}, {0, ř[R, R]-ř[24, R], ř[0, R]-ř[24, R]}},
PlotRange->{{-24, 24}, {1, -16}}, Frame->True, PlotTheme->"Classic", PlotStyle->{Blue, Thick},
 AspectRatio->17/48, ImageSize->516, ImagePadding->18],
Graphics[{
{Cyan, Opacity[0.1], Rectangle[{-R, 1}, {R, -19}]},
{Magenta, Opacity[0.1], Rectangle[{-R, 1}, {-24, -19}]},
{Magenta, Opacity[0.1], Rectangle[{R, 1}, {24, -19}]}}]]},
{Show[
Plot[ṫ[я, R], {я, -24, 24}, GridLines->{{-R, R}, {0, 1, ṫ[R, R], ṫ[0, R]}},
PlotRange->{{-24, 24}, {1.1, -0.6}}, Frame->True, PlotTheme->"Classic", PlotStyle->{Red, Thick},
 AspectRatio->17/48/2, ImageSize->516, ImagePadding->18],
Graphics[{
{Cyan, Opacity[0.1], Rectangle[{-R, -1}, {R, 18}]},
{Magenta, Opacity[0.1], Rectangle[{-R, -1}, {-24, 18}]},
{Magenta, Opacity[0.1], Rectangle[{R, -1}, {24, 18}]}}]]
}, {Grid[{{
s["   r(M)"->N[R]], s["R(M)"->N[Ȓ[R,R]]]
},{"               ", " "}},
Alignment->Left]}}, Alignment->Left];

(* Plot A *)
Do[Print[Quiet[Rasterize[plot[Max[R, 2]]]]], {R, 2, 4, 1/5}]

(* Plot B *)
Do[Print[Quiet[Rasterize[Grid[{{
Show[
Plot[{5}, {x, -4.8, +4.8},
Frame->True, PlotRange->{{-4.8, +4.8}, {-4.8, +4.8}},
Frame->True, PlotTheme->"Classic", AspectRatio->1,
ImageSize->516, ImagePadding->18],
Graphics[{Table[Circle[{0, 0}, г[n, R]], {n, 1/5, 20, 1/5}]}],
Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {R, 8}]}],
Graphics[{White, Opacity[0.7], Disk[{0, 0}, R]}],
Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, R]}],
Graphics[{Blue, Thick, Dashed, Opacity[0.9], Circle[{0, 0}, R]}]
]}, {Grid[{{
s["   r(M)"->N[R]], s["R(M)"->N[Ȓ[R,R]]]
},{"               ", " "}},
Alignment->Left]
}}, Alignment->Left]
]]], {R, 2, 4, 1/5}]
Bild
Симон Тыран @ vk || wikipedia || stackexchange || wolfram

Benutzeravatar
Yukterez
Administrator
Beiträge: 194
Registriert: Mi 21. Okt 2015, 02:16

Flammsches Paraboloid für die Kerr Metrik

Beitragvon Yukterez » Di 27. Feb 2018, 02:19

Bild
Wie oben, nur mit der grr-Komponente der Kerr-Metrik in Boyer-Lindquist Koordinaten:



Kartesischer Radius auf der äquatorialen Ebene:



Einbettung in die Ř,w-Fläche:



Irreduzible Masse:



Die Spriralen werden durch die Pfade radial einfallender Photonen symbolisiert. Äquatorialer Querschnitt, kartesische Projektion; alle a/M von 0..0.99999, kombinierte und separate Ansicht:

Bild

Das transversale Grid zeigt Schalen die zueinander einen physikalischen Abstand von ΔR=0.5Gℳ/c² haben, wobei ℳ hier für die irreduzible Masse steht (der äquatoriale Ereignishorizont liegt in dieser Maßeinheit immer bei 2Gℳ/c²). Das radiale Grid kann anstatt durch einfallende Photonen auch durch Freifaller aus dem Unendlichen gezeichnet werden. Dafür wird v0 im folgenden ausklappbaren Code von 1 auf ж0 geändert (Vergleich siehe hier). Zoom auf das Thorn'sche Limit von a/M=0.998: klick, Code:
Bild
Bild
1D (r) auf 2D (r,w) Einbettung, Anzeige im Programmfenster:

Code: Alles auswählen

(* Syntax: Wolfram || kerr.yukterez.net || flamm.yukterez.net *)

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;

A=a;                                             (* pseudosphärisch: A=0, kartesisch: A=a *)
crd=1;                                       (* Boyer-Lindquist: crd=1, Kerr-Schild crd=2 *)
dsp=1;
tmax=120;                                                                    (* Eigenzeit *)
Tmax=120;                                                              (* Koordinatenzeit *)
TMax=Min[Tmax, т[plunge-1*^-4]]; tMax=Min[tmax, plunge];              (* Integrationsende *)

r0=Sqrt[2]10;                                                              (* Startradius *)
θ0=Pi/2;                                                                   (* Breitengrad *)
φ0=0;                                                                       (* Längengrad *)
a=998/1000;                                                              (* Spinparameter *)
μ=-0;                                                        (* Baryon: μ=-1, Photon: μ=0 *)

v0=1;                                                           (* Anfangsgeschwindigkeit *)
α0=-Pi/2;                                                    (* vertikaler Abschusswinkel *)
ψ0=0;                                                           (* Bahninklinationswinkel *)

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[δ \[CapitalXi]/\[Chi]]/j[v0]+Lz ω0;           (* Energie und Drehimpulskomponenten *)
Lz=vφ0 Ы/j[v0];
pθ0=vθ0 Sqrt[\[CapitalXi]]/j[v0];
pr0=vr0 Sqrt[(\[CapitalXi]/δ)/j[v0]^2];
Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];        (* Carter Konstante *)
k=Q+Lz^2+a^2 (ε^2+μ);                                                         (* Carter k *)

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<=\[Pi]&&δ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 - *)

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];

horizons[A_, mesh_, w1_, w2_]:=Show[
ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 \[Pi]}, {θ, 0, \[Pi]},
Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 \[Pi]}, {θ, 0, \[Pi]},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 \[Pi]}, {θ, 0, \[Pi]},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 \[Pi]}, {θ, 0, \[Pi]},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
BLKS:=Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];

j[v_]:=Sqrt[1+μ v^2];                                                    (* Lorentzfaktor *)
mirr=Sqrt[(Sqrt[1-a^2]+1)/2];                                        (* irreduzible Masse *)
я=Sqrt[\[CapitalChi]/Σ]Sin[θ[τ]];                                (* axialer Umfangsradius *)
яi[τ_]:=Sqrt[\[CapitalChi]i[τ]/Σi[τ]]Sin[Θ[τ]];
Ы=Sqrt[\[Chi]/\[CapitalXi]]Sin[θ0];
Σ=r[τ]^2+a^2 Cos[θ[τ]]^2;                                    (* poloidialer Umfangsradius *)
Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
\[CapitalXi]=r0^2+a^2 Cos[θ0]^2;
\[CapitalDelta]=r[τ]^2-2r[τ]+a^2;
\[CapitalDelta]i[τ_]:=R[τ]^2-2R[τ]+a^2;
δ=r0^2-2r0+a^2;
\[CapitalChi]=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 \[CapitalDelta];
\[CapitalChi]i[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 \[CapitalDelta]i[τ];
\[Chi]=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
ц=2r[τ]/Σ; ц0=2r0/\[CapitalXi];

т[τ_]:=Evaluate[t[τ]/.sol][[1]];                        (* Koordinatenzeit nach Eigenzeit *)
д[\[Xi]_] :=Quiet[\[CapitalXi] /.FindRoot[т[\[CapitalXi]]-\[Xi], {\[CapitalXi], 0}]];
                                                        (* Eigenzeit nach Koordinatenzeit *)
T :=Quiet[д[tk]];                           
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*)
ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
γ[τ_]:=ю[τ];                                                                 (* totale ZD *)
R[τ_]:=Evaluate[r[τ]/.sol][[1]];                                (* Boyer-Lindquist Radius *)
\[CapitalPhi][τ_]:=Evaluate[φ[τ]/.sol][[1]];                               
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
ß[τ_]:=Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ]];

\[FinalSigma][τ_]:=Sqrt[\[CapitalChi]i[τ]/\[CapitalDelta]i[τ]/Σi[τ]]; \[FinalSigma]0=Sqrt[\[Chi]/δ/\[CapitalXi]];
                                                                        (* gravitative ZD *)
ω[τ_]:=2R[τ] a/\[CapitalChi]i[τ]; ω0=2r0 a/\[Chi]; ωd=2r[τ] a/\[CapitalChi];
                                                  (* Frame Dragging Winkelgeschwindigkeit *)
Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];            (* Frame Dragging beobachtete Geschwindigkeit *)
й[τ_]:=ω[τ] яi[τ] \[FinalSigma][τ]; й0=ω0 Ы \[FinalSigma]0;
                                                 (* Frame Dragging lokale Geschwindigkeit *)

ж[τ_]:=Sqrt[\[FinalSigma][τ]^2-1]/\[FinalSigma][τ]; ж0=Sqrt[\[FinalSigma]0^2-1]/\[FinalSigma]0;
                                                                 (* Fluchtgeschwindigkeit *)
     
vd[τ_]:=Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2-
        r[τ]^4(ε-Lz ωd)^2+\[CapitalDelta](Σ+a^2 Sin[θ[τ]]^2 (ε-
        Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+
        a^2 Sin[θ[τ]]^2 \[CapitalDelta]](ε - 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[τ]/\[CapitalDelta]i[τ] R'[τ];
sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
epot[τ_]:=ε+μ-ekin[τ];                                             (* potentielle Energie *)
ekin[τ_]:=If[μ==0, \[FinalSigma][τ], 1/Sqrt[1-v[τ]^2]-1];           (* kinetische Energie *)

dp= Y^Y; n0[z_]:=Chop[N[z]];

(* Boyer-Lindquist-Koordinaten *)

pr2τ[τ_]:=1/(Σ \[CapitalDelta]) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] \[CapitalDelta] μ+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)/(Σ \[CapitalDelta]),
t[0]==0,
r'[τ]==(pr[τ] \[CapitalDelta])/Σ,
r[0]==r0,
θ'[τ]==0,
θ[0]==θ0,
φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ \[CapitalDelta]),
φ[0]==φ0,
pr'[τ]==pr2τ[τ],
pr[0]==pr0,
pθ'[τ]==0,
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 δ)/\[CapitalXi]; dθ=pθ0/\[CapitalXi]; dφ=(2a r0 ε+(\[CapitalXi]-2r0)Lz Csc[θ0]^2)/(\[CapitalXi] δ)+dr a/δ; d\[CapitalPhi]=If[θ0==0, 0, If[θ0==\[Pi], 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\[CapitalPhi]))/(-1+ц0)+\[Sqrt]((1/((-1+ц0)^2 \[CapitalXi]))(\[CapitalXi] (dr^2+(-1+ц0) (μ-\[CapitalXi] dθ^2))+Sin[θ1]^2 d\[CapitalPhi] (-2a \[CapitalXi] dr-(-1+ц0) \[Chi] d\[CapitalPhi]+ц0^2 a^2 \[CapitalXi] Sin[θ1]^2 d\[CapitalPhi]))), θ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\[CapitalPhi],
φ[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];

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])}];

x[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]+n Pi/180]/.sol][[1]];
y[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]+n Pi/180]/.sol][[1]];
z[τ_]:=0;

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]];

xKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
yKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];

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_}, \[Beta]_]:={x Cos[\[Beta]]+z Sin[\[Beta]], y, z Cos[\[Beta]]-x Sin[\[Beta]]};
xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};

VP={r0, r0, r0};                                                     (* Perspektive x,y,z *)
d1=1;                                                                    (* Schweiflänge *)
plp=50;                                                            (* Flächenplot Details *)
w1l=0; w2l=0; w1r=0; w2r=0;                                          (* Startperspektiven *)
Mrec=100; mrec=10;                                       (* Parametric Plot Subdivisionen *)
imgsize=380;                                                                 (* Bildgröße *)

s[text_]:=Style[text, FontSize->font]; font=11;                            (* Anzeigestil *)

display[T_]:=Grid[{
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c\.b3"], s[dp]},
{s[" т coord"], " = ", s[n0[tk]], s["GM/c\.b3"], s[dp]},
{s[" γ total"], " = ", s[n0[γ[T]-If[crd==1, 0, -2 R'[T] R[T]/\[CapitalDelta]i[T]]]], s["dt/dτ"], s[dp]},
{s[" \[FinalSigma] gravt"], " = ", s[n0[\[FinalSigma][T]]], s["dt/dτ"], s[dp]},
{s[" r coord"], " = ", s[n0[R[T]]], s["GM/c\.b2"], s[dp]},
{s[" φ longd"], " = ", s[n0[\[CapitalPhi][T] 180/\[Pi]]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[T] 180/\[Pi]]], s["deg"], s[dp]},
{s[" M irred"], " = ", s[n0[mirr]], s["M"], s[dp]},
{s[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c\.b2"], s[dp]},
{s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c\.b2"], s[dp]},
{s[" E kinet"], " = ", s[n0[ekin[T]]], s["mc\.b2"], s[dp]},
{s[" E poten"], " = ", s[n0[epot[T]]], s["mc\.b2"], s[dp]},
{s[" E total"], " = ", s[n0[ε]], s["mc\.b2"], s[dp]},
{s[" CarterQ"], " = ", s[N[Q]], s["GMm/c"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM\.b2/c"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c\.b2"], s[dp]},
{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c\.b2"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c\.b2"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c\.b2"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c\.b2"], s[dp]},
{s[" ω fdrag"], " = ", s[n0[ω[T]]], s["c\.b3/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[й[T]]], s["c"], s[dp]},
{s[" v propr"], " = ", s["Infinity"], 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}];

tk=55;
setj={"GlobalAdaptive","MaxErrorIncreases"->100,Method->"GaussKronrodRule"};

Σj[r_]:=r^2+a^2 Cos[θ0]^2;
\[CapitalDelta]j[r_]:=r^2-2 r+a^2;
\[CapitalChi]j[r_]:=(r^2+a^2)^2-a^2 Sin[θ0]^2 \[CapitalDelta]j[r];

grr[r_]:=Σj[r]/\[CapitalDelta]j[r];

Ȓj[я_]:=Sqrt[NIntegrate[Sqrt[Abs[grr[r]]],{r,0,Abs[я]},Method->setj,MaxRecursion->100]^2+a^2];
гj[я_]:=x/.FindRoot[Ȓj[x]-я,{x,1}];

grpl=Graphics[{Table[{Magenta,Circle[{0,0},гj[n]]},{n,1/2,40,1/2}]}]
pppl=ParametricPlot[

Table[{x[tt,n], y[tt,n]}, {n,10,360,10}],
{tt, 0, T},
PlotStyle-> {Thickness[0.003], Blue},
PlotPoints-> Automatic,
MaxRecursion-> mrec,
PlotRange->{{-15,15},{-15,15}},PlotTheme->"Classic",AspectRatio->1,Axes->False]
plot1a[{xx_, yy_, zz_, tk_, w1_, w2_,PR_}]:=                                     (* Animation *)
Show[

ParametricPlot[

{0, 0},
{ttt, 0,1},
PlotStyle-> {Thickness[0.003], Blue},
PlotPoints-> Automatic,
MaxRecursion-> mrec,
Frame->True,PlotRange->{{-PR,PR},{-PR,PR}},Frame->True,PlotTheme->"Classic",AspectRatio->1,ImageSize->516,ImagePadding->18,Axes->False],

pppl,
grpl,

Graphics[{White,Disk[{0,0},Sqrt[rA^2+a^2]]}],
Graphics[{Magenta,Opacity[0.1],Annulus[{0,0},{Sqrt[rA^2+a^2],25}]}],Graphics[{White,Opacity[0.7],Disk[{0,0},Sqrt[rA^2+a^2]]}],Graphics[{Cyan,Opacity[0.1],Disk[{0,0},Sqrt[rA^2+a^2]]}],
Graphics[{Magenta,Circle[{0,0},Sqrt[rA^2+a^2]]}]
];

Do[
Print[Rasterize[Grid[{{
plot1a[{0, 0, Infinity, tk, w1r, w2r,PR}]
}}, Alignment->Left]]],
{PR,10,2,-1/20}]
Flamm's paraboloid, rubber-band-model for the Kerr-metric
Bild
Симон Тыран @ vk || wikipedia || stackexchange || wolfram

Benutzeravatar
Yukterez
Administrator
Beiträge: 194
Registriert: Mi 21. Okt 2015, 02:16

Isometrische Einbettung der Kerr Metrik

Beitragvon Yukterez » Fr 2. Mär 2018, 06:13

Bild
Animation des 3D Paraboloids auf der äquatorialen θ=π/2 Ebene für alle Spinparameter; für eine alternative Darstellung ohne einlaufende Lichtstrahlen siehe Bardeen '72, Fig. 2. a/M läuft von 0 bis 1, Längeneinheiten wie oben in Gℳ/c² mit G=ℳ=c=1 wobei ℳ die irreduzible Masse ist. Die metrischen Abstände bleiben da ℳ im Gegensatz zu M von a unabhängig ist für alle a die gleichen, und der Horizont bleibt auf konstantem Ř: r=M+√(M²-a²)→Ř=√(r²+a²)=√(x²+y²)→Ř=2ℳ:

Bild

Um die Animation auf ¼ des Tempos zu verlangsamen hier klicken. Bilder: Yukterez (Simon Tyran), CC BY-SA 4.0. Code:
Bild
Bild
2D (x,y) auf 3D (x,y,w) Einbettung, Export als Bildsequenz in den Dokumente-Ordner:

Code: Alles auswählen

(* Yukterez, version 02.03.2018 || syntax: Wolfram Mathematica || kerr.yukterez.net || flamm.yukterez.net *)
 
ClearAll["Global`*"]
Quiet[Do[
 
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;
 
a; A=a;
crd=1; dsp=1;
tmax=If[a>99/100, 500, 120]; Tmax=tmax;
TMax=Min[Tmax,т[plunge-1*^-4]];tMax=Min[tmax,plunge];
 
r0=Sqrt[2]10 mirr;
θ0=Pi/2;
φ0=0;
μ=-0;
v0=1;
α0=-Pi/2;
ψ0=0;
 
vr0=v0 Sin[α0];
vθ0=v0 Cos[α0] Sin[ψ0];
vφ0=v0 Cos[α0] Cos[ψ0];
x0BL[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];
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;
Lz=vφ0 Ы/j[v0];
pθ0=vθ0 Sqrt[Ξ]/j[v0];
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2,θ1->θ0];
k=Q+Lz^2+a^2 (ε^2+μ);
 
rE=1+Sqrt[1-a^2 Cos[θ]^2]; 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]; 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]; 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]; RI[A_,w1_,w2_]:=Xyz[xyZ[{Sqrt[rI^2+A^2] Sin[θ] Cos[φ],Sqrt[rI^2+A^2] Sin[θ] Sin[φ],rI Cos[θ]},w1],w2];
 
j[v_]:=Sqrt[1+μ v^2];
mirr=Sqrt[(Sqrt[1-a^2]+1)/2];
я=Sqrt[Χ/Σ]Sin[θ[τ]];
яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
Ы=Sqrt[χ/Ξ]Sin[θ0];
Σ=r[τ]^2+a^2 Cos[θ[τ]]^2;
Σ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]];
д[ξ_]:=Quiet[Ξ/.FindRoot[т[Ξ]-ξ,{Ξ,0}]];
T:=Quiet[д[tk]];
 
ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
γ[τ_]:=ю[τ];
R[τ_]:=Evaluate[r[τ]/.sol][[1]];
Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
ß[τ_]:=Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2]/ю[τ]];
 
ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]];ς0=Sqrt[χ/δ/Ξ];
ω[τ_]:=2R[τ] a/Χi[τ];ω0=2r0 a/χ;ωd=2r[τ] a/Χ;
Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];
й[τ_]:=ω[τ] яi[τ] ς[τ];й0=ω0 Ы ς0;
ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ];ж0=Sqrt[ς0^2-1]/ς0;
 
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]]];
dst[τ_]:=Evaluate[str[τ]/.sol][[1]];
pΘ[τ_]:=Evaluate[pθ[τ]/.sol][[1]];pΘks[τ_]:=Σi[τ] Θ'[τ];
pR[τ_]:=Evaluate[pr[τ]/.sol][[1]];pRks[τ_]:=Σi[τ]/Δi[τ] R'[τ];
sh[τ_]:=Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
epot[τ_]:=ε+μ-ekin[τ];
ekin[τ_]:=If[μ==0,ς[τ],1/Sqrt[1-v[τ]^2]-1];
dp=Y^Y;n0[z_]:=Chop[N[z]];
 
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+μ));
 
DGL={t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),t[0]==0,r'[τ]==(pr[τ] Δ)/Σ,r[0]==r0,θ'[τ]==0,θ[0]==θ0,φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),φ[0]==φ0,pr'[τ]==pr2τ[τ],pr[0]==pr0,pθ'[τ]==0,pθ[0]==pθ0,str'[τ]==vd[τ]/Max[1*^-16,Abs[Sqrt[1-vd[τ]^2]]],str[0]==0,vlt'[τ]==vd[τ],vlt[0]==0};
 
dr=(pr0 δ)/Ξ;dθ=pθ0/Ξ;dφ=(2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ;dΦ=If[θ0==0,0,If[θ0==\[Pi],0,dφ]];
φk=φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}];

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])}];
 
x[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]+n Pi/180]/.sol][[1]];
y[τ_,n_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]+n Pi/180]/.sol][[1]];
z[τ_]:=0;
 
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]];
 
xKS[τ_]:=Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
yKS[τ_]:=Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
 
XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2];XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];
Xyz[{x_,y_,z_},α_]:={x Cos[α]-y Sin[α],x Sin[α]+y Cos[α],z};
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[ψ]};
 
d1=1;
plp=50;
w1l=0;
w2l=0;
w1r=0;
w2r=0;
Mrec=300;
mrec=10;
imgsize=482;
 
s[text_]:=Style[text,FontSize->font];font=11;
 
tk=TMax;
setj={"GlobalAdaptive","MaxErrorIncreases"->100,Method->"GaussKronrodRule"};
 
Σj[r_]:=r^2+a^2 Cos[θ0]^2;
Δj[r_]:=r^2-2 r+a^2;
Χj[r_]:=(r^2+a^2)^2-a^2 Sin[θ0]^2 Δj[r];
 
grr[r_]:=Σj[r]/Δj[r];
 
Ȓj[я_]:=Sqrt[NIntegrate[Sqrt[Abs[grr[r]]],{r,0,Abs[я]},Method->setj,MaxRecursion->100]^2+a^2];
гj[я_]:=Re[x/.FindRoot[Ȓj[x]-я,{x,3}]];
 
Ȓ[я_]:=NIntegrate[Sqrt[Abs[grr[r]]],{r,0,Abs[я]},Method->setj,MaxRecursion->100];
Ř[я_]:=NIntegrate[r/Sqrt[r^2+a^2] Sqrt[Abs[(grr[r]-1)]],{r,Abs[я],r0},Method->setj,MaxRecursion->100];
ř[я_]:=Ř[0]-Ř[я];
г[я_]:=x/.FindRoot[Ȓ[x]-я,{x,1}];
z[я_]:=ř[я]-ř[r0];
 
rint=Evaluate[Interpolation[Table[{tint,R[д[tint]]},{tint,0,tk,1}]]];
pint=Evaluate[Interpolation[Table[{tint,Φ[д[tint]]},{tint,0,tk,1}]]];
zint=Evaluate[Interpolation[Table[{tint,z[rint[tint]]},{tint,0,tk,1}]]];
xint[τ_,n_]:=Sqrt[rint[τ]^2+a^2] Sin[θ0] Cos[pint[τ]+n Pi/180];
yint[τ_,n_]:=Sqrt[rint[τ]^2+a^2] Sin[θ0] Sin[pint[τ]+n Pi/180];
 
plot3D[PR_,w2_]:=
Rasterize[Grid[{{
Show[
ParametricPlot3D[Table[xYz[{xint[tt,n]/mirr, yint[tt,n]/mirr,zint[tt]/mirr},w2],{n,10,360,10}],{tt,0,tk},
PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic",
PlotRange->{{-PR,PR},{-PR,PR},{-PR,PR}},ImageSize->imgsize,ImagePadding->1,ViewPoint->{Infinity,0,0}],
ParametricPlot3D[xYz[{2 Sin[uuu], 2 Cos[uuu],zint[tk]/mirr-0.005},w2],{uuu,0,2Pi},
PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic"],
If[a == 1,ParametricPlot3D[xYz[{2 Sin[uuu], 2 Cos[uuu], zint[tk]/mirr - uuu/100}, w2],
{uuu, 0, 2000}, PlotStyle -> {Thickness[0.0025]}, PlotPoints -> Automatic, PlotTheme -> "Classic",
ColorFunction -> Function[{x, y, z}, Darker[Hue[2/3, 1, 1, z]]]], {}],
ParametricPlot3D[xYz[{r0/mirr Sin[uuu], r0 /mirr Cos[uuu],zint[0]/mirr},w2],{uuu,0,2Pi},
PlotStyle->{Thickness[0.0025],Darker[Blue]},PlotPoints->Automatic,PlotTheme->"Classic"]
]},{s["   a"->N[a]]}},Alignment->Left]];
 
Export["A-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,0]];
Export["B-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/8]];
Export["C-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/4]];
Export["D-KerrParaboloid"<>ToString[100a]<>".png",plot3D[16,Pi/2]];
 
ClearAll["Global`*"];
ClearAll["Local`*"],
{a,1,0,-1/5}]];
 
(* Bildsequenz wird ins Hauptverzeichnis exportiert *)
Bild
Dieser Beitrag ist ein Unterkapitel von kerr.yukterez.net.
Bild
Симон Тыран @ vk || wikipedia || stackexchange || wolfram


Zurück zu „Yukterez Notizblock“

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 1 Gast