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