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)
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
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}]
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:
And the graph I am looking to obtain is
