User friendly means “easy to use”. We can all agree on that, I think. And it’s what we mean when we say it. So far, so good.

Except, not. When people say “easy to use”, what they REALLY mean is “easy to learn”, or “easy to use without learning”.

From What User friendly means

I am trying to make my blog user friendly. Sometimes I think I succeed. But not always. We are following giants. It is not easy. Sometimes we need to overcome our fears.

I am studying the paper and the code:

Elena Celledoni, Antonella Zanna, et al. “

Algorithm 903: FRB–Fortran routines for the exact computation of free rigid body motions.” ACM Transactions on Mathematical Software, Volume 37 Issue 2, April 2010, Article No. 23, doi:10.1145/1731022.1731033

I am trying to understand it and to reproduce the results. If I am able to reproduce the results, then there is a good chance that I really understand what I am studying, and then I can share my understanding here, on the blog. Not all I am able to understand, but what to do? **Sometimes I am completely lost**.

Here is one such case. I am quoting one paragraph from Standing on the shoulders of giants – Reboot

I am going to use the following paper:

Elena Celledoni, Antonella Zanna, et al.

“Algorithm 903: FRB–Fortran routines for the exact computation of free rigid body motions.”ACM Transactions on Mathematical Software, Volume 37 Issue 2, April 2010, Article No. 23, doi:10.1145/1731022.1731033.The paper is only seven years old. And it has all what we need. As in Standing on the shoulders of giants we will split the attitude matrix into a product: but this time we will do it somewhat smarter. Celledoni and Zanna denote the angular momentum vector in the body frame using letter and they assume that is normalized: They assume, as we do, and use to denote the kinetic energy. Let us first look at the following part of their paper (which is available here):

In the same post we introduced the matrix that transforms vector with components into In Attitudes and Behavior we wrote it as (and I am replacing with etc.):

(1)

We have the formula unit mapping quaternions into rotation matrices (see Chiromancy in the rotation group).

If

(2)

(3)

If we use the formula (21) from the paper (see image above), then

(4)

Now I am using Mathematica to calculate

We obtain Two details will catch the attention of the Reader:

In the Mathematica code I use, for instance, while in the formula for I have But for the unit quaternion, with , these are the same! Similarly the Mathematica code gives as the component of the matrix

while in we have simply . But again, with the Reader will easily check that these are the same!

So it seems that we are progressing with our understanding of the language and the science of giants! Success. Therefore we must not stop! We must continue!

After the case Celledoni and Zanna move to the case of :

So, if we understand it correctly, for we should replace matrix with the new one, constructed the same way as before, but with the third axis being replaced by the first. The quaternion that implements such a rotation, according to the formula (23) from the paper, should be given by

(5)

I calculate for this quaternion, and I get …. mess. Certainly not the matrix that acting on would give So, **something is wrong**! I am trying to compare the formulas in the paper with those in the Fortran code. But then I am getting confused even more. For the case the code has these lines:

APPP = 1.0D0 + EM30

APP = SQRT(APPP)

ADUE=SQRT(2.0D0)

Q2(1) = APP/ADUE

Q2(4) = -0.0D0

Q2(2) = -Q2(1)*EM20/APPP

Q2(3) = Q2(1)*EM10/APPP

We recognize that Q2(1) is our in Eqs. (\label{eq{wxyz3}), Q2(4) is our EM10,EM20,EM30 is probably our If so, then Q2(2) is **minus** our and Q2(3) is **minus** our **I have no idea why there are minuses, but there are there**!

What about the case ? The code has these lines:

APPP = 1.0D0 + EM10

APP = SQRT(APPP)

ADUE=SQRT(2.0D0)

Q2(1) = APP/ADUE

Q2(2) = 0.0D0

Q2(3) = -Q2(1)*EM30/APPP

Q2(4) = Q2(1)*EM20/APPP

Again we recognize Q2(1) as our we recognize Q2(2) as But now the signs in Q2(2) and Q2(3) **agree** with those in the paper. Something weird is going on. So, as an experiment I change the signs in formula (23):

(6)

So that, consistently, the two signs in the paper are opposite to those in the code. Now gives what it is supposed to do: rotation matrix that transforms into

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

Should I write to the authors? My thinking goes into the direction: it would be of no use.

The paper is 7 years old. In fact the preprint tells me that the work has been essentially finished in 2008. Now the authors are working on a different project, they have other things to do than worrying about signs in a formula. Perhaps the formula is wrong, so what? Who cares? No one is reading paper anyway. Scientists write papers in order to document that they were busy with something, to show that they are experts, to get new grants, to be quoted. Anyway intelligent reader should be able to fix typos. It should even be advised to make some errors. Fixing someone’s errors takes our time, true. But, on the other hand, what a wonderful opportunity to think and to verify that we really understand the subject! Finding errors in pear reviewed papers gives us reassurance that referees in these papers, as a rule, do not pay attention to the question as to whether a given paper is correct or incorrect. It should simply **give a strong impression that it is correct.**. That is all.

True, some authors take mistakes in their published works seriously. I like to quote this particular very good example: ERRATA TO ANALYSIS, MANIFOLDS AND PHYSICS. PART I, the monograph by Yvonne Choquet-Bruhat and Cecile DeWitt-Morette, is 11 pages long. With all kinds of errors corrected, including correction to wrong proofs. That is life.

So, let us continue. We take what is good. We correct what can be corrected. And we move forward.

I was going to make this post about the case – to give the complete algorithm. But somehow I got deflected, I stumbled. I will recover in the next post.

How is the quaternion p(t) produced.

Given a vector m=(m1(t),m2(t),m3(t)), from Standing on the shoulders of giants -Reboot

cos(theta)=m3(t) and sin(theta)=mp(t) assuming you are rotating the vector up to the z axis. The Rodrigues formula is simple to apply.

But quaternions use 1/2 angles so would need etc. Don’t see how to convert it at the moment to produce the p(t) in the paper.

I figured it out. Half angle formulas.

etc.

It was still a bit tricky though. One mystery removed!

->

I ma ->

I am

->

In eq.(3) R21 right parethesis lost

Thank you. Fixed.

Correct me is my intrepertation is incorrect.

The quaternion is the encoding of the axis and angle of rotation. Changing the sign 0f say (W,X,-Y,0) to (W,-X,Y,0) changes the dirction vector by 180 degrees so the rotation rotates the vector in the opposite direction. If found(-W,-X,Y,0) rotates on opposite direction and works corectly. Have not tested othe x axis yet. Have just derived the formula on paper so far and was unclear which way the direction vector should be pointing. Some of these formulas are easier to simplify be hand than get Maple too.

Yes, for unit quaternions, the conjugated quaternion (with negative imaginary part) determines opposite (inverse) rotation.