orbitUtils.hpp

Go to the documentation of this file.

/** \file orbitUtils.hpp
 * \author Jared R. Males
 * \brief Utilities related to orbits.
 * \ingroup astrofiles
 *
 */
 
#ifndef mx_astro_orbitUtils_hpp
#define mx_astro_orbitUtils_hpp
 
#include "../math/constants.hpp"
#include "kepler.hpp"
 
namespace mx
{
namespace astro
{
 
/// Calculate the period of an orbit, given the masses and semi-major axis.
/**
 * \tparam units specifies the unit system and precision, see \ref astrounits.
 *
 * \returns the period in the specified units
 *
 * \ingroup orbits
 */
template <typename units>
typename units::realT orbitPeriod( typename units::realT m1, ///< [in] mass of one object.
                                   typename units::realT m2, ///< [in] mass of the other object (can be 0).
                                   typename units::realT a   ///< [in] the semi-major axis of the orbit.
)
{
    return math::two_pi<typename units::realT>() * sqrt( ( a * a * a ) / ( constants::G<units>() * ( m1 + m2 ) ) );
}
 
/// Calculate the semi-major axis of an orbit, given the masses and Period, using SI units.
/**
 * \tparam units specifies the unit system and precision, see \ref astrounits.
 *
 * \returns the semi-major axis in the specified units
 *
 * \ingroup orbits
 */
template <typename units>
typename units::realT orbitSemiMaj( typename units::realT m1, ///< [in] mass of one object.
                                    typename units::realT m2, ///< [in] mass of the other object (can be 0).
                                    typename units::realT P   ///< [in] the period of the orbit.
)
{
    typedef typename units::realT realT;
 
    constexpr realT two = static_cast<realT>( 2 ); // Just for convenience.
 
    return pow( pow( P / ( math::two_pi<realT>() ), two ) * constants::G<units>() * ( m1 + m2 ), math::third<realT>() );
}
 
/// Calculate the mean anomaly at time t, given the time of pericenter passage t0 and period P.
/**
 * \ingroup orbits
 */
template <typename realT>
realT orbitMeanAnol( realT t, realT t0, realT P )
{
    return math::two_pi<realT>() * ( t - t0 ) / P;
}
 
/// Calculate the separation and true anomaly for an orbit at the specified mean anomaly.
/**
 * \tparam realT is the type used for arithmetic
 *
 * \returns the numer of iterations in solving Kepler's equation, if successful.
 * \returns -1 if itmax is exceeded.
 *
 * \ingroup orbits
 */
template <typename realT>
long orbitElements( realT &r,                 ///< [out] is the orbital separation
                    realT &f,                 ///< [out] is the true anomaly
                    realT E,                  ///< [out] is the eccentric anomaly, provided since it is free
                    realT D,                  ///< [out] is the error after the last iteration of the Kepler solution.
                    realT e,                  ///< [in] is the eccentricity
                    realT M,                  ///< [in] is the mean anomaly
                    realT a,                  ///< [in] is the semi-major axis
                    realT tol = KEPLER_TOL,   ///< [in] is the desired tolerance
                    long itmax = KEPLER_ITMAX ///< [in] is the maximum number of iterations
)
{
    long its;
 
    if( e < 0.0 )
    {
        std::cerr << "e < 0 in rf_elements\n";
        return -1;
    }
 
    if( e == 1.0 )
    {
        std::cerr << "e = 1 in rf_elements, which is not currently handled.\n";
        return -1;
    }
 
    its = solve_kepler( E, D, e, M, tol, itmax );
 
    if( e > 1.0 )
    {
        f = 2.0 * std::atan( std::sqrt( ( e + 1.0 ) / ( e - 1.0 ) ) * std::tanh( E / 2.0 ) );
        r = a * ( 1.0 - e * std::cosh( E ) );
    }
    else if( e < 1.0 )
    {
        if( e == 0.0 )
        {
            f = M;
            r = a;
        }
        else
        {
            f = 2.0 * std::atan( std::sqrt( ( 1.0 + e ) / ( 1.0 - e ) ) * std::tan( E / 2.0 ) );
            r = a * ( 1.0 - e * std::cos( E ) );
        }
    }
 
    return its;
}
 
/// Get the orbital phase at true anomaly f. Calculates the cos(alfa) where alfa is the orbital phase.
/** Only calculates cos(alfa) to save an operation when possible.
 *
 * \tparam realT is the type used for arithmetic
 *
 * \returns the cosine of the phase angle
 *
 * \ingroup orbits
 */
template <typename realT>
realT orbitPhaseCosine( realT f,  ///< [in] is the true anomoaly at which to calculate alfa
                        realT w,  ///< [in] the argument of pericenter of the orbit
                        realT inc ///< [in] is the inclinatin of the orbit
)
{
    return sin( f + w ) * sin( inc );
}
 
/// Get the lambertian phase function at orbital phase specified as cos(alfa), where alfa is the phase angle.
/** Uses cos(alfa) to save an operation after calculating orbit phase.
 *
 * \tparam realT is the type used for arithmetic
 *
 * \returns the lambert phase function at alpha
 *
 * \ingroup orbits
 */
template <typename realT>
realT orbitLambertPhase( realT cos_alf /**< [in] the cosine of the phase angle*/ )
{
    realT alf = std::acos( cos_alf );
    return ( sin( alf ) + ( math::pi<realT>() - alf ) * cos_alf ) / math::pi<realT>();
}
 
/// Calculate the radial velocity of an orbiting body inclined relative to an observer.
/** \returns the radial velocity of the body with mass m1 as it orbits the barycenter of its orbit w.r.t. m2.
 *
 * \ingroup orbits
 */
template <typename units>
typename units::realT orbitRV( typename units::realT m1,  ///< [in] the mass of the body for which the RV is computed
                               typename units::realT m2,  ///< [in] the mass of the other body in the orbiting pair
                               typename units::realT inc, ///< [in] the inclination of the orbit of m2.
                               typename units::realT a,   ///< [in] the semi-major axis of the orbit
                               typename units::realT e,   ///< [in] the eccentricity of the orbit
                               typename units::realT w,   ///< [in] the argument of pericenter
                               typename units::realT f ///< [in] the true anomaly at which the RV is to be calculated.
)
{
    typedef typename units::realT realT;
 
    realT msini = m2 * sin( inc );
 
    realT rv = (msini)*std::sqrt( ( constants::G<units>() / ( m1 + m2 ) ) / ( a * ( 1 - e * e ) ) );
 
    rv *= ( std::cos( w + f ) + e * std::cos( w ) );
 
    return rv;
}
 
/// Calculate various quantities of an orbit given the keplerian elements and a vector of times.
/** Only those quantities with non-null pointers are actually calculated.
 *
 * The orbital parameters with dimensions, (a, P, and t0) must be in consistent units.
 *
 * \tparam vectorT a type whose elements are accessed with the [] operator.  No other requirements are placed on this
 * type, so it could be a native pointer (i.e. double *). \tparam realT the type in which to perform calculations.  Does
 * not need to match the storage type of vectorT.
 *
 * \retval -1 on error calculating the orbit (from rf_elements)
 * \retval 0 on success.
 *
 * \ingroup orbits
 */
template <typename vectorT, typename realT>
int orbitCartesianWork(
    vectorT *x, ///< [out] [optional] If not NULL, will be the projected x positions of the orbit.  Must be at least as
                ///< long as t.
    vectorT *y, ///< [out] [optional] If not NULL, will be the projected y positions of the orbit.  Must be at least as
                ///< long as t.
    vectorT *z, ///< [out] [optional] If not NULL, will be the projected z positions of the orbit.  Must be at least as
                ///< long as t.
    vectorT *r, ///< [out] [optional] If not NULL, will be the separation of the orbit.  Must be at least as long as t.
    vectorT *rProj, ///< [out] [optional] If not NULL, will be the projected separation of the orbit.  Must be at least
                    ///< as long as t.
    vectorT
        *f, ///< [out] [optional] If not NULL, will be the true anomaly of the orbit.  Must be at least as long as t.
    vectorT *cos_alf, ///< [out] [optional] If not NULL, will be the phase angle of the orbit.  Must be at least as long
                      ///< as t.
    vectorT
        *phi, ///< [out] [optional] If not NULL, will be the Lambertian phase function.  Must be at least as long as t.
    const vectorT &t, ///< [in] the times at which to calculate the projected positions
    const size_t N,   ///< [in] the number of points contained in t, and allocated in nx, ny, nz.  Avoids requiring a
                      ///< size() member (i.e. if vectorT is native pointers).
    realT a,          ///< [in] the semi-major axis of the orbit.
    realT P,          ///< [in] the orbital period.
    realT e,          ///< [in] the eccentricity of the orbit.
    realT t0,         ///< [in] the time of pericenter passage of the orbit.
    realT i,          ///< [in] the inclination of the orbit [radians].
    realT w,          ///< [in] the argument of pericenter of the orbit [radians].
    realT W           ///< [in] the longitude of the ascending node of the orbit [radians].
)
{
    long rv;
    realT _M, _r, _f, _E, _D;
 
    // Just do these calcs once
    realT cos_i = cos( i );
    realT sin_i = sin( i );
    realT cos_W = cos( W );
    realT sin_W = sin( W );
 
    realT cos_wf, sin_wf;
 
    for( size_t j = 0; j < N; j++ )
    {
        _M = orbitMeanAnol( t[j], t0, P );
 
        rv = orbitElements( _r, _f, _E, _D, e, _M, a );
 
        if( rv < 0 )
            return -1;
 
        // one calc
        cos_wf = cos( w + _f );
        sin_wf = sin( w + _f );
 
        // Assumption is that these will be optimized away if NULL
        if( x )
            ( *x )[j] = _r * ( cos_W * cos_wf - sin_W * sin_wf * cos_i );
        if( y )
            ( *y )[j] = _r * ( sin_W * cos_wf + cos_W * sin_wf * cos_i );
        if( z )
            ( *z )[j] = _r * sin_wf * sin_i;
        if( r )
            ( *r )[j] = _r;
        if( rProj )
        {
            if( x && y )
            {
                ( *rProj )[j] = sqrt( ( *x )[j] * ( *x )[j] + ( *y )[j] * ( *y )[j] );
            }
            else if( x && !y )
            {
                realT _y = _r * ( sin_W * cos_wf + cos_W * sin_wf * cos_i );
                ( *rProj )[j] = sqrt( ( *x )[j] * ( *x )[j] + _y * _y );
            }
            else if( !x && y )
            {
                realT _x = _r * ( cos_W * cos_wf - sin_W * sin_wf * cos_i );
                ( *rProj )[j] = sqrt( _x * _x + ( *y )[j] * ( *y )[j] );
            }
            else
            {
                realT _x = _r * ( cos_W * cos_wf - sin_W * sin_wf * cos_i );
                realT _y = _r * ( sin_W * cos_wf + cos_W * sin_wf * cos_i );
                ( *rProj )[j] = sqrt( _x * _x + _y * _y );
            }
        }
 
        if( f )
            ( *f )[j] = _f;
 
        if( cos_alf )
            ( *cos_alf )[j] = sin_wf * sin_i;
 
        if( phi )
            ( *phi )[j] = orbitLambertPhase<realT>( sin_wf * sin_i );
    }
    return 0;
}
 
/// Calculate the cartesian x-y position and the Lambert phase function of an orbit given Keplerian elements and a
/// vector of times.
/**
 * \tparam vectorT a type whose elements are accessed with the [] operator.
 * \tparam arithT the type in which to perform arithmetic
 *
 * \retval -1 on error calculating the orbit (from orbitCartesianWork)
 * \retval 0 on success.
 *
 * \ingroup orbits
 */
template <typename vectorT, typename arithT>
int orbitCartesian2DPhi(
    vectorT &x,       ///< [out] the projected x positions of the orbit.  Must be at least as long as t.
    vectorT &y,       ///< [out] the projected y positions of the orbit.  Must be at least as long as t.
    vectorT &r,       ///< [out] the separation of the orbit.  Must be at least as long as t.
    vectorT &rProj,   ///< [out] the projected separation of the orbit.  Must be at least as long as t.
    vectorT &phi,     ///< [out] the Lambert phase function.  Must be at least as long as t.
    const vectorT &t, ///< [in] the times at which to calculate the projected positions
    const size_t N,   ///< [in] the number of points contained in N, and allocated in nx, ny, nz
    arithT a,         ///< [in] the semi-major axis of the orbit
    arithT P,         ///< [in] the orbital period
    arithT e,         ///< [in] the eccentricity of the orbit
    arithT t0,        ///< [in] the time of pericenter passage of the orbit
    arithT i,         ///< [in] the inclination of the orbit
    arithT w,         ///< [in] the argument of pericenter of the orbit
    arithT W          ///< [in] the longitude of the ascending node of the orbit.
)
{
    return orbitCartesianWork(
        &x, &y, (vectorT *)0, &r, &rProj, (vectorT *)0, (vectorT *)0, &phi, t, N, a, P, e, t0, i, w, W );
}
 
} // namespace astro
} // namespace mx
 
#endif // mx_astro_orbitUtils_hpp