## Gravity of Disks and Rings

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

### Gravity of Disks and Rings

This is the english version.   Für die deutsche Version geht es hier entlang.
Framework: Newton's law; Keywords: inverse square law, calculation, differential equation, formula, equations of motion

1) Constant density

Newtonian gravitational Potential V of a thin disc with mass M, radius я and constant surface density ρ=M/π/я², aligned on the z=0-plane:

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

Semianalytical solution in terms of elliptic integrals:

Gravitational acceleration vector in the vicinity of the disk:

Equation of motion, vector components:

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

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

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

The orbit velocity v results from the centrifugal force:

G is the gravitational constant. Required variables:

Definitions of the required functions and elliptical integrals:

For 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.

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:

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:

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

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:

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

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

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

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

Another inclined orbit around the same disk as above:

Closed polar orbit around the same disk as above:

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:

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

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

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

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

Gravity tunnel through a ring with M=1:

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

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

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

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

Video in Full HD: click here. To view the full code click here, and here for an other example.

Remarks:
• To 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.
• Plot 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).
• In 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.
• The 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.

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

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

### Gravity of Disks and Rings

2) Arbitrary density functions

Gravitational potential:

Solution with elliptic integral:

Equations of motion:

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

Exponentialy decreasing density function for plot :

Derivation of the solution:

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

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

Explanation:
• Out[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.

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

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

### Gravity of Disks and Rings

3) Gravity of a disk with a halo

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:

Green: density of the disk, orange: density of the halo, blue: R=x, violet: R=z

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

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

### Gravity of Disks and Rings

4) Orbits with random parameters

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

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

Orbits of population Ⅱ stars can be chaotic; for more harmonic orbits see animation Ⅰ to Ⅻ.

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

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

### Gravity of Disks and Rings

5) Simulator, Plot & Solver Codes

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

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

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

### Gravity of Disks and Rings

6) References and tutorials

Subchapter

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

by Simon Tyran, Vienna @ youtube || rumble || odysee || minds || wikipedia || stackexchange || License: CC-BY 4 ▣ If images don't load: [ctrl]+[F5]

Zurück zu „Yukterez Notepad“

### Wer ist online?

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