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….
So, here is the code:
The solution of the Euler’s equations has two trajectories. For one with we take
(1)
while for the other one, with we take
(2)
Let
(3)
The solution of the Euler’s equations is given by
(4)
With constants and defined as
(5)
(6)
we set the phase variable as
(7)
where the Jacobi function is defined as
Let
(8)
(9)
Then the attitude matrix is given by
(10)
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 rather than with the signs in formulas for 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 that depends on the trajectory we are taking. The constant 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 I use Mathematica to calculate the transposed final
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 is fixed, and equal ! 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 I hope I am right here….