CR3BPDynamics.hpp

/******************************************************************************
* Copyright (c) ATTX LLC 2024. All Rights Reserved.
*
* This software and associated documentation (the "Software") are the 
* proprietary and confidential information of ATTX, LLC. The Software is 
* furnished under a license agreement between ATTX and the user organization 
* and may be used or copied only in accordance with the terms of the agreement.
* Refer to 'license/attx_license.adoc' for standard license terms.
*
* EXPORT CONTROL NOTICE: THIS SOFTWARE MAY INCLUDE CONTENT CONTROLLED UNDER THE
* INTERNATIONAL TRAFFIC IN ARMS REGULATIONS (ITAR) OR THE EXPORT ADMINISTRATION 
* REGULATIONS (EAR99). No part of the Software may be used, reproduced, or 
* transmitted in any form or by any means, for any purpose, without the express 
* written permission of ATTX, LLC.
******************************************************************************/
/*
Circular Restricted Three Body Dynamics header file
---------------------------------------------------
 
Author: Alex Reynolds
*/
#ifndef CR3BPUTILS_CR3BP_DYNAMICS_HPP
#define CR3BPUTILS_CR3BP_DYNAMICS_HPP
 
#include "utils/Rates.hpp"
#include "core/CartesianVector.hpp"
#include "core/Matrix.hpp"
 
/**
 * @brief   Class to propagate CR3BP dynamics in characteristic units
 * 
 * This class is used to propagate circular restricted 3 body problem
 * dynamics in characteristic units. It sets internal parameters on
 * the basis of the setters for mu and l, which are used in the propagation.
 */
namespace modelspace {
    const unsigned int NUM_CR3BP_STATES = 6;
 
    class CR3BPDynamics : public clockwerk::Rates<double, NUM_CR3BP_STATES> {
    public:
        CR3BPDynamics() : Rates() {}
 
        /// @brief   Function to calculate frame rates of change for a frame
        /// @param   time    The time input for the calculation. Uses current time
        /// @param   state   The current state vector in the calculation
        /// @param   out_rates   Implicit return of rates of change in states
        int calculateRates(double time,
                           const std::array<double, NUM_CR3BP_STATES> &state,
                           std::array<double, NUM_CR3BP_STATES> &out_rates) {
            // Pre-calculate some parameters to make life easier down the line
            double one_minus_mu_star = 1.0 - mu_star;
            double x_plus_mu_star = state[0] + mu_star;
            double x_plus_one_minus_mu_star = state[0] - 1.0 + mu_star;
            double y_squared = state[1]*state[1];
            double z_squared = state[2]*state[2];
            double r1 = std::sqrt(x_plus_mu_star*x_plus_mu_star + y_squared + z_squared);
            double r2 = std::sqrt(x_plus_one_minus_mu_star*x_plus_one_minus_mu_star + y_squared + z_squared);
            double one_over_r1_cubed = 1.0/(r1*r1*r1);
            double one_over_r2_cubed = 1.0/(r2*r2*r2);
 
            // Now calculate our state vector rate of change
            out_rates[0] = state[3];
            out_rates[1] = state[4];
            out_rates[2] = state[5];
            out_rates[3] = -one_minus_mu_star*x_plus_mu_star*one_over_r1_cubed - mu_star*(state[0] - one_minus_mu_star)*one_over_r2_cubed + 2*state[4] + state[0];
            out_rates[4] = -one_minus_mu_star*state[1]*one_over_r1_cubed - mu_star*state[1]*one_over_r2_cubed - 2.0*state[3] + state[1];
            out_rates[5] = -one_minus_mu_star*state[2]*one_over_r1_cubed - mu_star*state[2]*one_over_r2_cubed;
 
            return NO_ERROR;
        }
 
        /// @brief This is the ratio of second body mass to total system mass
        double mu_star;
    private:
 
    };
 
}
 
#endif