1
$\begingroup$

In this article, they solve the following integral $$ \mathcal{I}_{o b s}\left(\nu_0\right)=\int_{\gamma_i} \frac{(A(r))^{3 / 2}}{r^2} \sqrt{\frac{1}{A(r)}+r^2\left(\frac{d \varphi}{d r}\right)^2} d r $$

where $$ A(r)=1-\frac{2 M r^2}{\left(g^2+r^2\right)^{3 / 2}} $$ $$ \left( \frac{d \varphi}{d r} \right)^2 = \frac{1}{r^4} \left( \frac{1}{b^2} - \frac{A(r)}{r^2} \right)^{-1} = \frac{1}{r^2} \left( \frac{b^2}{r^2 - A(r) b ^2} \right) $$ so $$ I_{obs} (\nu_0) = \int_{\gamma_i} \frac{\left( A(r) \right)^{3/2}}{r^2} \sqrt{ \frac{1}{A(r)} + \frac{b^2}{r^2 - A (r) b^2} } d r $$

and when evaluating the path of the photon, they obtain a result as a function of b that when graphed is (for $g = 0$ and $M = 1$ we have the curve in black)

enter image description here

So I wrote the following in wolfram mathematica

ClearAll;

M = 1;

g = 0;

A[r_] := 1 - (2 M r^2)/(g^2 + r^2)^(3/2);

Arm[r_] := (A[r])^(3/2)/r^2 Sqrt[1/A[r] + (b^2/(r^2 - A[r] b^2))];

Simplify[ToRadicals[Integrate[Arm[r], r, Assumptions -> r > 0 && b > 0]]]

we have

enter image description here

I named the above $F \left[ r, b \right]$$ and evaluated the integral separately as follows

FF[b_] := F[100, b] - F[2.4, b]

This is because when I try to evaluate it when doing the integral, it does not give me results and they keep executing. When graphing, I get a part of the graph I want to get

Plot[FF[b], {b, 1, 10}, PlotRange -> {0, 1}]

enter image description here

Could you give me any suggestions on how to optimize the routine to be able to evaluate and graph the result on the line where it is integrated? Any suggestion or observation is welcome.

Update In the case of falling spherical accretion flow, the intensity is (eq. 24) $$ I_{o b s} \left(v_0^{i} \right) \propto \int_{\gamma_i} \frac{g_i^3 k_t d r}{r^2\left|k_r\right|} $$

$$ g_i=\left(u_e^t+\left(\frac{\mathcal{K}_r}{\mathcal{K}_t}\right) u_e^r\right)^{-1} $$

$$\mathcal{K}_t=\frac{1}{b}, \quad \frac{\mathcal{K}_r}{\mathcal{K}_t}= \pm \sqrt{B(r)\left(\frac{1}{A(r)}-\frac{b^2}{r^2}\right)}$$

$$ B(r) = 1/A(r) $$

$$ u_e^t=A(r)^{-1}, \quad u_e^r=-\sqrt{\frac{1-A(r)}{A(r) B(r)}}, \quad u_e^\theta=u_e^{\varphi}=0 $$

In Mathematica, I wrote:

