Attitude matrix for m<1

Posted on by arkajad

The last post was about The final answer for the Universe in which m>1. Yet, as I see it today, it was neither complete nor final. It was a very bad approximation. It may be also that our universe is also a very bad approximation to the one that has been intended. But that is another story. I think I will return to this problem in my blogging, but now I will concentrate on the algorithm for the rotation matrix. At the same time I will make a comment how the algorithm from the last post, for the quaternions, should be improved. The fact is that I was making wrong assumptions in my mind. Making wrong assumptions, without being conscious of the fact that one is making such assumptions, that is very dangerous. Lot of bad things happens around us for this reason….

Plane crash
The result of wrong assumptions

So, here is the code:

The case m<1.
The solution \mathbf{L}(t) of the Euler’s equations has two trajectories. For one with L_3(t)>0 we take

(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}},\\ \end{eqnarray*}

while for the other one, with L_3(t)<0 we take

(2)   \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}}.\\ \end{eqnarray*}

Let

(3)   \begin{eqnarray*} 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 is given by

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

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

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

we set the phase variable \psi(t) as

(7)   \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).
Let

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

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

Then the attitude matrix Q(t) is given by

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

If you compare this with the formulas I gave yesterday, you will find the following important difference: now I am playing with signs in the constants A_1,A_3 rather than with the signs in formulas for L_1(t),L_3(t). It took me whole day to realize that this way is a better way. I was not realizing before that there is one sign in the formula for \psi(t) that depends on the trajectory we are taking. The constant A_2/A_3 does this job now.

In a couple of days I will also change the formulas in the quaternionic algorithm to make them universal as well.

Now, I check the modified example3. The modified initial data in init_example3.dat file are as follows:
1.0000000000D0
10.0000000000D0
-0.544332842491675D0
0.729131780907662D0
-0.414811526666455D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
1.000000000000000D0
1.012686988782515D0
3.306237422473038D0
5
They are the same as in the quaternionic case discussed yesterday. The initial attitude matrix is the identity matrix. We get the same t_0. I use Mathematica to calculate the transposed final Q

Transpose[Inverse[Q[t0]].Q[t0 + 10.0]]
{{0.0656754, 0.550118, -0.832501}, {0.995487, 0.0211483,
0.0925082}, {0.0684964, -0.834819, -0.546246}}

We compare it with the Fortran output from 903 code:

0.0656754297
0.55011776
-0.832500564
0.995487311
0.0211482921
0.0925081752
0.0684963551
-0.834819262
-0.546246326

And it seems that our code is doing the same job as the Fortran code from the giants!
Yet, as I said at the beginning, it took me the whole days to figure it out why for some initial data I was in agreement, but for another initial data there was a disagreement.

There is one little plus of the code above. The coefficient in front of \arctan is fixed, and equal 1.! Zanna and Celledoni in their paper and in their code have complicated formulas for this coefficient. It would take some work to verify that their complicated formula simplifies to 1. I hope I am right here….

4 thoughts on “Attitude matrix for m<1”

Leave a Reply