Gravitation einer Scheibe (yukterez.net Backup)

Das ist die deutschsprachige Version. Bild

English version: click here. Bild

Framework: Gravitationsgesetz nach Newton; Keywords: Abstandsquadrat, Rechnung, Differentialgleichung, Formel, Bewegungsgleichungen

1) Konstante Dichte

Neutonisches gravitatives Potential V einer auf die z=0-Ebene ausgerichteten Scheibe mit Masse M, Radius я und Flächendichte ρ=M/π/я²:

$Bild$

Polares und äquatoriales Potential mit r=0 als Funktion von z bzw. mit z=0 als Funktion von r:

$Bild$

Semianalytische Lösung als elliptisches Integral:

$Bild$

Schwerefeld der Scheibe:

$Bild$

Vektorkomponenten der Bewegungsgleichung:

$Bild$

$Bild$

$Bild$
Bild

Sitzt zusätzlich eine Kugel mit Masse Ḿ und Radius ʁ im Zentrum addieren sich die betreffenden Terme hinzu. Potential der Kugel:

$Bild$

Hier wird von einer durchlässigen Kugel mit konstanter Dichte ausgegangen. Funktion für die innere Restmasse:

$Bild$

Zusätzliche Komponenten des Beschleunigungsvektors in Anwesenheit der Kugel:

$Bild$

$Bild$
Bild

Die Umlaufgeschwindigkeit v ergibt sich über die Zentrifugalkraft:

$Bild$
Bild

G ist die Gravitationskonstante. Benötigte Variablen:

$Bild$

$Bild$
Bild

Definitionen der benötigten Funktionen und elliptischen Integrale:

$Bild$

$Bild$

$Bild$

$Bild$
Bild

Für einen Annulus mit Innen- und Außenradius я1 und я2 wird eine Scheibe mit я=я1 von einer Scheibe mit я=я2 subtrahiert, bzw. wird statt von 0 bis я von я1 bis я2 integriert. Für einen Vergleich der semianalytischen Lösung mit einer Brute Force n-Body Simulation klick hier.
Bild

⒜ 3D-Simulator Code für Orbits um einen kugelförmigen Planeten mit Ring oder Scheibe; die Bahn kann entweder bis zum Aufprall auf die Oberfläche des Planeten, oder auch durch denselben hindurch berechnet werden, wobei im letzteren Fall eine konstante Dichte und keine Reibung angenommen werden (so als bestünde der Planet aus dunkler Materie oder als fiele der Testpartikel durch einen Gravitationstunnel):

