Gravity of Disks and Rings

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

Gravity of Disks and Rings

Beitragvon Yukterez » Do 23. Apr 2020, 19:38

Bild
Bild This is the english version.   Bild Für die deutsche Version geht es hier entlang.Bild
Bild
Framework: Newton's law; Keywords: inverse square law, calculation, differential equation, formula, equations of motion
Bild
Bild
Bild1) Constant density
Bild
BildNewtonian gravitational Potential V of a thin disc with mass M, radius я and constant surface density ρ=M/π/я², aligned on the z=0-plane:

Bild

Polar and equatorial potential with r=0 as a function of z respective with z=0 as a function of r:

Bild

Semianalytical solution in terms of elliptic integrals:

Bild

Gravitational acceleration vector in the vicinity of the disk:

Bild

Equation of motion, vector components:

Bild

Bild

Bild
Bild

If there is also a ball with mass Ḿ and radius ʁ in the center, the relevant terms add. Potential of a homogenous sphere:

Bild

Here we assume that the ball is permeable and of constant density, so the function for its inner mass is:

Bild

Additional components of the acceleration vector in the vicinity of the ball:

Bild

Bild
Bild

The orbit velocity v results from the centrifugal force:

Bild
Bild

G is the gravitational constant. Required variables:

Bild

Bild
Bild

Definitions of the required functions and elliptical integrals:

Bild

Bild

Bild

Bild
Bild
BildFor an annular ring with inner and outer radius я1 and я2, a disk with я=я1 is subtracted from a disk with я=я2, respectively the integral from 0 to я is done from я1 to я2 instead. For a comparison of the semianalytic solution with a brute force n-body simulation click here.
Bild
Bild  Comparison of acceleration g (upper plot), orbit velocity v (middle) and gravitational potential Φ (lower plot); the gray curves show g, v and Φ in the field of a homogenous sphere with Ḿ=ʁ=1, the blue in the field of a disk with M=я=1 at R=x, y=z=0 and the violet at R=z, x=y=0. The horizontal axis is the distance of the test particle to the center, R:

Bild

  Plot for the combined field of a saturnoid planet with a central sphere of Ḿ=1, ʁ=я3=10 and a ring with M=1, я1=15, я2=20:

Bild

 Sample calculation for an infinite sheet with я=∞, ρ=1 (Solution: g=2πGρ=constant):

Bild

 Beside the disk, g is higher than above or below it, while with a sphere or a point mass is is the same in all directions. M=1, R=2:

Bild

 Field vectors and contours for different values of constant g in the gravitational field of a disk with M=1, я=20:

Bild

 Vector- and contour-plot with an annular disk with M=1, я1=15, я2=20 (obviously the shell-theorem only holds for spheres):

Bild

 Vector- and contour-plot with a ring like in the example above plus a spherical mass with Ḿ=1 and ʁ=10 in the center:

Bild

 Inclined orbit around a thin disk with mass M=1 and radius я=20 (for the initial conditions click on the images):

Bild

 Another inclined orbit around the same disk as above:

Bild

 Closed polar orbit around the same disk as above:

Bild

 Polar orbit around a sphere with Ḿ=1, ʁ=я3=10 with an annular ring of mass M=1/2 and an inner/outer-radius of я1=15, я2=20:

Bild

 Closed polar orbit around a ring planet with the same properties as above:

Bild

 Circular equatorial orbit, sphere and ring with the same properties as in the previous example:

Bild

 Inclined orbit, sphere: Ḿ=1/2, ʁ=я3=10, ring: M=1, я1=15, я2=20:

Bild

Gravity tunnel through a saturnoid and permeable planet with Ḿ=1, ʁ=10 with a ring of mass M=1/2:

Bild

 Gravity tunnel through a ring with M=1:

Bild

 Gravity tunnel with inclination, ring like in the last example:

Bild

 Star-shaped orbit around an annular ring with M=1, я1=15, я2=20:

Bild

 Orbit inside a permeable disk of M=1, я=20 (z=0+ε, z"=0), Initial velocity v⊥=√(GM/я):

Bild

Display, columns [1|2|3|4]→[position|velocity|acceleration|distribution]

Bild

Video in Full HD: click here. To view the full code click here, and here for an other example.
Bild
BildRemarks:
  • BildTo reduce the CPU time, the integral for V can be changed from ∫{...}d[θ=0..2π] to 2∫{...}d[θ=0..π] due to the symmetry of the disk if the density function is isotrope. With constant density the semianalytical solution is recommended, since it is significantly faster than the purely numerical integral.
  • BildPlot shows that the circular orbit velocity inside a homogenous disk grows stronger with increasing radius than it would inside a homogenous sphere. In row 2 the blue v-curve for R=r, z=0 shows the circular orbit velocity in the equatorial plane, while the violet v-curve shows the required horizontal velocity for z"=0 (due to the asymmetry of the body and it's field, at low heights above the disk's center closed polar orbits are not always possible).
  • BildIn plot g diverges on the inner and outer edges of the annular disk; so to say, the disk gets squeezed into a ring. This is one of the reasons why protoplanetary disks can collapse into planets. The divergence on the exact edges can be avoided by keeping a minimum distance on the scale of the mean distance between the particles that form the disk, so the test particle does not plunge onto the infinitely thin disk with finite area-density, but infinite volume-density. To study orbits in the purely equatorial plane that reach into the disk, one sets z'=z"=0, z=0+ε whereby ε in the context of a galaxy would be the mean distance between the stars, or in the context of saturn's rings the radius of the particles that form the annular disk.
  • BildThe 2nd line in Plot draws the surfaces of constant gravity which define the sea level if the disk was covered with water and the water's mass was small compared to the mass of the disk.
Bild
by Simon Tyran, Vienna @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild

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

Gravity of Disks and Rings

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

Bild
Bild2) Arbitrary density functions
Bild
BildGravitational potential:

Bild

Solution with elliptic integral:

Bild

Equations of motion:

Bild

Bild

Bild

In plot the density function decreases linearly with increasing disc radius:

Bild

Exponentialy decreasing density function for plot :

Bild
Bild
Bild  Derivation of the solution:

Bild

 Surface density, gravitational acceleration, circular velocity and potential around a disk whose density ρ decreases linearly:

Bild

 Plot for a disk with exponential density drop (model for the distribution of luminous matter in a disk galaxy):

Bild
Bild
BildExplanation:
  • BildOut[7] shows the density function of the disk (in the center at r=0 the density is ρ=3/π, and at the edge ar r=я=1 it decreases to ρ=0).
  • Out[8] shows the gravitational acceleration in the field of the disk (the density in the center of plot is 3 times higher than in plot with constant density, therefore the acceleration in the center is also 3 times higher than in that example).
  • Out[9] shows the orbital velocity, whereby the same restrictions apply to the violet curve as in the example above (see here).
  • Out[10] shows the absolute value of the gravitational potential above, beside and inside the disk.
  • Out[11] is the test calculation to make sure the gravitational acceleration converges to the value of a point mass at far distances.
Bild
by Simon Tyran, Vienna @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild

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

Gravity of Disks and Rings

Beitragvon Yukterez » Mi 13. Mai 2020, 03:10

Bild
Bild3) Gravity of a disk with a halo
Bild
Bild Parameters like in the previous example, but with a halo: disk with exponentially decreasing density (M=1, ρ=ρ0·e-10r/я2) embedded in a spherical halo (Ḿ=4), color coding like above:

Bild

Green: density of the disk, orange: density of the halo, blue: R=x, violet: R=z
Bild
by Simon Tyran, Vienna @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild

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

Gravity of Disks and Rings

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

Bild
Bild4) Orbits with random parameters
Bild
BildⅩⅢ  Disk: я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

ⅩⅣ  Disk: я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 of population Ⅱ stars can be chaotic; for more harmonic orbits see animation Ⅰ to Ⅻ.
Bild
by Simon Tyran, Vienna @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild

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

Gravity of Disks and Rings

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

Bild
Bild5) Simulator, Plot & Solver Codes
Bild
Bild  3D-simulator code for orbits around a spherical planet surrounded by a ring or a disk; the trajectory can be shown until the test particle hits the surface of the planet, or the test particle can fly through the planet, in which case a constant density and no friction is assumed (as if the planet was composed of dark matter or the test particle flew through a gravitational tunnel):

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 for the purely horizontal or vertical gravitational acceleration, orbital velocity and 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}]

  Vector and contour plot of the gravitative field:

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,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],
If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];

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

  The codes are set to natural units (G=1) by default. Constants is SI-units (kg, m, sec):

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

  Simulator code for orbits around or inside a disk with variable density and a central bulge or 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 for acceleration, velocity and potential with variable density:

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

  Fieldlines and contours, version for variable disk density:

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

  Rotation animation for galaxies, use in combination with the 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}]

Comparision with the brute force calculation with 10 million point masses: click here
Bild
by Simon Tyran, Vienna @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild

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

Gravity of Disks and Rings

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

Bild
Bild6) References and tutorials
Bild
Bild
BildSubchapter
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 @ minds || gab || wikipedia || stackexchange || License: CC-BY 4. If images don't load: [ctrl]+[F5]Bild


Zurück zu „Yukterez Notepad“

Wer ist online?

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