# Plot solution of an integral with parameter

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 • Oerhaps it helps to do the integral in 2 pieces, one up to the cusp and then starting agin from there. Commented Nov 2, 2023 at 8:55 ## 1 Answer 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]  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}]  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]  • Wow! Thank you very much for all the clarification. I don't have much experience with Mathematica. Again, thank you very much! Commented Nov 2, 2023 at 13:41 • You are welcome! Commented Nov 2, 2023 at 14:53 • 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 Commented Nov 3, 2023 at 6:27
• @Soliton-104 Please, see Update 1 to my answer. Commented Nov 4, 2023 at 7:34
• @Soliton-104 If so, then see Update 2 to my answer. :) Commented Nov 16, 2023 at 8:02