Definite integral $\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$ In general relativity, null geodesics (in the unbounded case) can be written under the following form :
$$\frac{d\varphi}{dr}=\frac{1}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$
with:


*

*$\left(r, \varphi\right)$ the polar coordinates of the photon

*$R_S$ the Schwarschild radius of the central object ($R_S\in\mathbb{R^{+}_{*}}$)

*$R_0$ the distance of closest approach ($R_0 > \frac{3\sqrt{3}}{2}R_S$)


Consequently, to compute the exact trajectory up to a radius $R$ (with $R > R_0$), one can evaluate:
$$\varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$
And now comes my question : is there an analytical formula (in terms of special functions for example) corresponding to this integral ? (I would like to compute this integral numerically and an expression in terms of special functions would help a lot).

There is not an analytic pretty closed form to this one, but we can still try and learn a tiny bit more about it.  However it is related to Elliptic integrals of the first kind.   Note, Jacobi Elliptic functions are inverses for the elliptic integrals.
If you would like to compute it numerically still, let me know, I did so and have codes in mathematica from using NIntegrate.   I also computed it using just Integrate in mathematica as well and have this code.
I hope this helps, writing your integral we have
$$\varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{\frac{R_0-R_S}{R_0^3}-\frac{1}{r^2}\left(1-\frac{R_{s}}{r}\right)}}$$
We define $a\equiv (R_0-R_s)/R^3_0$, $b\equiv R_s$ and write
$$
\varphi\left(R\right)=\int_{R_0}^{R}\frac{dr}{r^2\sqrt{a-\frac{1}{r^2}\left(1-\frac{b}{r}\right)}}$$
Now we factor our some constants under square root and do algebra to obtain
$$
\varphi(r)=a^{-1/2}\int_{R_0}^R \frac{dr}{r^2\sqrt{1-\frac{1}{r^2}\left(a^{-1}-\frac{b}{ar}\right)}}=a^{-1/2}\int_{R_0}^R  \frac{dr}{r \sqrt{r^2+\frac{b}{ar}-a^{-1}}}.
$$
Now we can briefly look at the indefinite integral of this
$$
\varphi(r)=a^{-1/2}\int  \frac{dr}{r \sqrt{r^2+\frac{b}{ar}-a^{-1}}}
$$
This integral does not have a closed form however can be written in a not so pretty     Elliptic Integral of the First Kind (Elliptic-F) Function. So to answer your question, it is expressible in terms of elliptical integrals.  It is something like $I\propto r F(\arcsin(...)) $. Mathematica could not compute much else unless you want to use NIntegrate.
Elliptic integrals of first kind in Jacobi's form is given by
$$
F(x,k)=\int_0^x \frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}},
$$
where the usual form you may know is 
$$
F(\zeta,k)=\int_0^\zeta \frac{d\theta}{\sqrt{1-k^2\sin^2\theta}},
$$ 
they are related by a change of variables $t=\sin \theta$ and k is related to the eccentricity  of the orbit. Sorry I can't find a closed form for your post, but I hope elliptic F functions will do.  This is definitely also why they do not put the solution in any book, because it is not so clean.
