Angular distribution of lines passing through two squares. Let's say I've got two squares with side length $d$ that are held parallel at a distance $m$ apart.  
Suppose that particles are randomly falling from above in such a way that the polar angle $\varphi$ of the trajectory of the particles has a probability distribution proportional to $\operatorname{cos}^2(\varphi)$, while the azimuthal angle is uniformly distributed.
A particle is admissible if it passes through both squares.


Can we find a closed form expression for the probability distribution of the polar angle of admissible particles?  (or at least a good numerical solution?)


My attempt at a stronger formulation:
First, let's normalize the problem by letting $\alpha=m/d$ and assuming each square has unit side length.  We may assume that the point at which the particle passes through Panel $1$ (the upper panel) is uniformly distributed, that is, there are a uniformly distributed variables $X_1$ and $Y_1$, each between $0$ and $1$, and the point of contact on the plane of panel $1$ is $$\mathfrak{X}_1=\left(\begin{array}{c}X_1\\Y_1\\ \alpha\end{array}\right).$$  In order to be admissible, a particle must pass through Panel $1$, and we are only interested in admissible particles, so this is our universe of particles.  
Let's find the point of contact $\mathfrak{X}_2$ on the plane of Panel $2$ (the lower panel).  We need $$\left(\begin{array}{c}X_1+t\operatorname{sin}(\varphi)\operatorname{cos}(\theta)\\ Y_1+t\operatorname{sin}(\varphi)\operatorname{sin}(\theta)\\ \alpha+t\operatorname{cos}(\varphi)\end{array}\right)=\left(\begin{array}{c}X_2\\Y_2\\0\end{array}\right)$$
so $t=-\alpha\operatorname{sec}(\varphi)$, from which we have $$\mathfrak{X}_2=\left(\begin{array}{c}X_1-\alpha\operatorname{tan}(\varphi)\operatorname{cos}(\theta)\\Y_1-\alpha\operatorname{tan}(\varphi)\operatorname{sin}(\theta)\\0\end{array}\right).$$
So, the question is, what is the probability distribution of $\varphi$ for the lines satisfying
$$\begin{array}{rcl}0 \leq& X_1-\alpha\operatorname{tan}(\varphi)\operatorname{cos}(\theta) &\leq  1 \\
0 \leq& Y_1-\alpha\operatorname{tan}(\varphi)\operatorname{sin}(\theta) &\leq 1 \end{array}
$$
given that $$\begin{array}{rl}X_1 \text{ and }Y_1& \text{are uniformly distributed from }0\text{ to }1 \\ \theta &\text{is uniformly distributed from }0\text{ to }2\pi \\ \varphi & \text{has probability distribution }\frac{4}{\pi}\operatorname{cos}^2\varphi\text{ from }\frac{\pi}{2}\text{ to }\pi\end{array}$$
?

EDIT 2:
So I ran a simulation, yielding the following histogram (displayed below with $3$ different bin widths).  I haven't fit a curve to this yet, but maybe it will give us a hint.  Here $\alpha$ is about $6$.

 A: 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)$:

A: 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.
