This came out a bit long, so first you may want to scroll down to
the algorithm, skipping the explanations and the justifications. I
have written a discrete dynamical model of the case when the
asteroid rotates freely in space, twisting and turning, with
precessing angular velocity (so there is angular acceleration),
without any thrusting boosters. As you will see, the dynamics is non-trivial and interesting enough.
I can suggest a sketch but you may have to try to figure out the
details. I will use standard $3$ dimensional vectors, orthogonal
and skew symmetric matrices but one can do the absolutely the same
thing using quaternions. I am not sure which one is
computationally cheaper or more convenient but maybe there is not
that much difference. Ether way, the philosophy is always the same
and conceptually everything looks the same regardless of the usage
of quaternions or $3$D vectors and matrices.
We will be working a lot with the set of all orthogonal matrices
(also known as the orthogonal group) $$\text{SO}(3) = \Big\{ \, U
\in \mathbb{R}^{3 \times 3} \,\, \big{|} \,\, U\, U^T = U^T U = I
\, \Big\}$$ where $U^T$ is the transpose matrix of $U$.
Furthermore, we will be working with the set of all skew-symmetric
matrices (also known as the Lie algebra of $\text{SO}(3)$)
$$\text{so}(3) = \Big\{ \, A \in \mathbb{R}^{3 \times 3} \,\, \big{|}
\,\, A^T = -\, A \, \Big\}$$
A matrix $A$ is skew-symmetric if and only if
$$A = \begin{pmatrix} \, 0 & z & -y \\
-z & 0 & x \\
y & -x & 0 \, \end{pmatrix}$$
Thus, for the sake of clarity, I introduce two functions:
$$\overrightarrow{V} = \text{vector}(V)$$ which turns a skew-symmetric matrix
$$V = \begin{pmatrix} \, 0 & V_3 & -V_2 \\
-V_3 & 0 & V_1 \\
V_2 & -V_1 & 0 \, \end{pmatrix}$$
into the vector column
$$\overrightarrow{V} = \begin{pmatrix} V_1 \\
V_2 \\
V_3 \end{pmatrix} $$
and its inverse function
$$V = \text{Matrix}(\overrightarrow{V})$$
which turns a vector column $\overrightarrow{V}$ into the
corresponding skew-symmetric matrix $V$ with. We also will use the
notation
$$|V| = |\overrightarrow{V}| = \sqrt{V_1^2 + V_2^2 + V_3^2}$$ which is the
length of the vector $\overrightarrow{V}$, and we can also think
of it as the magnitude $|V|$ of the skew-symmetric matrix $V$.
The shape of the asteroid and the distribution of its mass is
encoded in a $3 \times 3$ matrix ${\bf \mathcal{I}}$, called the
inertia tensor of the asteroid.
First, you have the coordinate system $S = O\,
\vec{e}_1\,\vec{e}_2\,\vec{e}_3$ fixed in space, according to
which you determine the positions of the point $P$ that fires upon
the asteroid and the position of the center of mass of the
asteroid $A$. Furthermore, you have another coordinate system
$S_A=A\, \vec{E}_1\,\vec{E}_2\,\vec{E}_3$, firmly attached to the
asteroid. Thus, the coordinate system $S_A$ rotates together with
the asteroid. Despite the fact that the system $S_A$ is
non-inertial, it is the right choice of system with respect to
which to describe the dynamics of the motion because in $S_A$ the
inertia tensor is constant, i.e. doesn't change with time, while
in $S$ it does, making things complicated.
Thus, the idea is to study the (discrete) time evolution (i.e.
evolution in small time-increments) of the quantities of motion in
$S_A$ and then transform them into the system $S$. If $\vec{X}$ is
the position vector of a point on the asteroid with respect to
system $S_A$, the position vector $\vec{x}(n)$ at time $n$ of the
same point with respect to $S$ can be calculate by applying the
following transformation
$$\vec{x}(n) = g(n) \,\vec{X} + {a}(n)$$ where $n=0, 1, 2, 3, ...$
is the discrete time, $g(n) \in \text{SO}(3)$ is a time-dependent
orthogonal matrix, representing rotation at time $n$, and
$\vec{a}(n)$ is the position of the center of mass $A$ of the
asteroid at time $n$. Furthermore, let $\vec{v}$ be a fixed free
vector with respect to $S$ (free vector is not a position vector).
Then the same free vector, observed in $S_A$ will change with time
$n$, so in $S_A$ vector $\vec{v}$ will be $\vec{V}(n)$ at time $n$
and the two are connected by
$$\vec{v} = g(n) \, \vec{V}(n)$$
Next, after one time step (i.e. one time increment) time changes
from $n$ to $n+1$ so we get
$$\vec{v} = g(n+1)\, \vec{V}(n+1)$$ Thus,
$$\vec{V}(n+1) = g(n+1)^{-1} \, \vec{v} = g(n+1)^{-1} \,g(n) \, \vec{V}(n)$$
Denote by $W(n) = g(n)^{-1}\, g(n+1)\, \in \, \text{SO}(3)$. Then
$$\vec{V}(n+1) = W(n)^{-1} \, \vec{V}(n) = W(n)^T \, \vec{V}(n)$$
because $W(n)^{-1} = W(n)^T$ for any orthogonal matrix. The
orthogonal operator $W(n)^T$ describes the one time step rotation
of the vector $\vec{V}(n)$ to the vector $\vec{V}(n+1)$ in the
coordinate system $S_A$. The matrix $W(n)$ is the discrete time
analogue of angular velocity. If one finds the matrices $W(n)$
then one can generate the matrices $g(n)$ by
$$g(n+1) = g(n) \, W(n)$$
If $h>0$ is the small time step (the time increment), then the
actual time is $t_{n} = n\, h$ and $t_{n+1} = t_{n} + h = (n+1)h$.
Then, during the time period $h$ between $t_n$ and $t_{n+1}$ the
rotation $W(n)$ is around an axis determined by a vector
$\vec{\Omega}(n)$ in $S_A$ for which $\Omega(n) =
\text{Matrix}(\vec{\Omega}(n))$ and so
\begin{align} W(n) = e^{h\, \Omega(n)} = {\bf 1} \, + \,
\frac{\sin\Big(h\,
\big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|}\,\, \Omega(n) \, +
\, \frac{1 - \cos\Big(h\,
\big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|^2}\,\, \Omega(n)^2\end{align}
The vector $\vec{\Omega}(n)$ is the angular velocity written in
the coordinate system $S_A$ and ${\bf 1}$ is the identity matrix.
With analogy of the time continuous case of Euler's differential
equations of motion, the latter are based on two laws of physics.
(i) The angular momentum $\vec{m}$ of the asteroid with
respect to $S$ is constant, i.e. it doesn't change with time. This
is the so called conservation of angular momentum law. Thus,
$$\vec{m} = g(n) \, \vec{M}(n)$$ where $\vec{M}(n)$ is the same
angular momentum but written with respect to the moving coordinate
system $S_A$, so $\vec{M}(n)$ must change with time $n$. But for
$n+1$
$$\vec{m} = g(n+1) \, \vec{M}(n+1)$$ Thus
$$g(n) \, \vec{M}(n) = \vec{m} = g(n+1) \, \vec{M}(n+1)$$ which
means that
$$\vec{M}(n+1) = g(n+1)^{-1}\, g(n) \, \vec{M}(n)$$ and hence
$$\vec{M}(n+1) = W(n)^T \, \vec{M}(n) = \Big( e^{-\,h\, \Omega(n)} \Big)\, \vec{M}(n)$$
where $g(n)^{-1}\, g(n+1) = W(n) = e^{h\, \Omega(n)}$ and $W(n)^T
= e^{-\, h\, \Omega(n)}$.
(ii) The inertia matrix $\mathcal{I}$ gives the link between
the angular velocity and the angular momentum
$$\vec{M}(n) = \mathcal{I} \, \vec{\Omega}(n)$$ so the discrete
time analogue of Euler's differential equations are the following
discrete time difference equations
$$\vec{M}(n+1) = \left( e^{-\,h\, \text{Matrix}\big(\mathcal{I}^{-1} \, \vec{M}(n)\big)} \right)\, \vec{M}(n) $$
$$M(n) = \text{Matrix}(\vec{M}(n)) \,\,\, \text{ and } \,\,\,
\vec{M}(n) = \text{Vector}({M}(n)) $$
Generating $\vec{M}(n)$ leads to $W(n) = e^{h\, \mathcal{I}^{-1}
\, {M}(n)}$ and consequently to
$$g(n) = g_0\, W(1)\, W(2) \, ... \, W(n)$$
If $Q$ is a point on the asteroid with coordinates (i.e. position
vector) $$\overrightarrow{Q} = \overrightarrow{AQ} = \begin{pmatrix}Q_1 \\
Q_2\\
Q_3 \end{pmatrix}$$
in the coordinate system $S_A$, the point $Q$ has coordinates
(i.e. position vector) that change with time
$$\vec{q}(n) = \overrightarrow{OQ}(n) = \begin{pmatrix} q_1(n) \\
q_2(n)\\
q_3(n) \end{pmatrix}$$
with respect to the
coordinate system $S$
and
$$\vec{q}(n) = g(n) \, \overrightarrow{Q} + \vec{a}(n) =
\Big( g_0\, W(1)\, W(2) \, ... \, W(n) \Big) \, \overrightarrow{Q}
+ \vec{a}(n)$$
Algorithm
1. Start with initial position of asteroid's center of
gravity $A$ given by a position vector $\vec{a}(0)$ with respect
to the inertial coordinate system $S = O\, \vec{e}_1\, \vec{e}_2\,
\vec{e}_3$.
2. Furthermore, fix an initial orientation of the asteroid
with respect to $S$, given by a choice of three vectors forming
the coordinate system $S_A = A\, \vec{E}_1\, \vec{E}_2\,
\vec{E}_3$ attached to the asteroid at time $n=0$. The latter
three vectors are written as $3$D vector columns with respect to
$S$.
3. Form the orthogonal matrix $g(0) = \Big(\vec{E}_1\,
\vec{E}_2\, \vec{E}_3\Big)^T \, \in \, \text{SO}(3)$.
4. Choose a point $Q$ on the asteroid, which is then given
by a position vector $\overrightarrow{Q}$ with respect to $S_A$.
5. Take a diagonal inertia matrix $\mathcal{I} =
\text{diag}(I_1, I_2, I_3)$.
6. Fix a time step $h>0$, say $h = \frac{1}{25}$ ($25$
frames per second is a reasonable choice to simulate smooth
continuous motion on the screen, but the rate of change of frames
can be chosen even faster if you want). Time is then $t = n \, h$
seconds, where $n = 0, 1, 2, ...$. Choose an integer $N$ so that
time runs from $t=0$ to $t = N\, h$.
7. Give the asteroid an initial angular momentum
$\vec{M}(0)$ written as a vector column with respect to $S_A$.
Alternatively, you can choose initial angular velocity
$\vec{\Omega}(0) = \mathcal{I}^{-1} \, \vec{M}(0)$.
8. Execute the following loop: For $n = 0\, ... \, N$ do
the following
\begin{align*}
&\vec{\Omega}(n) = \mathcal{I}^{-1} \, \vec{M}(n) \\
&{\Omega}(n) = \text{Matrix}\Big( \vec{\Omega}(n) \Big) \\
&W(n) = e^{h\, {\Omega}(n)} = {\bf 1} \, + \, \frac{\sin\Big(h\,
\big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|}\,\, \Omega(n) \, +
\, \frac{1 - \cos\Big(h\,
\big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|^2}\,\Omega(n)^2\\
&\vec{M}(n+1) = W(n)^T\, \vec{M}(n)\\
&g(n+1) = g(n) \, W(n)
\end{align*}
The result is a sequence of matrices $g(n) \, : \, n = 0\, ... \,
N$ and so the motion of point $Q$ with respect to the inertial
non-moving frame $S$ is
$$\vec{q}(n) = g(n) \, \overrightarrow{Q} + \vec{a}(n) =
g(n) \, \overrightarrow{Q} + \vec{a}(0) + nh\, \vec{v}_A $$ where
I have assumed that the center of mass $A$ of the asteroid moves
uniformly with constant velocity $\vec{v}_A$ along a straight-line
trajectory.
I will be very interested to get some feed back from you whether
you have been able to realize the simulation of this kind of
motion, regardless of the involvement of shooting projectiles at
point $Q$ from the asteroid. So. please let me know whether it
works.
