[#1325] Added relative velocity atmo in cannonballDrag

[carlo98]

• Tickets addressed: Issue Add atmosphere-relative velocity and density correction support to Drag effectors #1325
• Review: By commit
• Merge strategy: Merge (no squash)

Description

Adds optional atmosphere-relative velocity support to cannonballDrag, allowing drag to be computed against the rotating atmosphere rather than the inertial frame.

Key implementation details reviewers should be aware of:

• Default behavior is unchanged: The new feature is opt-in. Existing users keep the current inertial-velocity behavior by default (useAtmosphereRelativeVelocity = false).

• Relative-velocity formulation: When enabled, the module computes v_rel = v_sc - (planetOmega_N x r_sc) before transforming into the site frame and evaluating the cannonball drag force and torque.

• Earth default, configurable for other bodies: planetOmega_N defaults to [0, 0, OMEGA_EARTH], so Earth cases work without extra configuration. Non-Earth cases can override this vector from Python via setPlanetOmega_N.

• Getter and setter methods: setUseAtmosphereRelativeVelocity / getUseAtmosphereRelativeVelocity and setPlanetOmega_N / getPlanetOmega_N are exposed to Python via SWIG.

Verification

• Added test_cannonballDragAtmoRelVelDisabled: confirms drag is unchanged when useAtmosphereRelativeVelocity = False and matches the analytical inertial-velocity cannonball drag formula.

• Added test_cannonballDragAtmoRelVelEnabled: confirms the computed force matches a reference when useAtmosphereRelativeVelocity = True, verifying that omega x r is correctly subtracted from the inertial velocity before evaluating drag.

• Extended the integration test test_orbit: parametrized over useRelVel = [False, True], orbitCase = ["LPO", "LTO"], and planetCase = ["earth", "mars"], comparing the module output to an independent reference implementation at every simulation step over a full orbital period.

Documentation

• Updated cannonballDrag.rst with the atmosphere-relative velocity formulation, including the v_rel equation, Earth-default note, non-Earth configuration note, and description of the detailed per-step behavior.

Future work

• Source planetOmega_N automatically from Spice kernels.

[juan-g-bonilla]

Thanks Carlo for the PR. I think it's a great step to undo the assumption of a stationary atmosphere. However, I think we can take this idea and make it more generic to support arbitrary wind models relatively easily.

The main functional change in this pr is:

const auto v_S = dcm_SN*Eigen::Map<const Eigen::Vector3d>(siteStatePayload.v_BN_N);

is now (paraphrasing a bit):

Eigen::Vector3d vRel_N = Eigen::Map<const Eigen::Vector3d>(siteStatePayload.v_BN_N);
if (this->useAtmosphereRelativeVelocity) {
    const Eigen::Map<const Eigen::Vector3d> r_N(siteStatePayload.r_BN_N);
    auto vWind_N = this->planetOmega_N.cross(r_N);
    vRel_N -= vWind_N;
}
const auto v_S = dcm_SN * vRel_N;

I think we could improve this module further by not baking in the assumption that we want a co-rotating atmosphere model of the local wind, and instead take that local wind vWind_N as an input message to the model (similarly to how the atmospheric properties are taken as an input to the model).

If the model takes in an arbitrary vWind_N through a message:

typedef struct {
    double wind_N[3];   //!< inertial wind velocity [m/s]
} WindMsgPayload;

then we could add a new, very simple model CorotatingAtmosphereWindModel that computes the wind as you describe in this PR:

#ifndef COROTATING_ATMOSPHERE_WIND_MODEL_H
#define COROTATING_ATMOSPHERE_WIND_MODEL_H

#include <Eigen/Dense>
#include "architecture/_GeneralModuleFiles/sys_model.h"
#include "architecture/messaging/messaging.h"
#include "architecture/msgPayloadDefC/WindMsgPayload.h"
#include "architecture/msgPayloadDefC/SCStatesMsgPayload.h"

/*!
 * @brief Computes atmosphere velocity assuming rigid corotation with planet.
 *
 * wind_N = omegaPlanet_N × r_BN_N
 */
class CorotatingAtmosphereWindModel : public SysModel {
public:
    CorotatingAtmosphereWindModel();
    ~CorotatingAtmosphereWindModel() = default;

    void Reset(uint64_t CurrentSimNanos) override;
    void UpdateState(uint64_t CurrentSimNanos) override;

    /** Output wind message */
    Message<WindMsgPayload> windOutMsg;

    /** Input spacecraft state message */
    ReadFunctor<SCStatesMsgPayload> scStateInMsg;

    /**
     * @brief Set the planetary rotation vector in the inertial frame [rad/s].
     * Defaults to Earth's rotation rate. Used when useAtmosphereRelativeVelocity is enabled.
     */
    void setPlanetOmega_N(const Eigen::Vector3d& omega);

    /** @brief Returns the planetary rotation vector in the inertial frame [rad/s]. */
    Eigen::Vector3d getPlanetOmega_N() const;

private:
    Eigen::Vector3d omegaPlanet_N;
};

#endif

#include "simulation/environment/CorotatingAtmosphereWindModel/CorotatingAtmosphereWindModel.h"

CorotatingAtmosphereWindModel::CorotatingAtmosphereWindModel()
{
    this->ModelTag = "CorotatingAtmosphereWindModel";

    // Default Earth rotation rate
    double earthRate = 7.2921159e-5;   // rad/s
    this->omegaPlanet_N = Eigen::Vector3d(0.0, 0.0, earthRate);
}

void CorotatingAtmosphereWindModel::Reset(uint64_t CurrentSimNanos)
{
    if (!this->scStateInMsg.isLinked()) {
        bskLogger.bskLog(BSK_ERROR, "CorotatingAtmosphereWindModel.scStateInMsg was not connected.");
    }
}

void CorotatingAtmosphereWindModel::UpdateState(uint64_t CurrentSimNanos)
{
    SCStatesMsgPayload scState = this->scStateInMsg();

    Eigen::Vector3d r_BN_N(
        scState.r_BN_N[0],
        scState.r_BN_N[1],
        scState.r_BN_N[2]
    );

    Eigen::Vector3d wind = this->omegaPlanet_N.cross(r_BN_N);

    WindMsgPayload windMsg;
    windMsg.wind_N[0] = wind[0];
    windMsg.wind_N[1] = wind[1];
    windMsg.wind_N[2] = wind[2];

    this->windOutMsg.write(&windMsg, this->moduleID, CurrentSimNanos);
}

void CorotatingAtmosphereWindModel::setPlanetOmega_N(const Eigen::Vector3d& omega)
{
    this->omegaPlanet_N = omega;
}

Eigen::Vector3d CorotatingAtmosphereWindModel::getPlanetOmega_N() const
{
    return this->omegaPlanet_N;
}

The clear benefit is that this opens the door for arbitrary wind models to work in our simulation.

NOTE: take the code above with a grain of salt, just a quick sketch of the final product. Feel free to improve as you see fit. For example, in keeping with our environmental models pattern, I would add a windBase base sys model that this specific model inherits from.

[schaubh]

Interesting idea @juan-g-bonilla , if we do all this math, and we have several modules that could use this info, I wonder if we should have a separate module that computes the relative wind that can be fed to drag cannon ball model, faceted model, etc? This way we are not duplicating all this functionality in each module?

[juan-g-bonilla]

Interesting idea @juan-g-bonilla , if we do all this math, and we have several modules that could use this info, I wonder if we should have a separate module that computes the relative wind that can be fed to drag cannon ball model, faceted model, etc? This way we are not duplicating all this functionality in each module?

Exactly. We should have a set of modules that compute the wind at some location in the atmosphere (like we have several modules that compute the density at some location) that all feed a wind vector to the drag models (cannonBall, facet, etc.). It would actually reduce repetition of math in PR #1326 and #1327 .

[carlo98]

Ok, so I can implement that new simple module (CorotatingAtmosphereWindModel) and windBase in a new PR and then update this one and #1326 . I'm going to convert #1326 and #1327 into drafts while working on the wind module. Does it sound good? @juan-g-bonilla @schaubh

Also, to update dragDynamicEffector do you want me to open another PR or use one of these?

[juan-g-bonilla]

Ok, so I can implement that new simple module (CorotatingAtmosphereWindModel) and windBase in a new PR and then update this one and #1326 . I'm going to convert #1326 and #1327 into drafts while working on the wind module. Does it sound good? @juan-g-bonilla @schaubh

Also, to update dragDynamicEffector do you want me to open another PR or use one of these?

Sounds good to me, thanks for the push!