Angular distribution of lines passing through two squares.

Solution 1:

Let $d=\alpha \tan (\varphi) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$ be the distance (projected over the $z$ axis) between the point points at which the ray intersect the square planes. (Disregard for now the restriction that the ray passes through the squares).

Then, the probability density of $d$ can obtained, by the usual transformation of variables, from the density of $\varphi$, as

$$P(d) \propto \frac{1}{(\alpha^2 +d^2)^2} \,\, (d\ge 0)$$

Now, instead of considering just one pair of unit squares, let's imagine that the planes are tiled by a grid of such squares, and let $E$ be the event that the rays passes through a pair or aligned squares. Then, we are insterested in computing

$$P(d|E) = \frac{P(E|d) P(d)}{P(E)}$$

The denominator is just a normalization constant, so we are left with computing $P(E|d)$ Because the azimuthal angle is uniform, fixing $d$ means considering a circle of that radius and random (uniform) center, with uniform distribution along its circumference. We want to compute the probability that a point of that circle falls inside the same square as its center. Which is the same of computing the fraction of its circumference that falls inside the unit square (conditioned on the center coordinates, and integrated uniformly over the square). Because of symmetry, it's enough to consider the west-south quarter of circle; lets compute the fraction of this quarter-of circumference that falls outside the unit square. We have two ranges of interest: for $0\le d<1$ we have to consider 5 zones, and for $1\le d <\sqrt{2}$ we have 2 zones.

Lets deal with $0\le d<1$, where $P(E|d) = 1 - (I_A+2 I_{Bx} +I_C+ I_D)$. When the center of the circle falls in zone A, : the quater of circumference fall all inside the square (so $I_A=0$). At zone D, it falls all outside, the fraction is 1, which integrated over the region gives a contribution of $I_D=\pi d^2 /4$

Zone C is more complex:

The relevant fraction is the ratio of the red arcs over $\pi/2$, which must be integrated over the zone. This is

$$ \frac{2}{\pi} \int_C acos(x/d)+acos(y/d) dx dy = \frac{4 d^2}{\pi} \int_0^1(1-\sqrt{1-t^2}) acos(t) \, dt = \frac{d^2(12-\pi^2)}{4 \pi}$$

Similarly:

$$I_{BX}= I_{BY}= \frac{2}{\pi} \int_d^1 \int_0^d acos(y/d) dy dx = \frac{2}{\pi} (1-d) d $$

So $$P(E|d) = g_1(d) \triangleq 1 - (I_A+2 I_{Bx} +I_C+ I_D) = 1 + \frac{d^2-4d}{\pi}$$.

(The range $1 < d < \sqrt{2}$, though we have less zones, it's more complicated [*] )

From that, we get the form of $P(\varphi,E)$ in the range $0\le d\le 1$ , or $0 \le \varphi \le atan (1/\alpha)$ :

$$P(\varphi|E) \propto \cos^2 \varphi \, (\pi - 4 \alpha \tan \varphi + \alpha^2 \tan^2 \varphi) = (\pi -\alpha^2) \cos^2 \varphi - 2 \alpha \cos(2 \varphi) + \alpha^2 $$

This (unless errors) is exact in that range, but to get the proportionality factor one would need to consider also the remaining (small) tail, on $1 < d < \sqrt{2}$.

[*] In the "tail zone", $P(E|d) = 1 - (I_E+I_{F} )$ with

$$ I_E = \frac{4}{\pi} \int_E acos(x/d) dy dx = \frac{4}{\pi} \int_{\sqrt{d^2-1}}^1 (1-\sqrt{d^2-x^2}) acos(x/d) dx$$

$$ I_F = \int_F 1 dy dx = \sqrt{d^2-1} + \int_{\sqrt{d^2-1}}^1 \sqrt{d^2-x^2}) dx$$ 

$$P(E|d)= g_2(d) \triangleq \frac{1}{2 \pi}\left[ -2\,{d}^{2}\,{\mathrm{asin}\left( \frac{\sqrt{{d}^{2}-1}}{d}\right) }^{2}+\left( \pi \,{d}^{2}-4\,{d}^{2}\,\mathrm{acos}\left( \frac{\sqrt{{d}^{2}-1}}{d}\right) \right) \,\mathrm{asin}\left( \frac{\sqrt{{d}^{2}-1}}{d}\right) +\sqrt{{d}^{2}-1}\,\left( 4\,\mathrm{acos}\left( \frac{\sqrt{{d}^{2}-1}}{d}\right) +4\,\mathrm{acos}\left( \frac{1}{d}\right) -2\,\pi +8\right) +\left( 2\,{\mathrm{asin}\left( \frac{1}{d}\right) }^{2}+\left( 4\,\mathrm{acos}\left( \frac{1}{d}\right) -\pi \right) \,\mathrm{asin}\left( \frac{1}{d}\right) -2\right) \,{d}^{2}-8\,\mathrm{acos}\left( \frac{1}{d}\right) +2\,\pi -4 \right]$$

(Sorry, Maxima cannot simplify this further). The integrals of $P(E|d)$ over the two regions are $1-5/(3 \pi) = 0.469483523$ and $0.00371748$ respectively. Here's a graph of $P(E|d)$:

Solution 2:

Use Bayes' theorem: $$ P(\phi\mid\text{admissible}) = \frac{P(\text{admissible}|\phi)P(\phi)}{P(\text{admissible})}. $$

Unless I am misreading your question (what do you do about rays that are parallel to an admissible ray? I am assuming here they are admissible, so only the direction or a ray matters), a ray parametrized by $$ (t \sin\phi \cos\theta, t \sin\phi \sin\theta, t \cos\phi), \qquad t\in\mathbb{R}, \quad 0\leq\phi\leq\frac\pi2, \qquad 0\leq\theta\leq2\pi $$ is admissible if at the point $m = z(t) = t\cos\phi$, we have $$ |t\sin\phi\cos\theta| \leq d, \quad\text{and}\quad |t\sin\phi\sin\theta| \leq d, $$ because if it passes through the two squares, the most it can move at most $d$ in either $x$- or $y$-direction as it moves a distance $m$ in $z$-direction. This condition is equivalent to $$ |\tan\phi\cos\theta| \leq d/m, \qquad |\tan\phi\sin\theta| \leq d/m = \alpha, $$ which is equivalent to $$ |\tan\phi|^2 \leq 2\alpha^2 = (\tan\phi_0)^2, \qquad \max(|\cos\theta|,|\sin\theta|) \leq \frac{\alpha}{|\tan\phi|}. $$

Finally, $ P(\phi) = \frac{4}{\pi}\cos^2\phi$, and, setting $\beta=\alpha/|\tan\phi|$, $P(\text{admissible}) = \text{const}$, and $$ P(\text{admissible}|\phi) = [\phi\leq\phi_0] P\left( \max(|\cos\theta|,|\sin\theta|) \leq \beta \right) \\ = [\phi\leq\phi_0] 4 P(\theta\geq\arccos\beta, \theta\leq\arcsin\beta\mid0\leq\theta\leq\pi/2) \\ = [\phi\leq\phi_0]([\beta\geq1] + \frac{2}{\pi}(\arcsin\beta-\arccos\beta)[\beta\leq1]). $$

Putting things together, $$ P(\phi\mid\text{admissible}) \propto \cos^2\phi \times [\phi\leq\phi_0]([\beta\geq1]+\frac{2}{\pi}[\beta\leq1](\arcsin\beta-\arccos\beta)). $$ The constant of proportionality is found by integrating the rhs over $0\leq\phi\leq\frac\pi2$ and setting it to $1$, but probably doesn't have a closed form.