Solution (source code)

= Solution
{c}

By spherical symmetry, choose the orbital plane to be $\theta=\pi/2$. The <geodesic Lagrangian> of the <Schwarzschild metric> is
$$
\mathcal L=\frac12\left(-Bc^2\dot t^2+B^{-1}\dot r^2+r^2\dot\phi^2\right),
\qquad B=1-\frac{2GM}{rc^2}.
$$
For a massive <test particle>, dots denote <proper time> derivatives and $2\mathcal L=-c^2$. The <cyclic coordinates> $t$ and $\phi$ give, through the <Euler-Lagrange equations>, two constants:
$$
B\dot t=\mathcal E,\qquad r^2\dot\phi=\ell.
$$
At infinity, $\mathcal E=\gamma_\infty=1/\sqrt{1-v_\infty^2/c^2}$ is the <Lorentz factor>. The constant $\ell$ is the <specific angular momentum>, and the energy per unit <rest mass> is $E_\infty/m=c^2\gamma_\infty$. The normalization of the <four-velocity> gives
$$
-c^2=-\frac{c^2\gamma_\infty^2}{B}
+\frac{\dot r^2}{B}+\frac{\ell^2}{r^2},
\qquad
\dot r^2=c^2\gamma_\infty^2-B\left(c^2+\frac{\ell^2}{r^2}\right).
$$
Use $dr/d\phi=\dot r\,r^2/\ell$ and define the <energy-normalized angular momentum for a Schwarzschild geodesic>
$$
J=\frac{\ell}{\gamma_\infty}=\frac{c^2L_\infty}{E_\infty}.
$$
Then
$$
\boxed{\left(\frac{dr}{d\phi}\right)^2
+r^2B\left(1+\frac{c^2r^2}{\gamma_\infty^2J^2}\right)
=\frac{c^2r^4}{J^2}.}
$$
For an incoming massive <test particle> with <impact parameter> $b$, $L_\infty=m\gamma_\infty v_\infty b$, so $J=bv_\infty$. In particular, $J$ is not the same as the <specific angular momentum> $\ell$ unless $\gamma_\infty=1$.

A <photon> has no <proper time> parametrization and no finite <Lorentz factor>. For it, use an <affine parameter> in the same <geodesic Lagrangian> and impose $2\mathcal L=0$. The ratio $J=c^2L_\infty/E_\infty=bc$ remains well defined, and the null result is
$$
\boxed{\left(\frac{dr}{d\phi}\right)^2+r^2B=\frac{r^4}{b^2}.}
$$
Equivalently this is the ultrarelativistic limit of the massive formula at fixed $J$, with $1/\gamma_\infty^2\to0$. The printed massive formula needs this interpretation before it is applied to a <photon>.