TidalPotential.h

Go to the documentation of this file.

#ifndef __UNIT_TESTS_TIDAL_POTENTIAL_H
#define __UNIT_TESTS_TIDAL_POTENTIAL_H

#include "EccentricOrbit.h"
#include <cmath>

namespace Evolve {
    class TidalPotential {
    private:
        EccentricOrbit __orbit;

        double
            __inclination,

            __arg_of_periapsis;
            
    public:
        TidalPotential(
            double primary_mass=Core::NaN,

            double secondary_mass=Core::NaN,

            double semimajor=Core::NaN,

            double eccentricity=Core::NaN,

            double inclination=Core::NaN,
            
            double arg_of_periapsis=Core::NaN
        ) :
            __orbit(primary_mass, secondary_mass, semimajor, eccentricity),
            __inclination(inclination),
            __arg_of_periapsis(arg_of_periapsis)
        {}

        double inclination() const {return __inclination;}

        double &inclination() {return __inclination;}

        double arg_of_periapsis() const {return __arg_of_periapsis;}

        double &arg_of_periapsis() {return __arg_of_periapsis;}

        const EccentricOrbit &orbit() const {return __orbit;}

        EccentricOrbit &orbit() {return __orbit;}

        template<class POSITION_TYPE>
            double operator()(
                const POSITION_TYPE &position,

                double time
            ) const;
    };

    template<class POSITION_TYPE>
        double TidalPotential::operator()(const POSITION_TYPE &position,
                                          double time) const
        {
            Eigen::Vector3d secondary_position = __orbit.secondary_position(
                2.0 * M_PI * time / __orbit.orbital_period()
            );


            //Rotate around L_hat to a coordinate system with z along L and y
            //along SxL
            double z_rotated_secondary_x = (
                secondary_position[0] * std::cos(__arg_of_periapsis)
                -
                secondary_position[1] * std::sin(__arg_of_periapsis)
            );
            double z_rotated_secondary_y = (
                secondary_position[0] * std::sin(__arg_of_periapsis)
                +
                secondary_position[1] * std::cos(__arg_of_periapsis)
            );

            //Rotated around SxL to the final coordinate system.
            Eigen::Vector3d transformed_secondary_position(
                (
                    z_rotated_secondary_x * std::cos(__inclination)
                    +
                    secondary_position[2] * std::sin(__inclination)
                ),
                z_rotated_secondary_y,
                (
                    -z_rotated_secondary_x * std::sin(__inclination)
                    +
                    secondary_position[2] * std::cos(__inclination)
                )
            );

            double center_to_secondary = transformed_secondary_position.norm();
            double position_to_secondary = (
                position
                -
                transformed_secondary_position
            ).norm();
            return (
                Core::AstroConst::G
                *
                __orbit.secondary_mass() * Core::AstroConst::solar_mass
                *
                (
                    position.dot(transformed_secondary_position)
                    /
                    std::pow(center_to_secondary, 3)
                    -
                    1.0 / position_to_secondary
                    +
                    1.0 / center_to_secondary
                ) / Core::AstroConst::solar_radius
            );
        }
} //End Evolve namespace

#endif