Clear[Int]; M = 1;
A[r_] := 1 - (2 M r^2)/(g^2 + r^2)^(3/2);
uet = 1/A[r];
uer = -Sqrt[1 - A[r]];
kt = 1/b;
kr = kt/A[r] Sqrt[1 - (A[r] b^2)/r^2];
gi = (uet + kr/kt uer)^-1;
in = (gi^3/r^2) (kt/kr) // Simplify ;
Int[bt_?NumericQ, gt_?NumericQ] := NIntegrate[in /. {b -> bt, g -> gt}, {r, 2, 100}, Method -> "LocalAdaptive", AccuracyGoal -> 10, PrecisionGoal -> 10]
g0 = {0, .475, .75};
Plot[Evaluate[Int[b, #] & /@ g0 // Re], {b, 0, 20}, PlotRange -> All, PlotStyle -> {Black, Green, Red}, AxesLabel -> {b, Subscript[I, obs]}]

However, I do not get the expected graph (Figure 5 in the paper) instead:

enter image description here

And the graph I am looking to obtain is

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Oerhaps it helps to do the integral in 2 pieces, one up to the cusp and then starting agin from there. $\endgroup$ Nov 2, 2023 at 8:55

1 Answer 1

2
$\begingroup$

Using $\frac {d\phi}{dr}$ and Ueff from the paper we have

M = 1; L = 1;

A[r_] := 1 - (2 M r^2)/(g^2 + r^2)^(3/2); Veff = 
 1/r^2 (1 - 2 M r^2/(g^2 + r^2)^(3/2)); dphidr2 = 
 L^2/r^4/(1/b^2 - Veff);

in = (A[r])^(3/2)/r^2 Sqrt[1/A[r] + r^2 dphidr2] // Simplify

For this expression it could be better to use numerical integration in a form

 Int[b_?NumericQ,g_?NumericQ] := 
     NIntegrate[((1 - (2 r^2)/(g^2 + r^2)^(3/2))^(3/2) Sqrt[
 1/(r^2 (1/b^2 - 1/r^2 + 2/(g^2 + r^2)^(3/2))) + 1/(
  1 - (2 r^2)/(g^2 + r^2)^(3/2))])/r^2, {r, 2, 100}, Method -> "LocalAdaptive"]

Visualization

g0 = {0, .475, .75}; Plot[
 Evaluate[Int[b, #] & /@ g0 // Re], {b, 0, 20}, 
 PlotStyle -> {Black, Green, Red}, 
 AxesLabel -> {"b", "\!\(\*SubscriptBox[\(I\), \(obs\)]\)"}, PlotRange -> {0, 1.2}, PlotPoints -> 100]

Figure 1

It looks similar but not exactly the same as in the paper. Maybe we should play with parameters like integration limits {r, 2, 100} and L.

Update 1 To reproduce Figure 4 from the paper we use code

Int[b_?NumericQ, g_?NumericQ] := 
 NIntegrate[((1 - (2 r^2)/(g^2 + r^2)^(3/2))^(3/2) Sqrt[
    1/(r^2 (1/b^2 - 1/r^2 + 2/(g^2 + r^2)^(3/2))) + 1/(
     1 - (2 r^2)/(g^2 + r^2)^(3/2))])/r^2, {r, 2, 100}, 
   Method -> "LocalAdaptive", AccuracyGoal -> 5, PrecisionGoal -> 4]//Quiet
g0 = {0, .475, .75}; ParallelTable[
 DensityPlot[
  Evaluate[Re[Int[Sqrt[x^2 + y^2], g]]], {x, -15, 15}, {y, -15, 15}, 
  PlotLabel -> Row[{"g/M = ", g}], PlotLegends -> Automatic, 
  ColorFunction -> "SunsetColors", PlotPoints -> 150, 
  PlotRange -> All], {g, g0}]

Figure 2

Update 2. To compute Figure 5 from the paper we use code

Clear[Int]; M = 1;
A[r_] := 1 - (2 M r^2)/(g^2 + r^2)^(3/2); B[r_] := 1/A[r];
uet = 1/A[r];
uer = -Sqrt[1 - A[r]];
kt = 1/b;
kr = -kt  Sqrt[B[r] (1/A[r] - b^2/r^2)];
gi = kt/(kt uet + kr uer);
in = (gi^3/r^4) (kt/(kr^2)^.5);
Int[bt_?NumericQ, gt_?NumericQ] := 
 NIntegrate[in /. {b -> bt, g -> gt}, {r, 2., 100}, 
  Method -> "LocalAdaptive", AccuracyGoal -> 8, PrecisionGoal -> 8]
g0 = {0, .475, .75};
Plot[Evaluate[Int[b, #] & /@ g0 // Abs], {b, 0, 20}, PlotRange -> All,
  PlotStyle -> {Black, Green, Red}, 
 AxesLabel -> {b, Subscript[I, obs]}, PlotPoints -> 150]

Figure3

$\endgroup$
11
  • $\begingroup$ Wow! Thank you very much for all the clarification. I don't have much experience with Mathematica. Again, thank you very much! $\endgroup$ Nov 2, 2023 at 13:41
  • $\begingroup$ You are welcome! $\endgroup$ Nov 2, 2023 at 14:53
  • $\begingroup$ hi, I'm trying to obtain the two-dimensional intensity map (Fig. 4 in the article) from the results you gave me, however, I still can't get it. Could you give me any suggestions if possible? Likewise, I found in this article (page 8) that the intensity is circularly symmetric, with the impact parameter $b$ of the radius, which satisfies $$b^2 = x^2 + y^2$$ but I can't think of how to implement it $\endgroup$ Nov 3, 2023 at 6:27
  • 1
    $\begingroup$ @Soliton-104 Please, see Update 1 to my answer. $\endgroup$ Nov 4, 2023 at 7:34
  • 1
    $\begingroup$ @Soliton-104 If so, then see Update 2 to my answer. :) $\endgroup$ Nov 16, 2023 at 8:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.