Attitude matrix and quaternion path for m>1

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,


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

The mysterious paths on the three-sphere

Wherever you are, whatever you do, there is a certain special direction that takes you out of the infinite labyrinth, and leads to the reincarnation cycle getting closer and closer to the ideal path. Ordinary people do not know about it. Warriors do know.

Beyond a certain point there is no return. This point has to be reached.
Franz Kafka

I am discussing geodesic lines on the three-sphere, geodesics of the left-invariant metric determined by asymmetric rigid body. There is a special, very special class of geodesics there. At every point there is a special direction. If you start the geodesic in this special direction – it has, at both ends, in the future and in the past – a limit cycle/circle. Enough esoteric talk. Let’s go to the math. The following math describes one such geodesic line. Other are obtained by left translations.

I_1,I_2,I_3 are moments of inertia of our rigid body, ordered as I_1<I_2<I_3.
We define

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1(I_3-I_2)}{I_2(I_3-I_1)}},\\ A_2&=&1,\\ A_3&=&\sqrt{\frac{I_3(I_2-I_1)}{I_2(I_3-I_1)}},\\ B&=&\frac{1}{I_2}\,\sqrt{\frac{(I_2-I_1)(I_3-I_2)}{I_1I_3}},\\ \delta&=&\frac{\sqrt{I_2(I_3-I_1)}-\sqrt{I_1(I_3-I_2)}}{\sqrt{I_3(I_2-I_1)}} \end{eqnarray*}


(2)   \begin{eqnarray*} l_1(t)&=& A_1\,\mathrm{sech }(B t),\\ l_2(t)&=&A_2\,\tanh( B t),\\ l_3(t)&=&A_3\,\mathrm{sech }(B t), \end{eqnarray*}

where \mathrm{sech }(x)=1/\cosh(x).

(3)   \begin{equation*} \psi(t)=\frac{t}{I_2}+2\arctan\left(\delta \tanh(B t/2)\right). \end{equation*}


