[F&E]: Time Conversion Issues for Legacy MPC Observations

[luigigisolfi]

Description:

Functions from the DateTime module do not work before a given date (1958-09-18 19:31:01).
This is not a problem for most applications, but it arises whenever we want to MPC-query objects like Eros, for which observations are available starting from 1893-10-29.

For context, I would like to change the following code:

                # convert JD to J2000 and UTC, convert deg to rad
                obs = (
                    obs.assign(
                        epochJ2000secondsTDB=lambda x: (
                            Time(x.epoch, format="jd", scale="utc").tdb.value
                            - 2451545.0
                        )
                        * 86400
                    )
                    .assign(RA=lambda x: np.radians(x.RA))
                    .assign(DEC=lambda x: np.radians(x.DEC))
                    .assign(epochUTC=lambda x: Time(x.epoch, format="jd").to_datetime())
                )

to the following (faster and more robust, leveraging tudatpy's DateTime)

                dt_objects = [DateTime.from_julian_day(jd) for jd in obs.epoch]

                # Convert DateTime objects to epochJ2000secondsTDB
                obs['epochJ2000secondsTDB'] = [
                    time_scale_converter.convert_time(
                        input_scale=time_representation.utc_scale,
                        output_scale=time_representation.tdb_scale,
                        input_value=dt_obj.epoch()
                    ) for dt_obj in dt_objects
                ]

                # Convert DateTime objects to Python datetime objects for epochUTC
                obs['epochUTC'] = [DateTime.to_python_datetime(dt_obj) for dt_obj in dt_objects]

However, this is not possible because it would break Eros's estimation example, since the current way to query an object from the MPC only allows date filtering AFTER querying all observations for the object itself. In other words, the mpc wrapper will ask for all dates and only then one could perform: batch.filter(start_date,end_date)

This breaks the workflow of the get_observations() function in mpc.py.

[DominicDirkx]

Could you add some context on why this does not work (which error do you het, etc.)

[luigigisolfi]

The current sofa wrapper handles conversions flawlessly but, unlike the previous (astropy) time conversion functions used in the mpc wrapper, it throws an error if a date is far in the past (in this case, very far = 1893-10-29, date of the first Eros observations):

An error occurred while retrieving observations of MPC: 433
Provided caledar date cannot properly give Delta AT1893 10 29

CORRESPONDING SOURCE CODE:

> Get number of leap seconds and check feasibility of calculation
>     double deltaAt;
>     int deltaAtReturn = iauDat( year, month, day, fractionOfDay, &deltaAt );
>     if( deltaAtReturn != 0 )
>     {
>         throw std::runtime_error( "Provided caledar date cannot properly give Delta AT" + std::to_string( year ) + " " +
>                                   std::to_string( month ) + " " + std::to_string( day ) );
>     }
> 
>     return deltaAt;
> }
> 

However, the MPC query by astroquery does not allow to write something like:

astroquery.MPC.get_observations('433', start_time, end_time)

but only

astroquery.MPC.get_observations('433')

which means we can only filter our dates after that call, writing something like

batch.filter(epoch_start, epoch_end)

However, to create the batch object in the way tudat wants it, time conversions are applied to the original query, before being able to filter. And, since SOFA does not allow conversion before a certain date, the error above is thrown and the

estimation_with_mpc.py script fails.

Moral of the story: as far as I understand, astropy uses lower accuracy transformations for dates in the past, but the coverage is still there. Our SOFA wrapper gives up instead. One quick solution would be to check dates and use SOFA and astropy according to the date. But this would slow down the mpc query.

Another option (but not a real solution) would be to create a batch from astroquery and allow filtering before time conversion (notice that this would be good for most of our applications, but not for analysis with very old observations).

[DominicDirkx]

There are values for Delta AT in Tudat that go back to 1620, but apparently they only get used in the time_scale_converter when going from TAI->UTC, not UTC->TAI. But, luckily this should be fairly straighforward to resolve.

I'll put some code snippets below rather than links to code (let me know if you'd prefer it differently). They're all from the TerrestrialTimeScaleConverter class definition in tudatpy/include/tudat/astro/earth_orientation/terrestrialTimeScaleConverter.h. The code below gets called when computing UTC from TAI:

//! Function to update the universal times (UT1 and UTC) in CurrentTimes member at requested precision
    template< typename TimeType >
    void calculateUniversalTimes( )
    {
        if( getCurrentTimeList< TimeType >( ).tai < utcIntroductionEpochInTai_ )
        {
            double approximateYear = getCurrentTimeList< TimeType >( ).tt / physical_constants::JULIAN_YEAR + 2000.0;
            try
            {
                getCurrentTimeList< TimeType >( ).ut1 =
                    getCurrentTimeList< TimeType >( ).tt - historicalDeltaTInterpolator_->interpolate( approximateYear );
            }
            catch( std::runtime_error& caughtException )
            {
                throw std::runtime_error( "Error in historical UT calculation.\nOriginal error: " + std::string( caughtException.what( ) ) );
            }


            getCurrentTimeList< TimeType >( ).utc = getCurrentTimeList< TimeType >( ).ut1;
        }
        else
        {
            getCurrentTimeList< TimeType >( ).utc = sofa_interface::convertTAItoUTC< TimeType >( getCurrentTimeList< TimeType >( ).tai );

            if( dailyUtcUt1CorrectionInterpolator_ != nullptr )
            {
                try
                {
                    getCurrentTimeList< TimeType >( ).ut1 = static_cast< TimeType >( dailyUtcUt1CorrectionInterpolator_->interpolate(
                                                                    getCurrentTimeList< TimeType >( ).utc ) ) + getCurrentTimeList< TimeType >( ).utc;
                }
                catch( std::runtime_error& caughtException )
                {
                    throw std::runtime_error( "Error in UT1 calculation.\nOriginal error: " + std::string( caughtException.what( ) ) );
                }


                getCurrentTimeList< TimeType >( ).ut1 += static_cast< TimeType >(
                        shortPeriodUt1CorrectionCalculator_->getCorrections( getCurrentTimeList< TimeType >( ).tt ) );
            }
        }
    }

if the requested epoch is outside the bounds where the SOFA Delta AT function is defined, we use interpolated values from the historicalDeltaTInterpolator_.

The conversions from different scales as 'starting point' are in the

template< typename TimeType >
    void updateTimes( const basic_astrodynamics::TimeScales inputScale,
                      const TimeType& inputTimeValue,
                      const Eigen::Vector3d& earthFixedPosition )

function.

In your case, when going from UTC to any other scale, this gets called:

    ```
case basic_astrodynamics::utc_scale:
            timesToUpdate.utc = inputTimeValue;
            timesToUpdate.tai = sofa_interface::convertUTCtoTAI< TimeType >( timesToUpdate.utc );

            timesToUpdate.tt = basic_astrodynamics::convertTAItoTT< TimeType >( timesToUpdate.tai );
            tdbMinusTt = static_cast< TimeType >( this->getTDBminusTT( timesToUpdate.tt, earthFixedPosition ) );
            timesToUpdate.tdb = timesToUpdate.tt + tdbMinusTt;

Your error occurs in the line

`timesToUpdate.tai = sofa_interface::convertUTCtoTAI< TimeType >( timesToUpdate.utc );`

this should be modified to a separate function (like the `calculateUniversalTimes`) that checks whether to use sofa or our interpolator for Delta AT.