astroDynamics.hpp

Go to the documentation of this file.

/** \file astroDynamics.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Various utilities for astrodynamics.
 * \ingroup astrofiles
 *
 */
 
//***********************************************************************//
// Copyright 2018 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef astroDynamics_hpp
#define astroDynamics_hpp
 
#include <cmath>
 
#include "sofa.hpp"
 
#include "../math/constants.hpp"
#include "../math/geo.hpp"
 
namespace mx
{
namespace astro
{
 
/// Breaks down a decimal day into hours, minutes and decimal point seconds.
/** Assumes 86400 seconds per day, ignoring leap seconds.
 *
 * \tparam realT the floating point type for calculations
 *
 * \retval 0 on success
 * \retval -1 on failure.
 */
template <typename realT>
int hrMinSecFromDay( int &Dy,    ///< [out] the integer day
                     int &hr,    ///< [out] the integer hour
                     int &min,   ///< [out] the integer minute
                     realT &sec, ///< [out] the integer second
                     realT day   ///< [in] the decimal day to be broken down
)
{
    Dy = (int)day;
    hr = ( day - (realT)Dy ) * 24.0;
    min = ( day - (realT)Dy - ( (realT)hr ) / 24.0 ) * 60. * 24.0;
    sec = ( day - (realT)Dy - ( (realT)hr ) / 24.0 - ( (realT)min ) / ( 60.0 * 24.0 ) ) * 3600.0 * 24.0;
 
    return 0;
}
 
/// Returns Greenwich Mean Sidereal Time for a given UTC time.
/** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
 *
 * \retval 0 on success
 */
template <typename realT>
int getGMST( realT &GMST, ///< [out] GMST in hours
             int Yr,      ///< [in] integer year
             int Mo,      ///< [in] integer month (1-12)
             int Dy,      ///< [in] integer day (1-31)
             int hr,      ///< [in] integer hour (0-23)
             int min,     ///< [in] integer minute (0-59)
             realT sec    ///< [in] decimal second (0-59)
)
{
    int rv;
    double utc0, utc1, tai0, tai1, tt0, tt1, ut10, ut11;
 
    // Get UTC time in SOFA 2-double format
    rv = sofa::iauDtf2d( "UTC", Yr, Mo, Dy, hr, min, static_cast<double>( sec ), &utc0, &utc1 );
    if( rv < 0 )
    {
        /// \todo handle SOFA error.
        return -1;
    }
 
    // Convert UTC to TAI
    rv = sofa::iauUtctai( utc0, utc1, &tai0, &tai1 );
    if( rv < 0 )
    {
        /// \todo handle SOFA error.
        return -1;
    }
 
    // Convert TAI to TT
    rv = sofa::iauTaitt( tai0, tai1, &tt0, &tt1 );
    if( rv < 0 )
    {
        /// \todo handle SOFA error.
        return -1;
    }
 
    // Convert UTC to UT1
    rv = sofa::iauUtcut1( utc0, utc1, -.4, &ut10, &ut11 ); // The current DUT1 - how to keep updated?
    if( rv < 0 )
    {
        /// \todo handle SOFA error.
        return -1;
    }
 
    GMST = sofa::iauGmst06( ut10, ut11, tt0, tt1 ) / ( math::two_pi<realT>() ) * static_cast<realT>( 24 );
 
    return 0;
}
 
/// Returns Greenwich Mean Sidereal Time for a given UTC time.
/** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
 * \param GMST [output] GMST in hours
 * \param Yr [input] integer year
 * \param Mo [input] integer month (1-12)
 * \param day [input] decimal day (1-31)
 * \retval 0 on sucess
 */
template <typename realT>
int getGMST( realT &GMST, ///< [out] GMST in hours
             int Yr,      ///< [in] integer year
             int Mo,      ///< [in] integer month (1-12)
             realT day    ///< [in] decimal day (1-31)
)
{
    int Dy, hr, min;
    double sec;
 
    hrMinSecFromDay<realT>( Dy, hr, min, sec, day );
    return getGMST( GMST, Yr, Mo, Dy, hr, min, sec );
}
 
/// Returns Local Mean Sidereal Time for a given UTC time and longitude.
/** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
 *
 * \retval 0 on success
 *
 */
template <typename realT>
int getLMST( realT &LMST, ///< [out] LMST in hours
             int Yr,      ///< [in] integer year
             int Mo,      ///< [in] integer month (1-12)
             int Dy,      ///< [in] integer day (1-31)
             int hr,      ///< [in] integer hour (0-23)
             int min,     ///< [in] integer minute (0-59)
             realT sec,   ///< [in] decimal second (0-59)
             realT lon    ///< [in] longitude in degrees, E+, W-
)
{
    int rv;
    realT GMST;
 
    rv = getGMST( GMST, Yr, Mo, Dy, hr, min, sec );
 
    LMST = GMST + lon / static_cast<realT>( 15 );
 
    LMST = fmod( LMST, 24 );
    if( LMST < 0.0 )
        LMST += 24.;
 
    return rv;
}
 
/// Returns Local Mean Sidereal Time for a given UTC time.
/** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
 *
 * \tparam realT the floating point type for all calculations
 *
 * \retval 0 on sucess
 */
template <typename realT>
int getLMST( realT &LMST, ///< [out] LMST in hours
             int Yr,      ///< [in] integer year
             int Mo,      ///< [in] integer month (1-12)
             realT day,   ///< [in] decimal day (1-31)
             realT lon    ///< [in] longitude in degrees, E+, W-
)
{
    int Dy, hr, min;
    realT sec;
 
    hrMinSecFromDay<realT>( Dy, hr, min, sec, day );
 
    return getLMST( LMST, Yr, Mo, Dy, hr, min, sec, lon );
}
 
/// Calculates the Azimuth and Elevation for a given Hour-Angle, Declination, and Latitude.
/** References: J. Meeus "Astronomical Algorithms", 1991; V. Pisacane "Fundamentals of Space Systems", 2nd ed., 2005.
 *
 * \tparam realT the real floating point type for calculations
 *
 * \returns 0 on success
 * \returns -1 on error
 *
 * \ingroup posastro
 */
template <typename realT>
int calcAzEl( realT &az, ///< [out] the calculated azimuth [degrees]
              realT &el, ///< [out] the calculate elevation [degrees]
              realT ha,  ///< [in] the hour ange [degrees]
              realT dec, ///< [in] the declination [degrees]
              realT lat  ///< [in] the latitude [degrees]
)
{
    realT ha_rad, dec_rad, lat_rad, az_rad, el_rad;
 
    ha_rad = ha * static_cast<realT>( 15 ) * math::pi<realT>() / static_cast<realT>( 180 );
    dec_rad = dec * math::pi<realT>() / static_cast<realT>( 180 );
    lat_rad = lat * math::pi<realT>() / static_cast<realT>( 180 );
 
    az_rad = atan2( sin( ha_rad ), cos( ha_rad ) * sin( lat_rad ) - tan( dec_rad ) * cos( lat_rad ) );
 
    el_rad = asin( sin( lat_rad ) * sin( dec_rad ) + cos( lat_rad ) * cos( dec_rad ) * cos( ha_rad ) );
 
    az = az_rad * static_cast<realT>( 180 ) / math::pi<realT>() + static_cast<realT>( 180 );
    el = el_rad * static_cast<realT>( 180 ) / math::pi<realT>();
 
    return 0;
}
 
/// Calculate the Parallactic Angle from the pre-calculated trig functions.  Result in radians.
/**
 * \returns the  parallactic angle in radians.
 *
 * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
 *
 * \ingroup posastro
 */
template <typename realT>
realT parAngRad( realT sinHA,  ///< [in] the sine of the target hour angle
                 realT cosHA,  ///< [in] the cosine of the target hour angle
                 realT sinDec, ///< [in] the sine of the target declination
                 realT cosDec, ///< [in] the cosine of the target declination
                 realT tanLat  ///< [in] the tangent of the observer latitude
)
{
    return atan2( sinHA, cosDec * tanLat - sinDec * cosHA );
}
 
/// Calculate the Parallactic Angle from the pre-calculated trig functions.  Result in degrees.
/**
 * \returns the  parallactic angle in degrees.
 *
 *
 *
 * \ingroup posastro
 */
template <typename realT>
realT parAngDeg( realT sinHA,  ///< [in] the sine of the target hour angle
                 realT cosHA,  ///< [in] the cosine of the target hour angle
                 realT sinDec, ///< [in] the sine of the target declination
                 realT cosDec, ///< [in] the cosine of the target declination
                 realT tanLat  ///< [in] the tangent of the observer latitude
)
{
    return math::rtod( parAngRad( sinHA, cosHA, sinDec, cosDec, tanLat ) );
}
 
/// Calculate the Parallactic Angle, with angles in radians
/**
 * \return the parallactic angle in radians.
 *
 * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
 *
 * \ingroup posastro
 */
template <typename realT>
realT parAngRad( realT ha,  ///< [in] the hour angle, in radians, negative to the east
                 realT dec, ///< [in] the object declination, in radians.
                 realT lat  ///< [in] the observer latitude, in radians.
)
{
    return parAngRad( sin( ha ), cos( ha ), sin( dec ), cos( dec ), tan( lat ) );
}
 
/// Calculate the Parallactic Angle, with angles in degrees
/**
 * \return the parallactic angle in degrees.
 *
 * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
 *
 * \ingroup posastro
 *
 * \note 2020-Jan-19: re-reordered arguments and added newflag to prevent compilation of current programs so the new
 * order is implemented.
 */
template <typename realT>
realT parAngDeg( realT ha,    ///< [in] the hour angle, in degrees, negative to the east
                 realT dec,   ///< [in] the object declination, in degrees.
                 realT lat,   ///< [in] the observer latitude, in degrees.
                 bool newflag ///< [in] [temporary] to force adaptation of new argument order.  Unused, so can be true
                              ///< or false. Added 2020-Jan-19.
)
{
    static_cast<void>( newflag ); // be unused
 
    return math::rtod( parAngRad( math::dtor( ha ), math::dtor( dec ), math::dtor( lat ) ) );
}
 
} // namespace astro
} // namespace mx
 
#endif // astroDynamics_hpp