SphericalHarmonicsGravityModel.h

/******************************************************************************
* Copyright (c) ATTX INC 2025. All Rights Reserved.
*
* This software and associated documentation (the "Software") are the 
* proprietary and confidential information of ATTX INC. 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 INC.
******************************************************************************/
/*
Solar Radiation Pressure model header file

Author: Sam Matez
*/
/*
Metadata for MS GUI:
imdata = {"displayname" : "Spherical Harmonic Gravity",
          "exclude" : False,
          "category" : "Dynamics"
}
aliases = {"mu" : "GM",
           "filename" : "Gravity File",
           "r_planet" : "Planet Equatorial Radius",
           "body_mass" : "Spacecraft Mass",
           "n" : "Degree",
           "m" : "Order",
           "pos_body_planet__pcpf" : "Position",
           "grav_force__f" : "Gravity Force"
}
*/

#ifndef MODELS_SPHERICAL_HARMONICS_GRAVITY_MODEL_H 
#define MODELS_SPHERICAL_HARMONICS_GRAVITY_MODEL_H

#include <string>

#include "simulation/Model.h"
#include "utils/sphericalharmonicsutils.h"
#include "locations.h"
#include "constants/planetdefaults.h"

namespace warptwin {

    /**
     * @brief   Spherical Harmonics Model
     * 
     * Implements spherical harmonic gravity up to an 18x18 field based on the implementation
     * in the 42 simulation: 
     * https://github.com/ericstoneking/42/blob/4be580d8defc968da97d937ce09ef107c07860cb/Kit/Source/envkit.c#L82
     * 
     * Author: Sam Matez
     * Email: sam.matez@attx.tech
     * 
    */     
    MODEL(SphericalHarmonicsGravityModel)
    public:
        // Model params
        //         NAME                     TYPE                    DEFAULT VALUE
        START_PARAMS
            /** The filename for the file that contains the coefficients for spherical
             *  harmonic calculation. These coefficients should be un-normalized, as the 
             *  model applies a Neumann normalization to the loaded coefficients. */
            SIGNAL(filename,                std::string,            warptwinDir() + "data/WGS84/WGS84.txt")
            /** This is the gravitational parameter of our parent planet. Defaults to
             * Earth's gravitational parameter for ease of use */
            SIGNAL(mu,                      double,                 warpos::earth_wgs84.mu)
            /** This is the radius of our parent planet. Defaults to Earth's radius, in meters for ease of use */
            SIGNAL(r_planet,                double,                 warpos::earth_wgs84.eq_radius) 
            /** The coefficient n which represents the variations in gravitational field in the latitudinal direction */
            SIGNAL(n,                       int,                    0)
            /** The coefficient m which represents the variations in gravitational field in the longitudinal direction.
             *  must be less than or equal to n */
            SIGNAL(m,                       int,                    0)
            /** This is the mass of the body, which will be multiplied by our force
             *  to get the force calculation... TODO to replace this with accel methods */
            SIGNAL(body_mass,               double,                 1.0)
        END_PARAMS

        // Model inputs
        //         NAME                     TYPE                    DEFAULT VALUE
        START_INPUTS
            /** The body position wrt the planet in a planet-centered, planet fixed frame (m) */
            SIGNAL(pos_body_planet__pcpf,   CartesianVector3,      CartesianVector3({0.0, 0.0, 0.0}))
        END_INPUTS

        // Model Outputs
        //         NAME                     TYPE                    DEFAULT VALUE
        START_OUTPUTS
            /** This is the acceleration due to gravity calculated in the same reference frame as pos_body__f */
            SIGNAL(grav_accel__f,           CartesianVector3,      CartesianVector3({0.0, 0.0, 0.0}))
            /** This is the force due to gravity calculated in the same reference frame as pos_body__f */
            SIGNAL(grav_force__f,           CartesianVector3,      CartesianVector3({0.0, 0.0, 0.0}))
        END_OUTPUTS

    protected:
        int16 start() override;
        int16 execute() override; 

        // Pre-calculated parameters
        double _K;

        // Variables to hold gravity coefficients
        double _c[NMAX_SPHERICAL_HARMONICS+1][NMAX_SPHERICAL_HARMONICS+1];
        double _s[NMAX_SPHERICAL_HARMONICS+1][NMAX_SPHERICAL_HARMONICS+1];

        // Local variable to hold force in the spherical frame
        CartesianVector3 _accel_grav__spherical;

        // Local variables for lla coordinates
        double _lat;
        double _lon;
        double _alt;

        // Local variable for magnitude of position vector
        double _r_pos;

        // Local variable to hold colatitude
        double _colatitude;

        // Local variables for rotation back to pcpf
        double _sth;
        double _cth;
        double _sph;
        double _cph;
        CartesianVector3 _fe;

    };

}

#endif