Searching for large pairs of adjacent gaps in Euclid's orchard

You are at the point $(0, 0)$ on the coordinate plane. There is a tree at each point with nonnegative integer coordinates, such as $(5, 6)$, $(0, 7)$, and $(9, 1)$. The trees are very thin, so that they only obstruct trees that are directly behind them. For example, the tree at $(2, 2)$ is not visible, as it is obscured by the tree at $(1, 1).$

Now, you can’t see infinitely far in this forest. You’re not sure exactly how far you can see, but it’s pretty dang far. As you look around, you can make out very narrow angular gaps between the trees. The largest gaps are near the positive $x$-axis and the positive $y$-axis. After those, the next largest pair of gaps is on either side of the tree at $(1, 1)$, $45$ degrees up from the $x$-axis.

Consider the view between $0$ and $45$ degrees up from the $x$-axis. The next largest pair of adjacent gaps in this range are on either side of what angle up from the $x$-axis?

First, let's note that our forest is akin to what's known as Euclid's orchard, which assumes that each tree is the three dimensional line segment from $(m, n, 0)$ to $(m, n, 1)$ for $(m,n) \in \mathbb{N}^2.$ Well, it turns out that if you as an observer at $(0,0)$ are viewing this orchard along the line $y = x,$ that either the tree whose base is at $(m,n)$ is blocked by another tree if $\textsf{gcd}(m,n) \gt 1$, or the height of the tree will appear to be at height $\frac{1}{m+n}$ if $\textsf{gcd}(m,n) = 1.$ Therefore, we can associate any visible tree at $(m,n)$ with the point $x = x(m,n) = \frac{m}{m+n} \in (0,1)$ and then its height will be given by Thomae's function defined by $$T(x) = \begin{cases} \frac{1}{q}, &\text{if $x = \frac{p}{q}$ for $p \in \mathbb{Z},$ $q \in \mathbb{N}$ and $\textsf{gcd}(p,q) = 1$;}\\ 0, &\text{otherwise.}\end{cases}$$ Since we cannot actually see infinitely far into Euclid's orchard, let's assume that we can only see a tree if its base is at $(m,n) \in \mathbb{N}^2$ and $m^2 + n^2 \leq R^2$ for some $R \gg 1.$ Therefore, since $R\sqrt{2} = \max \{ x+y \mid x^2 + y^2 \leq R^2 \},$ the corresponding approximated Thomae's function is $$T_R(x) = \begin{cases} \frac{1}{q}, &\text{if $x = \frac{p}{q}$ for $p \in \mathbb{Z},$ $q \in \mathbb{N}$, $\textsf{gcd}(p,q) = 1$ and $q \leq R\sqrt{2}$;}\\ 0, &\text{otherwise.}\end{cases}.$$

Now Thomae's function is fascinating with lots of pathology (e.g., integrable, nowhere differentiable, continuous at each irrational, discountinuous at each rational, etc.). However, here we are interested in using adjacent discontinuities of the approximate Thomae's function but measuring the distances between them in some completely different way, since ultimately we are measuring the angle between the $x$-axis and the point $(m,n),$ which has much more to do with $\frac{n}{m}$ than $\frac{m}{m+n}.$ So, though I would love to have this problem yield a deep dive into Thomae's function, instead, let's just brute force throw a computer at this problem and see what comes out, since

We need to just implement the following pseudocode procedure: (1) enumerate all of the visible trees, for instance, $$\textsf{trees} = \left[ (m,n) \mid n \in \{1, \dots, R\}, m \in \{n, \dots, R \}, \textsf{gcd}(m,n) = 1, m^2 + n^2 \leq R^2 \right];$$ (2) convert each tree to its bearing (angle above the $x$-axis), that is $$\textsf{bearings} = \left[ \textsf{arctan2(n,m)} \mid (m,n) \in \textsf{trees} \right];$$ (3) order the values of $\textsf{bearings}$ in ascending order and let $\textsf{gaps}$ be the corresponding gaps between adjacent trees; and (4) take an $\textsf{argsort}$ of $\textsf{gaps}$ to pull out the location of the largest gaps.

Below is a Python implementation of this pseudocode:

import numpy as np

def get_trees(R):
    return [ (m, n) for n in range(R) for m in range(n, R) if np.gcd(m,n) == 1 and m*m + n*n <= R*R ]

def get_tree_bearings(R):
    bearings = [ np.arctan2(n, m) for (m, n) in get_trees(R) ]
    return sorted(bearings), np.argsort(bearings)

def get_big_gap_bearings(R,k=7):
    trees = get_trees(R)
    bearings, orders = get_tree_bearings(R)
    gaps = [ x - y for x, y in zip(bearings[1:], bearings[:-1])]
    big_gaps = np.argsort(gaps)[-2*k:]
    return [(trees[orders[i]], trees[orders[i+1]]) for i in big_gaps]

For the Classic problem we notice that for $k=2$ and various values of $R = 10, 20, 50, 100, 200, 500, 1000,$ that we confirm that the largest gap always is next to the $x$-axis, followed by the gap nearest to the tree at $(1,1),$ then the next largest pair of adjacent gaps is on either side of the tree at $(2,1),$ that is, with bearing $$\theta^* = \tan^{-1} \left(\frac{1}{2}\right) \approx 26.565051177\dots^\circ.$$ This interestingly, though still mysteriously since I'm not totally sold on the ability to directly infer much of anything from Thomae's function, is equivalent to the second tallest peak of Thomae's function at $x^* = 2/3.$