Posting the question yesterday has focussed my thoughts on this problem, and I think that I've now come up an answer. In terms of the rotation matrix $\mathbf{R}$ that I defined in the question, I think that the best definition of the azimuth $\phi$ (i.e. the angle through which the phone must be rotated so that it's pointing north) is given by
$$
\phi=\tan^{-1}\left( \frac{E_{y}-N_{x}}{E_{x}+N_{y}}\right) \text{ .}
$$
I give a full explanation below, but one of my main concerns is to address what people expect when using a smartphone. If a smartphone is lying flat on a table and there's a compass arrow pointing north in the direction of the longest edge of the phone as shown in the following image
then if the phone is rotated from its flat horizontal position into a vertical position by rotating along the short edge of the phone, the compass directions and hence the azimuth $\phi$ shouldn't change. Of course, with the phone upright the arrow that was pointing North when the phone was flat is now actually pointing away from the earth and straight into space, but people interpret that upward pointing arrow as pointing straight ahead in a horizontal direction. I've tested the above formula for $\phi$ in my android devices and as far as I can tell it works well. (NB: I've posted the relevant piece of the code that I'm using here on stackoverflow.com.)
In devising my answer, I use three angles which I call azimuth $\phi$, pitch $\theta$ and pitch axis $\psi$
$$
\begin{array}
[c]{lll}%
\text{azimuth} & =\phi & =\text{rotation about }\mathbf{G}\\
\text{pitch} & =\theta & =\text{rotation about axis at angle }\psi\text{ in
the horizontal }\mathbf{E}\text{-}\mathbf{N}\text{ plane}\\
\text{pitch axis} & =\psi & =\text{angle in horizontal plane for pitch
rotation}
\end{array}
$$
Whereas the three Euler angles correspond to three rotations, the $\phi$, $\theta$ and $\psi$ here only correspond to two rotations, namely
$$
\mathbf{R}=\mathbf{R}_{\mathbf{G}}\left( \phi\right) \mathbf{R}
_{\mathbf{EN}}\left( \theta|\psi\right)
$$
where
$$
\begin{align*}
\mathbf{R}_{\mathbf{G}}\left( \phi\right) & \mathbf{=}\left[
\begin{array}
[c]{ccc}%
\cos\phi & \sin\phi & 0\\
-\sin\phi & \cos\phi & 0\\
0 & 0 & 1
\end{array}
\right] \\
\mathbf{R}_{\mathbf{EN}}\left( \theta|\psi\right) & \mathbf{=}\left[
\begin{array}
[c]{ccc}
\cos^{2}\psi+\sin^{2}\psi\cos\theta & -\sin\psi\cos\psi\left( 1-\cos
\theta\right) & \sin\psi\sin\theta\\
-\sin\psi\cos\psi\left( 1-\cos\theta\right) & \sin^{2}\psi+\cos^{2}\psi
\cos\theta & \cos\psi\sin\theta\\
-\sin\psi\sin\theta & -\cos\psi\sin\theta & \cos\theta
\end{array}
\right]
\end{align*}
$$
so that doing the matrix multiplication results in a matrix $\mathbf{R}$ which is given by
$$
{\scriptsize
\left[
\begin{array}
[c]{ccc}%
\cos\psi\cos\left( \phi+\psi\right) +\sin\psi\sin\left( \phi+\psi\right)
\cos\theta & -\sin\psi\cos\left( \phi+\psi\right) +\cos\psi\sin\left(
\phi+\psi\right) \cos\theta & \sin\left( \phi+\psi\right) \sin\theta\\
-\cos\psi\sin\left( \phi+\psi\right) +\sin\psi\cos\left( \phi+\psi\right)
\cos\theta & \sin\psi\sin\left( \phi+\psi\right) +\cos\psi\cos\left(
\phi+\psi\right) \cos\theta & \cos\left( \phi+\psi\right) \sin\theta\\
-\sin\psi\sin\theta & -\cos\psi\sin\theta & \cos\theta
\end{array}
\right]}
$$
The idea of these coordinates is that there's always an Euler angle type of co-ordinate system where the roll angle is zero. If the roll angle is always zero it can be discarded, but then one needs an additional angle to specify the axis for the pitch rotation. Using Euler angles to describe the rotation of a smartphone out of the horizontal $\mathbf{E}$- $\mathbf{N}$ plane, the pitch corresponds to a rotation out of the plane along one particular axis, the roll corresponds to a rotation out of the plane along a second perpendicular axis, and rotations along other axes are described by a combination of pitch and roll. With the $\left( \phi,\theta,\psi\right) $ coordinates presented here there's just one angle $\theta$ to describe rotation out of the plane, but then the angle $\psi$ is required to describe which axis to rotate about. In the context of smartphones, the pitch angle $\theta$ just represents the phone's tilt out of the horizontal plane, irrespective of which axis it has been tilted along.
Thinking like this, any 3D rotation can be constructed from a rotation $\mathbf{R}_{\mathbf{G}}\left( \phi\right) $ about the vertical axis $\mathbf{G}$, and a second rotation $\mathbf{R}_{\mathbf{EN}}\left( \theta|\psi\right) $ to take the device of the plane perpendicular to $\mathbf{G}$. With $\mathbf{R}$ specified above as a function of $\phi$, $\theta$ and $\psi$, it is straightforward to show that
$$
\begin{align*}
\cos\phi\left( 1+\cos\theta\right) & =E_{x}+N_{y}\\
\sin\phi\left( 1+\cos\theta\right) & =E_{y}-N_{x}
\end{align*}
$$
and hence the $\tan^{-1}$ formula for $\phi$ given above. It is worth noting
that
$$
\begin{align*}
\sin\left( \phi+2\psi\right) \left( 1-\cos\theta\right) & =-E_{y}-N_{x}\\
\cos\left( \phi+2\psi\right) \left( 1-\cos\theta\right) & =E_{x}-N_{y}\text{ .}
\end{align*}
$$
Also note that if $\psi\rightarrow\psi+\pi$ and $\theta\rightarrow-\theta$ then $\mathbf{R}$ is unchanged. Whereas the Euler angle system is degenerate when its $\theta=\pm\frac{\pi}{2}$, this $\left( \phi,\theta,\psi\right) $ co-ordinate system is degenerate at $\theta=0$ and $\theta=\pi$. The $\theta=0$ degeneracy occurs when the device is lying flat on the table facing upwards, and in that situation no rotation out of the horizontal $\mathbf{E}$-$\mathbf{N}$ plane is required so $\psi$ is undefined. However, $\phi$ is well defined, so if $\phi$ is the goal then the fact that $\psi$ is undefined does not matter.
The $\theta=\pi$ degeneracy occurs when the device is lying face down on the table, when $E_x+N_y=E_y-N_x=0$, and this is a more interesting case. In this situation, $\phi$ is undefined, although $\phi+2\psi$ is defined. The purpose of this $\left(\phi,\theta,\psi\right) $ co-ordinate system is to capture the idea that people expect the compass direction (i.e. the azimuth $\phi$) to be unchanged when rotating a phone from a horizontal $\theta=0$ position to a vertical $\theta=\frac{\pi}{2}$ position. However, if the rotation continues and $\theta$ increases beyond $\frac{\pi}{2}$, then as the phone approaches the horizontal $\theta=\pi$ upside down position, a compass arrow that had been pointing north when the device was horizontal will now be pointing south if it's direction remains fixed on the smartphone display. So as the $\theta=\pi$ degeneracy is approached, the smartphone is no longer behaving in an appropriate way. One way of dealing with this would be to switch co-ordinate systems at some point between $\frac{\pi}{2}$ and $\pi$ which could alter the direction of the arrow by $\pi$ and make the smartphone behave in an appropriate way again. However, it is worth noting that almost all the time when people are using their phones $\left\vert \theta\right\vert\le\frac{\pi}{2}$, so the inappropriate behaviour near $\theta=\pi$ when the phone is upside down should not be a problem in practice.
I've searched the internet and had a look at the article Shuster, M., "A Survey of attitude representations", Journal of the Astronautical Sciences 41(4):1993 that was suggested in another answer, but I haven't been able to find any references to the $\left( \phi,\theta,\psi\right) $ co-ordinate system defined here. In such a mature area as 3D rotations, it seems unlikely that no one has defined rotations in this way before? However, the emergence of smartphones has provided a new application for rotation matrices, so perhaps there haven't previously been applications where this type of co-ordinate system has been appropriate?
Let me cite a theorem by Ribando (Measuring Solid Angles Beyond Dimension Three, Discrete Comput Geom 36:479–487 (2006)), which is a rediscovery of a result of Aomoto (Analytic structure of Schlafli function, Nagoya Math. J. 68:1-16 (1977)) as described in the paper by Beck, Robins and Sam (Positivity Theorems for Solid-Angle Polynomials, Beitrage zur Algebra und Geometrie, Vol. 51, No. 2, 493-507 (2010)). I cite Ribando's result as its statement more to my taste:
Let $\Omega \subseteq \Bbb{R}^n$ be a solid-angle spanned by unit vectors $\lbrace v_1 , \dots , v_n \rbrace$, let $V$ be the matrix whose ith column is $v_i$ , and let $\alpha _{ij} = v_i \cdot v_j$ as above. Let $T_{\alpha}$ be the following infinite multivariable Taylor series: $$T_{\alpha} = \dfrac{det \ V}{(4 \pi )^{n/2}} \sum _{a \in \Bbb{N}^{{n \choose 2}}} \left[ \dfrac{(-2)^{\sum _{i < j} a_{ij}}}{ \Pi _{i<j} a_{ij}!} \Pi _{i} \Gamma \left( \dfrac{1 + \sum _{m \neq i} a_{im}}{2} \right) \right] \alpha^{a}$$ The series $T_{\alpha}$ agrees with the normalized measure of solid-angle $\Omega$ whenever $T_{\alpha}$ converges.
Best Answer
Hint:
In the point $E$ are concurrent four sides of the pyramid $\overline{EA}$, $\overline{EB}$, $\overline{EC}$, $\overline{ED}$.
If you want find the angles between any two of them you have to find the vectors parallel to the sides, e.g. $$ \vec {EA}=(-x,-4-y,4)^T \qquad \vec {EB}=(-x,4-y,4)^T $$ than the angle $\theta$ between them is given can be found by means of the dot product:
$$ \theta=\arccos\left(\frac{\vec{EA}\cdot \vec{EB}}{|\vec {EA}||\vec{EB}|} \right) $$