Right off the bat, a reality check shows something amiss here.
Given ANY axis U = (xi, yj, zk) and a rotational angle alpha, the w (i.e., real) value of the resulting quaternion should be cos(alpha/2). But if alpha is 45 degrees, then w = cos(pi/8.0) should be around 0.924... and not somewhere about 0.963..., as you seem to expect when 'correctly' calculating ccc=quat([.3 .1 .6], 45).
Your error is that the axis rotation quaternion given axis U and rotation angle alpha assumes that U is already normalized (has length 1.0). So your formula should (in an ideal world, without any round-off error) be:
norm(U)*sin(alpha/2) + cos(alpha/2).
where norm(U) is U, 'normalized' to have length 1. What you have coded up is instead
norm(U*sin(alpha/2) + cos(alpha/2)).
As a test, try using your existing code to evaluate ccc = quat([300, 100, 600], 45). It will be much more 'off' than your example (where U is 'relatively' close to being normalized).
That should make sense: the quaternion that rotates 45 degrees around the axis (.3, .1, .6) should be exactly the same quaternion that rotates 45 degrees around the axis (300, 100, 600). In both caes, it's really the same axis: the latter is just a 'longer version' of the former.
Since we are not in an idealized world lacking round-off error, what you want to do is first, normalize the axis of rotation U; then construct the quaternion, and then normalize the quaternion:
norm(norm(U)*sin(alpha/2) + cos(alpha/2))
BTW, while it's quite readable, I don't recognize the syntax of the language your are working in. What language is it?
Okay, so you have an East-North-Down coordinate system as the fixed frame, and you have the phone's Up-Left-Out (of the screen) set of axes, and the quaternion is given in terms of the fixed coordinate system's basis vectors. I'll call the fixed axes $\hat x, \hat y, \hat z$ and the phone's axes $\hat x', \hat y', \hat z'$.
If you're using the phone as a camera, then I imagine it's the phone's Out-of-screen direction that is most like the direction the camera is facing (or rather, the negative of this). So the vector $-\hat z'$ should be the direction that the camera is facing. You should then project this vector onto the $xy$-plane by zeroing out its $\hat z$ component: $P_{xy}(\hat z') = \hat z' - \hat z (\hat z \cdot \hat z')$. You can then renormalize the result and find the angle in the $xy$-plane by basic trig. This should give you the compass direction that the camera is pointing (as an angle relative to east).
Angle off the horizon is easier. Take $-\hat z' \cdot \hat z = \cos \theta$. $\theta$ is then the angle off the vertical. $\theta - \pi/2$ should then give the angle you want: if $\theta = \pi/2$, then the camera is level. If it's larger, then the camera is tilted upwards. That checks out.
I admit, I'm not sure why you want to know if the device is horizontal to the ground; it's not clear to me how this relates to the direction of the camera, and so I'm not sure what answer to give for that part.
Best Answer
The simplest way to obtain the relative orientation is the integrating of kinematic equations Quaternion kinematics for the error-state KF (formula 107). All the explanations about quaternions are in the book.
Gyroscope measures angular velocity $\omega$, so the relative orientation you can evaluate by integrating (in real-time) the kinematic equation $\dot{q}=\frac{1}{2}q\circ\omega$, where normalized quaternion $q$ defines orientation of the body-frame relative to the initial frame. The disadvantage of this method is the result diverges with time (because of integration errors and precision of gyroscope).
There is a better approach uses other sensors.
If you want to represent relative orientation as a sequence of rotations around 3 axis you should learn a bit about the Euler angles. Actually the second angle $\beta$ should be always in range $-\frac{\pi}{2}..\frac{\pi}{2}$ or in $0..\pi$.