Yukterez

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