Nika Oman, september 2006

Kerrova metrika in izračun geodetk svetlobe

<< Graphics`Arrow`
Remove[MultipleListPlot]
<< Graphics`MultipleListPlot`

Kerrova metrika in inverz ter nastavitev koordinat:

Clear[g, invg]
g[0, 
          0] = -(1 - ((2  M r[λ])/Σ));
g[0, 1] = 0; g[0, 2] = 0;
g[0, 
          3] = -((2  
                a M r[λ] Sin[ϑ[λ]]^2)/Σ);
g[1, 
          1] = (Σ/Δ);
g[2, 2] = Σ;
g[3, 0] = g[0, 3];
g[3, 
          3] = (r[λ]^2 + 
              (a^2) + ((2  (a^2) M r[λ] Sin[ϑ[λ]]^2)/Σ)) 
          Sin[ϑ[λ]]^2;
g[1, 0] = 0; g[1, 2] = 0; g[1, 3] = 0;
g[2, 0] = 0; g[2, 1] = 0; g[2, 3] = 0;
g[3, 1] = 0; 
    g[3, 2] = 
      0;


invg[0, 0] = 
        g[3, 3]/(-g[0, 3] g[3, 0] + g[0, 0] g[3, 3]);
invg[0, 3] = 
        g[0, 3]/(g[0, 3] g[3, 0] - g[0, 0] g[3, 3]);
invg[1, 1] = (1/g)[1, 1];
invg[2, 2] = (1/g)[2, 2];
invg[3, 0] = invg[0, 3];
invg[3, 3] = 
        g[0, 0]/(-g[0, 3] g[3, 0] + g[0, 0] g[3, 3]);
invg[0, 1] = 0; invg[0, 2] = 0;
invg[1, 0] = 0; invg[1, 2] = 0; 
    invg[1, 3] = 0;
invg[2, 0] = 0; invg[2, 1] = 0; 
    invg[2, 3] = 0;
invg[3, 1] = 0; 
    invg[3, 2] = 
      0;


Σ = 
        r[λ]^2 + (a^2) 
            Cos[ϑ[λ]]^2;
Δ = 
        r[λ]^2 - 2  M r[λ] + (a^2);

    
xx[0] = t[λ];
xx[1] = r[λ];
xx[2] = ϑ[λ];
xx[3] = φ[λ];

Funkcija za izračun Christoffelovega simbola (Subscript[Γ, l m]^k) (k,l,m=0,...,3):

Clear[Christoffel];
Christoffel[k_, l_, 
          m_] := Christoffel[k, l, m] = Christoffel[k, m, l] = 
            (1/2)*Sum[
                invg[n, 
                    k]*(D[g[l, n], xx[m]] + D[g[m, n], xx[l]] - 
                      D[g[l, m], xx[n]]), {n, 0, 3}];

Fukcija, ki vrne i-to enačbo v sistemu enačb geodetke (i=0,...,3):

Clear[EnacbaGeodetke, M, a]; 
    EnacbaGeodetke[k_] := EnacbaGeodetke[k] = 
        D[xx[k], {λ, 2}] + 
            Sum[Christoffel[k, l, m]*D[xx[l], λ]*
                D[xx[m], λ], {l, 0, 3}, {m, 0, 3}] == 
          0;

Določitev pogojev za ekvatorialno ravnino in izračun začetnega pogoja za čas:

