Yukterez

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,
  70. ImagePadding->1],
  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