Trajectory of a photon around a Schwarzschild black hole?

Consider a photon coming from the infinity in a unbounded orbit to a Schwarzschild black hole (Schwarzschild radius $r_{s}$) (see this for illustration). Its impact parameter is $b$ and its distance of closest approach is $r_{0}$ with $$b^2=\frac{r_{0}^{3}}{r_{0}-r_{s}}$$.

Then its trajectory in polar coordinates is defined by :

$$\frac{d\varphi}{dr} = \frac{1}{r^2\sqrt{\frac{1}{b^2}-\left(1-\frac{r_s}{r}\right)\frac{1}{r^2}}}$$

Consequently : $$\varphi\left(r\right) = \int_{r_{0}}^{r} \frac{dp}{p^2\sqrt{\frac{1}{b^2}-\left(1-\frac{r_s}{p}\right)\frac{1}{p^2}}}$$

and one can compute the total deviation using : $\Delta\varphi = 2\times\left(\lim_{r\to+\infty}\varphi\left(r\right)-\frac{\pi}{2}\right)$

But my question is : how can I plot/draw the trajectory using the integral expression of $\varphi\left(r\right)$ ?

Because if I compute : $$f\left(r\right) = 2\times\left(\int_{r_{0}}^{r} \frac{dp}{p^2\sqrt{\frac{1}{b^2}-\left(1-\frac{r_s}{p}\right)\frac{1}{p^2}}}-\frac{\pi}{2}\right)$$ I obtain $f\left(r_{0}\right) = -\pi$, and then $f$ increases up to zero, crosses zero, and tends to its positive value at infinity $\Delta\varphi$. It does not make sense for me and I do not understand how to compute the trajectory from that.

If I compute : $$g\left(r\right) = \int_{r_{0}}^{r} \frac{dp}{p^2\sqrt{\frac{1}{b^2}-\left(1-\frac{r_s}{p}\right)\frac{1}{p^2}}}$$ it starts from $0$, and increase up to $\frac{\pi}{2}+\frac{\Delta\varphi}{2}$.

I would like to compute the trajectory in the $\left(x, y\right)$ plane, so how to use the values of $f\left(r\right)$ or $g\left(r\right)$ to compute the function $y\left(x\right)$ ?

• BTW, it's convenient to work with the parameter $u=\frac{r_s}r$, rather than $r$ (and use units such that $r_s=1$). Also, that integral is an elliptic integral of the 1st kind, which can be rapidly computed to high precision using Carlson's algorithm. See Numerical computation of real or complex elliptic integrals, Bille C. Carlson (1994). Commented Jul 1, 2022 at 9:56

In fact, the problem was that computing the deflection is not very intuitive. So the trajectory in polar coordinates is : $$\varphi\left(r\right) = \int_{r_{0}}^{r} \frac{dp}{p^2\sqrt{\frac{1}{b^2}-\left(1-\frac{r_s}{p}\right)\frac{1}{p^2}}}$$ and in cartesian coordinates it is nothing else than :
• $x=r\cos\left(\varphi\left(r\right)\right)$
• $y=r\sin\left(\varphi\left(r\right)\right)$
and it represents a photon starting from $\left(x_0, y_0\right) = \left(r_{0}, 0\right)$ and going up.
But what is non-intuitive (a drawing helps a lot), it that the half deflection $\alpha/2$ is in fact :
$$\frac{\alpha}{2} = \varphi\left(r\right)-\cos^{-1}\left(\frac{r_{0}}{r}\right)$$