EarthModel.hpp

Go to the documentation of this file.

// Vulcan Earth Model
// Ellipsoid parameters and Earth rotation models for coordinate transformations
#pragma once
 
#include <vulcan/core/Constants.hpp>
 
#include <cmath>
#include <memory>
 
namespace vulcan {
 
// =============================================================================
// EarthModel - Ellipsoid Parameters
// =============================================================================

struct EarthModel {
    double a;     
    double f;     
    double omega; 
    double mu;    
 
    // Derived quantities (computed at construction)
    double b;        
    double e2;       
    double e_prime2; 

    constexpr EarthModel(double a_, double f_, double omega_, double mu_)
        : a(a_), f(f_), omega(omega_), mu(mu_), b(a_ * (1.0 - f_)),
          e2(2.0 * f_ - f_ * f_),
          e_prime2((2.0 * f_ - f_ * f_) / (1.0 - (2.0 * f_ - f_ * f_))) {}

    static constexpr EarthModel WGS84() {
        return EarthModel(constants::wgs84::a,     // 6378137.0 m
                          constants::wgs84::f,     // 1/298.257223563
                          constants::wgs84::omega, // 7.292115e-5 rad/s
                          constants::wgs84::mu     // 3.986004418e14 m³/s²
        );
    }

    static constexpr EarthModel Spherical() {
        return EarthModel(constants::earth::R_mean, // 6371008.8 m
                          0.0, // No flattening (perfect sphere)
                          constants::earth::omega, // 7.292115e-5 rad/s
                          constants::earth::mu     // 3.986004418e14 m³/s²
        );
    }
};
 
// =============================================================================
// EarthRotationModel - Abstract Interface
// =============================================================================

struct EarthRotationModel {
    virtual ~EarthRotationModel() = default;

    [[nodiscard]] virtual double gmst(double t_seconds) const = 0;

    [[nodiscard]] virtual double ecef_to_eci_angle(double t_seconds) const {
        return gmst(t_seconds);
    }
};
 
// =============================================================================
// ConstantOmegaRotation - Simple Rotation Model (Phase 3a)
// =============================================================================

struct ConstantOmegaRotation : EarthRotationModel {
    double omega;  
    double theta0; 
 
    constexpr ConstantOmegaRotation(double omega_, double theta0_ = 0.0)
        : omega(omega_), theta0(theta0_) {}
 
    [[nodiscard]] double gmst(double t_seconds) const override {
        return theta0 + omega * t_seconds;
    }

    static constexpr ConstantOmegaRotation from_wgs84(double theta0 = 0.0) {
        return ConstantOmegaRotation(constants::wgs84::omega, theta0);
    }

    static constexpr ConstantOmegaRotation from_model(const EarthModel &model,
                                                      double theta0 = 0.0) {
        return ConstantOmegaRotation(model.omega, theta0);
    }
};
 
// =============================================================================
// GMSTRotation - GMST-Based Rotation Model (Optional Enhancement)
// =============================================================================

struct GMSTRotation : EarthRotationModel {
    double jd_epoch; 
 
    explicit constexpr GMSTRotation(double jd_epoch_ = 2451545.0) // J2000.0
        : jd_epoch(jd_epoch_) {}
 
    [[nodiscard]] double gmst(double t_seconds) const override {
        // Convert seconds to Julian centuries from J2000.0
        constexpr double sec_per_century = 86400.0 * 36525.0;
        const double T = t_seconds / sec_per_century;
 
        // IAU 1982 GMST formula (in seconds of time)
        const double gmst_sec = 67310.54841 +
                                (876600.0 * 3600.0 + 8640184.812866) * T +
                                0.093104 * T * T - 6.2e-6 * T * T * T;
 
        // Convert seconds of time to radians
        // 1 second of time = 15 arcsec = 15/3600 deg = π/(43200) rad
        constexpr double sec_to_rad = constants::angle::pi / 43200.0;
        return std::fmod(gmst_sec * sec_to_rad, 2.0 * constants::angle::pi);
    }
};
 
} // namespace vulcan