main.cpp

// main.cpp
// Driver for the distlink MOID library (Baluev & Mikryukov, 2018-2020).
//
// Compile:
//   g++ -std=c++20 main.cpp distlink.cpp -O3 -march=native -mfpmath=sse -o moid
// Do not add -ffast-math: the authors warn it introduces the numerical
// errors this library exists to control.

#include <iostream>
#include <iomanip>
#include <numbers>
#include "distlink.h"

int main()
{
    const double deg = std::numbers::pi_v<double> / 180.0;

    // Orbital elements: a in AU, e dimensionless, i/w/Om in radians.
    // distlink.h does not state units in its comments, but this matches
    // every published use of the algorithm, so degrees are converted here.

    // Earth, approximate J2000 mean elements.
    COrbitData<double> earth(
        1.00000011,       // a  (AU)
        0.01671022,       // e
        0.00005 * deg,    // i
        102.94719 * deg,  // w  (argument of pericenter)
        -11.26064 * deg   // Om (longitude of ascending node)
    );

    // Illustrative asteroid orbit. Replace with real elements pulled from
    // the JPL Small-Body Database for actual work.
    COrbitData<double> asteroid(
        1.4583,           // a  (AU)
        0.2226,           // e
        3.3312 * deg,     // i
        126.401 * deg,    // w
        204.061 * deg     // Om
    );

    double max_root_error, min_root_error, max_anom_error;
    detect_suitable_options<double>(max_root_error, min_root_error, max_anom_error);

    // Step 1: the fast, exact method (root-finding on the degree-16 polynomial).
    SMOIDResult<double> result =
        MOID_fast<double>(earth, asteroid, max_root_error, min_root_error);

    // Step 2 (per distlink.h): if unreliable, try the orbits in the other order.
    if (!result.good) {
        result = MOID_fast<double>(asteroid, earth, max_root_error, min_root_error);
    }

    // Step 3: if still unreliable, retry at higher precision (long double).
    if (!result.good) {
        COrbitData<long double> earth_ld(earth);
        COrbitData<long double> asteroid_ld(asteroid);

        long double ld_max_root_error, ld_min_root_error;
        [[maybe_unused]] long double ld_max_anom_error;
        detect_suitable_options<long double>(ld_max_root_error, ld_min_root_error, ld_max_anom_error);

        SMOIDResult<long double> result_ld =
            MOID_fast<long double>(earth_ld, asteroid_ld, ld_max_root_error, ld_min_root_error);

        if (result_ld.good) {
            std::cout << std::setprecision(15);
            std::cout << "MOID (AU, long double): " << result_ld.distance << "\n";
            std::cout << "reliable result:        yes\n";
            return 0;
        }
    }

    // Step 4: last resort, one-dimensional scan-and-refine search.
    if (!result.good) {
        const unsigned int densities[] = {1000, 30, 3, 4, 0};
        double max_dist_error = 1e-10;
        result = MOID_direct_search<double>(
            earth, asteroid, densities, max_dist_error, max_anom_error);
    }

    std::cout << std::setprecision(12);
    std::cout << "MOID (AU):      " << result.distance << "\n";
    std::cout << "distance error: " << result.distance_error << "\n";
    std::cout << "u1 (rad):       " << result.u1 << "\n";
    std::cout << "u2 (rad):       " << result.u2 << "\n";
    std::cout << "reliable:       " << (result.good ? "yes" : "no") << "\n";

    return 0;
}