Clear[a, M, tz, dtz];
ekvator = {ϑ[λ] -> (π/2), ϑ'[λ] -> 0};
tz = 0;
dtz := (((-g[0, 3] dφz - (Sqrt[dφ(z^2) g[0, 3]^2 - g[0, 0]*(g[1, 1] (drz^2) + g[2, 2] dϑ(z^2) + g[3, 3] dφ(z^2))]))/g)[0, 0] /. 
                r[λ] -> 
                  rz /. ϑ[λ] -> ϑz /. φ[λ] -> φz);

Izračun zunanjega horizonta:

outHorizon := M + (Sqrt[(M^2) - (a^2)]);

Dolčitev  numeričnih metod za reševanje sistema enačb geodetke in definicija funkcij za izris horizonta in geodetke (normalizacija na dt):

metFS := {"FixedStep", 
          Method -> "ExplicitEuler"};
\(\(metEvent := {EventLocator, "Event" -> \((\((Abs[outHorizon] /. 
                    aM)\) - 
                r[λ])\), "EventLocationMethod" -> "LinearInterpolation", "Method" -> 
            metFS};\)

\)
Horizon := {GrayLevel[0.85], 
          Disk[{0, 0}, Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. aM], GrayLevel[0], 
          Circle[{0, 0}, Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. aM], 
          GrayLevel[0], 
          Line[{{Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. aM, 
                0}, {-Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. aM, 0}}], 
          Line[{{0, 
                Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. 
                  aM}, {0, -Abs[(Sqrt[(outHorizon^2) + (a^2)])] /. aM}}], 
          GrayLevel[
            0]};


narisi[
          options___] := (
{{rz, ϑz, φz}, {drz, dϑz, 
                dφz}} = {{rz, ϑz, φz}, (1/dtz) {drz, dϑz, 
                  dφz}}; 
geodetka = 
            NDSolve[Join[
Table[
                      EnacbaGeodetke[ii], {ii, 0, 3}], 
{ 
                      t[0] == tz, 
                      r[0] == 
                        rz, ϑ[
                          0] == ϑz, φ[
                          0] == φz}, {(t^′)[
                          0] == dtz, (r^′)[0] == 
                        drz, (ϑ^′)[0] == 
                        dϑz, (φ^′)[
                          0] == dφz}
] /. 
                  aM /. ekvator, {t, 
                r, φ, ϑ}, {λ, 0, 
                lambdaMax}, 
Method -> 
                metEvent
, 
              StartingStepSize -> startStep, MaxSteps -> 50000]; 
          ParametricPlot[{(Sqrt[r[λ]^2 + (a^2)])*
                    Cos[φ[λ]]*
                    Sin[ϑ[λ]], (Sqrt[r[λ]^2 + (a^2)])*Sin[φ[λ]]*Sin[ϑ[λ]]} /. aM /. 
              geodetka, {λ, 0, 
              geodetka[[1, 1, 2, 1, 1, 
                2]]}, 
            options, 
PlotPoints -> 200, 
            AspectRatio -> Automatic, 
            Prolog -> 
              Horizon]
);

DF = 
        DisplayFunction -> $DisplayFunction;
DI = 
        DisplayFunction -> 
          Identity;

rotation = {Circle[{0, 0}, 0.3, {(π/8), 2  π - (π/8)}], 
          Arrow[0.4 {Cos[-(π/8)] - 0.35, Sin[-(π/8)]}, 
            0.4 {Cos[(π/8)] - 0.25, Sin[(π/8)]}, 
            HeadLength -> 0.03, HeadWidth -> 0.3, 
            HeadCenter -> 
              0.8]};

Nastavitev začetnih vrednosti (a, M, {rz,ϑz,φz}, zač. hitrosti (dr,dϑ,dφ) z ustrezno parametrizacijo):

\(\(\(
\)\(\(aM = {a -> 1, 
          M -> 1};\)

    startStep = 0.02; \ lambdaMax = 50;


    \({rz, ϑz, φz} = {7, π\/2, 
          0};\)

    \(hitrosti = {{0, 0, 1}, {\(-1\), 0, 
            Tan[Pi\/3] rz\/\(rz\^2 + a\^2\)}, {\(-1\), 0, 
            Tan[Pi\/4] rz\/\(rz\^2 + a\^2\)}, {\(-1\), 0, 
            Tan[Pi\/6] rz\/\(rz\^2 + a\^2\)}, {\(-1\), 0, 
            Tan[Pi\/8] rz\/\(rz\^2 + a\^2\)}, {\(-1\), 0, 0}, {\(-1\), 
            0, \(-Tan[Pi\/8]\) rz\/\(rz\^2 + a\^2\)}, {\(-1\), 
            0, \(-Tan[Pi\/6]\) rz\/\(rz\^2 + a\^2\)}, {\(-1\), 
            0, \(-Tan[Pi\/4]\) rz\/\(rz\^2 + a\^2\)}, {\(-1\), 
            0, \(-Tan[Pi\/3]\) rz\/\(rz\^2 + a\^2\)}, {0, 
            0, \(-1\)}};\)


    \)\)\)

Izračun in izris slike:

thc = AbsoluteThickness[1];
\(\(colors = {{GrayLevel[0], thc}};\)
\)
\(\(slika = 
        Table[{drz, dϑz, dφz} = 
            hitrosti[[ii]] /. aM; 
          narisi[PlotStyle -> {colors[[
                  Mod[ii, Length[colors], 1]]]}, DI], {ii, 
            1, Length[hitrosti]}];\)
\)
\(Show[slika, PlotRange -> {{\(-9\), 9}, {\(-9\), 9}}, 
      Frame -> True, Axes -> False, 
      FrameLabel -> {"x", "y"}, 
Epilog -> {
Arrow[{rz\ \ Cos[φz] + 1.2, 
              rz\ \ Sin[φz]}, {rz\ \ Cos[φz] + 0.2, 
              rz\ Sin[φz]}, 
            HeadLength -> 
              0.01], 
Horizon, 
Text[\*""(Subscript[T, z])"", 
{rz\ \ Cos[φz] + 
                0.7, rz\ \ Sin[φz] + 0.5}, {\(-1\), 0}, 
            TextStyle -> \((FontSize -> 14)\)], 
          Text[\*""(Subscript[R, H])"", \((Abs[\@\(outHorizon\^2 + a\^2\)] /. 
                  aM)\)*{\@3\/2, 1\/2}, {1, 1}], 
          Text[\*""(a/M)="" <> ToString[a /. aM], {\(-rz\), rz}, 
            TextStyle -> \((FontSize -> 14)\)], rotation}, 
      Prolog -> {}, 

DF]\)

Fotonske sfere

Clear[a, M, ϑ, r]

Pogoji za fotonsko sfero v ekvatorialni ravnini:

CondPhotS = {ϑ[λ] -> (π/2), 
ϑ'[λ] -> 
            0, 
r'[λ] -> 
            0, 
r''[λ] -> 
            0, 
r[λ] -> 
            r
};

Rešitev za Ω iz metrike:

solΩ = 
        FullSimplify[
          Solve[g[0, 0] + 2 g[0, 3] Ω + 
                  g[3, 3] (Ω^2) == 0, Ω] /. 
            CondPhotS, {Ω ∈ Reals, r > 0, 
            r ∈ Reals}];

Ω1 = Ω /. 
          solΩ[[1]];
Ω2 = Ω /. 
          solΩ[[2]];

Nastavitev parametrov a in M:

aM = {a -> 1, M -> 1};

Enačbe za izračun radijev fotonskih sfer (za obe rešitvi Ω):

eqn1 = 
        Simplify[(EnacbaGeodetke[1] /. 
                    CondPhotS /. t'[λ] -> (φ'[λ]/Ω)) /. φ'[λ] -> 1 /. 
            solΩ[[1]], 
          Assumptions -> {a > 0, r > 0, M > 0}];

