  The numerical integration of particle trajectories in curved spacetimes is
fundamental for obtaining realistic models of the particle dynamics around
massive compact objects such as black holes and neutron stars. Generalized
algorithms capable of handling generic metrics are required for studies of both
standard spacetimes (Schwarzschild and Kerr metrics) and non-standard
spacetimes (e.g. Schwarzschild metric plus non-classical perturbations or
multiple black hole metrics). The most commonly employed explicit numerical
schemes (e.g. Runge-Kutta) are incapable of producing highly accurate results
at critical points, e.g. in the regions close to the event horizon where
gravity causes extreme curvature of the spacetime, at an acceptable
computational cost. Here, we describe a generalized algorithm for the numerical
integration of time-like (massive particles) and null (photons) geodesics in
any given 3+1 split spacetime. We introduce a new, exactly energy-conserving
implicit integration scheme based on the preservation of the underlying
Hamiltonian, and we compare its properties with a standard fourth-order
Runge-Kutta explicit scheme and an implicit midpoint scheme. We test the
numerical performance of the three schemes against analytic solutions of
particle and photon orbits in Schwarzschild and Kerr spacetimes. We also prove
the versatility of our framework in handling more exotic metrics such as
Morris-Thorne wormholes and quantum-perturbed Schwarzschild black holes. The
generalized approach is also discussed in the perspective of future extensions
to more complex particle dynamics, e.g. the addition of the Lorentz force
acting on charged particles, which allows for test particle diagnostics in
GRMHD simulations.
