On the distribution of angles between geodesic rays associated with hyperbolic lattice points

For every two points $z_0,z_1$ in the upper half-plane, consider all elements $\gamma$ in the principal congruence group $\Gamma(N)$, acting on the upper half-plane by fractional linear transformations, such that the hyperbolic distance between $z_1$ and $\gamma z_0$ is at most $R>0$. We study the distribution of angles between the geodesic rays $[z_1,\gamma z_0]$ as $R\to \infty$, proving that the limiting distribution exists independently of $N$ and explicitly computing it. When $z_1=z_0$ this is found to be the uniform distribution on the interval $[-\pi/2,\pi/2]$.