In here there are functions that allow to plot geodesics in the upper half plane (namely, half circles perpendicular to the x axis and vertical lines). Are there similar functions that do the same thing but only in the fundamental domain of $SL_2(\mathbb{Z})$, so that the geodesic is folded up into the domain $|x|<\frac{1}{2}$ and $x^2+y^2>1$?
What I would like to have is something like the picture in page 5 of this paper.