From f75c5243c23fe6c177b2151f8b55e1952b05a4b7 Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Wed, 3 Jul 2024 12:17:07 +0200 Subject: [PATCH 1/8] fix: EulerZYX -> EulerXYZ on moordyn::Euler2Quat() --- source/Misc.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/Misc.hpp b/source/Misc.hpp index 2514ca72..54832b3e 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -262,9 +262,9 @@ inline quaternion Euler2Quat(const vec3& angles) { using AngleAxis = Eigen::AngleAxis; - quaternion q = AngleAxis(angles.x(), vec3::UnitX()) * + quaternion q = AngleAxis(angles.z(), vec3::UnitZ()) * AngleAxis(angles.y(), vec3::UnitY()) * - AngleAxis(angles.z(), vec3::UnitZ()); + AngleAxis(angles.x(), vec3::UnitX()); return q; } From c88755ea34dbe261863d9674a35a487f79e855b3 Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Thu, 4 Jul 2024 06:39:28 +0200 Subject: [PATCH 2/8] fix: accelerations of Coupled/fixed bodies/rods when there are several isntances of them --- source/Body.hpp | 16 +++++++++++++--- source/MoorDyn2.cpp | 44 ++++++++++++-------------------------------- source/MoorDyn2.hpp | 6 ------ source/Rod.hpp | 12 ++++++++++-- 4 files changed, 35 insertions(+), 43 deletions(-) diff --git a/source/Body.hpp b/source/Body.hpp index 26b00c20..456ffcfe 100644 --- a/source/Body.hpp +++ b/source/Body.hpp @@ -248,7 +248,17 @@ class Body final : public io::IO, public SuperCFL * @throw moordyn::invalid_value_error If the body is of type * moordyn::Body::FREE */ - void initializeUnfreeBody(vec6 r = vec6::Zero(), vec6 rd = vec6::Zero(), vec6 rdd = vec6::Zero()); + void initializeUnfreeBody(vec6 r = vec6::Zero(), + vec6 rd = vec6::Zero(), + vec6 rdd = vec6::Zero()); + + /** @brief Get the last setted velocity for an unfree body + * + * For free bodies the behaviour is undetermined + * + * @return The velocity (6 dof) + */ + inline vec6 getUnfreeVel() const { return rd_ves; } /** @brief Initialize the FREE point state * @return The 6-dof position (first) and the 6-dof velocity (second) @@ -349,11 +359,11 @@ class Body final : public io::IO, public SuperCFL * boundary conditions (body kinematics) for the proceeding time steps * @param r The input position * @param rd The input velocity - * @param rdd The input velocity + * @param rdd The input acceleration * @throw moordyn::invalid_value_error If the body is not of type * moordyn::Body::COUPLED or moordyn::Body::FIXED */ - void initiateStep(vec6 r_in, vec6 rd_in, vec6 rdd_in); + void initiateStep(vec6 r, vec6 rd, vec6 rdd); /** @brief Sets the kinematics based on the position and velocity of the * fairlead. diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp index ca14d00d..a6a89a37 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -433,14 +433,10 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) << x[ix + 1] << ", " << x[ix + 2] << "..." << endl; // BUG: These conversions will not be needed in the future - vec6 r, rd, rdd; + vec6 r, rd; if (BodyList[l]->type == Body::COUPLED){ moordyn::array2vec6(x + ix, r); moordyn::array2vec6(xd + ix, rd); - // determine acceleration - rdd = (rd - rd_b) / dtM0; - rd_b = rd; - ix += 6; } else { // for pinned body 3 entries will be taken @@ -449,13 +445,9 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) r(Eigen::seqN(0, 3)) = r3; moordyn::array2vec(xd + ix, rd3); rd(Eigen::seqN(0, 3)) = rd3; - // determine acceleration - rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_b) / dtM0; - rd3_b = rd3; - ix += 3; } - BodyList[l]->initializeUnfreeBody(r, rd, rdd); + BodyList[l]->initializeUnfreeBody(r, rd, vec6::Zero()); } for (auto l : CpldRodIs) { @@ -466,10 +458,6 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) // for cantilevered rods 6 entries will be taken moordyn::array2vec6(x + ix, r); moordyn::array2vec6(xd + ix, rd); - // determine acceleration - rdd = (rd - rd_r) / dtM0; - rd_r = rd; - ix += 6; } else { // for pinned rods 3 entries will be taken @@ -478,13 +466,9 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) r(Eigen::seqN(0, 3)) = r3; moordyn::array2vec(xd + ix, rd3); rd(Eigen::seqN(0, 3)) = rd3; - // determine acceleration - rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_r) / dtM0; - rd3_r = rd3; - ix += 3; } - RodList[l]->initiateStep(r, rd, rdd); + RodList[l]->initiateStep(r, rd, vec6::Zero()); RodList[l]->updateFairlead(0.0); // call this just to set up the output file header RodList[l]->initialize(); @@ -655,13 +639,12 @@ moordyn::MoorDyn::Step(const double* x, for (auto l : CpldBodyIs) { // BUG: These conversions will not be needed in the future vec6 r, rd, rdd; + const vec6 rd_b = BodyList[l]->getUnfreeVel(); if (BodyList[l]->type == Body::COUPLED){ moordyn::array2vec6(x + ix, r); moordyn::array2vec6(xd + ix, rd); // determine acceleration - rdd = (rd - rd_b) / dtM0; - rd_b = rd; - + rdd = (rd - rd_b) / dt; ix += 6; } else { // for pinned body 3 entries will be taken @@ -671,23 +654,21 @@ moordyn::MoorDyn::Step(const double* x, moordyn::array2vec(xd + ix, rd3); rd(Eigen::seqN(0, 3)) = rd3; // determine acceleration - rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_b) / dt; - rd3_b = rd3; - + rdd(Eigen::seqN(0, 3)) = (rd3 - rd_b.head<3>()) / dt; ix += 3; } - BodyList[l]->initiateStep(r, rd, rdd); // acceleration required for inertial terms + // acceleration required for inertial terms + BodyList[l]->initiateStep(r, rd, rdd); } for (auto l : CpldRodIs) { vec6 r, rd, rdd; + const vec6 rd_r = RodList[l]->getUnfreeVel(); if (RodList[l]->type == Rod::COUPLED) { // for cantilevered rods 6 entries will be taken moordyn::array2vec6(x + ix, r); moordyn::array2vec6(xd + ix, rd); // determine acceleration rdd = (rd - rd_r) / dt; - rd_r = rd; - ix += 6; } else { // for pinned rods 3 entries will be taken @@ -697,12 +678,11 @@ moordyn::MoorDyn::Step(const double* x, moordyn::array2vec(xd + ix, rd3); rd(Eigen::seqN(0, 3)) = rd3; // determine acceleration - rdd(Eigen::seqN(0, 3)) = (rd3 - rd3_r) / dt; - rd3_r = rd3; - + rdd(Eigen::seqN(0, 3)) = (rd3 - rd_r.head<3>()) / dt; ix += 3; } - RodList[l]->initiateStep(r, rd, rdd); // acceleration required for inertial terms + // acceleration required for inertial terms + RodList[l]->initiateStep(r, rd, rdd); } for (auto l : CpldPointIs) { vec r, rd; diff --git a/source/MoorDyn2.hpp b/source/MoorDyn2.hpp index 923e153c..599a8234 100644 --- a/source/MoorDyn2.hpp +++ b/source/MoorDyn2.hpp @@ -627,12 +627,6 @@ class MoorDyn final : public io::IO /// (if using env.WaveKin=1) unsigned int npW; - /// Previous time step velocity for calculating coupled point acceleration (first time step assumed stationary) - vec6 rd_b = vec6::Zero(); // body - vec6 rd_r = vec6::Zero(); // coupled rod - vec3 rd3_r = vec3::Zero(); // coupled pinned rod - vec3 rd3_b = vec3::Zero(); // coupled pinned body - /// main output file ofstream outfileMain; diff --git a/source/Rod.hpp b/source/Rod.hpp index d426183f..0107914c 100644 --- a/source/Rod.hpp +++ b/source/Rod.hpp @@ -473,13 +473,21 @@ class Rod final : public io::IO, public SuperCFL * boundary conditions (rod kinematics) for the proceeding time steps * @param r The input position * @param rd The input velocity - * @param rdd The input velocity + * @param rdd The input acceleration * @throw moordyn::invalid_value_error If the rod is not of type * moordyn::Rod::COUPLED or moordyn::Rod::CPLDPIN * @note If the rod is of type moordyn::Rod::CPLDPIN, then just 3 components * of @p r and @p rd are considered */ - void initiateStep(vec6 r_in, vec6 rd_in, vec6 rdd_in); + void initiateStep(vec6 r, vec6 rd, vec6 rdd); + + /** @brief Get the last setted velocity for an unfree rod + * + * For free rods the behaviour is undetermined + * + * @return The velocity (6 dof) + */ + inline vec6 getUnfreeVel() const { return rd_ves; } /** @brief Sets the kinematics * From 6d7b60475c1a8fe8313a399574e8be9c3d05fabf Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Thu, 4 Jul 2024 08:28:28 +0200 Subject: [PATCH 3/8] fix: Odd treatment was meant for indexes from 1 to 3, not 0 to 2, and the matrix indexes were transposed --- source/Misc.hpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/source/Misc.hpp b/source/Misc.hpp index 54832b3e..b12c3034 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -190,9 +190,9 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // s2s3 -c2s1s3 + c1c3 -c2c1s3 - s1c3 // -s2c3 c2s1c3 + c1s3 c2c1c3 - s1s3 // Note: s2 is always positive. - Scalar s2 = Eigen::numext::hypot(coeff(j, i), coeff(k, i)); - if (odd) { - res[0] = atan2(coeff(j, i), coeff(k, i)); + Scalar s2 = Eigen::numext::hypot(coeff(i, j), coeff(i, k)); + if (!odd) { + res[0] = atan2(coeff(i, j), coeff(i, k)); // s2 is always positive, so res[1] will be within the canonical [0, // pi] range res[1] = atan2(s2, coeff(i, i)); @@ -211,7 +211,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // arguments to atan2, while the calculation of the third angle does // not need special adjustment since it uses the adjusted res[0] as // the input and produces a correct result. - res[0] = atan2(-coeff(j, i), -coeff(k, i)); + res[0] = atan2(-coeff(i, j), -coeff(i, k)); res[1] = -atan2(s2, coeff(i, i)); } // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first @@ -225,8 +225,8 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 Scalar s1 = sin(res[0]); Scalar c1 = cos(res[0]); - res[2] = atan2(c1 * coeff(j, k) - s1 * coeff(k, k), - c1 * coeff(j, j) - s1 * coeff(k, j)); + res[2] = atan2(c1 * coeff(k, j) - s1 * coeff(k, k), + c1 * coeff(j, j) - s1 * coeff(j, k)); } else { // Tait-Bryan angles (all three axes are different; typically used for // yaw-pitch-roll calculations). The i, j, k indices enable addressing @@ -236,17 +236,17 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // c2c3 s2s1c3 - c1s3 s2c1c3 + s1s3 // c2s3 s2s1s3 + c1c3 s2c1s3 - s1c3 // -s2 c2s1 c2c1 - res[0] = atan2(coeff(j, k), coeff(k, k)); - Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(i, j)); + res[0] = atan2(coeff(k, j), coeff(k, k)); + Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(j, i)); // c2 is always positive, so the following atan2 will always return a // result in the correct canonical middle angle range [-pi/2, pi/2] - res[1] = atan2(-coeff(i, k), c2); + res[1] = atan2(-coeff(k, i), c2); Scalar s1 = sin(res[0]); Scalar c1 = cos(res[0]); - res[2] = atan2(s1 * coeff(k, i) - c1 * coeff(j, i), - c1 * coeff(j, j) - s1 * coeff(k, j)); + res[2] = atan2(s1 * coeff(i, k) - c1 * coeff(i, j), + c1 * coeff(j, j) - s1 * coeff(j, k)); } - if (!odd) { + if (odd) { res = -res; } return res; @@ -255,7 +255,8 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) inline vec3 Quat2Euler(const quaternion& q) { - return canonicalEulerAngles(q, 0, 1, 2); // 0, 1, 2 correspond to axes leading to XYZ rotation + // 0, 1, 2 correspond to axes leading to XYZ rotation + return canonicalEulerAngles(q, 0, 1, 2); } inline quaternion From b96f74f90c9ceda4d4b81c2e9b73c8f425f2707d Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Thu, 4 Jul 2024 08:58:45 +0200 Subject: [PATCH 4/8] test: Test the rotations --- .github/workflows/python-wrapper.yml | 2 +- tests/Mooring/rotations/nwu.txt | 44 +++++++++++ tests/test_generator.py | 2 - tests/test_minimal.py | 3 +- tests/test_rotations.py | 108 +++++++++++++++++++++++++++ 5 files changed, 154 insertions(+), 5 deletions(-) create mode 100644 tests/Mooring/rotations/nwu.txt create mode 100644 tests/test_rotations.py diff --git a/.github/workflows/python-wrapper.yml b/.github/workflows/python-wrapper.yml index da136713..ed80d237 100644 --- a/.github/workflows/python-wrapper.yml +++ b/.github/workflows/python-wrapper.yml @@ -31,7 +31,7 @@ jobs: id: setup-python - name: Install Python dependencies - run: pip install --upgrade build pytest + run: pip install --upgrade build pytest numpy scipy - name: Install VTK run: | diff --git a/tests/Mooring/rotations/nwu.txt b/tests/Mooring/rotations/nwu.txt new file mode 100644 index 00000000..6ab64ee9 --- /dev/null +++ b/tests/Mooring/rotations/nwu.txt @@ -0,0 +1,44 @@ +--------------------- MoorDyn Input File ------------------------------------------------------- +Testing the rotations assuming the NWU angles criteria, i.e. the bow is +pointing towards the x axis, the portside towards the y axis and the z axis +upwards. +To check the rotations we set 3 points pointing on the positive x, y and z + +The rotations are thus consistent with the EulerXYZ criteria: + + - Roll: rotates around x axis, in such a way the starboard side (-y dir) goes + down + - Pitch: rotates around y axis, in such a way the bow (+x dir) goes down + - Yaw: rotates around y axis, in such a way the bow (+x dir) goes to the + portside + +The rotations can be contatenated on a straight forward way, i.e. you can get +the current state and just compute the velocity as the difference. e.g. + +double r[6], rd[6]; +const auto [r, tmp] = MoorDyn_GetBodyState(body, r, rd); +const double r_next[6] = {x, y, z, roll, pitch, yaw}; +for (unsigned int i = 0; i < 6; i++) { + rd[i] = (r_next[i] - r[i]) / dt; +} +MoorDyn_Step(system, r, rd, t, dt); +---------------------------- BODIES ----------------------------------------------------- +ID Attachment X0 Y0 Z0 r0 p0 y0 Mass CG* I* Volume CdA* Ca +(#) (-) (m) (m) (m) (deg) (deg) (deg) (kg) (m) (kg-m^2) (m^3) (m^2) (-) +1 coupled 0 0 0 0 0 0 0 0 0 0 0 0 +----------------------- POINTS ---------------------------------------------- +Node Type X Y Z M V CdA CA +(-) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Body1 1 0 0 0 0 0 0 +2 Body1 0 1 0 0 0 0 0 +3 Body1 0 0 1 0 0 0 0 +-------------------------- SOLVER OPTIONS--------------------------------------------------- +0 writeLog - Write a log file +0.0 g - No gravity +1.0 dtM - time step to use in mooring integration +3.0e6 kb - bottom stiffness +3.0e5 cb - bottom damping +100 WtrDpth - water depth +0.0 TmaxIC - threshold for IC convergence +euler tScheme - Just one step per time step +--------------------------- need this line ------------------------------------------------- diff --git a/tests/test_generator.py b/tests/test_generator.py index 2589449d..f6d3bc79 100644 --- a/tests/test_generator.py +++ b/tests/test_generator.py @@ -1,6 +1,4 @@ import sys -print('HERE') -print(sys.path) from unittest import TestCase, main as unittest_main import moordyn diff --git a/tests/test_minimal.py b/tests/test_minimal.py index 7d668109..fe610a0f 100644 --- a/tests/test_minimal.py +++ b/tests/test_minimal.py @@ -1,5 +1,4 @@ import sys -print(sys.path) from unittest import TestCase, main as unittest_main import os import tempfile @@ -45,6 +44,7 @@ def test_one_iter(self): cwd = os.getcwd() os.chdir(tmp_folder) system = moordyn.Create() + os.chdir(cwd) x = [] for i in range(4, 7): point = moordyn.GetPoint(system, i) @@ -55,7 +55,6 @@ def test_one_iter(self): self.assertEqual(moordyn.NCoupledDOF(system), 9, "Wrong number of coupled DOFs") self.assertEqual(moordyn.GetNumberLines(system), 3) - os.chdir(cwd) v[0] = 0.1 forces = moordyn.Step(system, x, v, 0.0, 0.5) print(forces) diff --git a/tests/test_rotations.py b/tests/test_rotations.py new file mode 100644 index 00000000..12c4a2ea --- /dev/null +++ b/tests/test_rotations.py @@ -0,0 +1,108 @@ +import sys +from unittest import TestCase, main as unittest_main +import os +import moordyn +from scipy.spatial.transform import Rotation as R +import numpy as np +from test_minimal import setup_case + + +ENUREF = [[1, 0, 0], + [0, 1, 0], + [0, 0, 1]] +ANGLE = 5 + + +def rotate_vecs(vecs, angles, seq='xyz', degrees=True): + r = R.from_euler(seq, angles, degrees=degrees) + return [r.apply(vec) for vec in vecs] + + +def compare_vecs(vecs, refs, tol=1e-5): + assert len(vecs) == len(refs) + for i in range(len(vecs)): + assert len(vecs[i]) == len(refs[i]) + for j in range(len(vecs[i])): + assert abs(vecs[i][j] - refs[i][j]) < 1e-5 + + +def abseuler2releuler(org, dst, seq='xyz', degrees=True): + r0 = R.from_euler(seq, org, degrees=degrees).inv() + r1 = R.from_euler(seq, dst, degrees=degrees) + return (r1 * r0).as_euler(seq, degrees=degrees) + + +def rotate_moordyn(system, roll, pitch, yaw, t, dt=1.0, degrees=True): + if degrees: + roll = np.radians(roll) + pitch = np.radians(pitch) + yaw = np.radians(yaw) + pts = [moordyn.GetPoint(system, i + 1) for i in range(3)] + body = moordyn.GetBody(system, 1) + r, _ = moordyn.GetBodyState(body) + u = [0] * 3 + [roll / dt, pitch / dt, yaw / dt] + moordyn.Step(system, r, u, t, dt) + t += dt + r = [moordyn.GetPointPos(pt) for pt in pts] + return r, t + + +class RotationTests(TestCase): + def setUp(self): + pass + + def test_nwu(self): + tmp_folder = setup_case("rotations/nwu.txt") + cwd = os.getcwd() + os.chdir(tmp_folder) + system = moordyn.Create() + os.chdir(cwd) + pts = [moordyn.GetPoint(system, i + 1) for i in range(3)] + r = [0] * 6 + u = [0] * 6 + moordyn.Init(system, r, u) + # Check the intial condition + compare_vecs([moordyn.GetPointPos(pt) for pt in pts], ENUREF) + t = 0 + # Roll test + r, t = rotate_moordyn(system, ANGLE, 0, 0, t) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='x')) + assert(r[1][2] > 0.0) # Portside up, i.e. starboard down + r, t = rotate_moordyn(system, -ANGLE, 0, 0, t) + compare_vecs(r, ENUREF) + # Pitch test + r, t = rotate_moordyn(system, 0, ANGLE, 0, t) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='y')) + assert(r[0][2] < 0.0) # Bow down + r, t = rotate_moordyn(system, 0, -ANGLE, 0, t) + compare_vecs(r, ENUREF) + # Yaw test + r, t = rotate_moordyn(system, 0, 0, ANGLE, t) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='z')) + assert(r[0][1] > 0.0) # Bow to portside + r, t = rotate_moordyn(system, 0, 0, -ANGLE, t) + compare_vecs(r, ENUREF) + + # Test the angles order. For that we are rotating it in all directions + # simultaneously + r, t = rotate_moordyn(system, ANGLE, ANGLE, ANGLE, t) + compare_vecs(r, rotate_vecs(ENUREF, [ANGLE, ANGLE, ANGLE])) + r, t = rotate_moordyn(system, -ANGLE, -ANGLE, -ANGLE, t) + compare_vecs(r, ENUREF) + + # Test the rotations concatenations. To do that we are again rotating on + # the 3 axis simultaneously. However, this time we rotate twice, so the + # final angle shall be 2*ANGLE in all directions + rotate_moordyn(system, ANGLE, ANGLE, ANGLE, t) + r, t = rotate_moordyn(system, ANGLE, ANGLE, ANGLE, t) + roll, pitch, yaw = abseuler2releuler([ANGLE, ANGLE, ANGLE], + [2 * ANGLE, 2 * ANGLE, 2 * ANGLE]) + compare_vecs(r, rotate_vecs(rotate_vecs(ENUREF, [ANGLE, ANGLE, ANGLE]), + [roll, pitch, yaw])) + compare_vecs(r, rotate_vecs(ENUREF, [2 * ANGLE, 2 * ANGLE, 2 * ANGLE])) + + moordyn.Close(system) + + +if __name__ == '__main__': + unittest_main() From 75f6d5f531130fb69e386f82d6cae72f9b16a0b1 Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Thu, 4 Jul 2024 11:49:41 +0200 Subject: [PATCH 5/8] test: Typo on the config file description --- tests/Mooring/rotations/nwu.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/Mooring/rotations/nwu.txt b/tests/Mooring/rotations/nwu.txt index 6ab64ee9..945f79ad 100644 --- a/tests/Mooring/rotations/nwu.txt +++ b/tests/Mooring/rotations/nwu.txt @@ -9,7 +9,7 @@ The rotations are thus consistent with the EulerXYZ criteria: - Roll: rotates around x axis, in such a way the starboard side (-y dir) goes down - Pitch: rotates around y axis, in such a way the bow (+x dir) goes down - - Yaw: rotates around y axis, in such a way the bow (+x dir) goes to the + - Yaw: rotates around z axis, in such a way the bow (+x dir) goes to the portside The rotations can be contatenated on a straight forward way, i.e. you can get From e1274d224a33f10b083ee55dfc067477c02341f6 Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Tue, 9 Jul 2024 06:55:55 +0200 Subject: [PATCH 6/8] fix: EulerXYZ intrinsic angles instead of extrinsic --- instrinsic_angles.patch | 119 ++++++++++++++++++++++++++++++++++++++++ source/Misc.hpp | 36 ++++++------ tests/test_rotations.py | 10 ++-- 3 files changed, 142 insertions(+), 23 deletions(-) create mode 100644 instrinsic_angles.patch diff --git a/instrinsic_angles.patch b/instrinsic_angles.patch new file mode 100644 index 00000000..c780b2e6 --- /dev/null +++ b/instrinsic_angles.patch @@ -0,0 +1,119 @@ +diff --git a/source/Misc.hpp b/source/Misc.hpp +index b12c303..77d6b33 100644 +--- a/source/Misc.hpp ++++ b/source/Misc.hpp +@@ -168,11 +168,11 @@ EqualRealNos(const real a1, const real a2) + return std::abs(a1 - a2) <= fraction * tol; + } + +- +-inline vec3 ++inline vec3 + canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + { +- // From issue #163: https://github.com/FloatingArrayDesign/MoorDyn/issues/163 ++ // From issue #163: ++ // https://github.com/FloatingArrayDesign/MoorDyn/issues/163 + mat3 coeff = quat.normalized().toRotationMatrix(); + vec3 res{}; + using Index = int; +@@ -190,9 +190,9 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + // s2s3 -c2s1s3 + c1c3 -c2c1s3 - s1c3 + // -s2c3 c2s1c3 + c1s3 c2c1c3 - s1s3 + // Note: s2 is always positive. +- Scalar s2 = Eigen::numext::hypot(coeff(i, j), coeff(i, k)); +- if (!odd) { +- res[0] = atan2(coeff(i, j), coeff(i, k)); ++ Scalar s2 = Eigen::numext::hypot(coeff(j, i), coeff(k, i)); ++ if (odd) { ++ res[0] = atan2(coeff(j, i), coeff(k, i)); + // s2 is always positive, so res[1] will be within the canonical [0, + // pi] range + res[1] = atan2(s2, coeff(i, i)); +@@ -211,7 +211,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + // arguments to atan2, while the calculation of the third angle does + // not need special adjustment since it uses the adjusted res[0] as + // the input and produces a correct result. +- res[0] = atan2(-coeff(i, j), -coeff(i, k)); ++ res[0] = atan2(-coeff(j, i), -coeff(k, i)); + res[1] = -atan2(s2, coeff(i, i)); + } + // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first +@@ -225,8 +225,8 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 + Scalar s1 = sin(res[0]); + Scalar c1 = cos(res[0]); +- res[2] = atan2(c1 * coeff(k, j) - s1 * coeff(k, k), +- c1 * coeff(j, j) - s1 * coeff(j, k)); ++ res[2] = atan2(c1 * coeff(j, k) - s1 * coeff(k, k), ++ c1 * coeff(j, j) - s1 * coeff(k, j)); + } else { + // Tait-Bryan angles (all three axes are different; typically used for + // yaw-pitch-roll calculations). The i, j, k indices enable addressing +@@ -236,17 +236,17 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + // c2c3 s2s1c3 - c1s3 s2c1c3 + s1s3 + // c2s3 s2s1s3 + c1c3 s2c1s3 - s1c3 + // -s2 c2s1 c2c1 +- res[0] = atan2(coeff(k, j), coeff(k, k)); +- Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(j, i)); ++ res[0] = atan2(coeff(j, k), coeff(k, k)); ++ Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(i, j)); + // c2 is always positive, so the following atan2 will always return a + // result in the correct canonical middle angle range [-pi/2, pi/2] +- res[1] = atan2(-coeff(k, i), c2); ++ res[1] = atan2(-coeff(i, k), c2); + Scalar s1 = sin(res[0]); + Scalar c1 = cos(res[0]); +- res[2] = atan2(s1 * coeff(i, k) - c1 * coeff(i, j), +- c1 * coeff(j, j) - s1 * coeff(j, k)); ++ res[2] = atan2(s1 * coeff(k, i) - c1 * coeff(j, i), ++ c1 * coeff(j, j) - s1 * coeff(k, j)); + } +- if (odd) { ++ if (!odd) { + res = -res; + } + return res; +@@ -255,7 +255,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) + inline vec3 + Quat2Euler(const quaternion& q) + { +- // 0, 1, 2 correspond to axes leading to XYZ rotation ++ // 0, 1, 2 correspond to axes leading to XYZ rotation + return canonicalEulerAngles(q, 0, 1, 2); + } + +@@ -263,9 +263,9 @@ inline quaternion + Euler2Quat(const vec3& angles) + { + using AngleAxis = Eigen::AngleAxis; +- quaternion q = AngleAxis(angles.z(), vec3::UnitZ()) * ++ quaternion q = AngleAxis(angles.x(), vec3::UnitX()) * + AngleAxis(angles.y(), vec3::UnitY()) * +- AngleAxis(angles.x(), vec3::UnitX()); ++ AngleAxis(angles.z(), vec3::UnitZ()); + return q; + } + +diff --git a/tests/test_rotations.py b/tests/test_rotations.py +index 12c4a2e..dfa4a6a 100644 +--- a/tests/test_rotations.py ++++ b/tests/test_rotations.py +@@ -13,7 +13,7 @@ ENUREF = [[1, 0, 0], + ANGLE = 5 + + +-def rotate_vecs(vecs, angles, seq='xyz', degrees=True): ++def rotate_vecs(vecs, angles, seq="XYZ", degrees=True): + r = R.from_euler(seq, angles, degrees=degrees) + return [r.apply(vec) for vec in vecs] + +@@ -26,7 +26,7 @@ def compare_vecs(vecs, refs, tol=1e-5): + assert abs(vecs[i][j] - refs[i][j]) < 1e-5 + + +-def abseuler2releuler(org, dst, seq='xyz', degrees=True): ++def abseuler2releuler(org, dst, seq="XYZ", degrees=True): + r0 = R.from_euler(seq, org, degrees=degrees).inv() + r1 = R.from_euler(seq, dst, degrees=degrees) + return (r1 * r0).as_euler(seq, degrees=degrees) diff --git a/source/Misc.hpp b/source/Misc.hpp index b12c3034..77d6b336 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -168,11 +168,11 @@ EqualRealNos(const real a1, const real a2) return std::abs(a1 - a2) <= fraction * tol; } - -inline vec3 +inline vec3 canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) { - // From issue #163: https://github.com/FloatingArrayDesign/MoorDyn/issues/163 + // From issue #163: + // https://github.com/FloatingArrayDesign/MoorDyn/issues/163 mat3 coeff = quat.normalized().toRotationMatrix(); vec3 res{}; using Index = int; @@ -190,9 +190,9 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // s2s3 -c2s1s3 + c1c3 -c2c1s3 - s1c3 // -s2c3 c2s1c3 + c1s3 c2c1c3 - s1s3 // Note: s2 is always positive. - Scalar s2 = Eigen::numext::hypot(coeff(i, j), coeff(i, k)); - if (!odd) { - res[0] = atan2(coeff(i, j), coeff(i, k)); + Scalar s2 = Eigen::numext::hypot(coeff(j, i), coeff(k, i)); + if (odd) { + res[0] = atan2(coeff(j, i), coeff(k, i)); // s2 is always positive, so res[1] will be within the canonical [0, // pi] range res[1] = atan2(s2, coeff(i, i)); @@ -211,7 +211,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // arguments to atan2, while the calculation of the third angle does // not need special adjustment since it uses the adjusted res[0] as // the input and produces a correct result. - res[0] = atan2(-coeff(i, j), -coeff(i, k)); + res[0] = atan2(-coeff(j, i), -coeff(k, i)); res[1] = -atan2(s2, coeff(i, i)); } // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first @@ -225,8 +225,8 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 Scalar s1 = sin(res[0]); Scalar c1 = cos(res[0]); - res[2] = atan2(c1 * coeff(k, j) - s1 * coeff(k, k), - c1 * coeff(j, j) - s1 * coeff(j, k)); + res[2] = atan2(c1 * coeff(j, k) - s1 * coeff(k, k), + c1 * coeff(j, j) - s1 * coeff(k, j)); } else { // Tait-Bryan angles (all three axes are different; typically used for // yaw-pitch-roll calculations). The i, j, k indices enable addressing @@ -236,17 +236,17 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) // c2c3 s2s1c3 - c1s3 s2c1c3 + s1s3 // c2s3 s2s1s3 + c1c3 s2c1s3 - s1c3 // -s2 c2s1 c2c1 - res[0] = atan2(coeff(k, j), coeff(k, k)); - Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(j, i)); + res[0] = atan2(coeff(j, k), coeff(k, k)); + Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(i, j)); // c2 is always positive, so the following atan2 will always return a // result in the correct canonical middle angle range [-pi/2, pi/2] - res[1] = atan2(-coeff(k, i), c2); + res[1] = atan2(-coeff(i, k), c2); Scalar s1 = sin(res[0]); Scalar c1 = cos(res[0]); - res[2] = atan2(s1 * coeff(i, k) - c1 * coeff(i, j), - c1 * coeff(j, j) - s1 * coeff(j, k)); + res[2] = atan2(s1 * coeff(k, i) - c1 * coeff(j, i), + c1 * coeff(j, j) - s1 * coeff(k, j)); } - if (odd) { + if (!odd) { res = -res; } return res; @@ -255,7 +255,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) inline vec3 Quat2Euler(const quaternion& q) { - // 0, 1, 2 correspond to axes leading to XYZ rotation + // 0, 1, 2 correspond to axes leading to XYZ rotation return canonicalEulerAngles(q, 0, 1, 2); } @@ -263,9 +263,9 @@ inline quaternion Euler2Quat(const vec3& angles) { using AngleAxis = Eigen::AngleAxis; - quaternion q = AngleAxis(angles.z(), vec3::UnitZ()) * + quaternion q = AngleAxis(angles.x(), vec3::UnitX()) * AngleAxis(angles.y(), vec3::UnitY()) * - AngleAxis(angles.x(), vec3::UnitX()); + AngleAxis(angles.z(), vec3::UnitZ()); return q; } diff --git a/tests/test_rotations.py b/tests/test_rotations.py index 12c4a2ea..92e330b5 100644 --- a/tests/test_rotations.py +++ b/tests/test_rotations.py @@ -13,7 +13,7 @@ ANGLE = 5 -def rotate_vecs(vecs, angles, seq='xyz', degrees=True): +def rotate_vecs(vecs, angles, seq="XYZ", degrees=True): r = R.from_euler(seq, angles, degrees=degrees) return [r.apply(vec) for vec in vecs] @@ -26,7 +26,7 @@ def compare_vecs(vecs, refs, tol=1e-5): assert abs(vecs[i][j] - refs[i][j]) < 1e-5 -def abseuler2releuler(org, dst, seq='xyz', degrees=True): +def abseuler2releuler(org, dst, seq="XYZ", degrees=True): r0 = R.from_euler(seq, org, degrees=degrees).inv() r1 = R.from_euler(seq, dst, degrees=degrees) return (r1 * r0).as_euler(seq, degrees=degrees) @@ -66,19 +66,19 @@ def test_nwu(self): t = 0 # Roll test r, t = rotate_moordyn(system, ANGLE, 0, 0, t) - compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='x')) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='X')) assert(r[1][2] > 0.0) # Portside up, i.e. starboard down r, t = rotate_moordyn(system, -ANGLE, 0, 0, t) compare_vecs(r, ENUREF) # Pitch test r, t = rotate_moordyn(system, 0, ANGLE, 0, t) - compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='y')) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='Y')) assert(r[0][2] < 0.0) # Bow down r, t = rotate_moordyn(system, 0, -ANGLE, 0, t) compare_vecs(r, ENUREF) # Yaw test r, t = rotate_moordyn(system, 0, 0, ANGLE, t) - compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='z')) + compare_vecs(r, rotate_vecs(ENUREF, ANGLE, seq='Z')) assert(r[0][1] > 0.0) # Bow to portside r, t = rotate_moordyn(system, 0, 0, -ANGLE, t) compare_vecs(r, ENUREF) From 49fac259ab9d93d5d87505bf1a4267de32a7c81c Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Tue, 9 Jul 2024 07:11:11 +0200 Subject: [PATCH 7/8] docs: Clarify that intrinsic angles are considered, and link to the external resources --- docs/angles.svg | 203 ---------------------------------------------- docs/features.rst | 49 +++-------- 2 files changed, 11 insertions(+), 241 deletions(-) delete mode 100644 docs/angles.svg diff --git a/docs/angles.svg b/docs/angles.svg deleted file mode 100644 index 9a5fcd18..00000000 --- a/docs/angles.svg +++ /dev/null @@ -1,203 +0,0 @@ - - - - - - - - - - - - - xyz - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/features.rst b/docs/features.rst index 2a4b038d..9deb43d4 100644 --- a/docs/features.rst +++ b/docs/features.rst @@ -29,44 +29,17 @@ The main difference between MoorDyn-C and MoorDyn-F is that MoorDyn-C uses quate Orientation of 6 DOF objects: ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Euler angles – MoorDyn-F -"""""""""""""""""""""""" - -In the following figure the 6DOF object orientation angles convention is depicted: - -.. figure:: angles.svg - :alt: Angles criteria schematic view - -The roll and yaw angles, :math:`\phi` and :math:`\psi`, follow the -right hand criteria, while the pitch angle, :math:`\theta`, follows the left -hand criteria. -This way the classic rotation matrices can be considered, - -.. math:: - \begin{alignat}{1} - R_x(\phi) &= \begin{bmatrix} - 1 & 0 & 0 \\ - 0 & \cos \phi & -\sin \phi \\[3pt] - 0 & \sin \phi & \cos \phi \\[3pt] - \end{bmatrix} \\[6pt] - R_y(\theta) &= \begin{bmatrix} - \cos \theta & 0 & \sin \theta \\[3pt] - 0 & 1 & 0 \\[3pt] - -\sin \theta & 0 & \cos \theta \\ - \end{bmatrix} \\[6pt] - R_z(\psi) &= \begin{bmatrix} - \cos \psi & -\sin \psi & 0 \\[3pt] - \sin \psi & \cos \psi & 0 \\[3pt] - 0 & 0 & 1 \\ - \end{bmatrix} - \end{alignat} - - -Quaternions – MoorDyn-C -""""""""""""""""""""""" - -The latest MoorDyn-C internally uses quaternions to describe the location and orientation of 6 DOF objects. Externally MoorDyn-C behaves the same as MoorDyn-F, using Euler angles for both inputs and outputs. Quaternions are a common alternative to Euler angles for describing orientations of 3D objects. -Further description of quaternions can be found in PR #90 in the MoorDyn repository, put together by Alex Kinley of Kelson Marine: https://github.com/FloatingArrayDesign/MoorDyn/pull/90#issue-1777700494 +MoorDyn-C, MoorDyn-F and `MoorPy `_ share the +same Intrinsic Euler-XYZ (Tait-Bryan) angles criteria to compute orientations. +You can learn more about this on +`Hall M. Generalized Quasi-Static Mooring System Modeling with Analytic Jacobians. Energies. 2024; 17(13):3155. https://doi.org/10.3390/en17133155 `_ + +However, while on MoorDyn-F this is handled by considering orientation +matrices, on MoorDyn-C quaternions are considered to describe the location and +orientation of 6 DOF objects. +Further description of quaternions can be found in PR #90 in the MoorDyn +repository, put together by Alex Kinley of Kelson Marine: +https://github.com/FloatingArrayDesign/MoorDyn/pull/90#issue-1777700494 References ---------- From bfe9187f8f933cb7206d23bcf344b0a95809703c Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos-Pita Date: Wed, 10 Jul 2024 06:42:57 +0200 Subject: [PATCH 8/8] fix: Drop the patch to move from extrinsic to intrinsic Euler angles --- instrinsic_angles.patch | 119 ---------------------------------------- 1 file changed, 119 deletions(-) delete mode 100644 instrinsic_angles.patch diff --git a/instrinsic_angles.patch b/instrinsic_angles.patch deleted file mode 100644 index c780b2e6..00000000 --- a/instrinsic_angles.patch +++ /dev/null @@ -1,119 +0,0 @@ -diff --git a/source/Misc.hpp b/source/Misc.hpp -index b12c303..77d6b33 100644 ---- a/source/Misc.hpp -+++ b/source/Misc.hpp -@@ -168,11 +168,11 @@ EqualRealNos(const real a1, const real a2) - return std::abs(a1 - a2) <= fraction * tol; - } - -- --inline vec3 -+inline vec3 - canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - { -- // From issue #163: https://github.com/FloatingArrayDesign/MoorDyn/issues/163 -+ // From issue #163: -+ // https://github.com/FloatingArrayDesign/MoorDyn/issues/163 - mat3 coeff = quat.normalized().toRotationMatrix(); - vec3 res{}; - using Index = int; -@@ -190,9 +190,9 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - // s2s3 -c2s1s3 + c1c3 -c2c1s3 - s1c3 - // -s2c3 c2s1c3 + c1s3 c2c1c3 - s1s3 - // Note: s2 is always positive. -- Scalar s2 = Eigen::numext::hypot(coeff(i, j), coeff(i, k)); -- if (!odd) { -- res[0] = atan2(coeff(i, j), coeff(i, k)); -+ Scalar s2 = Eigen::numext::hypot(coeff(j, i), coeff(k, i)); -+ if (odd) { -+ res[0] = atan2(coeff(j, i), coeff(k, i)); - // s2 is always positive, so res[1] will be within the canonical [0, - // pi] range - res[1] = atan2(s2, coeff(i, i)); -@@ -211,7 +211,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - // arguments to atan2, while the calculation of the third angle does - // not need special adjustment since it uses the adjusted res[0] as - // the input and produces a correct result. -- res[0] = atan2(-coeff(i, j), -coeff(i, k)); -+ res[0] = atan2(-coeff(j, i), -coeff(k, i)); - res[1] = -atan2(s2, coeff(i, i)); - } - // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first -@@ -225,8 +225,8 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 - Scalar s1 = sin(res[0]); - Scalar c1 = cos(res[0]); -- res[2] = atan2(c1 * coeff(k, j) - s1 * coeff(k, k), -- c1 * coeff(j, j) - s1 * coeff(j, k)); -+ res[2] = atan2(c1 * coeff(j, k) - s1 * coeff(k, k), -+ c1 * coeff(j, j) - s1 * coeff(k, j)); - } else { - // Tait-Bryan angles (all three axes are different; typically used for - // yaw-pitch-roll calculations). The i, j, k indices enable addressing -@@ -236,17 +236,17 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - // c2c3 s2s1c3 - c1s3 s2c1c3 + s1s3 - // c2s3 s2s1s3 + c1c3 s2c1s3 - s1c3 - // -s2 c2s1 c2c1 -- res[0] = atan2(coeff(k, j), coeff(k, k)); -- Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(j, i)); -+ res[0] = atan2(coeff(j, k), coeff(k, k)); -+ Scalar c2 = Eigen::numext::hypot(coeff(i, i), coeff(i, j)); - // c2 is always positive, so the following atan2 will always return a - // result in the correct canonical middle angle range [-pi/2, pi/2] -- res[1] = atan2(-coeff(k, i), c2); -+ res[1] = atan2(-coeff(i, k), c2); - Scalar s1 = sin(res[0]); - Scalar c1 = cos(res[0]); -- res[2] = atan2(s1 * coeff(i, k) - c1 * coeff(i, j), -- c1 * coeff(j, j) - s1 * coeff(j, k)); -+ res[2] = atan2(s1 * coeff(k, i) - c1 * coeff(j, i), -+ c1 * coeff(j, j) - s1 * coeff(k, j)); - } -- if (odd) { -+ if (!odd) { - res = -res; - } - return res; -@@ -255,7 +255,7 @@ canonicalEulerAngles(const quaternion& quat, int a0, int a1, int a2) - inline vec3 - Quat2Euler(const quaternion& q) - { -- // 0, 1, 2 correspond to axes leading to XYZ rotation -+ // 0, 1, 2 correspond to axes leading to XYZ rotation - return canonicalEulerAngles(q, 0, 1, 2); - } - -@@ -263,9 +263,9 @@ inline quaternion - Euler2Quat(const vec3& angles) - { - using AngleAxis = Eigen::AngleAxis; -- quaternion q = AngleAxis(angles.z(), vec3::UnitZ()) * -+ quaternion q = AngleAxis(angles.x(), vec3::UnitX()) * - AngleAxis(angles.y(), vec3::UnitY()) * -- AngleAxis(angles.x(), vec3::UnitX()); -+ AngleAxis(angles.z(), vec3::UnitZ()); - return q; - } - -diff --git a/tests/test_rotations.py b/tests/test_rotations.py -index 12c4a2e..dfa4a6a 100644 ---- a/tests/test_rotations.py -+++ b/tests/test_rotations.py -@@ -13,7 +13,7 @@ ENUREF = [[1, 0, 0], - ANGLE = 5 - - --def rotate_vecs(vecs, angles, seq='xyz', degrees=True): -+def rotate_vecs(vecs, angles, seq="XYZ", degrees=True): - r = R.from_euler(seq, angles, degrees=degrees) - return [r.apply(vec) for vec in vecs] - -@@ -26,7 +26,7 @@ def compare_vecs(vecs, refs, tol=1e-5): - assert abs(vecs[i][j] - refs[i][j]) < 1e-5 - - --def abseuler2releuler(org, dst, seq='xyz', degrees=True): -+def abseuler2releuler(org, dst, seq="XYZ", degrees=True): - r0 = R.from_euler(seq, org, degrees=degrees).inv() - r1 = R.from_euler(seq, dst, degrees=degrees) - return (r1 * r0).as_euler(seq, degrees=degrees)