Foreword: the equation of the photon trajectory with Schwarzschild metric \(\frac{du}{d\varphi}=\pm\sqrt{R_su^3-u^2+\frac{1}{b^2}}\) allows to calculate \(\frac{d\varphi}{du}\), and \(\varphi(u)\) may then be calculated analytically by an elliptic integral of \(\frac{d\varphi}{du}\) when there is a pericentre, which is not always the case, depending on the value of \(b\). Moreover, the case of a photon emitted with a radial coordinate very close to \(R_s\) cannot be handled by this method.
By writing \((\frac{du}{d\varphi})^2=R_su^3-u^2+\frac{1}{b^2}\), the differential equation can be processed
without the \(\pm\) root singularity.
The derivation of \((\frac{du}{d\varphi})^2\) actually gives:
\(2\frac{du}{d\varphi}\frac{d^2u}{d\varphi ^2}=3R_su^2\frac{du}{d\varphi}-2u\frac{du}{d\varphi}\) that is \(\frac{d^2u}{d\varphi^2}={3\over 2}R_su^2-u\).
This differential equation may then be solved numerically by double integration
with the initial conditions and overcome the limits of the elliptic integral of \(\frac{d\varphi}{du}\) discussed earlier.
Double integration – 4th-order Runge-Kutta method
The equation with second derivative seen above can be numerically integrated as follows by a 4th-order Runge-Kutta method:
1) with \(k_1=\frac{d^2u}{d\varphi^2}(u(\varphi))\), \(k_2=\frac{d^2u}{d\varphi^2}(u(\varphi)+\frac{\delta \varphi}{2}\frac{du}{d\varphi})\), \(k_3=\frac{d^2u}{d\varphi^2}(u(\varphi)+\frac{\delta \varphi}{2}\frac{du}{d\varphi}+(\frac{\delta \varphi}{2})^2\ k_1)\) and \(k_4=\frac{d^2u}{d\varphi ^2}(u(\varphi)+\delta \varphi\ \frac{du}{d\varphi}+\frac{(\delta \varphi)^2}{2}\ k_2)\)
the value of \(\frac{du}{d\varphi}\) at the next step is:
\(\frac{du}{d\varphi}(\varphi+\delta \varphi)=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{6}(k_1+2k_2+2k_3+k_4)\), and
2) with \(l_1=\frac{du}{d\varphi}(u(\varphi))\), \(l_2=\frac{du}{d\varphi}(u(\varphi)+\frac{\delta \varphi}{2}l_1)\), \(l_3=\frac{du}{d\varphi}(u(\varphi)+\frac{\delta \varphi}{2}l_2)\) and \(l_4=\frac{du}{d\varphi}(u(\varphi)+\delta \varphi\ l_3)\)
the value of \(u\) at the next step is:
\(u(\varphi+\delta \varphi)=u(\varphi)+\frac{\delta \varphi}{6}(l_1+2l_2+2l_3+l_4)\)
which is written with \(l_2=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{2}k_1\), \(l_3=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{2}k_2\) and \(l_4=\frac{du}{d\varphi}(\varphi)+\delta \varphi\ k_3\)
and after development:
\(u(\varphi+\delta \varphi)=u(\varphi)+\delta \varphi\frac{du}{d\varphi}(\varphi)+\frac{(\delta \varphi)^2}{6}(k_1+k_2+k_3)\).
The calculation step size \(\delta \varphi\) can be chosen to be either constant or adaptive, depending on the required accuracy for \(r\) (\(={1\over u}\)).
The initial conditions of \(\frac{du}{d\varphi}\) et de \(u\) are found as follows, from the emission coordinates (\(r_{em}\) , \(\varphi_{em}\)) of the photon (\(u_{em}={1\over r_{em}}\)), and two potentiel cases for \(\frac{du}{d\varphi}\):
1) incoming photon (\(\frac{du}{dt}>0)\) and \(\frac{d\varphi}{dt}>0\) or outgoing photon (\(\frac{du}{dt}<0)\) and \(\frac{d\varphi}{dt}<0\): \(\frac{du}{d\varphi}=+\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\),
2) outgoing photon (\(\frac{du}{dt}<0)\) and \(\frac{d\varphi}{dt}>0\) or incoming photon (\(\frac{du}{dt}>0)\) and \(\frac{d\varphi}{dt}<0\): \(\frac{du}{d\varphi}=-\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\).
If the photon arrives from infinity, \(u_{em}=0\), the photon is incoming (\(\frac{du}{dt}>0\)) and the two potential cases become:
1) \(\frac{du}{d\varphi}=+{1\over b}\),
2) \(\frac{du}{d\varphi}=-{1\over b}\).