eqn2 = 
        Simplify[(EnacbaGeodetke[1] /. 
                    CondPhotS /. t'[λ] -> (φ'[λ]/Ω)) /. φ'[λ] -> 1 /. 
            solΩ[[2]], 
          Assumptions -> {a > 0, r > 0, M > 0}];

Funkcija prebere ustrezno (največjo realno) rešitev:

MaxSol[sols_] := 
        Max[Select[r /. sols, (# ∈ Reals) &]];

Vrnemo par (a,\(R\_f\%-\)):

{a /. aM, MaxSol[Solve[eqn1 /. aM, r]]}

Vrnemo par (a,\(R\_f\%+\)):

{a /. aM, MaxSol[Solve[eqn2 /. aM, r]]}

Counterrotating = {{0, 3}, {0.1, 3.113348589668769}, {0.2, 
            3.222813948079725}, {0.3, 3.3288545395174336}, {0.4, 
            3.4318433713717136}, {0.5, 3.5320888862379562}, {0.6, 
            3.6298496971319207}, {0.7, 3.7253452507617824}, {0.8, 
            3.818763716895554}, {0.9, 3.910267939103037}, {1, 
            4}};
Corotating = {{0, 3}, {0.1, 2.8821937284436094}, {0.2, 
            2.759193090270021}, {0.3, 2.6300259776020876}, {0.4, 
            2.4933625736569014}, {0.5, 2.3472963553338606}, {0.6, 
            2.188914019763633}, {0.7, 2.0133341235972897}, {0.8, 
            1.8110859802363666}, {0.9, 1.5578546274233829}, {1, 
            1}};

MultipleListPlot[{Counterrotating, Corotating},PlotJoined -> True,AxesLabel->{(a/M),r},PlotLegend->{"\!\(\*
StyleBox[\(R\_f\%-\),
FontSize->14]\)","\!\(\*
StyleBox[\(R\_f\%+\),
FontSize->14]\)"},LegendSize -> {0.4, 0.3},LegendPosition -> {0.7, 0}]