Code:: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* ||||| Mathematica Syntax | yukterez.net | Planet, Ring & Disk Simulator | Version 2 ||||||| *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}}; mt2 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20}; mt3 = {"EquationSimplification"-> "Residual"}; mt0 = Automatic; mta = mt1; wp = MachinePrecision; T = 2000; (* Simulationsdauer *) G = 1; (* Gravitationskonstante *) M = 1/2; (* Masse der Scheibe *) Ḿ = 1; (* Masse der zentralen Kugel *) я1 = 15; (* Scheibeninnenradius *) я2 = 20; (* Scheibenaußenradius *) я3 = 10; (* Kugelradius *) x0 = 25; (* Startposition x *) y0 = 0; (* Startposition y *) z0 = 1/100000; (* Startposition z *) v0 = Sqrt[vx^2+vy^2+vz^2]; (* Anfangsgeschwindigkeit *) vx = 0; (* Anfangsgeschwindigkeit x *) vy = 0; (* Anfangsgeschwindigkeit y *) vz = 900/1019 Sqrt[G (M+Ḿ)/x0]; (* Anfangsgeschwindigkeit z *) ρ = M/(π я2^2-π я1^2); (* Flächendichte der Scheibe *) m = If[я3==0, Ḿ, Ḿ Sqrt[x[t]^2+y[t]^2+z[t]^2]^3/я3]; (* innere Kugelrestmasse *) r[t_] := Sqrt[x[t]^2+y[t]^2]; (* zylindrischer Radius *) R[t_] := Sqrt[x[t]^2+y[t]^2+z[t]^2]; (* sphärischer Radius *) k[t_, я_] := Sqrt[4r[t] я/(R[t]^2+я^2+2r[t] я)]; (* Funktionen *) α[t_, я_] := Sqrt[4 r[t] я/(я+r[t])^2]; φ[t_, я_] := ArcSin[z[t]/(R[t]^2+я^2-2r[t] я)]; ax[t_, я_] := (-4G ρ x[t] Sqrt[я]/(k[t, я] r[t]^(3/2)))((1- k[t, я]^2/2) EllipticK[k[t, я]^2]-EllipticE[k[t, я]^2]); ay[t_, я_] := (-4G ρ y[t] Sqrt[я]/(k[t, я] r[t]^(3/2)))((1- k[t, я]^2/2) EllipticK[k[t, я]^2]-EllipticE[k[t, я]^2]); az[t_, я_] := If[z0==0 && vz==0, 0, 2G ρ z[t]/Abs[z[t]](-π+(R[t]^2+я^2+2r[t]я)^(-1/2)((R[t]-я)Sqrt[(R[t]-r[t])/(R[t]+ r[t])] EllipticPi[2r[t]/(R[t]+r[t]), k[t, я]^2]+(R[t]+я)Sqrt[(R[t]+r[t])/(R[t]- r[t])] EllipticPi[-2r[t]/(R[t]-r[t]), k[t, я]^2]))]; (* Scheibenfeld *) g[χ_] := -G Min[m, Ḿ] χ[t]/Sqrt[(x[t]^2+y[t]^2+z[t]^2)^3]; (* Kugelfeld *) V[t_]:=2G ρ(z[t](π/2+π/2 Sign[я-r[t]])-Sqrt[R[t]^2+я^2+2r[t]я] EllipticE[k[t]^2]-(я^2- r[t]^2)/Sqrt[R[t]^2+я^2+2r[t] я] EllipticK[k[t]^2]-(я-r[t])/(я+ r[t]) (z[t]^2)/Sqrt[R[t]^2+я^2+2r[t] я]EllipticPi[α[t]^2,k[t]^2]); (* Potential *) Z0=z0; If[N[Z0]==0.0&&N[vz]=!=0.0, (z0=1/1*^6;), (z0=Z0;)]; sol=NDSolve[{ (* Differentialgleichung *) x''[t]==ax[t, я2]-If[я1==0, 0, ax[t, я1]]+g[x], (* Beschleunigung x *) y''[t]==ay[t, я2]-If[я1==0, 0, ay[t, я1]]+g[y], (* Beschleunigung y *) z''[t]==az[t, я2]-If[я1==0, 0, az[t, я1]]+g[z], (* Beschleunigung z *) x'[0]==vx, y'[0]==vy, z'[0]==vz, x[0]==x0, y[0]==y0, z[0]==z0}, {x, y, z}, {t, 0, T}, WorkingPrecision-> wp, MaxSteps-> Infinity, Method-> mta, InterpolationOrder-> All, StepMonitor :> (laststep=plunge; plunge=t; stepsize=plunge-laststep;), Method->{"EventLocator", "Event" :> (If[stepsize<1*^-4, 0, 1])}]; X[t_]:=x[t]/.sol[[1]]; Y[t_]:=y[t]/.sol[[1]]; Z[t_]:=z[t]/.sol[[1]]; XYZ[t_]:=Sqrt[X[t]^2+Y[t]^2+Z[t]^2]; 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[ψ]}; φ[t_] := ArcTan[Y[t], X[t]]; (* Breitengrad *) θ[t_] := ArcCos[Z[t]/XYZ[t]]; (* Längengrad *) cd = If[Ḿ<M, 2/5, 4/5]; ck = If[Ḿ<M, 4/5, 2/5]; (* Farbfunktion *) annulus3dF[color_: GrayLevel[cd], o : OptionsPattern[]]:= (* Scheibe *) Graphics3D[{Glow[color], Opacity[0.8], EdgeForm[None], Black, ChartElementData["CylindricalSector3D"][{##}, 1]}, o, Boxed->False] &; dp= Style[\!$\*SuperscriptBox[\(Y$,$Y$]\), White]; n0[z_] := Chop[Re[N[Simplify[z]]]]; s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11]; (* Anzeigestil *) PR = 40; (* Plot Range *) d1 = 50; (* Schweiflänge *) Plp = Max[100, Round[tp/2]]; (* Kurven Details *) Mrec = 5000; (* Rekursionslimit *) mrec = 15; (* Parametric Plot Subdivisionen *) imgsize = 380; (* Bildgröße *) display[tp_] := Grid[{{ (* numerisches Display *) Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" t"], " = ", s[n0[tp]], s[dp]}, {s[" R"], " = ", s[n0[XYZ[tp]]], s[dp]}, {s[" θ"], " = ", s[n0[θ[tp]]], s[dp]}, {s[" φ"], " = ", s[n0[φ[tp]]], s[dp]}, {s[" x"], " = ", s[n0[X[tp]]], s[dp]}, {s[" y"], " = ", s[n0[Y[tp]]], s[dp]}, {s[" z"], " = ", s[n0[Z[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" vt"], " = ", s[n0[Sqrt[X'[tp]^2+Y'[tp]^2+Z'[tp]^2]]], s[dp]}, {s[" vR"], " = ", s[n0[XYZ'[tp]]], s[dp]}, {s[" vθ"], " = ", s[n0[XYZ[tp] θ'[tp]]], s[dp]}, {s[" vφ"], " = ", s[n0[XYZ[tp] φ'[tp]]], s[dp]}, {s[" vx"], " = ", s[n0[X'[tp]]], s[dp]}, {s[" vy"], " = ", s[n0[Y'[tp]]], s[dp]}, {s[" vz"], " = ", s[n0[Z'[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" at"], " = ", s[n0[Sqrt[X''[tp]^2+Y''[tp]^2+Z''[tp]^2]]], s[dp]}, {s[" aR"], " = ", s[n0[XYZ''[tp]]], s[dp]}, {s[" aθ"], " = ", s[n0[XYZ[tp] θ''[tp]]], s[dp]}, {s[" aφ"], " = ", s[n0[XYZ[tp] φ''[tp]]], s[dp]}, {s[" ax"], " = ", s[n0[X''[tp]]], s[dp]}, {s[" ay"], " = ", s[n0[Y''[tp]]], s[dp]}, {s[" az"], " = ", s[n0[Z''[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" M"], " = ", s[n0[M]], s[dp]}, {s[" Ḿ"], " = ", s[n0[Ḿ]], s[dp]}, {s[" я1"], " = ", s[n0[я1]], s[dp]}, {s[" я2"], " = ", s[n0[я2]], s[dp]}, {s[" я3"], " = ", s[n0[я3]], s[dp]}, {}, {s[" "], " ", s["yukterez.net"], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}]} }, Alignment-> Left, Spacings-> {0, 0}]; plot1b[{xx_, yy_, zz_, tp_}] := (* Animation *) Show[ Graphics3D[{Glow[GrayLevel[ck]], Black, Opacity[1], Sphere[{0, 0, 0}, я3]}, ImageSize-> imgsize, PlotRange-> { {-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR}, {-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR}, {-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR} }, SphericalRegion->False, ImagePadding-> 1], (* Kugel *) annulus3dF[][{0, 2 Pi}, {я1, я2}, {-PR/150, PR/150}], (* Scheibe *) Graphics3D[{ (* Testpartikel *) {PointSize[0.015], Red, Point[ {X[tp], Y[tp], Z[tp]}]}}, ImageSize-> imgsize, PlotRange-> PR, SphericalRegion->False, ImagePadding-> 1], If[tp==0, {}, (* Schweif *) Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ {X[tt], Y[tt], Z[tt]}, {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp}, PlotStyle-> {Thickness[PR/6000]}, 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]]], If[tp==0, {}, (* Trajektorie *) Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ {X[tt], Y[tt], Z[tt]}, {tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]}, PlotStyle-> {Thickness[PR/10000], Opacity[1], Darker@Gray}, PlotPoints-> Plp, MaxRecursion-> mrec]]], ViewPoint-> {xx, yy, zz}]; Quiet[Do[ Print[Rasterize[Grid[{{ Grid[{{ plot1b[{0, -Infinity, 0, tp}], plot1b[{0, 0, +Infinity, tp}] }}, Alignment->Left]}, {display[tp]}}, Alignment->Left]]], {tp, 0, plunge, plunge/2}]] (* Export als HTML Dokument *) (* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *) (* Export direkt als Bildsequenz *) (* Quiet[Do[Export["Y:\\export\\dateiname" <> ToString[tp] <> ".png", Rasterize[... *)

⒝ Solver für die horizontale und vertikale Fallbeschleunigung, Umlaufgeschwindigkeit und Potential:

Code:: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* ||||| Mathematica Syntax | yukterez.net | Ball, Ring, Disk g,v,Φ-Solver | Version 2 ||||||| *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) G = 1; (* Gravitationskonstante *) M = 1; (* Masse der Scheibe *) Ḿ = 0; (* Masse der zentralen Kugel *) я1 = 0; (* Scheibeninnenradius *) я2 = 1; (* Scheibenaußenradius *) я3 = 1; (* Kugelradius *) ρ = M/(π я2^2-π я1^2); (* Flächendichte der Scheibe *) m[R_, Ḿ_] := If[я3==0, Ḿ, Ḿ R^3/я3^3]; (* innere Kugelrestmasse *) gř[r_, я_] := (2 Sqrt[я] G ρ (-EllipticE[(4 я r)/(я^2+2 я r+r^2)]+EllipticK[(4 я r)/(я^2+2 я r+ r^2)] (1-(2 я r)/(я^2+2 я r+r^2))))/(Sqrt[r] Sqrt[(я r)/(я^2+2 я r+r^2)]); (* g(r) der Scheibe *) gž[z_, я_] := 2 G ρ(π-(π/2 (-я+z)+1/2 π (я+z))/Sqrt[я^2+z^2]); (* g(z) der Scheibe *) gk[R_, Ḿ_] := G Min[m[R, Ḿ], Ḿ]/R^2; (* g(R) der Kugel *) gr[r_, я1_, я2_] := gř[r, я2]-If[я1==0, 0, gř[r, я1]]+gk[r, Ḿ]; (* g(r) total *) gz[z_, я1_, я2_] := gž[z, я2]-If[я1==0, 0, gž[z, я1]]+gk[z, Ḿ]; (* g(z) total *) U[r_, z_, я_] := 2 G ρ (-Sqrt[я^2+r^2+2 я r+z^2] EllipticE[(4 я r)/(я^2+r^2+ 2 я r+z^2)]+((-я^2+r^2) EllipticK[(4 я r)/(я^2+r^2+2 я r+z^2)])/Sqrt[я^2+ r^2+2 я r+z^2]+((-я+r) z^2 EllipticPi[(4 я r)/(я+ r)^2, (4 я r)/(я^2+r^2+2 я r+z^2)])/((я+r) Sqrt[я^2+r^2+ 2 я r+z^2])+1/2 π z (1+Sign[я-r])); V[r_, z_, я1_, я2_] := U[r, z, я2]-If[я1==0, 0, U[r, z, я1]]; (* Potential der Scheibe *) W[R_, Ḿ_] := Integrate[-G Min[m[j, Ḿ], Ḿ]/j^2, {j, R, Infinity}]; (* Potential der Kugel *) Plot[{gk[R, M], gr[R, я1, я2], gz[R, я1, я2]}, {R, 0, 2 я2}, (* Plot g *) GridLines->{{я1, я2, я3}, {1}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Gray, Cyan, Magenta}, PlotRange->{All, All}] Plot[{Sqrt[gk[R, M] R], Sqrt[gr[R, я1, я2]R], Sqrt[gz[R, я1, я2]R]}, {R, 0, 2 я2}, (* Plot v *) GridLines->{{я1, я2, я3}, {1}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Gray, Cyan, Magenta}, PlotRange->{All, All}] Plot[{-W[R, M], -V[R, 0, я1, я2], -V[0, R, я1, я2]}, {R, 0, 2 я2}, (* Plot Φ *) GridLines->{{я1, я2, я3}, {1}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Gray, Cyan, Magenta}, PlotRange->{All, All}]

⒞ Vektor- und Konturplot des gravitativen Felds:

Code:: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* ||||| Mathematica Syntax | yukterez.net | Planet, Ring, Disk Fieldlines | Version 2 ||||||| *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) G = 1; (* Gravitationskonstante *) M = 1; (* Masse der Scheibe *) Ḿ = 1; (* Masse der zentralen Kugel *) я1 = 15; (* Scheibeninnenradius *) я2 = 20; (* Scheibenaußenradius *) я3 = 10; (* Kugelradius *) PR = 25; (* Plot Range *) IS = 400; (* Bildgröße *) m[x_, y_, z_] := If[я3==0, Ḿ, Ḿ Sqrt[x^2+y^2+z^2]^3/я3^3] (* innere Kugelrestmasse *) ε = 1/1000; g[{x_, y_, z_}]:=-{ If[Abs[x]<ε, 0, ((2 Sqrt[я2] G M x (-EllipticE[(4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+ 2 я2 Sqrt[x^2+y^2]+z^2)]+(1-(2 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+ z^2)) EllipticK[(4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+z^2)]))/((-я1^2 π+ я2^2 π) (x^2+y^2)^(3/4) Sqrt[(я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+z^2)]))- ((2 Sqrt[я1] G M x (-EllipticE[(4 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+ 2 я1 Sqrt[x^2+y^2]+z^2)]+(1-(2 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+ z^2)) EllipticK[(4 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)]))/((-я1^2 π+ я2^2 π) (x^2+y^2)^(3/4) Sqrt[(я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)]))+ If[Ḿ==0, 0, G Min[m[x, y, z], Ḿ] x/Sqrt[(x^2+y^2+z^2)^3]]], (* x Beschleunigung *) If[Abs[y]<ε, 0, ((2 Sqrt[я2] G M y (-EllipticE[(4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+ 2 я2 Sqrt[x^2+y^2]+z^2)]+(1-(2 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+ z^2)) EllipticK[(4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+z^2)]))/((-я1^2 π+ я2^2 π) (x^2+y^2)^(3/4) Sqrt[(я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+z^2)]))- ((2 Sqrt[я1] G M y (-EllipticE[(4 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+ y^2]+z^2)]+(1-(2 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+ z^2)) EllipticK[(4 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)]))/((-я1^2 π+ я2^2 π) (x^2+y^2)^(3/4) Sqrt[(я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)]))+ If[Ḿ==0, 0, G Min[m[x, y, z], Ḿ] y/Sqrt[(x^2+y^2+z^2)^3]]], (* y Beschleunigung *) If[Abs[z]<ε, 0, -1/((-я1^2 π+я2^2 π) Abs[z]) 2 G M z ((-π+((я2+Sqrt[x^2+y^2+z^2]) ((Sqrt[x^2+ y^2]+Sqrt[x^2+y^2+z^2])/(-Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2]))^(1/2) EllipticPi[-((2 Sqrt[x^2+ y^2])/(-Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2])), (4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+ y^2]+z^2)]+(-я2+Sqrt[x^2+y^2+z^2]) Sqrt[(-Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2])/(Sqrt[x^2+y^2]+ Sqrt[x^2+y^2+z^2])] EllipticPi[(2 Sqrt[x^2+y^2])/(Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2]), (4 я2 Sqrt[x^2+y^2])/(я2^2+x^2+y^2+2 я2 Sqrt[x^2+y^2]+z^2)])/(Sqrt[я2^2+x^2+y^2+2 я2 Sqrt[x^2+ y^2]+z^2]))-(-π+((я1+Sqrt[x^2+y^2+z^2]) Sqrt[(Sqrt[x^2+y^2]+Sqrt[x^2+y^2+ z^2])/(-Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2])] EllipticPi[-((2 Sqrt[x^2+y^2])/(-Sqrt[x^2+y^2]+ Sqrt[x^2+y^2+z^2])), (4 я1 Sqrt[x^2+y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)]+(-я1+Sqrt[x^2+ y^2+z^2]) Sqrt[(-Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2])/(Sqrt[x^2+y^2]+Sqrt[x^2+y^2+ z^2])] EllipticPi[(2 Sqrt[x^2+y^2])/(Sqrt[x^2+y^2]+Sqrt[x^2+y^2+z^2]), (4 я1 Sqrt[x^2+ y^2])/(я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2)])/(Sqrt[я1^2+x^2+y^2+2 я1 Sqrt[x^2+y^2]+z^2])))+ If[Ḿ==0, 0, G Min[m[x, y, z], Ḿ] z/Sqrt[(x^2+y^2+z^2)^3]]]}; (* z Beschleunigung *) annulus3dF[color_: GrayLevel[2/5], o : OptionsPattern[]]:= (* Scheibe *) Graphics3D[{Glow[color], Opacity[0.3], EdgeForm[None], Black, ChartElementData["CylindricalSector3D"][{##}, 1]}, o, Boxed->False] &; Show[ (* 3D Vektorplot *) VectorPlot3D[{g[{x, y, z}][[1]], g[{x, y, z}][[2]], g[{x, y, z}][[3]]}, {x, -PR, PR}, {y, -PR, PR}, {z, -PR, PR}, ImageSize->IS, VectorPoints->Fine], Graphics3D[{Glow[GrayLevel[2/5]], Black, Opacity[0.3], Sphere[{0, 0, 0}, я3]}, PlotRange-> PR, SphericalRegion->False, ImagePadding-> 1], annulus3dF[][{0, 2 π}, {я1, я2}, {-PR/150, PR/150}]] vcp1 = Show[ (* x,z-Ebene *) VectorPlot[{g[{x, 0, z}][[1]], g[{x, 0, z}][[3]]}, {x, -PR, PR}, {z, -PR, PR}, ImageSize->IS, VectorPoints->Fine, PlotRange->PR], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]]; vcp2 = Show[ (* x,y-Ebene *) VectorPlot[{g[{x, y, 0}][[1]], g[{x, y, 0}][[2]]}, {x, -PR, PR}, {y, -PR, PR}, ImageSize->IS, VectorPoints->Fine, PlotRange->PR], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Annulus[{0, 0}, {я1, я2}]}]]; Grid[{{vcp1, vcp2}}] (* 2D Vektorplot *) ctp1 = Show[Table[ (* x,z-Ebene *) ContourPlot[{Norm@g[{x, 0, z}]==n}, {x, -PR, PR}, {z, -PR, PR}, ImageSize->IS, MaxRecursion->8, PlotPoints->50, ClippingStyle->Automatic, PlotRange->PR], {n, 0.001, 0.01, 0.001}], (* Plot Range und Intervall für die Linien konstanter Gravitation *) Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]]; ctp2 = Show[Table[ (* x,z-Ebene *) ContourPlot[{Norm@g[{x, y, 0}]==n}, {x, -PR, PR}, {y, -PR, PR}, ImageSize->IS, MaxRecursion->8, PlotPoints->50, ClippingStyle->Automatic, PlotRange->PR], {n, 0.001, 0.01, 0.001}], (* Plot Range und Intervall für die Linien konstanter Gravitation *) Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}], Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Annulus[{0, 0}, {я1, я2}]}]]; Grid[{{ctp1, ctp2}}] (* Konturplot *)

⒟ Die Codes sind standardmäßig auf natürliche Einheiten (G=1) voreingestellt. Konstanten in SI-Einheiten (kg, m, sek):

Code:: G = 667384/10^16; (* Gravitationskonstante *) c = 299792458; (* Lichtgeschwindigkeit *) Au = 149597870700; (* Astronomische Einheit *) Mpc = 30856775777948584200000; (* Megaparsec *) Kpc = Mpc/1000; (* Kiloparsec *) dy = 24*3600; (* Tag *) yr = 36525*dy/100; (* Julianisches Jahr *)

⑴ Vergleich von Fallbeschleunigung g (oberer Plot), Umlaufgeschwindigkeit v (mittlerer Plot) und Potential Φ (unterer Plot); die grauen Kurven zeigen g, v und Φ im Feld einer Kugel mit Ḿ=ʁ=1, die blauen im Feld einer Scheibe mit M=я=1 bei R=x, y=z=0 und die violetten bei R=z, x=y=0. Die horizontale Achse bezeichnet den Abstand des Testpartikels vom Zentrum, R:

⑵ Plot für das kombinierte Feld eines Saturnoiden mit zentraler Kugel mit Ḿ=1, ʁ=я3=10 und Ring mit M=1, я1=15, я2=20:

Bild

⑶ Proberechnung für eine unendlich ausgedehnte Scheibe mit я=∞, ρ=1 (Lösung: g=2πGρ=konstant):

Bild

⑷ Neben der Scheibe ist g höher als über derselben, und bei einer Kugel oder Punktmasse in allen Richtungen gleich stark. M=1, R=2:

Bild

⑸ Zeile 1: Vektoren des Schwerefelds einer Scheibe mit M=1, я=20. Zeile 2: Konturen konstanter Schwere; links: y=0, rechts: z=0:

Bild

⑹ Feld einer annulären Scheibe mit M=1, я1=15, я2=20; (das Schalentheorem gilt wie man sieht nur bei Kugeln, nicht bei Scheiben):

Bild

⑺ Vektor- & Kontur-Plot mit einem Ring wie oben und zusätzlich einer Kugelmasse mit Ḿ=1 und ʁ=10 im Zentrum:

Bild

Ⅰ Geneigter Orbit um eine dünne Scheibe mit Masse M=1 und radius я=20 (für die Startbedingungen klick auf die Bilder):

Ⅱ Ein anderer geneigter Orbit um die gleiche Scheibe wie im oberen Beispiel:

Ⅲ Geschlossener polarer Orbit um eine Scheibe wie die obere:

Ⅳ Polarer Orbit um eine Kugel mit Ḿ=1, ʁ=я3=10 mit einem Ring mit M=1/2 und einem Innen/Außen-Radius von я1=15, я2=20:

Ⅴ Geschlossener polarer Orbit um einen Ringplaneten wie im oberen Beispiel:

Ⅵ Äquatorialer Kreisorbit, Kugel und Ring mit den gleichen Eigenschaften wie im vorigen Beispiel:

Ⅶ Inklinierter Orbit, Kugel : Ḿ=1/2, ʁ=я3=10; Ring: M=1, я1=15, я2=20:

Ⅷ Gravitationstunnel durch einen saturnoiden und durchlässigen Planeten mit Ḿ=1, ʁ=10 und einem Ring mit M=1/2:

Ⅸ Gravitationstunnel durch einen Ring mit M=1:

Ⅹ Gravitationstunnel mit Inklination, Ring wie im oberen Beispiel:

Ⅺ Sternförmiger Orbit um einen Ring bzw. Annulus mit M=1, я1=15, я2=20:

Ⅻ Orbit innerhalb einer durchlässigen Scheibe mit M=1, я=20 (z=0+ε, z"=0), Anfangsgeschwindigkeit v⊥=√(GM/я):

Video in Full HD: klick. Für die Anzeige des gesamten Codes klick hier, und für ein anderes Beispiel hier.

Anmerkungen:

Im Integral für V kann sofern die Dichtefunktion isotrop ist aufgrund der Scheibensymmetrie von ∫{...}d[θ=0..2π] auf 2∫{...}d[θ=0..π] umgestellt werden um die Rechenzeit am Computer zu beschleunigen. Bei konstanter Dichte empfiehlt sich die semianalytische Lösung, da diese um Größenordnungen schneller ist als das numerische Integral.
Plot ⑴ zeigt dass die Umlaufgeschwindigkeit innerhalb einer homogenen Scheibe mit steigendem Radius stärker steigt als das in einer homogenen Kugel der Fall wäre. Die blaue v-Kurve für R=r, z=0 zeigt die benötigte Kreisgeschwindigkeit in der äquatorialen Ebene, während die violette v-Kurve für R=z, r=0 die benötigte Geschwindigkeit für z"=0 zeigt (aufgrund der Asymmetrie sind bei zu geringer Höhe über dem Zentrum der Scheibe keine geschlossenen polaren Orbits möglich).
Plot ⑵ zeigt dass g an den Rändern der annulären Scheibe divergiert; die Scheibe wird sozusagen zu einem Ring zusammengedrückt. Das führt unter anderem dazu dass protoplanetare Scheiben sich zu Planeten verdichten können. Die Divergenz exakt am Rand lässt sich verhindern indem ein Mindestabstand in der Größe des mittleren Abstands der Teilchen aus denen sich die Scheibe zusammensetzt zur selben eingehalten wird, damit der Testpartikel nicht auf die unendlich dünne Fläche mit zwar endlicher Flächendichte, aber unendlicher Volumendichte knallt. Um Orbits in der rein äquatorialen Ebene die in die Scheibe hineinreichen zu studieren wird z'=z"=0, z=0+ε gesetzt wobei ε beispielsweise bei einer Galaxie sinnvollerweise der mittlere Abstand der Sterne untereinander wäre, oder beim Saturnring der Radius der Steinchen aus denen sich der Ring zusammensetzt.

2) Beliebige Dichtefunktionen

Gravitatives Potential:

$Bild$

Lösung mit elliptischem Integral:

$Bild$

Bewegungsgleichungen:

$Bild$

$Bild$

$Bild$

Mit steigendem Scheibenradius linear abnehmende Dichtefunktion für Plot ⑼:

$Bild$

Exponentiell abnehmende Dichtefunktion für Plot ⑽:

$Bild$
Bild

⒠ Orbit Simulator Code für eine Scheibe mit variabler Dichte und zentralem Bulge und Halo:

Code:: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* | Mathematica Syntax | yukterez.net | Simulator f. disks w. density function | Version 3 | *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}}; mt2 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20}; mt3 = {"EquationSimplification"-> "Residual"}; mt0 = Automatic; mta = mt0; wp = MachinePrecision; T = 500; (* Simulationsdauer *) G = 1; (* Gravitationskonstante *) ρ0 = E/450/(E-2)/π; (* Referenzdichte *) я1 = 0; (* Scheibeninnenradius *) я2 = 15; (* Scheibenaußenradius *) я3 = 20; (* Halo Radius *) я4 = 5; (* Bulge Radius *) ρ[r_] := ρ0 Exp[-r/я2]; (* Dichtefunktion *) M = Integrate[ρ[r] 2π r, {r, я1, я2}]; (* Masse der Scheibe *) Ḿ = 3/2; (* Halo Masse *) Ṃ = 1/5; (* Bulge Masse *) {"M"->N[M], "Ḿ"->N[Ḿ], "Ṃ"->N[Ṃ]} x0 = 25; (* Startposition x *) y0 = 0; (* Startposition y *) z0 = 1/1000; (* Startposition z *) v0 = Sqrt[vx^2+vy^2+vz^2]; (* Anfangsgeschwindigkeit *) vx = 0; (* Anfangsgeschwindigkeit x *) vy = 1/10; (* Anfangsgeschwindigkeit y *) vz = 1/10; (* Anfangsgeschwindigkeit z *) m = If[я3==0, Ḿ, Ḿ Sqrt[x[t]^2+y[t]^2+z[t]^2]^3/я3]; (* innere Bulge Restmasse *) ṃ = If[я4==0, Ṃ, Ṃ Sqrt[x[t]^2+y[t]^2+z[t]^2]^3/я4]; (* innere Halo Restmasse *) r[t_] := Sqrt[x[t]^2+y[t]^2]; (* zylindrischer Radius *) R[t_] := Sqrt[x[t]^2+y[t]^2+z[t]^2]; (* sphärischer Radius *) gx[x_, y_, z_] := -NIntegrate[(2 G r x (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+ y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))])/(r^2+x^2+y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)+ EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))]) ρ[r])/((x^2+ y^2) Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2]), {r, я1, я2}] (* x Beschleunigung *) gy[x_, y_, z_] := -NIntegrate[(2 G r y (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+ y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))])/(r^2+x^2+y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)+ EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))]) ρ[r])/((x^2+ y^2) Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2]), {r, я1, я2}] (* y Beschleunigung *) gz[x_, y_, z_] := -NIntegrate[(4 G r z EllipticE[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2- 2 Sqrt[r^2 (x^2+y^2)]+z^2))] ρ[r])/(Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2] (r^2+x^2+ y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)), {r, я1, я2}, Exclusions->z==0] (* z Beschleunigung *) gh[χ_] := -G Min[m, Ḿ] χ[t]/Sqrt[(x[t]^2+y[t]^2+z[t]^2)^3]; (* Halo Feld *) gb[χ_] := -G Min[ṃ, Ṃ] χ[t]/Sqrt[(x[t]^2+y[t]^2+z[t]^2)^3]; (* Bulge Feld *) Z0=z0; If[N[Z0]==0.0&&N[vz]=!=0.0, (z0=1/1*^6;), (z0=Z0;)]; sol=Quiet@NDSolve[{ (* Differentialgleichung *) x''[t]==gx[x[t], y[t], z[t]]+gh[x]+gb[x], (* Beschleunigung x *) y''[t]==gy[x[t], y[t], z[t]]+gh[y]+gb[y], (* Beschleunigung y *) z''[t]==gz[x[t], y[t], z[t]]+gh[z]+gb[z], (* Beschleunigung z *) x'[0]==vx, y'[0]==vy, z'[0]==vz, x[0]==x0, y[0]==y0, z[0]==z0}, {x, y, z}, {t, 0, T}, WorkingPrecision-> wp, MaxSteps-> Infinity, Method-> mta, InterpolationOrder-> All, StepMonitor :> (laststep=plunge; plunge=t; stepsize=plunge-laststep;), Method->{"EventLocator", "Event" :> (If[stepsize<1*^-4, 0, 1])}]; X[t_]:=x[t]/.sol[[1]]; Y[t_]:=y[t]/.sol[[1]]; Z[t_]:=z[t]/.sol[[1]]; XYZ[t_]:=Sqrt[X[t]^2+Y[t]^2+Z[t]^2]; 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[ψ]}; φ[t_] := ArcTan[Y[t], X[t]]; (* Breitengrad *) θ[t_] := ArcCos[Z[t]/XYZ[t]]; (* Längengrad *) cd = 4/5; ck = 3/5; (* Farbfunktion *) annulus3dF[color_: GrayLevel[cd], o : OptionsPattern[]]:= (* Scheibe *) Graphics3D[{Glow[color], Opacity[0.8], EdgeForm[None], Black, ChartElementData["CylindricalSector3D"][{##}, 1]}, o, Boxed->False] &; dp= Style[\!$\*SuperscriptBox[\(Y$,$Y$]\), White]; n0[z_] := Chop[Re[N[Simplify[Chop[z]]]]]; s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11]; (* Anzeigestil *) PR = 40; (* Plot Range *) d1 = 50; (* Schweiflänge *) Plp = Max[100, Round[tp/2]]; (* Kurven Details *) Mrec = 5000; (* Rekursionslimit *) mrec = 15; (* Parametric Plot Subdivisionen *) imgsize = 380; (* Bildgröße *) display[tp_] := Grid[{{ (* numerisches Display *) Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" t"], " = ", s[n0[tp]], s[dp]}, {s[" R"], " = ", s[n0[XYZ[tp]]], s[dp]}, {s[" θ"], " = ", s[n0[θ[tp]]], s[dp]}, {s[" φ"], " = ", s[n0[φ[tp]]], s[dp]}, {s[" x"], " = ", s[n0[X[tp]]], s[dp]}, {s[" y"], " = ", s[n0[Y[tp]]], s[dp]}, {s[" z"], " = ", s[n0[Z[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" vt"], " = ", s[n0[Sqrt[X'[tp]^2+Y'[tp]^2+Z'[tp]^2]]], s[dp]}, {s[" vR"], " = ", s[n0[XYZ'[tp]]], s[dp]}, {s[" vθ"], " = ", s[n0[XYZ[tp] θ'[tp]]], s[dp]}, {s[" vφ"], " = ", s[n0[XYZ[tp] φ'[tp]]], s[dp]}, {s[" vx"], " = ", s[n0[X'[tp]]], s[dp]}, {s[" vy"], " = ", s[n0[Y'[tp]]], s[dp]}, {s[" vz"], " = ", s[n0[Z'[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" at"], " = ", s[n0[Sqrt[X''[tp]^2+Y''[tp]^2+Z''[tp]^2]]], s[dp]}, {s[" aR"], " = ", s[n0[XYZ''[tp]]], s[dp]}, {s[" aθ"], " = ", s[n0[XYZ[tp] θ''[tp]]], s[dp]}, {s[" aφ"], " = ", s[n0[XYZ[tp] φ''[tp]]], s[dp]}, {s[" ax"], " = ", s[n0[X''[tp]]], s[dp]}, {s[" ay"], " = ", s[n0[Y''[tp]]], s[dp]}, {s[" az"], " = ", s[n0[Z''[tp]]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}], Grid[{ {s[" "], " ", s[" "], s[dp]}, {s[" M"], " = ", s[n0[M]], s[dp]}, {s[" Ḿ"], " = ", s[n0[Ḿ]], s[dp]}, {s[" Ṃ"], " = ", s[n0[Ṃ]], s[dp]}, {s[" я1"], " = ", s[n0[я1]], s[dp]}, {s[" я2"], " = ", s[n0[я2]], s[dp]}, {s[" я3"], " = ", s[n0[я3]], s[dp]}, {s[" я4"], " = ", s[n0[я4]], s[dp]} }, Alignment-> Left, Spacings-> {0, 0}]} }, Alignment-> Left, Spacings-> {0, 0}]; plot1b[{xx_, yy_, zz_, tp_}] := (* Animation *) Show[ Graphics3D[{Glow[GrayLevel[ck]], Black, Opacity[0.1], Sphere[{0, 0, 0}, я3]}, ImageSize-> imgsize, PlotRange-> { {-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR}, {-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR}, {-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR} }, SphericalRegion->False, ImagePadding-> 1], (* Halo *) Graphics3D[{Glow[GrayLevel[ck]], Black, Opacity[0.9], Sphere[{0, 0, 0}, я4]}], (* Bulge *) annulus3dF[][{0, 2 Pi}, {я1, я2}, {-PR/150, PR/150}], (* Scheibe *) Graphics3D[{ (* Testpartikel *) {PointSize[0.015], Red, Point[ {X[tp], Y[tp], Z[tp]}]}}, ImageSize-> imgsize, PlotRange-> PR, SphericalRegion->False, ImagePadding-> 1], If[tp==0, {}, (* Schweif *) Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ {X[tt], Y[tt], Z[tt]}, {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp}, PlotStyle-> {Thickness[PR/6000]}, 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]]], If[tp==0, {}, (* Trajektorie *) Block[{$RecursionLimit = Mrec}, ParametricPlot3D[ {X[tt], Y[tt], Z[tt]}, {tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]}, PlotStyle-> {Thickness[PR/10000], Opacity[1], Darker@Gray}, PlotPoints-> Plp, MaxRecursion-> mrec]]], ViewPoint-> {xx, yy, zz}]; Quiet[Do[ Print[Rasterize[Grid[{{ Grid[{{ plot1b[{0, -Infinity, 0, tp}], plot1b[{0, 0, +Infinity, tp}] }}, Alignment->Left]}, {display[tp]}}, Alignment->Left]]], {tp, 0, plunge, plunge/2}]] (* Export als HTML Dokument *) (* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *) (* Export direkt als Bildsequenz *) (* Quiet[Do[Export["Y:\\export\\dateiname" <> ToString[tp] <> ".png", Rasterize[... *)

⒡ Solver für die horizontale und vertikale Fallbeschleunigung, Umlaufgeschwindigkeit und Potential:

Code:: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* | Mathematica Syntax | yukterez.net | Disk w. arbitrary density, ρ g v φ Plot | Version 3 | *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) G = 1; (* Gravitationskonstante *) ρ0 = 16; (* Referenzdichte der Scheibe *) Ḿ = 0; (* Masse des Halo *) Ṃ = 0; (* Masse des Bulge *) M = Integrate[ρ[r] 2π r, {r, я1, я2}]; (* Masse der Scheibe *) я1 = 0; (* Scheibeninnenradius *) я2 = 1; (* Scheibenaußenradius *) я3 = 0; (* Halo Radius *) я4 = 0; (* Bulge Radius *) ε = 1/1000; (* Epsilon *) ρ[r_] := ρ0 Exp[-10r/я2]; (* Dichtefunktion der Scheibe *) Ρ[r_] := If[я3==0, 0, If[r>я3, 0, Ḿ/(4/3 π я3^3)]]; (* Dichtefunktion des Halo *) p[r_] := If[я4==0, 0, If[r>я4, 0, Ṃ/(4/3 π я4^3)]]; (* Dichtefunktion des Bulge *) m[R_] := If[я3==0, Ḿ, Ḿ R^3/я3^3]; (* innere Halo Restmasse *) ṃ[R_] := If[я4==0, Ṃ, Ṃ R^3/я4^3]; (* innere Bulge Restmasse *) {"M"->N[M], "Ḿ"->N[Ḿ], "Ṃ"->N[Ṃ]} V[x_, y_, z_] := NIntegrate[(4 G r EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+ y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))] ρ[r])/Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2], {r, 0, я2}] W[R_] := Integrate[G Min[m[j], Ḿ]/j^2, {j, R, Infinity}]; (* Potential des Halo *) U[R_] := Integrate[G Min[ṃ[j], Ṃ]/j^2, {j, R, Infinity}]; (* Potential des Bulge *) gh[R_] := G Min[m[R], Ḿ]/R^2; (* g(R) des Halo *) gb[R_] := G Min[ṃ[R], Ṃ]/R^2; (* g(R) des Bulge *) gk[R_] := gh[R]+gb[R]; gx[x_, y_, z_] := NIntegrate[(2 G r x (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+ y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))])/(r^2+x^2+y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)+ EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))]) ρ[r])/((x^2+ y^2) Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2]), {r, я1, я2}]+ x gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* x Beschleunigung *) gy[x_, y_, z_] := NIntegrate[(2 G r y (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+ y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))])/(r^2+x^2+y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)+ EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2))]) ρ[r])/((x^2+ y^2) Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2]), {r, я1, я2}]+ y gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* y Beschleunigung *) gz[x_, y_, z_] := NIntegrate[(4 G r z EllipticE[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2- 2 Sqrt[r^2 (x^2+y^2)]+z^2))] ρ[r])/(Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2] (r^2+x^2+ y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)), {r, я1, я2}, Exclusions->z==0]+ z gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* z Beschleunigung *) Plot[{Max[0, ρ[R]], Ρ[R], p[R]}, {R, 0, 2 я2}, GridLines->{{я1, я2, я3, я4}, {}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Green, Orange, Red}, PlotRange->{All, All}] (* Plot Dichte *) Plot[{gx[R, 0, ε], gz[0, 0, R]}, {R, 0, 2 я2}, GridLines->{{я1, я2, я3, я4}, {}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Beschleunigung *) Plot[{Sqrt[gx[R, 0, ε]R], Sqrt[gz[0, 0, R]R]}, {R, 0, 2 я2}, GridLines->{{я1, я2, я3, я4}, {}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Geschwindigkeit *) Plot[{V[R, 0, ε]+W[R]+U[R], V[0, 0, R]+W[R]+U[R]}, {R, 0, 2 я2}, GridLines->{{я1, я2, я3, я4}, {}}, Frame->True, ImageSize->640, AspectRatio->1/3, PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Potential *) Block[{R=100.0}, {gx[R, 0, 0], gz[0, 0, R], G M/R^2}] (* Proberechnung *)

⒢ Rotationsanimation für Galaxien, anzuwenden in Kombination mit dem Solver:

Code:: ω[r_]:=Sqrt[gx[r, 0, ε]r]/r; ωi=FunctionInterpolation[ω[r], {r, 0.1, 1}]; Manipulate[Show[ Table[Block[{r = k/10, j = k*4}, Table[Graphics[{PointSize[0.011], Gray, Point[ {r Sin[ωi[r] t + n], r Cos[ωi[r] t + n]} ]}], {n, 0, 2 Pi - Pi/j, Pi/j}]], {k, 1, 10, 1}], Table[Block[{r = k/10, j = k}, Table[Graphics[{PointSize[0.014], Black, Point[ {r Sin[ωi[r] t + n], r Cos[ωi[r] t + n]} ]}], {n, 0, 2 Pi - Pi/j, Pi/j}]], {k, 1, 10, 1}], PlotRange -> 1.2, Frame -> True], {t, 0, 1}]

⑻ Herleitung der Lösung:

Bild

⑼ Flächendichte, Fallbeschleunigung, Umlaufgeschwindigkeit und Potential bei einer Scheibe deren Dichte zum Rand hin linear abnimmt:

Bild

⑽ Plot für eine Scheibe mit exponentiellem Dichteabfall (Modell für die Verteilung leuchtender Materie in einer Scheibengalaxie):

Bild

Beschreibung:

Out[7] zeigt die Dichtefunktion der Scheibe (in der Mitte bei r=0 beträgt ρ=3/π, und bis zum Rand bei r=я=1 sinkt die Dichte auf ρ=0 ab).
Out[8] zeigt die gravitative Fallbeschleunigung im Feld der Scheibe (die Dichte in der Mitte ist in Plot ⑼ 3 mal höher als in Plot ⑴, weshalb auch die Fallbeschleunigung in der Mitte 3 mal höher ist als im Beispiel mit konstanter Dichte).
Out[9] zeigt die benötigte Umlaufgeschwindigkeit, wobei für die violette Kurve die selben Einschränkungen gelten wie hier.
Out[10] zeigt den Betrag des gravitativen Potentials über, neben und innerhalb der Scheibe.
Out[11] ist die Proberechnung dass die Fallbeschleunigung in weiter Entfernung auf den Wert wie bei einer Punktmasse zukonvergiert.

3) Gravitation einer Scheibe mit Halo

⑾ Parameter wie im letzten Beispiel, aber mit zusätzlichem Halo: Scheibe mit exponentiell abnehmender Dichte (M=1, ρ=ρ0·e-10r/я2) eingebettet in einen sphärischen Halo (Ḿ=4), Farbkodierung wie oben:

Grün: Dichte der Scheibe, Orange: Dichte des Halo, Blau: R=x, Violett: R=z

4) Orbits mit zufällig gewählten Parametern

ⅩⅢ Scheibe: я1=0, я2=15, M=1, ρ=ρ0-ρ0·r/я2, ρ0=1/75/π=0.004244; Halo: я3=20, Ḿ=1/2; Bulge: я4=5, Ṃ=2:

ⅩⅣ Scheibe: я1=0, я2=15, M=1, ρ=ρ0·e-r/я2, ρ0=e/450/(e-2)/π=0.0026769; Halo: я3=20, Ḿ=3/2; Bulge: я4=5, Ṃ=1/5:

Orbits von Population Ⅱ Sternen können chaotisch sein; für harmonischere Orbits siehe Animation Ⅰ bis Ⅻ.

Referenzen und Tutorials: Bild

Van Biezen: Gravity with mass distribution / The milky way galactic disc (Youtube playlists)
Alberti & Vidal: Dynamics of particles in a gravitational field of an annulus (Pdf)
Lass & Blitzer: Gravitational potential due to uniform disks & rings (Pdf)
Weiss: Certain aspects of the gravitational field of a disk (Pdf)
Krogh, Ng & Snyder: Gravitational field of a disk (Pdf)

images and animations by Simon Tyran, Vienna (Yukterez) - reuse permitted under the Creative Commons License CC BY-SA 4.0