The Euler XYZ angles were wrong

[sanguinariojoe]

There were several typos on the Euler angles treatment. Now I think is is quite consistent. I even added a test to demonstrate it.

@AlexWKinley I think you shall try this before continuing your debugging of the centripetal forces on several axes

[AlexWKinley]

I haven't dug into this yet, but my very first thought is that it's surprising that there's a need to change the canonicalEulerAngles function. Because that function isn't of my own doing, and is a function that's actually included in more recent version of Eigen. https://gitlab.com/libeigen/eigen/-/blame/master/Eigen/src/Geometry/EulerAngles.h?ref_type=heads#L45 You can also see the tests for that function here: https://gitlab.com/libeigen/eigen/-/blob/master/test/geo_eulerangles.cpp

I haven't looked enough to tell what the effect of the changes really is. My gut instinct is that if we're going to make changes to that function, there should a least be some comments explaining what the differences are and why.

I'll try and dig into this PR some more later today, because I'm not surprised that there are inconsistencies in how we angle euler angles, and it would be good to sort those out.

[sanguinariojoe]

Uhm, I thought it was part of the unsupported Eigen.

Anyway... Paradoxically, what I have done on that function is just follow what they say they are doing on the comments... I also checked that the comments are right, BTW.

We are now nailing the results obtained with scipy, which really makes me quite confident. Maybe I will try to dig on Eigen itself next week

On Fri, 5 Jul 2024, 15:01 AlexWKinley, @.***> wrote:

I haven't dug into this yet, but my very first thought is that it's surprising that there's a need to change the canonicalEulerAngles function. Because that function isn't of my own doing, and is a function that's actually included in more recent version of Eigen.

https://gitlab.com/libeigen/eigen/-/blame/master/Eigen/src/Geometry/EulerAngles.h?ref_type=heads#L45 You can also see the tests for that function here: https://gitlab.com/libeigen/eigen/-/blob/master/test/geo_eulerangles.cpp

I haven't looked enough to tell what the effect of the changes really is. My gut instinct is that if we're going to make changes to that function, there should a least be some comments explaining what the differences are and why.

I'll try and dig into this PR some more later today, because I'm not surprised that there are inconsistencies in how we angle euler angles, and it would be good to sort those out.

— Reply to this email directly, view it on GitHub https://github.com/FloatingArrayDesign/MoorDyn/pull/233#issuecomment-2210838882, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAMXKKFKUMIDZEVAX4ZLCPTZK2KLBAVCNFSM6AAAAABKK5QU3GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMJQHAZTQOBYGI . You are receiving this because you were assigned.Message ID: @.***>

[AlexWKinley]

Super quick thought: I wonder if part of the difference is intrinsic vs extrinsic euler angles? From scipy docs for from_euler:

Specifies sequence of axes for rotations. Up to 3 characters belonging to the set {‘X’, ‘Y’, ‘Z’} for intrinsic rotations, or {‘x’, ‘y’, ‘z’} for extrinsic rotations. Extrinsic and intrinsic rotations cannot be mixed in one function call.

I'm not actually 100% sure which convention MoorDyn intends to use, but that may be a source of difference, especially if it seems like the solution is flipping the order XYZ -> ZYX.

[sanguinariojoe]

Moordyn definitely wants to go extrinsic. But on the test I am using intrinsic on scipy. That's why I have the abs2releuler function (i do not remember the name. Look at the test, i really think it is pretty self descriptive

I am not with the computer anymore. Do I let you dig on the problem. Good luck and let me know please!

On Fri, 5 Jul 2024, 16:23 AlexWKinley, @.***> wrote:

Super quick thought: I wonder if part of the difference is intrinsic vs extrinsic euler angles? From scipy docs for from_euler:

Specifies sequence of axes for rotations. Up to 3 characters belonging to the set {‘X’, ‘Y’, ‘Z’} for intrinsic rotations, or {‘x’, ‘y’, ‘z’} for extrinsic rotations. Extrinsic and intrinsic rotations cannot be mixed in one function call.

I'm not actually 100% sure which convention MoorDyn intends to use, but that may be a source of difference, especially if it seems like the solution is flipping the order XYZ -> ZYX.

