I've been working on a different approach to the problem of finding the average area of a cube's projection after seeing 3Blue1Brown's video a few months ago.
For a cube of side length $s$, I am getting the following for the average area:
$$\bar{A}=s^2\cdot\frac{\pi+6}{2\pi}\approx 1.45s^2$$
I'm aware this should be $3/2$, but want to figure out why my answer is wrong--in particular, I want to see if it's an error in computation or methodology.
First, let's consider a cube of side length $s$ centered at the origin with each face orthogonal to one of the three axes. The coordinates of the 8 vertices are given below:
$$\left(\pm\frac{s}{2},\pm\frac{s}{2},\pm\frac{s}{2}\right)$$
Let's describe each vertex as a vector pointing from the origin to said vertex as outlined below:
$$\vec v_1=\frac{s}{2}\langle 1, 1, 1\rangle$$ $$\vec v_2=\frac{s}{2}\langle -1, 1, 1\rangle$$ $$\vec v_3=\frac{s}{2}\langle 1, -1, 1\rangle$$ $$\vec v_4=\frac{s}{2}\langle 1, 1, -1\rangle$$ $$\vec v_5=\frac{s}{2}\langle -1, -1, 1\rangle$$ $$\vec v_6=\frac{s}{2}\langle -1, 1, -1\rangle$$ $$\vec v_7=\frac{s}{2}\langle 1, -1, -1\rangle$$ $$\vec v_8=\frac{s}{2}\langle -1, -1, -1\rangle$$
The idea is to use rotation matrices to rotate the cube by an azimuthal angle of $\theta\in[0,\pi/2)$ and zenithal angle $\phi\in[0,\pi/2)$. The idea is that each unique orientation of the cube is decribed by a unique $(\theta',\phi')$
Let
$$R_\theta = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix}$$
$$R_\phi = \begin{bmatrix} \cos\phi & -\sin\phi & 0 \\ \sin\phi & \cos\phi & 0 \\ 0 & 0 & 1 \end{bmatrix}$$
Let $\vec p_i$ denote the projection of $R_\theta R_\phi\vec v_i$ onto the $yz$ plane. Then,
$$\vec p_1=\frac{s}{2}\begin{bmatrix} 0 \\ \sin\phi + \cos\phi \\ -\sin\theta\cos\phi + \sin\theta\sin\phi + \cos\phi \end{bmatrix}$$
$$\vec p_8=\frac{s}{2}\begin{bmatrix} 0 \\ -\sin\phi - \cos\phi \\ \sin\theta\cos\phi - \sin\theta\sin\phi - \cos\phi \end{bmatrix}$$
$$\vec p_2=\frac{s}{2}\begin{bmatrix} 0 \\ -\sin\phi + \cos\phi \\ \sin\theta\cos\phi + \sin\theta\sin\phi + \cos\phi \end{bmatrix}$$
$$\vec p_7=\frac{s}{2}\begin{bmatrix} 0 \\ \sin\phi - \cos\phi \\ -\sin\theta\cos\phi - \sin\theta\sin\phi - \cos\phi \end{bmatrix}$$
$$\vec p_3=\frac{s}{2}\begin{bmatrix} 0 \\ \sin\phi - \cos\phi \\ -\sin\theta\cos\phi - \sin\theta\sin\phi + \cos\phi \end{bmatrix}$$
$$\vec p_6=\frac{s}{2}\begin{bmatrix} 0 \\ -\sin\phi + \cos\phi \\ \sin\theta\cos\phi + \sin\theta\sin\phi - \cos\phi \end{bmatrix}$$
$$\vec p_4=\frac{s}{2}\begin{bmatrix} 0 \\ \sin\phi + \cos\phi \\ -\sin\theta\cos\phi + \sin\theta\sin\phi - \cos\phi \end{bmatrix}$$
$$\vec p_5=\frac{s}{2}\begin{bmatrix} 0 \\ -\sin\phi - \cos\phi \\ \sin\theta\cos\phi - \sin\theta\sin\phi + \cos\phi \end{bmatrix}$$
Note the edges $\bar e_{73},\bar e_{74},\bar e_{78}$ are interior to the cube's projection, and we may calculate the area of each of the 3 parallelograms that combine to form the projection by using determinants as follows:
$$\left|\det\begin{bmatrix}\vec p_7-\vec p_3 && \vec p_7-\vec p_4\end{bmatrix}\right|=s^2\cos^2\phi$$
$$\left|\det\begin{bmatrix}\vec p_7-\vec p_3 && \vec p_7-\vec p_8\end{bmatrix}\right|=s^2\sin\phi\cos\phi$$
$$\left|\det\begin{bmatrix}\vec p_7-\vec p_4 && \vec p_7-\vec p_8\end{bmatrix}\right|=s^2\sin\theta$$
Then setting up a double integral for $\bar A$:
$$\bar A=\frac{1}{\pi/2}\cdot\frac{1}{\pi/2}\int_0^{\pi/2}\int_0^{\pi/2}s^2(\cos^2\phi+\sin\phi\cos\phi+\sin\theta)\;d\theta d\phi$$
$$=s^2\cdot\frac{\pi+6}{2\pi}$$
I ran a computer program from this methodology that seemed to converge to the same answer, so where did I go wrong?