# Disk with Density Function, Fieldlines

May 18th, 2020
24
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
2. (* ||||| Mathematica Syntax | yukterez.net | Variable Density Fieldlines | Version 2 ||||||||| *)
3. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
4.
5. G = 1; (* Gravitationskonstante *)
6. ρ0 = 3/400/π; (* Referenzdichte der Scheibe *)
7. Ḿ = 1; (* Masse des Halo *)
8. Ṃ = 1; (* Masse des Bulge *)
9. M = Integrate[ρ[r] 2π r, {r, я1, я2}]; (* Masse der Scheibe *)
10.
11. я1 = 0; (* Scheibeninnenradius *)
12. я2 = 20; (* Scheibenaußenradius *)
13. я3 = 20; (* Halo Radius *)
14. я4 = 5; (* Bulge Radius *)
15.
16. ε = 1/1000; (* Epsilon *)
17. set = {"GlobalAdaptive", "MaxErrorIncreases"->100, Method->"GaussKronrodRule"};
18. n = 10; (* Integrationsregel und Rekursionstiefe *)
19.
20. ρ[r_] := ρ0-ρ0 r/я2; (* Dichtefunktion der Scheibe *)
21. Ρ[r_] := If[я3==0, 0, If[r>я3, 0, Ḿ/(4/3 π я3^3)]]; (* Dichtefunktion des Halo *)
22. p[r_] := If[я4==0, 0, If[r>я4, 0, Ṃ/(4/3 π я4^3)]]; (* Dichtefunktion des Bulge *)
23.
24. m[R_] := If[я3==0, Ḿ, Ḿ R^3/я3^3]; (* innere Halo Restmasse *)
25. ṃ[R_] := If[я4==0, Ṃ, Ṃ R^3/я4^3]; (* innere Bulge Restmasse *)
26.
27. PR = 25; (* Plot Range *)
28. IS = 400; (* Bildgröße *)
29.
30. V[x_, y_, z_] := NIntegrate[(4 G r EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+
31. 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}]
32. W[R_] := Integrate[G Min[m[j], Ḿ]/j^2, {j, R, Infinity}]; (* Potential des Halo *)
33. U[R_] := Integrate[G Min[ṃ[j], Ṃ]/j^2, {j, R, Infinity}]; (* Potential des Bulge *)
34.
35. gh[R_] := G Min[m[R], Ḿ]/R^2; (* g(R) des Halo *)
36. gb[R_] := G Min[ṃ[R], Ṃ]/R^2; (* g(R) des Bulge *)
37. gk[R_] := gh[R]+gb[R];
38.
39. g[{x_, y_, z_}]:=-{
40.
41. If[Abs[x]<ε, 0, NIntegrate[(2 G r x (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+
42. 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)+
43. 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+
44. 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]+
45. x gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]], (* x Beschleunigung *)
46.
47. If[Abs[y]<ε, 0, NIntegrate[(2 G r y (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+
48. 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)+
49. 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+
50. 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]+
51. y gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]], (* y Beschleunigung *)
52.
53. If[Abs[z]<ε, 0, NIntegrate[(4 G r z EllipticE[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-
54. 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+
55. y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)), {r, я1, я2}, Exclusions->z==0, Method->set, MaxRecursion->n]+
56. z gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]]}; (* z Beschleunigung *)
57.
58. annulus3dF[color_: GrayLevel[2/5], o : OptionsPattern[]]:= (* Scheibe *)
59. Graphics3D[{Glow[color], Opacity[0.3], EdgeForm[None], Black,
60. ChartElementData["CylindricalSector3D"][{##}, 1]}, o,
61. Boxed->False] &;
62.
63. Quiet@Show[ (* 3D Vektorplot *)
64. VectorPlot3D[{g[{x, y, z}][[1]], g[{x, y, z}][[2]], g[{x, y, z}][[3]]},
65. {x, -PR, PR}, {y, -PR, PR}, {z, -PR, PR},
66. ImageSize->IS, VectorPoints->Fine],
67. Graphics3D[{Glow[GrayLevel[2/5]], Black, Opacity[0.3], Sphere[{0, 0, 0}, я3]},
68. PlotRange->PR,
69. SphericalRegion->False,
71. Graphics3D[{Glow[GrayLevel[2/5]], Black, Opacity[0.3], Sphere[{0, 0, 0}, я4]}],
72. annulus3dF[][{0, 2 π}, {я1, я2}, {-PR/150, PR/150}]]
73.
74. vcp1 = Quiet@Show[ (* x,z-Ebene *)
75. VectorPlot[{g[{x, ε, z}][[1]], g[{x, ε, z}][[3]]}, {x, -PR, PR}, {z, -PR, PR},
76. ImageSize->IS, VectorPoints->31, VectorScale->0.05, PlotRange->PR],
77. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
78. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
79. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]];
80.
81. vcp2 = Quiet@Show[ (* x,y-Ebene *)
82. VectorPlot[{g[{x, y, ε}][[1]], g[{x, y, ε}][[2]]}, {x, -PR, PR}, {y, -PR, PR},
83. ImageSize->IS, VectorPoints->31, VectorScale->0.05, PlotRange->PR],
84. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
85. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
86. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3],
87. If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];
88.
89. Grid[{{vcp1, vcp2}}] (* 2D Vektorplot *)
90.
91. ctp1 = Quiet@Show[Table[ (* x,z-Ebene *)
92. ContourPlot[{Norm@g[{x, ε, z}]==n}, {x, -PR, PR}, {z, -PR, PR},
93. ImageSize->IS, MaxRecursion->8, PlotPoints->50, ClippingStyle->Automatic, PlotRange->PR],
94. {n, 0.001, 0.01, 0.001}], (* Plot Range und Intervall für die Linien konstanter Gravitation *)
95. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
96. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
97. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Thickness[0.01], Line[{{-я2, 0}, {я2, 0}}]}]];
98.
99. ctp2 = Quiet@Show[Table[ (* x,z-Ebene *)
100. ContourPlot[{Norm@g[{x, y, ε}]==n}, {x, -PR, PR}, {y, -PR, PR},
101. ImageSize->IS, MaxRecursion->8, PlotPoints->50, ClippingStyle->Automatic, PlotRange->PR],
102. {n, 0.001, 0.01, 0.001}], (* Plot Range und Intervall für die Linien konstanter Gravitation *)
103. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я3]}],
104. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3], Disk[{0, 0}, я4]}],
105. Graphics[{Glow[GrayLevel[2/5]], Opacity[0.3],
106. If[я1==0, Disk[{0, 0}, я2], Annulus[{0, 0}, {я1, я2}]]}]];
107.
108. Grid[{{ctp1, ctp2}}] (* Konturplot *)
RAW Paste Data