  The radiation observed in quasars and active galactic nuclei is mainly
produced by a relativistic plasma orbiting close to the black hole event
horizon, where strong gravitational effects are relevant. The observational
data of such systems can be compared with theoretical models to infer the black
hole and plasma properties. In the comparison process, ray tracing algorithms
are essential to computing the trajectories followed by the photons from the
source to our telescopes. In this paper, we present OSIRIS: a new stable
FORTRAN code capable of efficiently computing null geodesics around compact
objects, including general relativistic effects such as gravitational lensing,
redshift, and relativistic boosting. The algorithm is based on the Hamiltonian
formulation and uses different integration schemes to evolve null geodesics
while tracking the error in the Hamiltonian constrain to ensure physical
results. We found from an error analysis that the integration schemes are all
stable, and the best one maintains an error below $10^{-11}$. Particularly, to
test the robustness and ability of the code to evolve geodesics in curved
spacetime, we compute the shadow and Einstein rings of a Kerr black hole with
different rotation parameters and obtain the image of a thin Keplerian
accretion disk around a Schwarzschild black hole. Although OSIRIS is
parallelized neither with MPI nor with CUDA, the computation times are of the
same order as those reported by other codes with these types of parallel
computing platforms.
