Yukterez

Disk with Density Function, Solver & Plotter

May 14th, 2020
20
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* | Mathematica Syntax | yukterez.net | Disk w. arbitrary density, ρ g v φ Plot | Version 3 | *)
  3. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. G = 1; (* Gravitationskonstante *)
  6. ρ0 = 16; (* Referenzdichte der Scheibe *)
  7. Ḿ = 0; (* Masse des Halo *)
  8. Ṃ = 0; (* Masse des Bulge *)
  9. M = Integrate[ρ[r] 2π r, {r, я1, я2}]; (* Masse der Scheibe *)
  10.  
  11. я1 = 0; (* Scheibeninnenradius *)
  12. я2 = 1; (* Scheibenaußenradius *)
  13. я3 = 0; (* Halo Radius *)
  14. я4 = 0; (* Bulge Radius *)
  15.  
  16. ε = 1/1000; (* Epsilon *)
  17.  
  18. ρ[r_] := ρ0 Exp[-10r/я2]; (* Dichtefunktion der Scheibe *)
  19. Ρ[r_] := If[я3==0, 0, If[r>я3, 0, Ḿ/(4/3 π я3^3)]]; (* Dichtefunktion des Halo *)
  20. p[r_] := If[я4==0, 0, If[r>я4, 0, Ṃ/(4/3 π я4^3)]]; (* Dichtefunktion des Bulge *)
  21.  
  22. m[R_] := If[я3==0, Ḿ, Ḿ R^3/я3^3]; (* innere Halo Restmasse *)
  23. ṃ[R_] := If[я4==0, Ṃ, Ṃ R^3/я4^3]; (* innere Bulge Restmasse *)
  24. {"M"->N[M], "Ḿ"->N[Ḿ], "Ṃ"->N[Ṃ]}
  25.  
  26. V[x_, y_, z_] := NIntegrate[(4 G r EllipticK[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+
  27. 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}]
  28. W[R_] := Integrate[G Min[m[j], Ḿ]/j^2, {j, R, Infinity}]; (* Potential des Halo *)
  29. U[R_] := Integrate[G Min[ṃ[j], Ṃ]/j^2, {j, R, Infinity}]; (* Potential des Bulge *)
  30.  
  31. gh[R_] := G Min[m[R], Ḿ]/R^2; (* g(R) des Halo *)
  32. gb[R_] := G Min[ṃ[R], Ṃ]/R^2; (* g(R) des Bulge *)
  33. gk[R_] := gh[R]+gb[R];
  34.  
  35. gx[x_, y_, z_] := NIntegrate[(2 G r x (((-r^2+x^2+y^2-z^2) EllipticE[-((4 Sqrt[r^2 (x^2+
  36. 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)+
  37. 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+
  38. y^2) Sqrt[r^2+x^2+y^2-2 Sqrt[r^2 (x^2+y^2)]+z^2]), {r, я1, я2}]+
  39. x gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* x Beschleunigung *)
  40.  
  41. gy[x_, y_, z_] := NIntegrate[(2 G r y (((-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}]+
  45. y gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* y Beschleunigung *)
  46.  
  47. gz[x_, y_, z_] := NIntegrate[(4 G r z EllipticE[-((4 Sqrt[r^2 (x^2+y^2)])/(r^2+x^2+y^2-
  48. 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+
  49. y^2+2 Sqrt[r^2 (x^2+y^2)]+z^2)), {r, я1, я2}, Exclusions->z==0]+
  50. z gk[Sqrt[x^2+y^2+z^2]]/Sqrt[x^2+y^2+z^2]; (* z Beschleunigung *)
  51.  
  52. Plot[{Max[0, ρ[R]], Ρ[R], p[R]}, {R, 0, 2 я2},
  53. GridLines->{{я1, я2, я3, я4}, {}},
  54. Frame->True, ImageSize->640, AspectRatio->1/3,
  55. PlotStyle->{Green, Orange, Red}, PlotRange->{All, All}] (* Plot Dichte *)
  56.  
  57. Plot[{gx[R, 0, ε], gz[0, 0, R]}, {R, 0, 2 я2},
  58. GridLines->{{я1, я2, я3, я4}, {}},
  59. Frame->True, ImageSize->640, AspectRatio->1/3,
  60. PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Beschleunigung *)
  61.  
  62. Plot[{Sqrt[gx[R, 0, ε]R], Sqrt[gz[0, 0, R]R]}, {R, 0, 2 я2},
  63. GridLines->{{я1, я2, я3, я4}, {}},
  64. Frame->True, ImageSize->640, AspectRatio->1/3,
  65. PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Geschwindigkeit *)
  66.  
  67. Plot[{V[R, 0, ε]+W[R]+U[R], V[0, 0, R]+W[R]+U[R]}, {R, 0, 2 я2},
  68. GridLines->{{я1, я2, я3, я4}, {}},
  69. Frame->True, ImageSize->640, AspectRatio->1/3,
  70. PlotStyle->{Cyan, Magenta}, PlotRange->{All, All}] (* Plot Potential *)
  71.  
  72. Block[{R=100.0}, {gx[R, 0, 0], gz[0, 0, R], G M/R^2}] (* Proberechnung *)
RAW Paste Data