/*
 * 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>

/// Offset to time data in the JPL DE header, in bytes.
static constexpr std::size_t jpl_de_offset_time = 0xA5C;

/// Offset to the first coefficient table in the JPL DE header, in bytes.
static constexpr std::size_t jpl_de_offset_table1 = 0xA88;

/// Offset to the DE version number in the JPL DE header, in bytes.
static constexpr std::size_t jpl_de_offset_denum = 0xB18;

/// Offset to the second coefficient table in the JPL DE header, in bytes.
static constexpr std::size_t jpl_de_offset_table2 = 0xB1C;

/// Offset to the third coefficient table in the JPL DE header, in bytes, if the constant limit has not been exceeded.
static constexpr std::size_t jpl_de_offset_table3 = 0xB28;

/// Mask to detect bytes in the most significant word of the JPL DE version number.
static constexpr std::int32_t jpl_de_denum_endian_mask = 0xFFFF0000;

/// Number of items in the first coefficient table.
static constexpr std::size_t jpl_de_table1_count = 12;

/// Number of items in the second coefficient table.
static constexpr std::size_t jpl_de_table2_count = 1;

/// Number of items in the third coefficient table.
static constexpr std::size_t jpl_de_table3_count = 2;

/// Maximum number of items in a JPL DE file.
static constexpr std::size_t jpl_de_max_item_count = jpl_de_table1_count + jpl_de_table2_count + jpl_de_table3_count;

/// Maximum number of constants in the first set of constant names.
static constexpr std::size_t jpl_de_constant_limit = 400;

/// Length of a constant name, in bytes.
static constexpr std::size_t jpl_de_constant_length = 6;

/// Enumerated IDs of the JPL DE items.
enum
{
	/// Mercury
	jpl_de_id_mercury,
	
	/// Venus
	jpl_de_id_venus,
	
	/// Earth-Moon barycenter
	jpl_de_id_embary,
	
	/// Mars
	jpl_de_id_mars,
	
	/// Jupiter
	jpl_de_id_jupiter,
	
	/// Saturn
	jpl_de_id_saturn,
	
	/// Uranus
	jpl_de_id_uranus,
	
	/// Neptune
	jpl_de_id_neptune,
	
	/// Pluto
	jpl_de_id_pluto,
	
	/// Moon
	jpl_de_id_moon,
	
	/// Sun
	jpl_de_id_sun,
	
	/// Earth nutation
	jpl_de_id_earth_nutation,
	
	/// Lunar mantle libration
	jpl_de_id_luma_libration,
	
	/// Lunar mantle angular velocity
	jpl_de_id_luma_angular_velocity,
	
	/// TT-TDB
	jpl_de_id_tt_tdb
};

/// Number of components for each JPL DE item.
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)
};

/**
 * Deserializes an ephemeris.
 *
 * @param[out] file Text file to serialize.
 * @param[in,out] ctx Deserialize context.
 *
 * @throw deserialize_error Read error.
 */
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;
}