— Reply to this email directly, view it on GitHub https://github.com/FloatingArrayDesign/MoorDyn/pull/233#issuecomment-2210971187, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAMXKKHTB2P6PIKU4UMDCJ3ZK2T7LAVCNFSM6AAAAABKK5QU3GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMJQHE3TCMJYG4 . You are receiving this because you were assigned.Message ID: @.***>

[AlexWKinley]

I've spent a little more time looking at this (and figuring out how to get the python wrapper to compile and install on windows), and I think it's basically equivalent to switching intrinsic to extrinsic euler angles.

If you undo the changes to Euler2Quat and canonicalEulerAngles, and then switch the python test to use intrinsic euler angles seq='xyz'-> seq='XYZ' in both rotate_vecs and abseuler2releuler, the test passes. I've provided a simple patch that does this

instrinsic_angles.patch

Looking at MoorDyn-F, I'm not actually sure which convention they're using https://github.com/OpenFAST/openfast/blob/6a7a543790f3cad4a65b87242a619ac5b34b4c0f/modules/nwtc-library/src/NWTC_Num.f90#L1708

The comment says:

This function creates a rotation matrix, M, from a 1-2-3 rotation
sequence of the 3 Euler angles, \f$\theta_x\f$, \f$\theta_y\f$, and \f$\theta_z\f$, in radians.
M represents a change of basis (from global to local coordinates; 
not a physical rotation of the body). It is the inverse of EulerExtract (nwtc_num::eulerextract).

\begin{align}
M & = R(\theta_z) R(\theta_y) R(\theta_x) \\
  & = \begin{bmatrix}  \cos(\theta_z) & \sin(\theta_z) & 0 \\
                        -\sin(\theta_z) & \cos(\theta_z) & 0 \\
                          0      &  0      & 1 \end{bmatrix}
        \begin{bmatrix}  \cos(\theta_y) & 0 & -\sin(\theta_y) \\
                               0 & 1 & 0        \\
                         \sin(\theta_y) & 0 & \cos(\theta_y)  \end{bmatrix}
        \begin{bmatrix}   1 &  0       & 0       \\
                          0 &  \cos(\theta_x) & \sin(\theta_x) \\
                          0 & -\sin(\theta_x) & \cos(\theta_x) \end{bmatrix} \\
  & = \begin{bmatrix}  
   \cos(\theta_y)\cos(\theta_z) &   \cos(\theta_x)\sin(\theta_z)+\sin(\theta_x)\sin(\theta_y)\cos(\theta_z) &
                                    \sin(\theta_x)\sin(\theta_z)-\cos(\theta_x)\sin(\theta_y)\cos(\theta_z) \\
   -\cos(\theta_y)\sin(\theta_z)  & \cos(\theta_x)\cos(\theta_z)-\sin(\theta_x)\sin(\theta_y)\sin(\theta_z) & 
                                    \sin(\theta_x)\cos(\theta_z)+\cos(\theta_x)\sin(\theta_y)\sin(\theta_z) \\
   \sin(\theta_y)                & -\sin(\theta_x)\cos(\theta_y) & \cos(\theta_x)\cos(\theta_y) \\
        \end{bmatrix}   
\end{align}

I'll admit I don't really understand what the difference is between representing a change of basis and representing the physical rotation of the body. I also don't understand they their single axis rotation matrices are the inverse/transpose of like what is conventionally meant? From wikipedia: https://en.wikipedia.org/wiki/Rotation_matrix#Basic_3D_rotations

As their end result they get a matrix that is the transpose of what is listed on wikipedia as the matrix for X,Y,Z euler angles. https://en.wikipedia.org/wiki/Euler_angles#Conversion_to_other_orientation_representations

I guess I'll have to defer to others for what the intended behavior here is. As well as how we'd want to handle the fact we're potentially breaking other people's existing code if it was relying on how we currently handle rotation. This change fixes the issues I mentioned in the other pr about bodies undergoing multiple axis of rotation simultaneously. But that's not surprising given that I believe that issue is related to the fact that when an object that is spinning around the global z axis rotates 90 degrees around the global x axis, it will then be spinning around the global y axis, because it's local z axis has moved. (Ignoring gyroscopic effects).