TidalPotentialExpansion.h

Go to the documentation of this file.

#ifndef __UNIT_TESTS_TIDAL_POTENTIAL_EXPANSION_H
#define __UNIT_TESTS_TIDAL_POTENTIAL_EXPANSION_H

#include "EccentricOrbit.h"

#include "../../Evolve/TidalPotentialTerms.h"

#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <cassert>

namespace Evolve {
    class TidalPotentialExpansion {
    private:
        TidalPotentialTerms __expansion_coef;

        double
            __primary_mass,

            __secondary_mass,

            __semimajor,
            
            __eccentricity,

            __inclination,

            __arg_of_periapsis;

        double tidal_term(
            int mprime,
            
            double radial_distance,

            double azimuthal_angle,

            double polar_angle,

            double orbital_phase
        ) const;
    public:
        TidalPotentialExpansion(
            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
        ) :
            __primary_mass(primary_mass),
            __secondary_mass(secondary_mass),
            __semimajor(semimajor),
            __eccentricity(eccentricity),
            __inclination(inclination),
            __arg_of_periapsis(arg_of_periapsis)
        {}

        double evaluate_spherical_coords(
            double radial_distance,

            double azimuthal_angle,

            double polar_angle,

            double time
        );

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

                double time
            );

        void set_eccentricity_order(unsigned e_order)
        {__expansion_coef.change_e_order(e_order);}
    }; //End TidalPotentialExpansion class.

    template<class POSITION_TYPE>
        double TidalPotentialExpansion::operator()(
            const POSITION_TYPE &position,
            double time
        )
        {
            double radial_distance = position.norm(),
                   polar_angle = (
                       radial_distance == 0
                       ? 0
                       : std::acos(position[2]/radial_distance)
                   ),
                   azimuthal_angle = std::atan2(position[1], position[0]);

            if(azimuthal_angle < 0)
                azimuthal_angle += 2.0 * M_PI;

            return evaluate_spherical_coords(
                radial_distance,
                azimuthal_angle,
                polar_angle,
                time
            );
        }
} //End Evolve namespace.
#endif