Gravitation von Scheiben und Ringen

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » So 19. Apr 2020, 07:10

Bild
Bild Das ist die deutschsprachige Version.   Bild English version: click here.Bild
Bild
Framework: Gravitationsgesetz nach Newton; Keywords: Abstandsquadrat, Rechnung, Differentialgleichung, Formel, Bewegungsgleichungen
Bild
Bild
Bild1) Konstante Dichte
Bild
BildNeutonisches 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
BildFü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
Bild  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:

Bild

  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):

Bild

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

Bild

 Geschlossener polarer Orbit um eine Scheibe wie die obere:

Bild

 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:

Bild

 Geschlossener polarer Orbit um einen Ringplaneten wie im oberen Beispiel:

Bild

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

Bild

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

Bild

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

Bild

 Gravitationstunnel durch einen Ring mit M=1:

Bild

 Gravitationstunnel mit Inklination, Ring wie im oberen Beispiel:

Bild

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

Bild

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

Bild

Display, Spalten [1|2|3|4]→[Position|Geschwindigkeit|Beschleunigung|Massenverteilung]

Bild

Video in Full HD: klick. Für die Anzeige des gesamten Codes klick hier, und für ein anderes Beispiel hier.
Bild
BildAnmerkungen:
  • BildIm 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.
  • BildPlot 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).
  • BildPlot 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.
  • BildDie zweite Zeile in Plot zeichnet die Flächen konstanter Gravitation entlang derer der Meeresspiegel verliefe sofern die Scheibe mit Wasser bedeckt wäre und die Masse des Wassers im Vergleich zur Scheibe gering wäre.
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » So 19. Apr 2020, 10:07

Bild
Bild2) Beliebige Dichtefunktionen
Bild
BildGravitatives 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
Bild  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
Bild
Beschreibung:
  • BildOut[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.
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » Di 5. Mai 2020, 21:00

Bild
Bild3) Gravitation einer Scheibe mit Halo
Bild
Bild 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:

Bild

Grün: Dichte der Scheibe, Orange: Dichte des Halo, Blau: R=x, Violett: R=z
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » Mi 13. Mai 2020, 06:06

Bild
Bild4) Orbits mit zufällig gewählten Parametern
Bild
BildⅩⅢ  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:

Bild

ⅩⅣ  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:

Bild

Orbits von Population Ⅱ Sternen können chaotisch sein; für harmonischere Orbits siehe Animation Ⅰ bis Ⅻ.
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » Sa 16. Mai 2020, 01:40

Bild
Bild5) Simulator, Plot & Solver Codes
Bild
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: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| 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: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| 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: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| 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 *)
ρ->M/(π я2^2-π я1^2)
ε = 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, ε, z}][[1]], g[{x, ε, z}][[3]]}, {x, -PR, PR}, {z, -PR, PR},
ImageSize->IS, VectorPoints->35, VectorScale->0.05, 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, ε}][[1]], g[{x, y, ε}][[2]]}, {x, -PR, PR}, {y, -PR, PR},
ImageSize->IS, VectorPoints->35, VectorScale->0.05, PlotRange->PR],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3],
If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];

Grid[{{vcp1, vcp2}}]                                                          (* 2D Vektorplot *)           

ctp1 = Show[ParallelTable[                                                        (* x,z-Ebene *)
ContourPlot[{Norm@g[{x, ε, 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[ParallelTable[                                                        (* x,y-Ebene *)
ContourPlot[{Norm@g[{x, y, ε}]==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],
If[я1==0, Disk[{0, 0}, я2], 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: Alles auswählen

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 *)

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

Code: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* | 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 bei variabler Dichte:

Code: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* | 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 *)

  Feldlinien und Konturplot, Version für variable Scheibendichte:

Code: Alles auswählen

(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* ||||| Mathematica Syntax | yukterez.net | Variable Density Fieldlines | Version 2 ||||||||| *)
(* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
G  = 1;                                                               (* Gravitationskonstante *)
ρ0 = 3/400/π;                                                    (* Referenzdichte der Scheibe *)
Ḿ = 1;                                                                       (* Masse des Halo *)
Ṃ = 1;                                                                      (* Masse des Bulge *)
M = Integrate[ρ[r] 2π r, {r, я1, я2}];                                    (* Masse der Scheibe *)

я1 = 0;                                                                 (* Scheibeninnenradius *)
я2 = 20;                                                                (* Scheibenaußenradius *)
я3 = 20;                                                                        (* Halo Radius *)
я4 = 5;                                                                        (* Bulge Radius *)

ε  = 1/1000;                                                                        (* Epsilon *)
set = {"GlobalAdaptive", "MaxErrorIncreases"->100, Method->"GaussKronrodRule"};
n   = 10;                                             (* Integrationsregel und Rekursionstiefe *)

ρ[r_] := ρ0-ρ0 r/я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 *)

PR = 25;                                                                         (* Plot Range *)
IS = 400;                                                                         (* Bildgröße *)

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];
 
g[{x_, y_, z_}]:=-{
 
If[Abs[x]<ε, 0, 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}, Method->set, MaxRecursion->n]+
x gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]],                (* x Beschleunigung *)
 
If[Abs[y]<ε, 0, 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}, Method->set, MaxRecursion->n]+
y gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]],                (* y Beschleunigung *)
 
If[Abs[z]<ε, 0, 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, Method->set, MaxRecursion->n]+
z gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]]};               (* 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] &;
 
Quiet@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],
Graphics3D[{Glow[GrayLevel[2/5]], Black, Opacity[0.3], Sphere[{0, 0, 0}, я4]}],
annulus3dF[][{0, 2 π}, {я1, я2}, {-PR/150, PR/150}]]

vcp1 = Quiet@Show[                                                                (* x,z-Ebene *)
VectorPlot[{g[{x, ε, z}][[1]], g[{x, ε, z}][[3]]}, {x, -PR, PR}, {z, -PR, PR},
ImageSize->IS, VectorPoints->35, VectorScale->0.05, PlotRange->PR],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]];

vcp2 = Quiet@Show[                                                                (* x,y-Ebene *)
VectorPlot[{g[{x, y, ε}][[1]], g[{x, y, ε}][[2]]}, {x, -PR, PR}, {y, -PR, PR},
ImageSize->IS, VectorPoints->35, VectorScale->0.05, PlotRange->PR],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3],
If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];

Grid[{{vcp1, vcp2}}]                                                          (* 2D Vektorplot *)           

ctp1 = Quiet@Show[ParallelTable[                                                  (* x,z-Ebene *)
ContourPlot[{Norm@g[{x, ε, 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], Disk[{0, 0}, я4]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]];

ctp2 = Quiet@Show[ParallelTable[                                                  (* x,z-Ebene *)
ContourPlot[{Norm@g[{x, y, ε}]==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], Disk[{0, 0}, я4]}],
Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3],
If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];

Grid[{{ctp1, ctp2}}]                                                             (* Konturplot *)

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

Code: Alles auswählen

ω[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}]

Kontrollrechnung mit 10 Millionen Punktmassen: hier entlang
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild

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

Gravitation von Scheiben und Ringen

Beitragvon Yukterez » Di 19. Mai 2020, 06:00

Bild
BildReferenzen und Tutorials
Bild
BildUnterkapitel
Bild
Bildimages and animations by Simon Tyran, Vienna (Yukterez) - reuse permitted under the Creative Commons License CC BY-SA 4.0
Bild
by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]Bild


Zurück zu „Yukterez Notizblock“

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 3 Gäste