Attitude matrix and quaternion path for m>1

Posted on by arkajad

I noticed that somehow I did not finish with the case m>1. So today, without further ado, I am posting the algorithm.
We have the body with I_1<I_2<I_3, and we are considering the case with d>1/I_2, where d is, as always, the ratio of the doubled kinetic energy to angular momentum squared.
Then we define

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1(dI_3-1)}{I_3-I_1}},\\ A_2&=&\sqrt{\frac{I_2 (1-dI_1)}{I_2-I_1}},\\ A_3&=&\sqrt{\frac{I_3 (1-dI_1)}{I_3-I_1}},\\ B&=&\sqrt{\frac{(dI_3-1)(I_2-I_1)}{I_1I_2I_3}},\\ \mu&=&\frac{1}{m}=\frac{(1-dI_1)(I_3-I_2)}{(dI_3-1)(I_2-I_1)}. \end{eqnarray*}

(2)   \begin{eqnarray*} L_1(t)&=&A_1\,\dn(Bt,\mu),\\ L_2(t)&=&A_2\,\sn(Bt,\mu),\\ L_3(t)&=&A_3\,\cn(Bt,\mu). \end{eqnarray*}

(3)   \begin{equation*} \alpha=\frac{I_3-I_1}{\sqrt{\frac{I_1(dI_3-1)(I_2-I_1)I_3}{I_2}}}, \end{equation*}

(4)   \begin{equation*} \nu=\frac{I_3-dI_1I_3}{I_1-dI_1I_3}. \end{equation*}

(5)   \begin{equation*} \psi(t)=\frac{t}I_3-\arctan\left((A_2/A_1)\mathrm{sd}(Bt,\mu) \right)+\alpha \Pi(\nu,\am(Bt,\mu),\mu),\ \end{equation*}

(6)   \begin{equation*} Q_1(t)=\begin{bmatrix}1-\frac{L_1(t)^2}{1+L_3(t)}&-\frac{L_1(t)L_2(t)}{1+L_3(t)}&-L_1(t)\\ -\frac{L_1(t)L_2(t)}{1+L_3(t)}&1-\frac{L_2(t)^2}{1+L_3(t)}&-L_2(t)\\L_1(t)&L_2(t)&L_3(t)\end{bmatrix}. \end{equation*}

(7)   \begin{equation*} Q_2(t)=\begin{bmatrix}\cos\psi(t)&-\sin\psi(t)&0\\\sin\psi(t)&\cos\psi(t)&0\\0&0&1 \end{bmatrix}. \end{equation*}

(8)   \begin{equation*} Q(t)=Q_2(t)Q_1(t). \end{equation*}

(9)   \begin{eqnarray*} q_0(t)&=&\sqrt{\frac{1+L_3(t)}{2}}\cos\frac{\psi(t)}{2},\\ q_1(t)&=&\frac{1}{\sqrt{2(1+L_3(t))}}\left(L_2(t)\cos\frac{\psi(t)}{2}+L_1(t)\sin\frac{\psi(t)}{2}\right),\\ q_2(t)&=&\frac{1}{\sqrt{2(1+L_3(t))}}\left(L_2(t)\sin\frac{\psi(t)}{2}-L_1(t)\cos\frac{\psi(t)}{2}\right),\\ q_3(t)&=&\sqrt{\frac{1+L_3(t)}{2}}\,\sin\frac{\psi(t)}{2},\\ q(t)&=&(q_0(t),q_1(t),q_2(t),q_3(t))=q_0(t)+\mathbf{i}\,q_1(t)+\mathbf{j}\,q_1(t)+\mathbf{k}\,q_3(t). \end{eqnarray*}

I use the above formulas to draw a stereographic projection of one particular path. So, I take I_1=1,I_2=2,I_3=3,d=0.5000001, and do the parametric plot of the curve \mathbf{r}(t) in \mathbf{R}^3,
with

    \[\mathbf{r}(t)=\left(\frac{q_1(t)}{1-q_0(t)},\frac{q_2(t)}{1-q_0(t)},\frac{q_3(t)}{1-q_0(t)}\right).\]

I show below two plots. One with t\in(-1000,1000), and one with t\in(-10000,10000). For this selected value of d, the time between consecutive flips, given by the formula

(10)   \begin{equation*}\tau =4\mathrm{EllipticK}(\mu)/B= 116.472.\end{equation*}

So, for t\in(-10000,10000) we have somewhat less than 200 flips, and the lines are getting rather densely packed in certain regions.
Notice that is just one geodesic line, geometrically speaking the straightest possible line in the geometry determined by the inertial properties of the body.

It is this “geometry”that will become the main subject of the future notes.

Geodesic line for t between -1000 and 1000

Geodesic line for t between -10000 and 10000

How can such a line be “straight”? Well, it is….

