Drawing random lines in a cylinder - How does the distribution of angles look like? We have a 3D cylinder with a radius of $r = 1$ and a height of $h = 2$. 

Now, we pick two random points on the edges of the cylinder - one for the top and one for the bottom circle. By connecting the two dots with a line what we can do next is calculate the angle between the line and the vertical z-axis.
We repeat this procedure $N$-times and collect all the random angles. 
What can we say for sure? 


*

*The minimum angle should be $0$°, the case where we pick two points laying on top of each other giving us a parallel line to the z-axis. 

*The maximum angle should be $45$° - which is decided for us by $r$ and $h$.



My question here is what distribution of angles we should get. Running an experiment will give me the following distribution:

Which is a bit counter intuitive - at least to me - I expected a uniform distribution to be honest. So, I do hope I did not make a mistake in my code, which is why I provide the experiment code below.
In any case I would be interested in whether this distribution makes sense and if it's possible to explain it.
One thing that surely is important: np.random.rand() draws uniform values $x \in [0,1]$ so the points on the edges (v and w) should be uniform too.
%matplotlib inline  # Jupyter Notebooks
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.cm as cm
import seaborn as sns
import numpy as np
import pandas as pd

def angle(a, b):    
    a_u = a / np.linalg.norm(a)
    b_u = b / np.linalg.norm(b)    
    c = np.clip(np.dot(a_u, b_u), -1, 1)    
    return np.degrees(np.arccos(c))

# Create the random points on the cylinder's edges
N = 5000
a = 2 * np.pi * np.random.rand(N) 
b = 2 * np.pi * np.random.rand(N)
v = np.vstack([np.cos(a), np.sin(a), 2*np.ones(N)]).T
w = np.vstack([np.cos(b), np.sin(b), np.zeros(N)]).T
df_v = pd.DataFrame(v, columns=['x', 'y', 'z'])
df_w = pd.DataFrame(w, columns=['x', 'y', 'z'])
# Get the connecting lines
df_u = df_v - df_w

# Calculate the angles of the lines w.r.t. the z-axis
e3 = np.asarray([0, 0, 1])
ds_angles = df_u.apply(lambda x: angle(x, e3), 1)

# Plotting the distribution of angles
ax = plt.figure(figsize=(8, 4.5)).gca()
sns.distplot(df_u.apply(lambda x: angle(x, e3), 1), ax=ax, kde=False)
plt.xlabel('Degree')
plt.show()

 A: Let $X$ denote the point in the circle below, and $Y$ denote the point in the circle above.
For any given $X$, with probability $1/2$ we would have that $Y$ lies in the semi-circle whose midpoint is directly above $X$.
This means with probability $1/2$, the angle lies between $0$ and $\arccos(2/\sqrt6)\simeq 35.26^{\circ}$.
Hence, the assumption that the distribution of the angles would be uniform is flawed, and about half the time the angles should lie between $35.26^{\circ}$ and $45^{\circ}$ (which is far from half of the possible angles).

