Tudat/Astrodynamics/MissionSegments/gravityAssist.cpp
//! Calculate deltaV of a gravity assist.
Calculate deltaV of a gravity assist.
Calculate deltaV of a gravity assist .
double gravityAssist( const double centralBodyGravitationalParameter,
                      const Eigen::Vector3d& centralBodyVelocity,
                      const Eigen::Vector3d& incomingVelocity,
                      const Eigen::Vector3d& outgoingVelocity,
                      const double smallestPeriapsisDistance,
                      const bool useEccentricityInsteadOfPericenter,
                      const double speedTolerance,
                      RootFinderPointer rootFinder )
{
    using basic_mathematics::UnivariateProxyPointer;
    using basic_mathematics::UnivariateProxy;
    const Eigen::Vector3d incomingHyperbolicExcessVelocity
            = incomingVelocity - centralBodyVelocity;
    const Eigen::Vector3d outgoingHyperbolicExcessVelocity
            = outgoingVelocity - centralBodyVelocity;
    const double absoluteIncomingExcessVelocity = incomingHyperbolicExcessVelocity.norm( );
    const double absoluteOutgoingExcessVelocity = outgoingHyperbolicExcessVelocity.norm( );
    double bendingAngle = linear_algebra::computeAngleBetweenVectors(
                            incomingHyperbolicExcessVelocity, outgoingHyperbolicExcessVelocity );
    const double maximumBendingAngle =
            std::asin( 1.0 / ( 1.0 + ( smallestPeriapsisDistance *
                                       absoluteIncomingExcessVelocity *
                                       absoluteIncomingExcessVelocity /
                                       centralBodyGravitationalParameter ) ) ) +
            std::asin( 1.0 / ( 1.0 + ( smallestPeriapsisDistance *
                                       absoluteOutgoingExcessVelocity *
                                       absoluteOutgoingExcessVelocity /
                                       centralBodyGravitationalParameter ) ) );
    double bendingEffectDeltaV = 0.0;
    double velocityEffectDeltaV = 0.0;
    if ( bendingAngle > maximumBendingAngle )
    {
        const double extraBendingAngle = bendingAngle - maximumBendingAngle;
        bendingEffectDeltaV = 2.0 * std::min( absoluteIncomingExcessVelocity,
                                              absoluteOutgoingExcessVelocity ) *
                                    std::sin( extraBendingAngle / 2.0 );
        const double pericenterRadius = smallestPeriapsisDistance;
        const double incomingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter /
                                             absoluteIncomingExcessVelocity /
                                             absoluteIncomingExcessVelocity;
        const double outgoingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter /
                                             absoluteOutgoingExcessVelocity /
                                             absoluteOutgoingExcessVelocity;
        const double incomingEccentricity = 1 - pericenterRadius / incomingSemiMajorAxis;
        const double outgoingEccentricity = 1 - pericenterRadius / outgoingSemiMajorAxis;
        const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity *
                std::sqrt( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) );
        const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity *
                std::sqrt( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) );
        velocityEffectDeltaV = std::fabs( incomingVelocityAtPeriapsis -
                                          outgoingVelocityAtPeriapsis );
    }
    else if ( ( std::fabs( absoluteIncomingExcessVelocity - absoluteOutgoingExcessVelocity )
                <= speedTolerance ) )
    {
    }
    else if ( useEccentricityInsteadOfPericenter )
    {
        const double incomingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter /
                                             absoluteIncomingExcessVelocity /
                                             absoluteIncomingExcessVelocity;
        const double outgoingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter /
                                             absoluteOutgoingExcessVelocity /
                                             absoluteOutgoingExcessVelocity;
        EccentricityFindingFunctions eccentricityFindingFunctions( incomingSemiMajorAxis,
                                                                   outgoingSemiMajorAxis,
                                                                   bendingAngle );
        UnivariateProxyPointer rootFunction = std::make_shared< UnivariateProxy >(
                    std::bind( &EccentricityFindingFunctions::
                                 computeIncomingEccentricityFunction,
                                 eccentricityFindingFunctions, std::placeholders::_1 ) );
        rootFunction->addBinding( -1, std::bind(
                                      &EccentricityFindingFunctions::
                                      computeFirstDerivativeIncomingEccentricityFunction,
                                      eccentricityFindingFunctions, std::placeholders::_1 ) );
        double incomingEccentricity = TUDAT_NAN;
        if ( ( absoluteOutgoingExcessVelocity / absoluteIncomingExcessVelocity ) < 100.0 )
        {
            try
            {
                incomingEccentricity = rootFinder->execute( rootFunction, 1.0 + 1.0e-2 );
            }
            catch(std::runtime_error)
            {
                root_finders::RootFinderPointer rootFinder_temp
                  = std::make_shared< root_finders::Bisection >( 1.0e-12, 1000 ) ;
                incomingEccentricity = rootFinder_temp->execute( rootFunction, 1.0 + 1.0e-2 );
            }
        }
        else
        {
            try
            {
                incomingEccentricity = rootFinder->execute( rootFunction, 1.0 + 1.0e-10 );
            }
            catch(std::runtime_error)
            {
                root_finders::RootFinderPointer rootFinder_temp
                  = std::make_shared< root_finders::Bisection >( 1.0e-12, 1000 ) ;
                incomingEccentricity = rootFinder_temp->execute( rootFunction, 1.0 + 1.0e-10 );
            }
        }
        const double outgoingEccentricity = 1.0 - ( incomingSemiMajorAxis /
                                                    outgoingSemiMajorAxis ) *
                                            ( 1.0 - incomingEccentricity );
        const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity *
                std::sqrt( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) );
        const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity *
                std::sqrt( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) );
        velocityEffectDeltaV = std::fabs( incomingVelocityAtPeriapsis -
                                          outgoingVelocityAtPeriapsis );
    }
    else
    {
         This is the absolute semi-major axis, because
        const double absoluteIncomingSemiMajorAxis = 1.0 * centralBodyGravitationalParameter /
                                                     absoluteIncomingExcessVelocity /
                                                     absoluteIncomingExcessVelocity;
        const double absoluteOutgoingSemiMajorAxis = 1.0 * centralBodyGravitationalParameter /
                                                     absoluteOutgoingExcessVelocity /
                                                     absoluteOutgoingExcessVelocity;
        PericenterFindingFunctions pericenterFindingFunctions( absoluteIncomingSemiMajorAxis,
                                                               absoluteOutgoingSemiMajorAxis,
                                                               bendingAngle);
        UnivariateProxyPointer rootFunction = std::make_shared< UnivariateProxy >(
                    std::bind( &PericenterFindingFunctions::computePericenterRadiusFunction,
                                 pericenterFindingFunctions, std::placeholders::_1 ) );
        rootFunction->addBinding( -1, std::bind( &PericenterFindingFunctions::
                                                   computeFirstDerivativePericenterRadiusFunction,
                                                   pericenterFindingFunctions, std::placeholders::_1 ) );
        const double pericenterRadius = rootFinder->execute( rootFunction,
                                                             smallestPeriapsisDistance );
        const double incomingEccentricity = 1.0 + pericenterRadius / absoluteIncomingSemiMajorAxis;
        const double outgoingEccentricity = 1.0 + pericenterRadius / absoluteOutgoingSemiMajorAxis;
        const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity *
                std::sqrt( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) );
        const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity *
                std::sqrt( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) );
        velocityEffectDeltaV = std::fabs( incomingVelocityAtPeriapsis -
                                          outgoingVelocityAtPeriapsis );
    }
    return bendingEffectDeltaV + velocityEffectDeltaV;
}
double gravityAssist ( const double centralBodyGravitationalParameter , const Eigen :: Vector3d & centralBodyVelocity , const Eigen :: Vector3d & incomingVelocity , const Eigen :: Vector3d & outgoingVelocity , const double smallestPeriapsisDistance , const bool useEccentricityInsteadOfPericenter , const double speedTolerance , RootFinderPointer rootFinder ) { using basic_mathematics :: UnivariateProxyPointer ; using basic_mathematics :: UnivariateProxy ; const Eigen :: Vector3d incomingHyperbolicExcessVelocity = incomingVelocity - centralBodyVelocity ; const Eigen :: Vector3d outgoingHyperbolicExcessVelocity = outgoingVelocity - centralBodyVelocity ; const double absoluteIncomingExcessVelocity = incomingHyperbolicExcessVelocity . norm ( ) ; const double absoluteOutgoingExcessVelocity = outgoingHyperbolicExcessVelocity . norm ( ) ; double bendingAngle = linear_algebra :: computeAngleBetweenVectors ( incomingHyperbolicExcessVelocity , outgoingHyperbolicExcessVelocity ) ; const double maximumBendingAngle = std :: asin ( 1.0 / ( 1.0 + ( smallestPeriapsisDistance * absoluteIncomingExcessVelocity * absoluteIncomingExcessVelocity / centralBodyGravitationalParameter ) ) ) + std :: asin ( 1.0 / ( 1.0 + ( smallestPeriapsisDistance * absoluteOutgoingExcessVelocity * absoluteOutgoingExcessVelocity / centralBodyGravitationalParameter ) ) ) ; double bendingEffectDeltaV = 0.0 ; double velocityEffectDeltaV = 0.0 ; if ( bendingAngle > maximumBendingAngle ) { const double extraBendingAngle = bendingAngle - maximumBendingAngle ; bendingEffectDeltaV = 2.0 * std :: min ( absoluteIncomingExcessVelocity , absoluteOutgoingExcessVelocity ) * std :: sin ( extraBendingAngle / 2.0 ) ; const double pericenterRadius = smallestPeriapsisDistance ; const double incomingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter / absoluteIncomingExcessVelocity / absoluteIncomingExcessVelocity ; const double outgoingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter / absoluteOutgoingExcessVelocity / absoluteOutgoingExcessVelocity ; const double incomingEccentricity = 1 - pericenterRadius / incomingSemiMajorAxis ; const double outgoingEccentricity = 1 - pericenterRadius / outgoingSemiMajorAxis ; const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity * std :: sqrt ( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) ) ; const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity * std :: sqrt ( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) ) ; velocityEffectDeltaV = std :: fabs ( incomingVelocityAtPeriapsis - outgoingVelocityAtPeriapsis ) ; } else if ( ( std :: fabs ( absoluteIncomingExcessVelocity - absoluteOutgoingExcessVelocity ) <= speedTolerance ) ) { } else if ( useEccentricityInsteadOfPericenter ) { const double incomingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter / absoluteIncomingExcessVelocity / absoluteIncomingExcessVelocity ; const double outgoingSemiMajorAxis = -1.0 * centralBodyGravitationalParameter / absoluteOutgoingExcessVelocity / absoluteOutgoingExcessVelocity ; EccentricityFindingFunctions eccentricityFindingFunctions ( incomingSemiMajorAxis , outgoingSemiMajorAxis , bendingAngle ) ; UnivariateProxyPointer rootFunction = std :: make_shared < UnivariateProxy > ( std :: bind ( & EccentricityFindingFunctions :: computeIncomingEccentricityFunction , eccentricityFindingFunctions , std :: placeholders :: _1 ) ) ; rootFunction -> addBinding ( -1 , std :: bind ( & EccentricityFindingFunctions :: computeFirstDerivativeIncomingEccentricityFunction , eccentricityFindingFunctions , std :: placeholders :: _1 ) ) ; double incomingEccentricity = TUDAT_NAN ; if ( ( absoluteOutgoingExcessVelocity / absoluteIncomingExcessVelocity ) < 100.0 ) { try { incomingEccentricity = rootFinder -> execute ( rootFunction , 1.0 + 1.0e-2 ) ; } catch ( std :: runtime_error ) { root_finders :: RootFinderPointer rootFinder_temp = std :: make_shared < root_finders :: Bisection > ( 1.0e-12 , 1000 ) ; incomingEccentricity = rootFinder_temp -> execute ( rootFunction , 1.0 + 1.0e-2 ) ; } } else { try { incomingEccentricity = rootFinder -> execute ( rootFunction , 1.0 + 1.0e-10 ) ; } catch ( std :: runtime_error ) { root_finders :: RootFinderPointer rootFinder_temp = std :: make_shared < root_finders :: Bisection > ( 1.0e-12 , 1000 ) ; incomingEccentricity = rootFinder_temp -> execute ( rootFunction , 1.0 + 1.0e-10 ) ; } } const double outgoingEccentricity = 1.0 - ( incomingSemiMajorAxis / outgoingSemiMajorAxis ) * ( 1.0 - incomingEccentricity ) ; const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity * std :: sqrt ( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) ) ; const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity * std :: sqrt ( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) ) ; velocityEffectDeltaV = std :: fabs ( incomingVelocityAtPeriapsis - outgoingVelocityAtPeriapsis ) ; } else { const double absoluteIncomingSemiMajorAxis = 1.0 * centralBodyGravitationalParameter / absoluteIncomingExcessVelocity / absoluteIncomingExcessVelocity ; const double absoluteOutgoingSemiMajorAxis = 1.0 * centralBodyGravitationalParameter / absoluteOutgoingExcessVelocity / absoluteOutgoingExcessVelocity ; PericenterFindingFunctions pericenterFindingFunctions ( absoluteIncomingSemiMajorAxis , absoluteOutgoingSemiMajorAxis , bendingAngle ) ; UnivariateProxyPointer rootFunction = std :: make_shared < UnivariateProxy > ( std :: bind ( & PericenterFindingFunctions :: computePericenterRadiusFunction , pericenterFindingFunctions , std :: placeholders :: _1 ) ) ; rootFunction -> addBinding ( -1 , std :: bind ( & PericenterFindingFunctions :: computeFirstDerivativePericenterRadiusFunction , pericenterFindingFunctions , std :: placeholders :: _1 ) ) ; const double pericenterRadius = rootFinder -> execute ( rootFunction , smallestPeriapsisDistance ) ; const double incomingEccentricity = 1.0 + pericenterRadius / absoluteIncomingSemiMajorAxis ; const double outgoingEccentricity = 1.0 + pericenterRadius / absoluteOutgoingSemiMajorAxis ; const double incomingVelocityAtPeriapsis = absoluteIncomingExcessVelocity * std :: sqrt ( ( incomingEccentricity + 1.0 ) / ( incomingEccentricity - 1.0 ) ) ; const double outgoingVelocityAtPeriapsis = absoluteOutgoingExcessVelocity * std :: sqrt ( ( outgoingEccentricity + 1.0 ) / ( outgoingEccentricity - 1.0 ) ) ; velocityEffectDeltaV = std :: fabs ( incomingVelocityAtPeriapsis - outgoingVelocityAtPeriapsis ) ; } return bendingEffectDeltaV + velocityEffectDeltaV ; }
Calculate deltaV of a gravity assist.
Calculate deltaV of a gravity assist .
// Compute incoming and outgoing hyperbolic excess velocity. // Compute absolute values of the hyperbolic excess velocities. // Compute bending angle. // Compute maximum achievable bending angle. // Initialize bending effect deltaV, which is zero, unless extra bending angle is required. // Initialize velocity effect delta V parameter. // Check if an additional bending angle is required. If so, the additional bending angle // maneuver has to be performed. Also the pericenter radius will be the minimum pericenter // radius to obtain the largest possible bending angle 'for free'. Hence no root finding is // required for this case. As noted above, this may not be ideal for all cases. (for cases // in which the excess velocities are relatively small) // Compute required extra bending angle that cannot be delivered by an unpowered swing-by. // Compute necessary delta-V due to bending-effect. // This means the pericenter radius is now equal to the smallest pericenter radius, to // ensure the largest possible bending angle. // Compute semi-major axis of hyperbolic legs. // Compute incoming hyperbolic leg eccentricity. // Compute outgoing hyperbolic leg eccentricity. // Compute incoming and outgoing velocities at periapsis. // Compute necessary delta-V due to velocity-effect. // In this case no maneuver has to be performed. Hence no iteration is performed, and the // delta V is simply kept at 0.0. // Here the required maneuver to patch the incoming and outgoing excess velocities is // calculated. In this implementation, the eccentricity will be used as iteration parameter. // Compute semi-major axis of hyperbolic legs. // Set the gravity assist function with the variables to perform root finder calculations. // Create an object containing the function of which we whish to obtain the root from. // Add the first derivative of the root function. // Initialize incoming eccentricity. // Set initial guess of the variable computed in Newton-Rapshon method. // In these cases the very low estimate (which is given under else) may in some cases // result in no convergence. Hence a higher value of 1.01 is necessary. This will not // result in 'going through' 1.0 as mentioned below, because the eccentricity in these // cases is always high! // This is set to a value that is close to 1.0. This is more robust than higher values, // because for those higher values Newton Raphson sometimes 'goes through' 1.0. This // results in NaN values for the derivative of the eccentricity finding function. // Compute outgoing hyperbolic leg eccentricity. // Compute incoming and outgoing velocities at periapsis. // Compute necessary delta-V due to velocity-effect. // Here the required maneuver to patch the incoming and outgoing excess velocities is // calculated. In this implementation, the pericenter radius will be used as iteration // parameter. // Compute semi-major axis of hyperbolic legs. This is the absolute semi-major axis, because // it will otherwisely result in the root of a negative function for various cases during // the rootfinding process. // Set the gravity assist function with the variables to perform root finder calculations. // Create an object containing the function of which we whish to obtain the root from. // Add the first derivative of the root function. // Set pericenter radius based on result of Newton-Raphson root-finding algorithm. // Compute incoming hyperbolic leg eccentricity. // Compute outgoing hyperbolic leg eccentricity. // Compute incoming and outgoing velocities at periapsis. // Compute necessary delta-V due to velocity-effect. // Compute and return the total delta-V.