command to expand quaternion multiplication

Is there are any command to expand quaternion multiplication as below STEP by STEP.
Below is equation from my text.
I'm self learner and for some reason I can't wrap my head around quaternion multiplication.
I really need to see step by step workout.
I remember seeing step by step workout in wolframalpha but that's for calculus integration.
Just wondering if the same can be done in mathematica.
$$$$\begin{split} w =& qvq^* \\ =& (q_0+\vec{q})(0+\vec{v})(q_0-\vec{q}) \\ =&(2q_0^2-1)\vec{v}+2(\vec{q}\cdot\vec{v})\vec{q} + 2q_0(\vec{q}\times\vec{v}) \end{split}$$$$

-
The last equation is incorrect. This becomes obvious when you consider the special case $q_0+\vec{q}=0$ and $\vec{v}\ne \vec{0}$: then $w$ must be $0$, but the last expression equals $-\vec{v}$. –  whuber Oct 26 '12 at 13:17
the book mention about using triple vector product $a\times(b\times c)=(a\cdot c)b-(a\cdot b)c$ to derive the final equation. –  kypronite Oct 26 '12 at 15:14
That's fine--but the answer is still obviously wrong, as my counterexample shows. –  whuber Oct 26 '12 at 16:29

You can create quaternions in Mathematica. Here is a limited implementation, sufficient for the kinds of calculations shown in the question. As in the question, it represents quaternions as ordered pairs of a real part and a (vector) quaternionic part and defines multiplication in terms of dot and cross products:

ClearAll[quaternion];
quaternion[q0_, q_] + quaternion[t0_, t_] ^= quaternion[q0 + t0, q + t];
quaternion[q0_, q_] ** quaternion[t0_, t_] ^= quaternion[q0 t0 - q.t, q0 t + q t0 + Cross[q, t]];
Conjugate[quaternion[q0_, q_]] ^= quaternion[q0, -q];

(* Some helpful procedures to let reals interact with quaternions: *)
quaternion[q0_] := quaternion[q0, {0, 0, 0}];
Re[quaternion[q0_, _]] ^= q0;
Im[quaternion[_, q_]] ^= q;


The question asks for the purely quaternionic part of $q v q^*$ when--apparently--$v$ is a pure quaternion (with zero real part).

Im[quaternion[q0, q] ** quaternion[0, v] ** Conjugate[quaternion[q0, q]]]


$q0\ (q0\ v+q\times v)+(q0\ v+q\times v)\times (-q)+q\ q\cdot v$

Now the working out can be traced:

Clear[q0, q, v];
traceView2[quaternion[q0, q] ** quaternion[0, v] ** Conjugate[quaternion[q0, q]]]


(This is an excerpt from deeper into the trace.)

-