Here's an attempt at calculating things more explicitly.
Let $X, Y$ be uniform in the circle $C=\left\{(x,y,0)\in\mathbb R^3\,|\,x^2+y^2=1\right\}$.
Your random line has direction $\big(Y+(0,0,2)\big) - X$.
With $W = Y-X$, $h=(0,0,2)$ and $Z$ the cosine of the angle between that line an the $z$-axis, we have
$$Z = \frac{\langle W+h, (0,0,1)\rangle}{\lVert W+h\rVert} = \frac2{\lVert W+h\rVert}.$$
This implies
$$Z^2 = \frac{4}{4 + {\lVert W\rVert}^2} = \frac1{1+{\left\lVert\frac{W}2\right\rVert}^2}.$$
If $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$, then $X\sim \big(\cos(\theta_1),\sin(\theta_1)\big)$ and $Y\sim \big(\cos(\theta_2),\sin(\theta_2)\big)$ so that $W\sim\big(\cos(\theta_1)-\cos(\theta_2),\sin(\theta_1)-\sin(\theta_2)\big)$.
Hence,
\begin{align}
{\left\lVert\frac{W}2\right\rVert}^2
&\sim
\frac14\left(
\underbrace{(\cos\theta_1)^2+(\cos\theta_2)^2}_1-2\cos\theta_1\cos\theta_2
+
\underbrace{(\sin\theta_1)^2+(\sin\theta_2)^2}_1-2\sin\theta_1\sin\theta_2
\right)
\\&=
\frac14\left(
2 - 2\cos\theta_1\cos\theta_2 - 2\sin\theta_1\sin\theta_2
\right)
\\&=
\frac12\big(
1 - (\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2)
\big)
\\&=
\frac12\big(
1 - \cos(\theta_1-\theta_2)
\big)
=
\frac12 - \frac{\cos(\theta_1-\theta_2)}2
.
\end{align}
Therefore
$$Z^2 \sim \frac{2}{3-\cos(\theta_1-\theta_2)}.$$
Observing that $Z\geq 0$, it follows that the angle $\alpha$ is distributed as
$$\alpha \sim \arccos\left(\sqrt{\frac{2}{3-\cos(\theta_1-\theta_2)}}\right),$$
where $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$.
From here on, we can calculate $\mathbb P(\alpha \leq x)$ for $0<x<\pi/4$ as some area in the square $[0,2\pi]^2$ contained in the $\theta_1$-$\theta_2$ plane.
Alternatively, we can observe that $\Theta = \theta_1-\theta_2$ has pdf with support in $[-2\pi,2\pi]$ given by
$$f_{\Theta}(z) =
\left\{\begin{array}{cccc}
\frac{2\pi+z}{4\pi^2}&&&-2\pi\leq z <0\\
\frac{2\pi-z}{4\pi^2}&&&0\leq z \leq 2\pi
\end{array}\right.
$$
Then, some algebraic manipulations show that
$$\mathbb P(\alpha \leq x)
=
\mathbb P\left(\cos(\Theta) \geq 3-\frac{2}{{\cos(x)}^2}\right)
.$$
Using the pdf above, we can calculate this explicity in terms of integrals that simplify to
$$\mathbb P(\alpha \leq x)
=
\frac{1}{\pi}\,\arccos\left(3-\frac{2}{{\cos(x)}^2}\right)
.$$
The density function $f_\alpha$ of $\alpha$ has support in $[0,\pi/4]$, and according to Wolfram we have (after simplifications considering its support) that it equals
$$f_{\alpha}(z) = \frac{2}{\pi}\,\frac{1}{\cos(z)\,\sqrt{\cos(2z)}},$$
which is way simpler than I expected, to be honest.
Observe that it's indeed a probability density function, in that it is non-negative and its integral is $1$.
Your histogram does somewhat resemble the area under $f_{\alpha}$, as you can verify here, so this looks more like a case of 'wrong intuition' than of 'wrong code'.

I must point out that this sheds light to yet another rather curious identity involving $\pi$:

$$\int_{\sqrt{2}/2}^1\,\,\,\frac{1}{u\,\sqrt{2u^2-1}\,\sqrt{1-u^2}}\,\,\,du = \frac\pi2$$
  $${}$$
  $$\int_{1/2}^1\,\,\,\frac{1}{u\,\sqrt{2u-1}\,\sqrt{1-u}}\,\,\,du = \pi$$

A: It's possible to find the closed form for the cdf. The angle is
$$A = \arccos \frac
 {((\cos v, \sin v, 1) - (\cos u, \sin u, -1)) \cdot (0, 0, 1)}
 {| (\cos v, \sin v, 1) - (\cos u, \sin u, -1) |} =
\arccos \sqrt {\frac 2 {3 - \cos(v - u)}}.$$
Switching to $\xi = v + u, \eta = v - u$, the cdf is
$$\operatorname P(A < x) = \int_0^{2\pi} \int_\eta^{4\pi - \eta}
 \frac 1 {(2 \pi)^2} [A < x] d\xi d\eta \overset * = \\
\int_0^\pi
 \frac 2 {(2 \pi)^2} (4\pi - 2\pi) [A < x] d\eta = \\
\frac 1 \pi \int_0^\pi [\eta \lt \arccos(3 - 2 \sec^2 x)] d\eta = \\
\frac 1 \pi \arccos(3 - 2 \sec^2 x),$$
where $[\text {condition}]$ is $1$ if the condition is true and $0$ otherwise. For $*$, $f$ is symmetric around $\pi$, therefore $\int_0^{2\pi} \eta f(\eta) d\eta = 2 \pi \int_0^\pi f(\eta) d\eta$.
