Day: March 6, 2017

The final answer for the Universe in which m<1

Posted on by arkajad

This is not an independent post or note or whatever. It is a continuation. In fact it is the continuation of Approaching the ultimate answer for m<1. We have approached enough. We are ready.

There was a moment’s expectant pause while panels slowly came to life on the front of the console. Lights flashed on and off experimentally and settled down into a businesslike pattern. A soft low hum came from the communication channel.
“Good morning,” said Deep Thought at last.
“Er … good morning, O Deep Thought,” said Loonquawl nervously, “do you have … er, that is …”
“An answer for you?” interrupted Deep Thought majestically. “Yes. I have.”
The two men shivered with expectancy. Their waiting had not been in vain.
“There really is one?” breathed Phouchg.
“There really is one,” confirmed Deep Thought.

Douglas Adams, The Hitchhiker’s Guide to the Galaxy

We did not get to the answer yesterday, but today is the day. We will be using the formulas from the last post. Here they are again, repeated for your convenience:

The case m<1.

(1)   \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

(2)   \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

(3)   \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

(4)   \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*}

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

we set the phase variable \psi(t) as

(6)   \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

(7)   \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*}

We consider Example Four with

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

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

m_3(0) is negative, the therefore we take option as in Eq. 3). We calculate t_0 that reproduces the initial angular momenta. The method is known from the previous posts:

(10)   \begin{equation*}am=\arctan(-m_1(0)/A_1,m_2(0)/A_2)=0.925158,\end{equation*}

(11)   \begin{equation*}t_0=\mathrm{EllipticF}(am,m)=3.17243.\end{equation*}

So we have t_0. At t=t_0 our angular momentum is the same as angular momentum in example 4 at t=0. We calculate q_0=q(t_0) from Eq'(7)

    \[q_0=W(t_0)+X(t_0)\,\mathbf{i}+Y(t_0)\,\mathbf{j}+Z(t_0)\,\mathbf{k},\]

(12)   \begin{eqnarray*} W(t_0)&=&0.435964,\\ X(t_0)&=&0.194343,\\ Y(t_0)&=&-0.0858987,\\ Z(t_0)&=-&0.874521 \end{eqnarray*}

thus

(13)   \begin{equation*}q_0=0.435964+0.194343\,\mathbf{i}-0.0858987\,\mathbf{j}-0.874521\,\mathbf{k}.\end{equation*}

We verify that we have indeed a unit quaternion:

    \[||q_0||^2=W(t_0)^2+X(t_0)^2+Y(t_0)^2+Z(t_0)^2=1.0.\]

Thus the inverse quaternion is the same as conjugated one:

    \[q_0^{-1}=0.435964-0.194343\,\mathbf{i}+0.0858987\,\mathbf{j}+0.874521\,\mathbf{k}.\]

In the Example 4 the initial quaternion, at t=0 is 1. Therefore in order for our solution to reproduce the initial data of the example we set

(14)   \begin{equation*}\tilde{q}(t)=q_0^{-1}\,q(t+t_0).\end{equation*}

The final time in Example 4 is t=10. Therefore we calculate

(15)   \begin{equation*}\tilde{q}(10)=q_0^{-1}q(10+t_0),\end{equation*}

where the multiplication is the quaternion multiplication.
Calculated with Mathematica the answer is

(16)   \begin{eqnarray*} W&=&=-0.36762,\\ X&=&-0.630629,\\ Y&=&-0.612723,\\ Z&=&0.302874. \end{eqnarray*}

The result from the Fortran code in the file out_example4.dat is

-0.3676198430772359
-0.6306293413288832
-0.6127232632258010
0.3028737154869889


It seems that therefore that we have obtained the ultimate answer. At least for quaternions. We need to get it also for rotation matrices. So the saga will continue.