Have done lot of googling on this and am overwhelmed by the number of formulas; not a math major, just a developer struggling to understand quaternion rotations:)
Given a vector of: 0, 0, 0.15
And a quaternion: 0.86671, -0.40654, -0.21285, 0.19555 (s, x, y, z)
Vector should equal: -0.07919, 0.09322, 0.08683
Any help to understand what formulas should be applied to get to this result is appreciated.
On
You are calculating $p\mathbf{x}p^{-1}$, the conjugation of $\mathbf{x}$ by $p$, where
$$ \mathbf{x}=0.15\mathbf{k}, \qquad p=0.86671-0.40654\mathbf{j}-0.21285\mathbf{j}+0.19555\mathbf{k}. $$
Since $p$ is a unit quaternion, $p^{-1}=\overline{p}=0.86671+0.40654\mathbf{j}+0.21285\mathbf{j}-0.19555\mathbf{k}$.
There are many packages in programming languages and computer algebra systems that let you do calculations like these with quaternions. I explain quaternions below.
A quaternion can be thought of as a scalar plus a 3D vector (also known as real and imaginary parts). The product of a scalar and a 3D vector is the usual scalar multiplication. The product of two vectors produces a quaternion with both scalar and vector components, given by (minus) the dot product and cross product respectively. In other words, for 3D vectors $\mathbf{u}$ and $\mathbf{v}$, their quaternion product is
$$ \mathbf{uv}=-\mathbf{u}\cdot\mathbf{v}+\mathbf{u}\times\mathbf{v} $$
A fun way to think about this is that each of $\mathbf{i},\mathbf{j},\mathbf{k}$ is a square root of $-1$, and that multiplying two of them in cyclic order yields the third (e.g. $\mathbf{ij}=\mathbf{k}$) or out of order is the opposite (e.g. $\mathbf{ji}=-\mathbf{k}$). As a good exercise, you can discover the only square roots of $+1$ are $\pm1$; the square roots of $-1$ are any unit vector; two quaternions commute ($pq=qp$) if and only if their vector parts are parallel (i.e. scalar multiples of each other); and two quaternions anticommute ($qp=-pq$) if and only if they are perpendicular vectors.
The norm $|a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k}|=\sqrt{a^2+b^2+c^2+d^2}$ is multiplicative, i.e. $|pq|=|p||q|$. The quaternion conjugate of $p=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k}$ is $\overline{p}=a-b\mathbf{i}-c\mathbf{j}-d\mathbf{k}$, i.e. the vector part is negated. This is useful in calculating the multiplicative inverse of a quaternion, $p^{-1}=\overline{p}/|p|^2$.
Because all unit vectors are square roots of $-1$, we have Euler's formula $\exp(\theta\mathbf{u})=\cos\theta+\sin\theta\,\mathbf{u}$, which means every quaternion has a polar form $p=r\exp(\theta\mathbf{u})$ where $r=|p|\ge0$ is the magnitude, $\mathbf{u}$ is the vector part of $p$ normalized to be a unit vector, and $0\le\theta\le\pi$ is a convex angle. You calculate the polar form just as you would for a complex number (after normalizing the vector part for $\mathbf{u}$). For nonreal quaternions the polar form is unique with $0<\theta<\pi$, otherwise for real nonzero quaternions, $\theta$ is $0$ or $\pi$ for positive/negative reals respectively, and $\mathbf{u}$ is an arbitrary unit vector.
Given a quaternion $p$, left-multiplication $L_p(x):=px$ and right-multiplication $R_p(x):=xp$ are linear transformations when considered as functions of $x$, if we treat the quaternions as a 4D real vector space. Because the norm is multiplicative, these are orthogonal transformations (aka linear isometries, or 4D rotations) for unit quaternions $p$ (i.e. with $|p|=1$).
If $\mathbf{u}$ is a unit vector, it can be extended to an orthonormal basis $\{\mathbf{u},\mathbf{v},\mathbf{w}\}$ with the same orientation as the standard basis $\{\mathbf{i},\mathbf{j},\mathbf{k}\}$ (i.e. $\mathbf{w}=\mathbf{u}\times\mathbf{v}$ according to the right-hand rule, just as $\mathbf{k}=\mathbf{i}\times\mathbf{j}$). Then $\{1,\mathbf{u},\mathbf{v},\mathbf{w}\}$ is an orthonormal basis for the quaternions as a real vector space. We can calculate the $4\times4$ matrix representing $L_p$ and $R_p$ with respect to this basis fairly easily; both are block-diagonal with two $2\times2$ blocks which are rotation matrices with angle $\theta$, except the second block of $R_p$ uses $-\theta$ instead of $\theta$.
The cumulative effect is that "conjugating" a quaternion by $p$, i.e. $(L_p\circ R_p^{-1})(x)=pxp^{=1}$ rotates in the $\mathbf{vw}$-plane by $2\theta$ and fixed the $1\mathbf{u}$-plane pointwise. In other words, if $\mathbf{x}$ is a 3D vector, then $p\mathbf{x}p^{-1}$ is a 3D vector which is the rotation of $\mathbf{x}$ around $\mathbf{u}$ by an angle of $2\theta$. The fact there is a "$2$" in this formula is a manifestation of spin. It implies every 3D rotation is represented by two antipodal unit quaternions $\pm p$.
On
A quaternion can be faithfully represented by a real $4\times 4$ matrix, i.e. $$\eqalign{ \def\bbR#1{{\mathbb R}^{#1}} \def\o{{\tt1}}\def\p{\partial} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\v#1{\begin{Bmatrix}#1\end{Bmatrix}} \def\qiq{\quad\implies\quad} \def\qif{\quad\iff\quad} \def\c#1{\color{red}{#1}} q &= \v{s\\x\\y\\z} \qif Q = \m{ s&-x&-y&-z \\ x&s&-z&y \\ y&z&s&-x \\ z&-y&x&s \\ } \\ }$$ Note that $q$ appears in the first column of $Q$ and the conjugate $\:\bar q\iff Q^T$
You can now use ordinary matrix multiplication to address your problem: $$\eqalign{ Q &= \m{ 0.86671 & 0.40654 & 0.21285 & -0.19555 \\ -0.40654 & 0.86671 & -0.19555 & -0.21285 \\ -0.21285 & 0.19555 & 0.86671 & 0.40654 \\ 0.19555 & 0.21285 & -0.40654 & 0.86671 \\ } \\ V &= \m{ 0 & 0 & 0 & -0.15 \\ 0 & 0 & -0.15 & 0 \\ 0 & 0.15 & 0 & 0 \\ 0.15 & 0 & 0 & 0 \\ } \\ V_{rot} &= QVQ^T = \m{ 0 & 0.07919 & -0.09322 & -0.08683 \\ \c{-0.07919} & 0 & -0.08683 & 0.09322 \\ \c{0.09322} & 0.08683 & 0 & 0.07919 \\ \c{0.08683} & -0.09322 & -0.07919 & 0 \\ } \\ }$$ The desired vector can be read from the $\c{\rm first}$ column of $V_{rot}$
A further advantage of the matrix representation becomes apparent when calculating quaternionic functions. Formulas for exp(q), log(q), sqrt(q), etc are seemingly bespoke one-offs, whereas the theory of matrix functions is very mature and is a standard feature of programming languages like Matlab and Julia.
Formula
The formula to calculatise would be $$\mathrm{v}_{\text{rotated}} = \mathrm{q} \cdot \mathrm{v} \cdot \overline{\mathrm{q}}$$ with your vector $$\mathrm{v} = \left[\begin{matrix}0\\ 0\\ 0.15\end{matrix}\right] = 0 \cdot \mathrm{i} + 0 \cdot \mathrm{j} + 0.15 \cdot \mathrm{k}$$ and your quaternion $$\mathrm{q} = 0.86671 - 0.40654 \cdot \mathrm{i} - 0.21285 \cdot \mathrm{j} + 0.19555 \cdot \mathrm{k}$$ and $$\overline{\mathrm{q}} = 0.86671 + 0.40654 \cdot \mathrm{i} + 0.21285 \cdot \mathrm{j} - 0.19555 \cdot \mathrm{k}$$ with $\mathrm{i}^{2} = \mathrm{j}^{2} = \mathrm{k}^{2} = \mathrm{i} \cdot \mathrm{j} \cdot \mathrm{k} = -1$.
But don't forget that quaternions are not commutative with respect to multiplication.
Formula with roationangle and rotationaxis
If you use the vector / quaternion $r$ as the axis of rotation, $\varphi$ as the rotation angle, $\mathrm{v}$ as the vector to be rotated and the rotated vector $\mathrm{v}_{\text{rotated}}$, then the formula applies: $$ \begin{align*} \mathrm{v}_{\text{rotated}} &= \mathrm{q} \cdot \mathrm{v} \cdot \overline{\mathrm{q}}\\ \mathrm{v}_{\text{rotated}} &= \left( \cos\left( \frac{\varphi}{2} \right) + \mathrm{u} \cdot \sin\left( \frac{\varphi}{2} \right) \right) \cdot \mathrm{v} \cdot \left( \cos\left( \frac{\varphi}{2} \right) - \mathrm{u} \cdot \sin\left( \frac{\varphi}{2} \right) \right)\\ \mathrm{v}_{\text{rotated}} &= \left( \cos\left( \frac{\varphi}{2} \right) + \frac{\mathrm{r}}{||\mathrm{r}||} \cdot \sin\left( \frac{\varphi}{2} \right) \right) \cdot \mathrm{v} \cdot \left( \cos\left( \frac{\varphi}{2} \right) - \frac{\mathrm{r}}{||\mathrm{r}||} \cdot \sin\left( \frac{\varphi}{2} \right) \right)\\ \end{align*} $$
Code
Since you are a developer, my code for calculating rotation using quaternions might tell you more then the formulas:
HTML
JavaScript
CSS
free to use for everyone ;)