# Kerr Newman Orbital Velocity & ISCO Solver

Jul 16th, 2019
122
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
2. (* || Orbital Velocity & ISCO Solver | yukterez.net | 16.07.2019 | Simon Tyran, Vienna || *)
3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
4.
5. ClearAll["Global`*"]; ClearAll["Local`*"];
6.
7. a = 0.7; (* Spinparameter *)
8. ℧ = 0.7; (* spezifische Ladung des schwarzen Lochs *)
9.
10. si = rP; (* untere Grenze, prograder Photonenkreis *)
11. sr = 10; (* obere Grenze *)
12. st = 0.02; (* Interpolationsintervall für den Plot *)
13.
14. Σ[я_] := я^2; (* Komponenten für die äquatoriale Ebene *)
15. Δ[я_] := я^2-2 я+a^2+℧^2;
16. Χ[я_] := (я^2+a^2)^2-a^2 Δ[я];
17. κ[я_] := a;
18.
19. rA = 1+Sqrt[1-a^2-℧^2]; (* Horizont *)
20. rE = 1+Sqrt[1-℧^2]; (* Ergosphäre *)
21. R0 = If[Element[rA, Reals], rA, 0]; (* Mindestradius *)
22.
23. rp = rf/.Solve[4 a^2 (rf-℧^2)-(rf^2-3 rf+2 ℧^2)^2 == 0 && rf >= R0, rf];
24. rP = Min[rp]; Rp = Max[rp]; (* prograder und retrograder Photonenkreis *)
25.
26. isco = (* innermost stable circular orbit *)
27. Quiet[RI/.NSolve[0 == RI (6 RI-RI^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (RI-℧^2)^(3/2) &&
28. RI >= R0, RI]];
29.
30. {"r horizon" -> N@rA, "r ergosphere" -> N@rE, "r isco" -> N@Min[isco],
31. "r photon pro" -> N@Min[rp], "r photon ret" -> N@Max[rp], "r disk" -> N@sr}
32.
33. j[v_] := Sqrt[1-v^2]; (* inverser Lorentzfaktor *)
34. Ы[я_] := Sqrt[Χ[я]/Σ[я]]; (* axialer Gyrationsradius *)
35. ωs[я_] := (a (2 я - ℧^2))/Χ[я]; (* Frame Dragging *)
36.
37. ε[я_] := Sqrt[Δ[я] Σ[я]/Χ[я]]/j[v]+Lz[я] ωs[я]; (* Gesamtenergie *)
38. Lz[я_] := v Ы[я]/j[v]; (* Bahndrehimpuls*)
39.
40. red[я_] := Quiet[Reduce[ (* Gleichungen *)
41. dt == (Lz[я] (-a (a^2+я^2)+Δ[я] κ[я])+ε[я] ((a^2+я^2)^2-Δ[я] κ[я]^2))/(Δ[я] Σ[я])
42. &&
43. 0 == ((a^2+(-2+я) я+℧^2) (16 a dt dΦ я (я-℧^2)+8 dt^2 я (-я+℧^2)+dΦ^2 я (8 я (-a^2+
44. я^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 я^6)
45. &&
46. dΦ == (Lz[я] (-a^2+Δ[я])+ε[я] (a (a^2+я^2)-Δ[я] κ[я]))/(Δ[я] Σ[я])
47. &&
48. v > 0,
49. {v, dΦ, dt},
50. Reals]];
51. (* Lösung nach Radius *)
52. sol[x_, я_] := If[Quiet@NumericQ[red[я][[x, 2]]], red[я][[x, 2]], 0]
53. (* Interpolationstabelle *)
54. vs = Interpolation[Table[{я, sol[1, я]}, {я, si, sr, st}]];
55. φs = Interpolation[Table[{я, sol[2, я]}, {я, si, sr, st}]];
56. ts = Interpolation[Table[{я, sol[3, я]}, {я, si, sr, st}]];
57. (* Plotfunktion *)
58. plot[func_, label_] := Plot[func, {я, rP, sr},
59. GridLines -> {{Min[rp], Max[rp], rA, si, Min[isco], Max[isco], rE, sr}, {}},
60. Frame -> True, ImagePadding -> {{40, 1}, {12, 1}}, ImageSize -> 340, PlotLabel -> label,
61. PlotRange->{{0, sr}, Automatic}]
62. (* Plots *)
63. plot[Sqrt[Χ[я]/Δ[я]/Σ[я]], "Gravitational time dilation: y=dt/dт, x=r"]
64. plot[ts[я], "Total time dilation: y=dt/dτ, x=r"]
65. plot[(a (2 я-℧^2))/((a^2+я^2)^2-a^2 (a^2-2 я+я^2+℧^2)), "Frame dragging: y=dφ/dт, x=r"]
66. plot[φs[я]/ts[я], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]
67. plot[φs[я], "Coordinate speed: y=dφ/dτ, x=r"]
68. plot[vs[я], "Local velocity: y=v=dl/dτ, x=r"]
69.
70. r0 = Min[isco]; (* Radialkoordinate *)
71.
72. "r isco" -> r0 "GM/c²" (* r isco *)
73. "dt/dτ " -> sol[3, r0] "sec/sec" (* totale ZD dt/dτ *)
74. "dt/dт " -> Sqrt[Χ[r0]/Δ[r0]/Σ[r0]] "sec/sec" (* gravitative ZD dt/dт *)
75. "dφ/dτ " -> sol[2, r0] "c³/G/M" (* Koordinatenschnelligkeit dφ/dτ *)
76. "dφ/dt " -> sol[2, r0]/sol[3, r0] "c³/G/M" (* shapiroverzögerte Winkelgeschw. dφ/dt *)
77. "v " -> sol[1, r0] "c" (* prograde Kreisbahngeschwindigkeit v lokal *)
RAW Paste Data