20 thoughts on “Attitude matrix and quaternion path for m>1”

  1. I noticed, on this 10000 picture, that there appears something similar to Saturn’s rings. Don’t know why. Perhaps it is a numerical artefact, but perhaps it is a real thing! A DISCOVERY?

      1. What i was plotting is just one trajectory. I mean in theory, because taking into account numerical approximations who knows what I was plotting in reality?
        In the meantime I have made the plot increasing time from 10 000 to 100 000. It came out completely dense without any Saturn’s rings structure. But numerically it is even less reliable.
        Yes, some trajectories are exceptional, but how easy it is to spot such is not an easy question to answer – except of the three uniform rotations about inertia axes.

          1. I did some tests. With no clear conclusions. For instance for d=0.5000002 there are groups of 2 instead of 5 as in the pisture in the post.

    2. Here is my Maple code. Result agree. Perhaps you have forgotten to take square root of mu?
      I1 := 1.0; I2 := 2.0; I3 := 3.0; d := .5000001; A1 := sqrt(I1*(I3*d-1)/(I3-I1)); A2 := sqrt(I2*(-I1*d+1)/(I2-I1)); A3 := sqrt(I3*(-I1*d+1)/(I3-I1)); B := sqrt((I3*d-1)*(I2-I1)/(I1*I2*I3)); mu := (-I1*d+1)*(I3-I2)/((I3*d-1)*(I2-I1)); kappa := sqrt(mu); nu := (-I1*I3*d+I3)/(-I1*I3*d+I1); alpha := (I3-I1)/sqrt(I1*(I3*d-1)*(I2-I1)*I3/I2); epi := proc (t::float, nu::float, k::float) local t2, n, dt, ep0, res; ep0 := EllipticPi(nu, k); t2 := EllipticK(k); n := floor(t/t2); dt := t-n*t2; if type(n, even) then res := Re(n*ep0+EllipticPi(JacobiSN(dt, k), nu, k)) else res := Re((n+1)*ep0-EllipticPi(JacobiSN(t2-dt, k), nu, k)) end if end proc;
      1.0
      2.0
      3.0
      0.5000001
      0.5000001500
      0.9999999000
      0.8660253172
      0.2886752212
      0.9999992000
      0.9999996000
      -2.999997600
      2.309400384
      proc(t::float, nu::float, k::float) … end;

      psi(t):=t/(I3)-arctan((A2/(A1))*JacobiSD(B*t,kappa))+alpha*epi(B*t,nu,kappa);
      t /A2 JacobiSD(B t, kappa)\
      t -> — – arctan|———————–|
      I3 \ A1 /

      + alpha epi(B t, nu, kappa)
      psi(7);
      3.095170587
      psi(179);
      98.01928537

  4. I get strange kicks in the path. http://imgur.com/RIoHWha . My quaternion for t =179 evaluates to
    Vector[column](4, [0.8808818313e-1, .1154559730, -.8474014431, .5107065531])
    Cheching the quaternion formula, it is the negative of what you have posted in the formulas.

  5. I have to improve my communication level.
    Lets take d=0.6 t=179
    Your quaternion formula (9) evaluates to
    Vector[column](4,
    [-.6452072055, .2607319986, .5511334836, -.4604110882])
    Stereographicially projected above evaluates to
    Vector[column](3, [.1584797329, .3349933562, -.2798499098])
    The quaternion I arrived at evaluates to
    Vector[column](4,
    [.6452072054, -.2607319985, -.5511334835, .4604110882])
    Stereographicially projected evaluates to
    Vector[column](3, [-.7348852696, -1.553395367, 1.297690075])
    Appears very different in 3D.Interestingly similar but different paths are produced.
    See the plots in link
    There are twice as many lines in the combined plot as the paths don’t match.
    http://imgur.com/a/Vhy9q

    1. Looking at your file I am finding:
      My derivation.
      ` This is from Algorithm 903 paper and post:- Giants that Stumble`;

      q1(t):= 1/(sqrt(2)) Vector(4,[-sqrt(1+L3(t)),-(L2(t))/(sqrt(1+ L3(t))),(L1(t))/(sqrt(1+L3(t))) ,0]);

      Where did you get this formula from? It is the minus of what is in the paper, and what is in my formula.

  7. Now I remember/know. In giants that stumble and the paper I didn’t understand where the quaternion came from. Figured out how to derive the individual ones to rotate to x,y & z axis. You had mentioned a problem with incorrect “-” signs causing a mess. I realised a negative quaternion can also be used to achieve the rotation.http://arkadiusz-jadczyk.eu/blog/2017/03/giants-stumble/#comment-671. I guess I never changed the signs back then copy pasted to recent worksheet.
    Lesson here just because the negative unit quaternion produces the same rotation, it does not produce the same sterographic projection.

    1. Yes. I have chosen stereographic projection from the pole +1. I could have chosen from -1. Then minus quaternion projected from +1 is negative of + quaternion projected from -1.

Leave a Reply