# Attitude matrix and quaternion path for m>1

[latexpage]
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_11/I_2$, where $d$ is, as always, the ratio of the doubled kinetic energy to angular momentum squared.
Then we define

\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}
\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}

\alpha=\frac{I_3-I_1}{\sqrt{\frac{I_1(dI_3-1)(I_2-I_1)I_3}{I_2}}},

\nu=\frac{I_3-dI_1I_3}{I_1-dI_1I_3}.

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

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}.

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

Q(t)=Q_2(t)Q_1(t).

\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
$$\tau =4\mathrm{EllipticK}(\mu)/B= 116.472.$$
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.

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. Bjab says:

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

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. Bjab says:

Have you made pictures (also with time 10 000) but with d slightly greater and d slightly less than already shown 0.5000001 ?

2. Bjab says:

I suspect that small variation in d should thighten or loosen the “Saturn’s rings”.

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. Bjab says:

So the variation (0.0000001) was too big.

2. Ronan says:

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

1. 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

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

3. Ronan says:

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

4. Ronan says:

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.

1. Is that for d=0.6?

“Cheching the quaternion formula, it is the negative of what you have posted in the formulas.”

Which formulas?

5. Ronan says:

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.
There are twice as many lines in the combined plot as the paths don’t match.
http://imgur.com/a/Vhy9q

1. That is, communicationwise, better. But I do not know which way you got your numbers?

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.

6. Ronan says:

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.