Yukterez

Kerr Newman Orbital Velocity & ISCO Plotter

Jul 17th, 2019
117
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* || Orbital Velocity & ISCO Plotter | yukterez.net | 16.7.2019 | Simon Tyran, Vienna || *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. ClearAll["Global`*"]; ClearAll["Local`*"];
  6.  
  7. a = a; (* Spinparameter *)
  8. ℧ = 0; (* spezifische Ladung des schwarzen Lochs *)
  9.  
  10. Σ[я_] := я^2; (* Komponenten für die äquatoriale Ebene *)
  11. Δ[я_] := я^2-2 я+a^2+℧^2;
  12. Χ[я_] := (я^2+a^2)^2-a^2 Δ[я];
  13. κ[я_] := a;
  14.  
  15. rA := 1+Sqrt[1-a^2-℧^2]; (* Horizont *)
  16. rE := 1+Sqrt[1-℧^2]; (* Ergosphäre *)
  17. R0 := If[Element[rA, Reals], rA, 0]; (* Mindestradius *)
  18.  
  19. rp := rf/.Solve[4 a^2 (rf-℧^2)-(rf^2-3 rf+2 ℧^2)^2 == 0 && rf >= R0, rf];
  20. rP := Min[rp]; Rp := Max[rp]; (* prograder und retrograder Photonenkreis *)
  21.  
  22. isco := (* innermost stable circular orbit *)
  23. 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) &&
  24. RI >= R0, RI]];
  25.  
  26. j[v_] := Sqrt[1-v^2]; (* inverser Lorentzfaktor *)
  27. Ы[я_] := Sqrt[Χ[я]/Σ[я]]; (* axialer Gyrationsradius *)
  28. ωs[я_] := (a (2 я - ℧^2))/Χ[я]; (* Frame Dragging *)
  29.  
  30. ε[я_] := Sqrt[Δ[я] Σ[я]/Χ[я]]/j[v]+Lz[я] ωs[я]; (* Gesamtenergie *)
  31. Lz[я_] := v Ы[я]/j[v]; (* Bahndrehimpuls*)
  32.  
  33. red[я_] := Quiet[Reduce[ (* Gleichungen *)
  34. dt == (Lz[я] (-a (a^2+я^2)+Δ[я] κ[я])+ε[я] ((a^2+я^2)^2-Δ[я] κ[я]^2))/(Δ[я] Σ[я])
  35. &&
  36. 0 == ((a^2+(-2+я) я+℧^2) (16 a dt dΦ я (я-℧^2)+8 dt^2 я (-я+℧^2)+dΦ^2 я (8 я (-a^2+
  37. я^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 я^6)
  38. &&
  39. dΦ == (Lz[я] (-a^2+Δ[я])+ε[я] (a (a^2+я^2)-Δ[я] κ[я]))/(Δ[я] Σ[я])
  40. &&
  41. v > 0,
  42. {v, dΦ, dt},
  43. Reals]];
  44. (* Lösung nach Radius *)
  45. sol[x_, я_] := If[Quiet@NumericQ[red[я][[x, 2]]], red[я][[x, 2]], 0]
  46.  
  47. ri := Min[isco]; (* gegebene Radialkoordinate *)
  48. vi := sol[1, ri]; (* dazugehörige Kreisbahngeschwindigkeit *)
  49.  
  50. Plot[rP, {a, 0, 1}] (* Photonenkreisradius *)
  51. Plot[ri, {a, 0, 1}] (* ISCO *)
  52. Plot[vi, {a, 0, 1}] (* lokale Kreisbahngeschwindigkeit am ISCO *)
RAW Paste Data