Attitude matrix and quaternion path for m>1

Posted on 2017/03/22 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….

Print 🖨 PDF 📄 eBook 📱
This entry was posted in Attitude matrix, Elliptic functions, Moments of inertia, Quaternions, Stereographic projection. Bookmark the permalink.

20 Responses to Attitude matrix and quaternion path for m>1

  1. arkajad says:
    2017/03/22 at 21:08

    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?

    • Bjab says:
      2017/03/23 at 09:40

      It is not surprising – for some d’s the path should be closed.

    • Bjab says:
      2017/03/23 at 10:08

      I mean closed and relatively short.

      • arkajad says:
        2017/03/23 at 11:27

        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.

  2. Ronan says:
    2017/03/26 at 21:37

    Can you post a couple of values for psi(t) as a check please. I get
    psi(7)=5.010163143
    psi(179)=163.1996033

    • arkajad says:
      2017/03/26 at 22:41

      For d=0.5000001 Mathematica gives me
      psi(7)=3.09517
      psi(179=98.0193

      I will try now to reproduce these values with Maple

    • arkajad says:
      2017/03/26 at 23:13

      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

  3. Ronan says:
    2017/03/27 at 01:23

    Thank you. My problem was in alpha. Things agrees now. Oh I ommited to say d=0.6

  4. Ronan says:
    2017/03/27 at 17:14

    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. Ronan says:
    2017/03/28 at 01:05

    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

    • arkajad says:
      2017/03/29 at 12:00

      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. Ronan says:
    2017/03/29 at 23:15

    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.

    • arkajad says:
      2017/03/30 at 08:26

      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