orbits.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once
 
#include "shambackends/math.hpp"
#include "shammath/matrix.hpp"
#include "shamunits/Constants.hpp"
#include "shamunits/UnitSystem.hpp"
 
namespace shamphys {
 
    inline auto get_binary_pair(double m1, double m2, double a, double e, double nu, double G) {
 
        double M = m1 + m2;
        // double mu = m1 * m2 / M;
 
        double r = a * (1 - e * e) / (1 + e * sycl::cos(nu));
 
        double x_orb = r * sycl::cos(nu);
        double y_orb = r * sycl::sin(nu);
 
        double h      = sycl::sqrt(G * M * a * (1 - e * e));
        double vx_orb = -G * M / h * sycl::sin(nu);
        double vy_orb = G * M / h * (e + sycl::cos(nu));
 
        f64_3 r_orb = {x_orb, y_orb, 0};
        f64_3 v_orb = {vx_orb, vy_orb, 0};
 
        f64_3 r1 = -m2 / M * r_orb;
        f64_3 r2 = m1 / M * r_orb;
 
        f64_3 v1 = -m2 / M * v_orb;
        f64_3 v2 = m1 / M * v_orb;
 
        return std::make_tuple(r1, r2, v1, v2);
    }
 
    auto get_binary_pair(
        double m1,
        double m2,
        double a,
        double e,
        double nu,
        const shamunits::UnitSystem<double> usys = shamunits::UnitSystem<double>{}) {
 
        double G = shamunits::Constants{usys}.G();
 
        return get_binary_pair(m1, m2, a, e, nu, G);
    }
 
    static f64_3x3 rotation_matrix(double roll, double pitch, double yaw) {
 
        // clang-format off
        // Rotation matrix around X-axis (roll)
        f64_3x3 Rx = {1, 0, 0,
                           0, sycl::cos(roll), -sycl::sin(roll),
                           0, sycl::sin(roll), sycl::cos(roll)};
 
        // Rotation matrix around Y-axis (pitch)
        f64_3x3 Ry = {sycl::cos(pitch), 0, sycl::sin(pitch),
                           0, 1, 0,
                           -sycl::sin(pitch), 0, sycl::cos(pitch)};
 
        // Rotation matrix around Z-axis (yaw)
        f64_3x3 Rz = {sycl::cos(yaw), -sycl::sin(yaw), 0,
                           sycl::sin(yaw), sycl::cos(yaw), 0,
                           0, 0, 1};
        // clang-format on
 
        f64_3x3 Ryx = {};
        f64_3x3 R   = {};
 
        // Combine the rotations (R = Rz * Ry * Rx)
        shammath::mat_prod(Ry.get_mdspan(), Rx.get_mdspan(), Ryx.get_mdspan());
        shammath::mat_prod(Rz.get_mdspan(), Ryx.get_mdspan(), R.get_mdspan());
 
        return R;
    }
 
    static f64_3 rotate_point(const f64_3 &point, double roll, double pitch, double yaw) {
        shammath::mat<f64, 3, 1> r = {point.x(), point.y(), point.z()};
 
        // Get the rotation matrix
        f64_3x3 R = rotation_matrix(roll, pitch, yaw);
 
        // Perform the rotation by multiplying the point with the rotation matrix
        shammath::mat<f64, 3, 1> rotated_point = {};
        shammath::mat_prod(R.get_mdspan(), r.get_mdspan(), rotated_point.get_mdspan());
 
        return {rotated_point.data[0], rotated_point.data[1], rotated_point.data[2]};
    }
 
    auto get_binary_rotated(
        double m1,
        double m2,
        double a,
        double e,
        double nu,
        double G,
        double roll,
        double pitch,
        double yaw) {
 
        auto [r1, r2, v1, v2] = get_binary_pair(m1, m2, a, e, nu, G);
 
        r1 = rotate_point({r1.x(), r1.y(), r1.z()}, roll, pitch, yaw);
        r2 = rotate_point({r2.x(), r2.y(), r2.z()}, roll, pitch, yaw);
        v1 = rotate_point({v1.x(), v1.y(), v1.z()}, roll, pitch, yaw);
        v2 = rotate_point({v2.x(), v2.y(), v2.z()}, roll, pitch, yaw);
 
        return std::make_tuple(r1, r2, v1, v2);
    }
 
    auto get_binary_rotated(
        double m1,
        double m2,
        double a,
        double e,
        double nu,
        const shamunits::UnitSystem<double> usys,
        double roll,
        double pitch,
        double yaw) {
        return get_binary_rotated(
            m1, m2, a, e, nu, shamunits::Constants{usys}.G(), roll, pitch, yaw);
    }
 
} // namespace shamphys