ephemeris.cpp

Go to the documentation of this file.

/*
 * Copyright (C) 2023  Christopher J. Howard
 *
 * This file is part of Antkeeper source code.
 *
 * Antkeeper source code 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.
 *
 * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#include <engine/physics/orbit/ephemeris.hpp>
#include <bit>
#include <cstdint>
#include <engine/resources/deserializer.hpp>
#include <engine/resources/resource-loader.hpp>
 
static constexpr std::size_t jpl_de_offset_time = 0xA5C;
 
static constexpr std::size_t jpl_de_offset_table1 = 0xA88;
 
static constexpr std::size_t jpl_de_offset_denum = 0xB18;
 
static constexpr std::size_t jpl_de_offset_table2 = 0xB1C;
 
static constexpr std::size_t jpl_de_offset_table3 = 0xB28;
 
static constexpr std::int32_t jpl_de_denum_endian_mask = 0xFFFF0000;
 
static constexpr std::size_t jpl_de_table1_count = 12;
 
static constexpr std::size_t jpl_de_table2_count = 1;
 
static constexpr std::size_t jpl_de_table3_count = 2;
 
static constexpr std::size_t jpl_de_max_item_count = jpl_de_table1_count + jpl_de_table2_count + jpl_de_table3_count;
 
static constexpr std::size_t jpl_de_constant_limit = 400;
 
static constexpr std::size_t jpl_de_constant_length = 6;
 
enum
{
    jpl_de_id_mercury,
    
    jpl_de_id_venus,
    
    jpl_de_id_embary,
    
    jpl_de_id_mars,
    
    jpl_de_id_jupiter,
    
    jpl_de_id_saturn,
    
    jpl_de_id_uranus,
    
    jpl_de_id_neptune,
    
    jpl_de_id_pluto,
    
    jpl_de_id_moon,
    
    jpl_de_id_sun,
    
    jpl_de_id_earth_nutation,
    
    jpl_de_id_luma_libration,
    
    jpl_de_id_luma_angular_velocity,
    
    jpl_de_id_tt_tdb
};
 
static constexpr std::uint8_t jpl_de_component_count[jpl_de_max_item_count] =
{
    3, // Mercury: x,y,z (km)
    3, // Venus: x,y,z (km)
    3, // Earth-Moon barycenter: x,y,z (km)
    3, // Mars: x,y,z (km)
    3, // Jupiter: x,y,z (km)
    3, // Saturn: x,y,z (km)
    3, // Uranus: x,y,z (km)
    3, // Neptune: x,y,z (km)
    3, // Pluto: x,y,z (km)
    3, // Moon: x,y,z (km)
    3, // Sun: x,y,z (km)
    2, // Earth nutation: d_psi,d_epsilon (radians)
    3, // Lunar mantle libration: phi,theta,ps (radians)
    3, // Lunar mantle angular velocity: omega_x,omega_y,omega_z (radians/day)
    1  // TT-TDB: t (seconds)
};
 
template <>
void deserializer<physics::orbit::ephemeris<double>>::deserialize(physics::orbit::ephemeris<double>& ephemeris, deserialize_context& ctx)
{
    ephemeris.trajectories.clear();
    
    // Init file reading function pointers
    std::size_t (deserialize_context::*read32)(std::byte*, std::size_t) = &deserialize_context::read32<std::endian::native>;
    std::size_t (deserialize_context::*read64)(std::byte*, std::size_t) = &deserialize_context::read64<std::endian::native>;
    
    // Read DE version number
    std::int32_t denum;
    ctx.seek(jpl_de_offset_denum);
    ctx.read8(reinterpret_cast<std::byte*>(&denum), sizeof(std::int32_t));
    
    // If file endianness does not match host endianness
    if (denum & jpl_de_denum_endian_mask)
    {
        // Use endian-swapping read functions
        if constexpr (std::endian::native == std::endian::little)
        {
            read32 = &deserialize_context::read32<std::endian::big>;
            read64 = &deserialize_context::read64<std::endian::big>;
        }
        else
        {
            read32 = &deserialize_context::read32<std::endian::little>;
            read64 = &deserialize_context::read64<std::endian::little>;
        }
    }
    
    // Read ephemeris time
    double ephemeris_time[3];
    ctx.seek(jpl_de_offset_time);
    std::invoke(read64, ctx, reinterpret_cast<std::byte*>(ephemeris_time), 3);
    
    // Make time relative to J2000 epoch
    const double epoch = 2451545.0;
    ephemeris_time[0] -= epoch;
    ephemeris_time[1] -= epoch;
    
    // Read number of constants
    std::int32_t constant_count;
    std::invoke(read32, ctx, reinterpret_cast<std::byte*>(&constant_count), 1);
    
    // Read first coefficient table
    std::int32_t coeff_table[jpl_de_max_item_count][3];
    ctx.seek(jpl_de_offset_table1);
    std::invoke(read32, ctx, reinterpret_cast<std::byte*>(coeff_table), jpl_de_table1_count * 3);
    
    // Read second coefficient table
    ctx.seek(jpl_de_offset_table2);
    std::invoke(read32, ctx, reinterpret_cast<std::byte*>(&coeff_table[jpl_de_table1_count][0]), jpl_de_table2_count * 3);
    
    // Seek past extra constant names
    if (constant_count > jpl_de_constant_limit)
    {
        ctx.seek(jpl_de_offset_table3 + (constant_count - jpl_de_constant_limit) * jpl_de_constant_length);
    }
    
    // Read third coefficient table
    std::invoke(read32, ctx, reinterpret_cast<std::byte*>(&coeff_table[jpl_de_table1_count + jpl_de_table2_count][0]), jpl_de_table3_count * 3);
    
    // Calculate number of coefficients per record
    std::int32_t record_coeff_count = 0;
    for (int i = 0; i < jpl_de_max_item_count; ++i)
    {
        std::int32_t coeff_count = coeff_table[i][0] + coeff_table[i][1] * coeff_table[i][2] * static_cast<std::int32_t>(jpl_de_component_count[i]) - 1;
        record_coeff_count = std::max(record_coeff_count, coeff_count);
    }
    
    // Calculate record size and record count
    std::size_t record_size = record_coeff_count * sizeof(double);
    std::size_t record_count = static_cast<std::size_t>((ephemeris_time[1] - ephemeris_time[0]) / ephemeris_time[2]);
    
    // Calculate coefficient strides
    std::size_t strides[11];
    for (int i = 0; i < 11; ++i)
    {
        strides[i] = coeff_table[i][2] * coeff_table[i][1] * 3;
    }
    
    // Resize ephemeris to accommodate items 0-10
    ephemeris.trajectories.resize(11);
    
    // Init trajectories
    for (int i = 0; i < 11; ++i)
    {
        auto& trajectory = ephemeris.trajectories[i];
        trajectory.t0 = ephemeris_time[0];
        trajectory.t1 = ephemeris_time[1];
        trajectory.dt = ephemeris_time[2] / static_cast<double>(coeff_table[i][2]);
        trajectory.n = coeff_table[i][1];
        trajectory.a.resize(record_count * strides[i]);
    }
    
    // Read coefficients
    for (std::size_t i = 0; i < record_count; ++i)
    {
        // Seek to coefficient record
        ctx.seek((i + 2) * record_size + 2 * sizeof(double));
        
        for (int j = 0; j < 11; ++j)
        {
            std::invoke(read64, ctx, reinterpret_cast<std::byte*>(&ephemeris.trajectories[j].a[i * strides[j]]), strides[j]);
        }
    }
}
 
template <>
std::unique_ptr<physics::orbit::ephemeris<double>> resource_loader<physics::orbit::ephemeris<double>>::load(::resource_manager& resource_manager, deserialize_context& ctx)
{
    auto resource = std::make_unique<physics::orbit::ephemeris<double>>();
    
    deserializer<physics::orbit::ephemeris<double>>().deserialize(*resource, ctx);
    
    return resource;
}