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 <#233 (comment)>, 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 <#233 (comment)>, 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).

[RyanDavies19]

@AlexWKinley and @sanguinariojoe I can hopefully offer some clarity on intent here. First off, that comment in the NWTC library needs to be updated. That function gives the intrinsic z, y', x'' rotation matrix, which is confusing based on that comment. In general MoorDyn uses the z-y'-x" intrinsic Tait–Bryan rotation sequence to match the OpenFAST convention (i.e. what is in MoorDynF).

Matt's new theory paper for MoorPy, NREL's quasi-static mooring model, can help shed some light on this (see section 3.4). Both MoorPy and MoorDyn share input file formats, object types and respective coupling, and handling of 6 DOF positioning: https://www.mdpi.com/1996-1073/17/13/3155.

Right now the rotation matrices for MoorDynC, MoorDynF, and MoorPy are all the same (see below). Because Eigen is a well validated package, to maintain consistencies across the codes, and because of Alex's point about breaking users existing code I don't think we need to make changes to how these are handled.

Rotation Matrices
MoorPy (note that in var names and comments, 3 refers to x and 1 refers to z):
https://github.com/NREL/MoorPy/blob/018e976d63f2b5d6c8e1ba26f87d5221590886e3/moorpy/helpers.py#L97

MoorDynC (note that in var names and comments, 1 refers to x and 3 refers to z):

MoorDyn/source/Misc.hpp

Line 230 in 0a9a8aa

// Tait-Bryan angles (all three axes are different; typically used for

MoorDynF:
https://github.com/OpenFAST/openfast/blob/6a7a543790f3cad4a65b87242a619ac5b34b4c0f/modules/nwtc-library/src/NWTC_Num.f90#L1682

[sanguinariojoe]

OK then!

I just left the test and accordingly modified the documentation. I suppose this is then harmless to be merged

[RyanDavies19]

Agreed @sanguinariojoe. I will merge later today unless you have objections @AlexWKinley

[AlexWKinley]

This PR looks good to me, except we should remove that instrinsic_angles.patch file. Once that's fixed, no objections.

[sanguinariojoe]

Fuck my shit! Tomorrow I am fixing that, cherry picking #235, squashing and merging

…

On Tue, 9 Jul 2024, 17:50 AlexWKinley, ***@***.***> wrote: This PR looks good to me, except we should remove that instrinsic_angles.patch file. Once that's fixed, no objections. — Reply to this email directly, view it on GitHub <#233 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMXKKCQFLIFG3G7XQVAPGTZLQBGBAVCNFSM6AAAAABKK5QU3GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMJYGA3DMMRRGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[sanguinariojoe]

I am finally not squashing... I am too lazy