Foreword: the equation of the photon trajectory with Schwarzschild’s metric \(\frac{du}{d\phi}=\pm\sqrt{R_su^3-u^2+\frac{1}{b^2}}\) allows to calculate \(\frac{d\phi}{du}\), and \(\phi(u)\) may then be calculated analytically by an elliptic integral of \(\frac{d\phi}{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\phi})^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\phi})^2\) actually gives:
\(2\frac{du}{d\phi}\frac{d^2u}{d\phi^2}=3R_su^2\frac{du}{d\phi}-2u\frac{du}{d\phi}\) that is \(\frac{d^2u}{d\phi^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\phi}{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\phi^2}(u(\phi))\), \(k_2=\frac{d^2u}{d\phi^2}(u(\phi)+\frac{\delta \phi}{2}\frac{du}{d\phi})\), \(k_3=\frac{d^2u}{d\phi^2}(u(\phi)+\frac{\delta \phi}{2}\frac{du}{d\phi}+(\frac{\delta \phi}{2})^2\ k_1)\) and \(k_4=\frac{d^2u}{d\phi^2}(u(\phi)+\delta \phi\ \frac{du}{d\phi}+\frac{(\delta \phi)^2}{2}\ k_2)\)
the value of \(\frac{du}{d\phi}\) at the next step is:
\(\frac{du}{d\phi}(\phi+\delta \phi)=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{6}(k_1+2k_2+2k_3+k_4)\), and
2) with \(l_1=\frac{du}{d\phi}(u(\phi))\), \(l_2=\frac{du}{d\phi}(u(\phi)+\frac{\delta \phi}{2}l_1)\), \(l_3=\frac{du}{d\phi}(u(\phi)+\frac{\delta \phi}{2}l_2)\) and \(l_4=\frac{du}{d\phi}(u(\phi)+\delta \phi\ l_3)\)
the value of \(u\) at the next step is:
\(u(\phi+\delta \phi)=u(\phi)+\frac{\delta \phi}{6}(l_1+2l_2+2l_3+l_4)\)
which is written with \(l_2=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{2}k_1\), \(l_3=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{2}k_2\) and \(l_4=\frac{du}{d\phi}(\phi)+\delta \phi\ k_3\)
and after development:
\(u(\phi+\delta \phi)=u(\phi)+\delta \phi\frac{du}{d\phi}(\phi)+\frac{(\delta \phi)^2}{6}(k_1+k_2+k_3)\).
The calculation step size \(\delta \phi\) 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\phi}\) et de \(u\) are found as follows, from the emission coordinates (\(r_{em}\) , \(\phi_{em}\)) of the photon (\(u_{em}={1\over r_{em}}\)), and two potentiel cases for \(\frac{du}{d\phi}\):
1) incoming photon (\(\frac{du}{dt}>0)\) and \(\frac{d\phi}{dt}>0\) or outgoing photon (\(\frac{du}{dt}<0)\) and \(\frac{d\phi}{dt}<0\): \(\frac{du}{d\phi}=+\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\)
2) outgoing photon (\(\frac{du}{dt}<0)\) and \(\frac{d\phi}{dt}>0\) or incoming photon (\(\frac{du}{dt}>0)\) and \(\frac{d\phi}{dt}<0\): \(\frac{du}{d\phi}=-\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\phi}=+{1\over b}\)
2) \(\frac{du}{d\phi}=-{1\over b}\)