(4)   \begin{eqnarray*} q_0(t)&=&\frac{\sqrt{1+l_1(t)}\cos\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_1(t)&=&\frac{\sqrt{1+l_1(t)}\sin\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_2(t)&=&\frac{l_3(t)\cos\frac{\psi(t)}{2}+l_2(t)\sin\frac{\psi(t)}{2}}{\sqrt{2(1+l_1(t))}},\\ q_3(t)&=&\frac{-l_2(t)\cos\frac{\psi(t)}{2}+l_3(t)\sin\frac{\psi(t)}{2}}{\sqrt{2(1+l_1(t))}}. \end{eqnarray*}

The story is this. The body can rotate uniformly about its middle axis forever. Either left or right. These are two circles in the rotation group that, using its double cover, topologically is S^3. Stereographic projection maps circles into circles. So they become two circles in \mathbf{R}^3. These circles:

But there is another possibility, when this uniform rotation along the middle axis is only asymptotic. It happens in the infinite past and in the infinite future. But only approximately with the real life, and excluding the short metamorphosis period. This is the trajectory described by the formulas above. Here is the trajectory for time t from t=-20000 to t=20000:

And here is this trajectory together with asymptotic circles. Here the red circle is for the bad past, the blue circle, for the good future. Under microscope it reveals rich structure – infinite mystery.

Approaching the ultimate answer for m<1

Far out in the uncharted backwaters of the unfashionable end of the Western Spiral arm of the Galaxy lies a small unregarded yellow sun.
Orbiting this at a distance of roughly ninety-eight million miles is an utterly insignificant little blue-green planet whose ape-descended life forms are so amazingly primitive that they still think digital watches are a pretty neat idea.
This planet has—or rather had—a problem, which was this: most of the people living on it were unhappy for pretty much of the time. Many solutions were suggested for this problem, but most of these were largely concerned with the movements of small green pieces of paper, which is odd because on the whole it wasn’t the small green pieces of paper that were unhappy.
And so the problem remained; lots of the people were mean, and most of them were miserable, even the ones with digital watches.

So sets Douglas Adams the scenery in his Hitchhiker’s Guide to the Galaxy. And he is right, neither green pieces of paper nor digital watches will make us happy. But digital computers – that’s another story. Digital computer, with some patience and determination made me happy today (or I made him happy, doesn’t mater which way), and I am going now to share my happiness with you. Because there is no true happiness which is not shared happiness.

And no, I am not going to calculate the answer to the Ultimate Question of Life, the Universe and Everything. Though we will be not very far away from obtaining such an answer. But please, don’t panic!

In the last post I have perversely asked this question:

So, my guess is that the formula (23) in the paper has wrong signs. What to do?
Should I write to the authors?

Now I admit, I was not completely sincere. I was teasing, I was playing with You, the Reader. Because the truth it that months ago I discovered wrong signs in the formula (20) of this paper. So, being young and naive, I wrote to prof. Zanna announcing my discovery. Here is the reply I got:

So, we should not consider the paper as being sacred, typos sometimes happen. But software – that is something different. Software is a serious matter. The software has been thoroughly tested, so if we test our software against something that has been thoroughly tested, that will be even better than through tests of our own code!

In The Hitchhiker’s Guide to the Galaxy the answer came as forty two.

“I checked it very thoroughly,” said the computer, “and that quite definitely is the answer.

We will also check very thoroughly. Our answer will not exactly be forty two. But we will come out with a quaternion, quaternio, quatro – it can explain 4. It remains to explain the last 2.

We will be testing example four from the code 903 – you need to read a couple of previous posts if you do not know what I am talking about. Alternatively, you can download the file, unzip it, go to the folder 903\Fortran77\Drivers, and contemplate the following three files there: driver_example4.f, init_example4.dat, out_example4.dat. All three with number four. That is what we will be doing now.

 Helsinky OH 52 - Emil ZÁTOPEK maraton
Emil Zátopek was a Czechoslovak long-distance runner best known for winning three gold medals at the 1952 Summer Olympics in Helsinki. He won gold in the 5,000 metres and 10,000 metres runs, but his final medal came when he decided at the last minute to compete in the first marathon of his life.

According to README.txt

driver_example3.f (example for FRB integration with many steps
of lenght h from 0 to Tfin, using
semi-exact methods)
To be linked with frb_step.f
driver_example4.f (Like driver_example3.f but uses quaternions)
To be linked with quat_step.f

Both examples have the body with the same moments of inertia, but initial momenta are different. So it is impossible to compare their outputs. We return to this issue later on, for now let us concentrate on example4.

The moments of inertia can be read from init_example4.dat. They are

(1)   \begin{eqnarray*} I_1&=&1.0,\\ I_2&=&1.012686988782515,\\ I_3&=&3.306237422473038 \end{eqnarray*}

So it is the same “future physics” body that we have already met in , when discussing example2.

Then comes initial angular momentum (which sometimes I call \mathbf{m} and sometimes \methbf{L}:

(2)   \begin{eqnarray*} m_1(0)&=&-0.544332842491675,\\ m_2(0)&=&0.729131780907662,\\ m_3(0)&=&-0.414811526666455. \end{eqnarray*}

As in Spinning Gulliver among giants we have

(3)   \begin{equation*} 1/I_2= 0.987472.\end{equation*}

We calculate d

    \[d =\frac{m_1(0)^2}{I_1}+\frac{m_2(0)^2}{I_2}+\frac{m_3(0)^2}{I_3}=0.873315.\]

So d<1/I_2, therefore we have m<1 case.
Now I am presenting the complete m<1 algorithm based on the method introduced in Standing on the shoulders of giants – Reboot.

The case m<1.

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

The solution \mathbf{L}(t) of the Euler’s equations has two trajectories: one with L_3(t)>0 is given by

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

while the other one, with L_3(t)<0 is given by

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

With constants \alpha and \nu defined as

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

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

we set the phase variable \psi(t) as

(9)   \begin{equation*} \psi(t)=\frac{t}I_1-\arctan\left((A_2/A_3)\mathrm{sd}(Bt,m)  \right)-\alpha \Pi(\nu,\am(Bt,m),m),\ \end{equation*}

where the Jacobi function \mathrm{sd} is defined as \mathrm{sd}(u,m)=\sn(u,m)/\dn(u,m).
Then the quaternionic attitude solution is given by q(t)=W(t)+\mathbf{i}X(t)+\mathbf{j}Y(t)+\mathbf{k}Z(t), with

(10)   \begin{eqnarray*} W(t)&=&\cos \frac{\psi(t)}{2}\sqrt{\frac{1+L_1(t)}{2}},\\ X(t)&=&\sin \frac{\psi(t)}{2}\sqrt{\frac{1+L_1(t)}{2}},\\ Y(t)&=&\frac{L_3(t)\cos \frac{\psi(t)}{2}+L_2(t)\sin \frac{\psi(t)}{2}}{\sqrt{2(1+L_1(t))}},\\ Z(t)&=&\frac{-L_2(t)\cos \frac{\psi(t)}{2}+L_3(t)\sin \frac{\psi(t)}{2}}{\sqrt{2(1+L_1(t))}} \end{eqnarray*}

My original intention was to verify the above with the 903 paper code using Mathematica, and, if there will be differences in the output, to fix the problems.
But now this Sunday is coming to the end, and I am getting tired. And before sleep I want to watch the second half of the movie “Master and Commander”. It fits this post:

Roger Ebert gave the film 4 stars out of 4, saying that “it achieves the epic without losing sight of the human”