#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/functional.h>
#include <array>
#include <fstream>
#include <limits>
#include <mutex>
#include <sstream>
#include <string>
#include "geometry.hpp"
#include "julian.hpp"
#include "math_utils.hpp"
#include "interpolation.hpp"
#include "solvers.hpp"
#include "coordinates.hpp"
#include "sidereal.hpp"
#include "daf.hpp"
#include "pck.hpp"
#include "light_time.hpp"
#include "evaluators.hpp"
#include "separation.hpp"
#include "events.hpp"
#include "eclipse_footprint.hpp"
#include "search_pool.hpp"
#include "lola.hpp"
#include "visibility.hpp"
#include "precession.hpp"
#include "harmograms.hpp"
#include "planetary_evaluator.hpp"
#include "physical_visibility_kernels.hpp"
#include "houses.hpp"
#include "primary_directions.hpp"
#include "marriage_enclosure.hpp"

namespace py = pybind11;
using namespace moira::native;

namespace {

struct NativeLeapSecondEntry {
    double jd_utc = 0.0;
    double tai_minus_utc = 0.0;
};

struct NativeNaifLsk {
    double delta_t_a = 32.184;
    double k = 1.657e-3;
    double eb = 1.671e-2;
    double m0 = 6.239996;
    double m1 = 1.99096871e-7;
    std::vector<NativeLeapSecondEntry> delta_at;
    bool loaded = false;
    std::string source_path;
};

NativeNaifLsk g_native_naif_lsk;
std::mutex g_native_naif_lsk_mutex;

constexpr double NATIVE_J2000 = 2451545.0;

std::string trim_copy(const std::string& value) {
    const size_t first = value.find_first_not_of(" \t\r\n");
    if (first == std::string::npos) {
        return "";
    }
    const size_t last = value.find_last_not_of(" \t\r\n");
    return value.substr(first, last - first + 1);
}

std::string replace_fortran_d(std::string value) {
    for (char& ch : value) {
        if (ch == 'D' || ch == 'd') {
            ch = 'E';
        }
    }
    return value;
}

int month_from_naif_abbrev(const std::string& month) {
    if (month == "JAN") return 1;
    if (month == "FEB") return 2;
    if (month == "MAR") return 3;
    if (month == "APR") return 4;
    if (month == "MAY") return 5;
    if (month == "JUN") return 6;
    if (month == "JUL") return 7;
    if (month == "AUG") return 8;
    if (month == "SEP") return 9;
    if (month == "OCT") return 10;
    if (month == "NOV") return 11;
    if (month == "DEC") return 12;
    throw std::runtime_error("Unsupported NAIF month token in LSK: " + month);
}

double jd_from_naif_date_token(const std::string& token) {
    if (token.size() < 2 || token[0] != '@') {
        throw std::runtime_error("Invalid NAIF date token in LSK: " + token);
    }

    std::string raw = token.substr(1);
    std::vector<std::string> parts;
    std::stringstream ss(raw);
    std::string item;
    while (std::getline(ss, item, '-')) {
        parts.push_back(trim_copy(item));
    }
    if (parts.size() != 3) {
        throw std::runtime_error("Unsupported NAIF date token format in LSK: " + token);
    }

    const int year = std::stoi(parts[0]);
    const int month = month_from_naif_abbrev(parts[1]);
    const int day = std::stoi(parts[2]);
    return julian_day(year, month, day, 0.0);
}

NativeNaifLsk parse_naif_lsk(const std::string& path) {
    std::ifstream input(path);
    if (!input) {
        throw std::runtime_error("Could not open NAIF LSK: " + path);
    }

    NativeNaifLsk parsed;
    parsed.source_path = path;

    auto parse_delta_at_row = [&parsed](const std::string& row_text) {
        const std::string row = trim_copy(row_text);
        if (row.empty()) {
            return;
        }
        const size_t comma_pos = row.find(',');
        if (comma_pos == std::string::npos) {
            throw std::runtime_error("Malformed DELTET/DELTA_AT row in LSK: " + row);
        }
        const double offset = std::stod(trim_copy(row.substr(0, comma_pos)));
        const std::string date_token = trim_copy(row.substr(comma_pos + 1));
        parsed.delta_at.push_back({jd_from_naif_date_token(date_token), offset});
    };

    std::string line;
    bool in_delta_at = false;
    while (std::getline(input, line)) {
        const std::string trimmed = trim_copy(line);
        if (trimmed.empty()) {
            continue;
        }

        if (in_delta_at) {
            const size_t close_pos = trimmed.find(')');
            const std::string row = close_pos == std::string::npos ? trimmed : trimmed.substr(0, close_pos);
            parse_delta_at_row(row);
            if (close_pos != std::string::npos) {
                in_delta_at = false;
            }
            continue;
        }

        if (trimmed.rfind("DELTET/DELTA_T_A", 0) == 0) {
            const size_t eq = trimmed.find('=');
            parsed.delta_t_a = std::stod(replace_fortran_d(trim_copy(trimmed.substr(eq + 1))));
            continue;
        }
        if (trimmed.rfind("DELTET/K", 0) == 0) {
            const size_t eq = trimmed.find('=');
            parsed.k = std::stod(replace_fortran_d(trim_copy(trimmed.substr(eq + 1))));
            continue;
        }
        if (trimmed.rfind("DELTET/EB", 0) == 0) {
            const size_t eq = trimmed.find('=');
            parsed.eb = std::stod(replace_fortran_d(trim_copy(trimmed.substr(eq + 1))));
            continue;
        }
        if (trimmed.rfind("DELTET/M", 0) == 0) {
            const size_t open = trimmed.find('(');
            const size_t close = trimmed.find(')');
            const std::string values = trimmed.substr(open + 1, close - open - 1);
            std::stringstream value_stream(values);
            std::string first;
            std::string second;
            value_stream >> first >> second;
            parsed.m0 = std::stod(replace_fortran_d(first));
            parsed.m1 = std::stod(replace_fortran_d(second));
            continue;
        }
        if (trimmed.rfind("DELTET/DELTA_AT", 0) == 0) {
            const size_t open = trimmed.find('(');
            if (open != std::string::npos) {
                const std::string inline_row = trim_copy(trimmed.substr(open + 1));
                if (!inline_row.empty()) {
                    parse_delta_at_row(inline_row);
                }
            }
            in_delta_at = true;
            continue;
        }
    }

    if (parsed.delta_at.empty()) {
        throw std::runtime_error("NAIF LSK did not provide any DELTET/DELTA_AT entries: " + path);
    }
    parsed.loaded = true;
    return parsed;
}

void load_naif_lsk_py(const std::string& path) {
    NativeNaifLsk parsed = parse_naif_lsk(path);
    std::lock_guard<std::mutex> lock(g_native_naif_lsk_mutex);
    g_native_naif_lsk = std::move(parsed);
}

double native_tai_minus_utc(double jd_utc) {
    if (!g_native_naif_lsk.loaded) {
        throw std::runtime_error("Native NAIF LSK is not loaded");
    }
    if (jd_utc < g_native_naif_lsk.delta_at.front().jd_utc) {
        throw std::runtime_error("Native NAIF LSK time admission does not support pre-1972 UTC epochs");
    }
    double result = g_native_naif_lsk.delta_at.front().tai_minus_utc;
    for (const NativeLeapSecondEntry& entry : g_native_naif_lsk.delta_at) {
        if (jd_utc >= entry.jd_utc) {
            result = entry.tai_minus_utc;
        } else {
            break;
        }
    }
    return result;
}

double jd_utc_to_et_seconds_past_j2000_py(double jd_utc) {
    std::lock_guard<std::mutex> lock(g_native_naif_lsk_mutex);
    const double utc_seconds_past_j2000 = (jd_utc - NATIVE_J2000) * 86400.0;
    const double delta_at = native_tai_minus_utc(jd_utc);

    double et = utc_seconds_past_j2000 + g_native_naif_lsk.delta_t_a + delta_at;
    for (int i = 0; i < 6; ++i) {
        const double mean_anomaly = g_native_naif_lsk.m0 + g_native_naif_lsk.m1 * et;
        const double eccentric_anomaly = mean_anomaly + g_native_naif_lsk.eb * std::sin(mean_anomaly);
        const double et_minus_utc = g_native_naif_lsk.delta_t_a + delta_at + g_native_naif_lsk.k * std::sin(eccentric_anomaly);
        et = utc_seconds_past_j2000 + et_minus_utc;
    }
    return et;
}

std::vector<double> load_double_vector(const py::handle& src_obj, const char* label) {
    py::sequence src = py::reinterpret_borrow<py::sequence>(src_obj);
    std::vector<double> out;
    out.reserve(static_cast<size_t>(py::len(src)));
    for (const py::handle& item : src) {
        out.push_back(py::cast<double>(item));
    }
    if (out.empty()) {
        throw std::runtime_error(std::string(label) + " cannot be empty");
    }
    return out;
}

std::vector<int32_t> load_int32_vector(const py::handle& src_obj, const char* label) {
    py::sequence src = py::reinterpret_borrow<py::sequence>(src_obj);
    std::vector<int32_t> out;
    out.reserve(static_cast<size_t>(py::len(src)));
    for (const py::handle& item : src) {
        out.push_back(py::cast<int32_t>(item));
    }
    if (out.empty()) {
        throw std::runtime_error(std::string(label) + " cannot be empty");
    }
    return out;
}

struct DenseMatrixView {
    std::vector<double> values;
    size_t rows = 0;
    size_t cols = 0;
};

DenseMatrixView load_double_matrix(const py::handle& src_obj, const char* label) {
    py::sequence rows_src = py::reinterpret_borrow<py::sequence>(src_obj);
    DenseMatrixView out;
    out.rows = static_cast<size_t>(py::len(rows_src));
    if (out.rows == 0) {
        throw std::runtime_error(std::string(label) + " cannot be empty");
    }

    py::sequence first_row = py::reinterpret_borrow<py::sequence>(rows_src[0]);
    out.cols = static_cast<size_t>(py::len(first_row));
    if (out.cols == 0) {
        throw std::runtime_error(std::string(label) + " rows cannot be empty");
    }

    out.values.reserve(out.rows * out.cols);
    for (const py::handle& row_obj : rows_src) {
        py::sequence row = py::reinterpret_borrow<py::sequence>(row_obj);
        if (static_cast<size_t>(py::len(row)) != out.cols) {
            throw std::runtime_error(std::string(label) + " rows must have consistent length");
        }
        for (const py::handle& item : row) {
            out.values.push_back(py::cast<double>(item));
        }
    }
    return out;
}

struct DenseTensor3View {
    std::vector<double> values;
    size_t dim0 = 0;
    size_t dim1 = 0;
    size_t dim2 = 0;
};

DenseTensor3View load_double_tensor3(const py::handle& src_obj, const char* label) {
    py::sequence outer = py::reinterpret_borrow<py::sequence>(src_obj);
    DenseTensor3View out;
    out.dim0 = static_cast<size_t>(py::len(outer));
    if (out.dim0 == 0) {
        throw std::runtime_error(std::string(label) + " cannot be empty");
    }

    py::sequence first_mid = py::reinterpret_borrow<py::sequence>(outer[0]);
    out.dim1 = static_cast<size_t>(py::len(first_mid));
    if (out.dim1 == 0) {
        throw std::runtime_error(std::string(label) + " inner dimensions cannot be empty");
    }

    py::sequence first_inner = py::reinterpret_borrow<py::sequence>(first_mid[0]);
    out.dim2 = static_cast<size_t>(py::len(first_inner));
    if (out.dim2 == 0) {
        throw std::runtime_error(std::string(label) + " innermost dimension cannot be empty");
    }

    out.values.reserve(out.dim0 * out.dim1 * out.dim2);
    for (const py::handle& mid_obj : outer) {
        py::sequence mid = py::reinterpret_borrow<py::sequence>(mid_obj);
        if (static_cast<size_t>(py::len(mid)) != out.dim1) {
            throw std::runtime_error(std::string(label) + " mid dimensions must have consistent length");
        }
        for (const py::handle& inner_obj : mid) {
            py::sequence inner = py::reinterpret_borrow<py::sequence>(inner_obj);
            if (static_cast<size_t>(py::len(inner)) != out.dim2) {
                throw std::runtime_error(std::string(label) + " inner dimensions must have consistent length");
            }
            for (const py::handle& item : inner) {
                out.values.push_back(py::cast<double>(item));
            }
        }
    }
    return out;
}

py::list matrix_to_py_list(const double* values, size_t rows, size_t cols) {
    py::list out;
    for (size_t i = 0; i < rows; ++i) {
        py::list row;
        for (size_t j = 0; j < cols; ++j) {
            row.append(values[i * cols + j]);
        }
        out.append(std::move(row));
    }
    return out;
}

NutationSeriesTerm load_native_nutation_term(const py::handle& row_obj) {
    py::sequence row = py::reinterpret_borrow<py::sequence>(row_obj);
    const size_t size = static_cast<size_t>(py::len(row));
    if (size < 3) {
        throw std::runtime_error("Nutation term rows must contain at least 3 entries");
    }

    NutationSeriesTerm term;
    term.c1 = py::cast<double>(row[0]);
    term.c2 = py::cast<double>(row[1]);
    term.argument_count = size - 2;
    if (term.argument_count > term.arguments.size()) {
        throw std::runtime_error("Nutation term rows exceed the supported 14 arguments");
    }
    for (size_t i = 0; i < term.argument_count; ++i) {
        term.arguments[i] = py::cast<int>(row[i + 2]);
    }
    return term;
}

void set_nutation_2000a_tables_py(
    const py::sequence& ls_terms_src,
    const py::sequence& pl_terms_src,
    size_t ls_j0_count,
    size_t pl_j0_count
) {
    std::vector<NutationSeriesTerm> ls_terms;
    std::vector<NutationSeriesTerm> pl_terms;
    ls_terms.reserve(static_cast<size_t>(py::len(ls_terms_src)));
    pl_terms.reserve(static_cast<size_t>(py::len(pl_terms_src)));

    for (const py::handle& row : ls_terms_src) {
        ls_terms.push_back(load_native_nutation_term(row));
    }
    for (const py::handle& row : pl_terms_src) {
        pl_terms.push_back(load_native_nutation_term(row));
    }

    register_nutation_2000r06_series(
        std::move(ls_terms),
        std::move(pl_terms),
        ls_j0_count,
        pl_j0_count
    );
}

py::tuple nutation_2000a_py(double jd_tt) {
    NutationResult result;
    {
        py::gil_scoped_release release;
        result = nutation_2000r06(jd_tt);
    }
    return py::make_tuple(result.longitude * RAD2DEG, result.obliquity * RAD2DEG);
}

py::list vector_to_py_list(const std::vector<double>& values) {
    py::list out;
    for (double v : values) out.append(v);
    return out;
}

py::tuple vector_to_py_tuple(const std::vector<double>& values) {
    py::tuple out(values.size());
    for (std::size_t index = 0; index < values.size(); ++index) {
        out[index] = values[index];
    }
    return out;
}

py::tuple physical_visibility_response_weights_py(
    const py::handle& band_wavelength_nm,
    const py::handle& band_differential_magnitude,
    const py::handle& spectral_bin_start_nm,
    double base_scotopic_to_photopic_ratio,
    const py::handle& base_photopic,
    const py::handle& base_scotopic
) {
    const std::vector<double> band_wavelengths = load_double_vector(
        band_wavelength_nm, "Band wavelengths"
    );
    const std::vector<double> band_differentials = load_double_vector(
        band_differential_magnitude, "Band differential magnitudes"
    );
    const std::vector<double> spectral_bins = load_double_vector(
        spectral_bin_start_nm, "Spectral bins"
    );
    const std::vector<double> photopic = load_double_vector(
        base_photopic, "Base photopic weights"
    );
    const std::vector<double> scotopic = load_double_vector(
        base_scotopic, "Base scotopic weights"
    );
    PhysicalVisibilityResponseWeights result;
    {
        py::gil_scoped_release release;
        result = resolve_physical_visibility_response_weights(
            band_wavelengths,
            band_differentials,
            spectral_bins,
            base_scotopic_to_photopic_ratio,
            photopic,
            scotopic
        );
    }
    return py::make_tuple(
        result.scotopic_to_photopic_ratio,
        vector_to_py_tuple(result.photopic),
        vector_to_py_tuple(result.scotopic)
    );
}

void load_mat3(const py::sequence& src, double out[3][3]) {
    if (py::len(src) != 3) {
        throw std::runtime_error("mat3 must have exactly 3 rows");
    }
    for (size_t i = 0; i < 3; ++i) {
        py::sequence row = py::reinterpret_borrow<py::sequence>(src[i]);
        if (py::len(row) != 3) {
            throw std::runtime_error("mat3 rows must have exactly 3 columns");
        }
        for (size_t j = 0; j < 3; ++j) {
            out[i][j] = py::cast<double>(row[j]);
        }
    }
}

void load_vec3(const py::sequence& src, double out[3]) {
    if (py::len(src) != 3) {
        throw std::runtime_error("vec3 must have exactly 3 components");
    }
    for (size_t i = 0; i < 3; ++i) {
        out[i] = py::cast<double>(src[i]);
    }
}

py::tuple mat3_to_py_tuple(const double src[3][3]) {
    return py::make_tuple(
        py::make_tuple(src[0][0], src[0][1], src[0][2]),
        py::make_tuple(src[1][0], src[1][1], src[1][2]),
        py::make_tuple(src[2][0], src[2][1], src[2][2])
    );
}

py::tuple vec3_to_py_tuple(const double src[3]) {
    return py::make_tuple(src[0], src[1], src[2]);
}

py::tuple rotation_matrix_multiply_py(const py::sequence& a_src, const py::sequence& b_src) {
    double a[3][3];
    double b[3][3];
    double out[3][3];
    load_mat3(a_src, a);
    load_mat3(b_src, b);

    for (size_t i = 0; i < 3; ++i) {
        for (size_t j = 0; j < 3; ++j) {
            out[i][j] = (
                a[i][0] * b[0][j] +
                a[i][1] * b[1][j] +
                a[i][2] * b[2][j]
            );
        }
    }
    return mat3_to_py_tuple(out);
}

py::tuple rotation_matrix_apply_py(const py::sequence& m_src, const py::sequence& v_src) {
    double m[3][3];
    double v[3];
    double out[3];
    load_mat3(m_src, m);
    load_vec3(v_src, v);

    for (size_t i = 0; i < 3; ++i) {
        out[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
    }
    return vec3_to_py_tuple(out);
}

py::tuple apply_aberration_velocity_py(const py::sequence& xyz_src, const py::sequence& velocity_src) {
    double xyz_vals[3];
    double velocity_vals[3];
    load_vec3(xyz_src, xyz_vals);
    load_vec3(velocity_src, velocity_vals);
    constexpr double c_km_per_day = 299792.458 * 86400.0;

    const double x = xyz_vals[0];
    const double y = xyz_vals[1];
    const double z = xyz_vals[2];
    const double dist = std::sqrt(x * x + y * y + z * z);
    if (dist < 1e-10) {
        return py::make_tuple(x, y, z);
    }

    const double ux = x / dist;
    const double uy = y / dist;
    const double uz = z / dist;

    const double bx = velocity_vals[0] / c_km_per_day;
    const double by = velocity_vals[1] / c_km_per_day;
    const double bz = velocity_vals[2] / c_km_per_day;

    const double beta2 = bx * bx + by * by + bz * bz;
    const double gamma = 1.0 / std::sqrt(1.0 - beta2);
    const double dot = ux * bx + uy * by + uz * bz;
    const double factor1 = 1.0 + dot / (1.0 + gamma);
    const double factor2 = gamma * (1.0 + dot);

    double ax = (ux + factor1 * bx) / factor2;
    double ay = (uy + factor1 * by) / factor2;
    double az = (uz + factor1 * bz) / factor2;

    const double scale = dist / std::sqrt(ax * ax + ay * ay + az * az);
    const double out[3] = {ax * scale, ay * scale, az * scale};
    return vec3_to_py_tuple(out);
}

py::tuple apply_frame_bias_py(const py::sequence& xyz_src) {
    double xyz_vals[3];
    load_vec3(xyz_src, xyz_vals);

    constexpr double arcsec2rad = 3.141592653589793238462643383279502884 / 648000.0;
    constexpr double dA_r = (-14.6 / 1000.0) * arcsec2rad;
    constexpr double xi0_r = (-16.6170 / 1000.0) * arcsec2rad;
    constexpr double de0_r = (-6.8192 / 1000.0) * arcsec2rad;

    const double x = xyz_vals[0];
    const double y = xyz_vals[1];
    const double z = xyz_vals[2];
    const double xb = x - de0_r * y + xi0_r * z;
    const double yb = de0_r * x + y - dA_r * z;
    const double zb = -xi0_r * x + dA_r * y + z;

    const double out[3] = {xb, yb, zb};
    return vec3_to_py_tuple(out);
}

/**
 * @brief 2D version: evaluates one record provided as a (coefficient, component) array.
 */
py::list spk_chebyshev_record_py(
    const py::sequence& coeff_record,
    double s
) {
    DenseMatrixView coeff_record_view = load_double_matrix(coeff_record, "Evaluator coefficients");
    const size_t component_count = coeff_record_view.rows;
    const size_t coefficient_count = coeff_record_view.cols;
    const auto* coeffs = coeff_record_view.values.data();

    std::vector<double> result(component_count);
    spk_chebyshev_record_inplace(
        coeffs, coefficient_count, component_count, s, result.data(), 
        1,                // coeff_stride = 1
        coefficient_count // component_stride = n
    );

    return vector_to_py_list(result);
}

py::tuple spk_chebyshev_record_with_derivative_py(
    const py::sequence& coeff_record,
    double s,
    double derivative_scale
) {
    DenseMatrixView coeff_record_view = load_double_matrix(coeff_record, "Evaluator coefficients");
    const size_t component_count = coeff_record_view.rows;
    const size_t coefficient_count = coeff_record_view.cols;
    const auto* coeffs = coeff_record_view.values.data();

    std::vector<double> values(component_count);
    std::vector<double> rates(component_count);

    spk_chebyshev_record_with_derivative_inplace(
        coeffs, coefficient_count, component_count, s, values.data(), rates.data(),
        1,                // coeff_stride = 1
        coefficient_count // component_stride = n
    );

    py::list v_out = vector_to_py_list(values);
    py::list r_out;
    for (double r : rates) r_out.append(r * derivative_scale);

    return py::make_tuple(v_out, r_out);
}

/**
 * @brief 3D version: evaluates one record from a (record, component, coefficient) series.
 */
py::list spk_chebyshev_series_record_py(
    const py::sequence& coefficients,
    int32_t record_index,
    double s
) {
    DenseTensor3View coeff_series_view = load_double_tensor3(coefficients, "Series evaluator coefficients");
    const size_t component_count = coeff_series_view.dim1;
    const size_t coefficient_count = coeff_series_view.dim2;
    const size_t record_stride = component_count * coefficient_count;
    if (record_index < 0 || static_cast<size_t>(record_index) >= coeff_series_view.dim0) {
        throw std::runtime_error("Series evaluator record index is out of range");
    }
    const auto* coeffs = coeff_series_view.values.data();
    std::vector<double> result(component_count);

    spk_chebyshev_record_inplace(
        coeffs + static_cast<size_t>(record_index) * record_stride,
        coefficient_count,
        component_count,
        s,
        result.data(),
        1,                // coeff_stride = 1 (contiguous)
        coefficient_count // component_stride = n
    );

    return vector_to_py_list(result);
}

py::tuple spk_chebyshev_series_record_with_derivative_py(
    const py::sequence& coefficients,
    int32_t record_index,
    double s,
    double derivative_scale
) {
    DenseTensor3View coeff_series_view = load_double_tensor3(coefficients, "Series evaluator coefficients");
    const size_t component_count = coeff_series_view.dim1;
    const size_t coefficient_count = coeff_series_view.dim2;
    const size_t record_stride = component_count * coefficient_count;
    if (record_index < 0 || static_cast<size_t>(record_index) >= coeff_series_view.dim0) {
        throw std::runtime_error("Series evaluator record index is out of range");
    }
    const auto* coeffs = coeff_series_view.values.data();
    std::vector<double> values(component_count);
    std::vector<double> rates(component_count);

    spk_chebyshev_record_with_derivative_inplace(
        coeffs + static_cast<size_t>(record_index) * record_stride,
        coefficient_count,
        component_count,
        s,
        values.data(),
        rates.data(),
        1,
        coefficient_count
    );

    py::list v_out = vector_to_py_list(values);
    py::list r_out;
    for (double r : rates) r_out.append(r * derivative_scale);

    return py::make_tuple(v_out, r_out);
}

py::tuple spk_chebyshev_series_bulk_evaluate_py(
    const py::sequence& coeff_series,
    const py::sequence& record_indices,
    const py::sequence& s_values,
    const py::sequence& derivative_scales,
    bool need_rates
) {
    DenseTensor3View coeff_series_view = load_double_tensor3(coeff_series, "Bulk evaluator coefficients");
    std::vector<int32_t> indices = load_int32_vector(record_indices, "Bulk evaluator record indices");
    std::vector<double> s_vec = load_double_vector(s_values, "Bulk evaluator s values");
    std::vector<double> ds_vec = load_double_vector(derivative_scales, "Bulk evaluator derivative scales");

    const size_t workload_size = indices.size();
    const size_t record_count = coeff_series_view.dim0;
    const size_t component_count = coeff_series_view.dim1;
    const size_t coefficient_count = coeff_series_view.dim2;
    const size_t record_stride = component_count * coefficient_count;
    if (s_vec.size() != workload_size || ds_vec.size() != workload_size) {
        throw std::runtime_error("Bulk evaluator workload arrays must have matching length");
    }

    const auto* coeffs = coeff_series_view.values.data();
    std::vector<double> values_buffer(workload_size * component_count, 0.0);
    auto* v_out_ptr = values_buffer.data();

    if (need_rates) {
        std::vector<double> rates_buffer(workload_size * component_count, 0.0);
        auto* r_out_ptr = rates_buffer.data();

        {
            py::gil_scoped_release release;
            for (size_t i = 0; i < workload_size; ++i) {
                const int32_t record_index = indices[i];
                if (record_index < 0 || static_cast<size_t>(record_index) >= record_count) {
                    throw std::runtime_error("Bulk evaluator record index is out of range");
                }
                spk_chebyshev_record_with_derivative_inplace(
                    coeffs + static_cast<size_t>(record_index) * record_stride,
                    coefficient_count, component_count, s_vec[i],
                    v_out_ptr + i * component_count,
                    r_out_ptr + i * component_count,
                    1, coefficient_count
                );

                const double ds = ds_vec[i];
                for (size_t j = 0; j < component_count; ++j) {
                    r_out_ptr[i * component_count + j] *= ds;
                }
            }
        }
        return py::make_tuple(
            matrix_to_py_list(values_buffer.data(), workload_size, component_count),
            matrix_to_py_list(rates_buffer.data(), workload_size, component_count)
        );
    } else {
        {
            py::gil_scoped_release release;
            for (size_t i = 0; i < workload_size; ++i) {
                const int32_t record_index = indices[i];
                if (record_index < 0 || static_cast<size_t>(record_index) >= record_count) {
                    throw std::runtime_error("Bulk evaluator record index is out of range");
                }
                spk_chebyshev_record_inplace(
                    coeffs + static_cast<size_t>(record_index) * record_stride,
                    coefficient_count, component_count, s_vec[i],
                    v_out_ptr + i * component_count,
                    1, coefficient_count
                );
            }
        }
        return py::make_tuple(
            matrix_to_py_list(values_buffer.data(), workload_size, component_count),
            py::none()
        );
    }
}

py::tuple daf_descriptor_to_py_tuple(const DafSummaryEntry& entry) {
    py::list descriptor;
    for (double value : entry.double_components) {
        descriptor.append(value);
    }
    for (int32_t value : entry.integer_components) {
        descriptor.append(value);
    }
    return py::tuple(descriptor);
}

py::tuple mat3_to_py_tuple(const Mat3& src) {
    return py::make_tuple(
        py::make_tuple(src[0][0], src[0][1], src[0][2]),
        py::make_tuple(src[1][0], src[1][1], src[1][2]),
        py::make_tuple(src[2][0], src[2][1], src[2][2])
    );
}

py::dict read_daf_catalog_py(const std::string& path) {
    DafCatalog catalog;
    {
        py::gil_scoped_release release;
        catalog = read_daf_catalog(path);
    }
    py::list summaries;
    for (const DafSummaryEntry& entry : catalog.summaries) {
        py::dict item;
        item["name"] = py::bytes(entry.name);
        item["descriptor"] = daf_descriptor_to_py_tuple(entry);
        summaries.append(std::move(item));
    }

    py::dict out;
    out["locidw"] = catalog.locidw;
    out["locfmt"] = catalog.locfmt;
    out["nd"] = catalog.nd;
    out["ni"] = catalog.ni;
    out["fward"] = catalog.fward;
    out["bward"] = catalog.bward;
    out["free"] = catalog.free;
    out["little_endian"] = catalog.little_endian;
    out["summaries"] = std::move(summaries);
    return out;
}

py::dict catalog_to_py_dict(const DafCatalog& catalog) {
    py::list summaries;
    for (const DafSummaryEntry& entry : catalog.summaries) {
        py::dict item;
        item["name"] = py::bytes(entry.name);
        item["descriptor"] = daf_descriptor_to_py_tuple(entry);
        summaries.append(std::move(item));
    }

    py::dict out;
    out["locidw"] = catalog.locidw;
    out["locfmt"] = catalog.locfmt;
    out["nd"] = catalog.nd;
    out["ni"] = catalog.ni;
    out["fward"] = catalog.fward;
    out["bward"] = catalog.bward;
    out["free"] = catalog.free;
    out["little_endian"] = catalog.little_endian;
    out["summaries"] = std::move(summaries);
    return out;
}

py::dict read_spk_chebyshev_segment_payload_py(
    const std::string& path, int32_t start_i, int32_t end_i, bool little_endian, int32_t data_type
) {
    SpkChebyshevSegmentPayload payload;
    {
        py::gil_scoped_release release;
        payload = read_spk_chebyshev_segment_payload(path, start_i, end_i, little_endian, data_type);
    }

    py::list records;
    for (size_t i = 0; i < payload.record_count; ++i) {
        py::list components;
        for (size_t j = 0; j < payload.component_count; ++j) {
            py::list coeffs;
            for (size_t k = 0; k < payload.coefficient_count; ++k) {
                coeffs.append(payload.coefficients[(i * payload.component_count + j) * payload.coefficient_count + k]);
            }
            components.append(py::tuple(coeffs));
        }
        records.append(py::tuple(components));
    }

    py::dict out_dict;
    out_dict["init"] = payload.init;
    out_dict["intlen"] = payload.intlen;
    out_dict["record_count"] = payload.record_count;
    out_dict["component_count"] = payload.component_count;
    out_dict["coefficient_count"] = payload.coefficient_count;
    out_dict["coefficients"] = py::tuple(records);
    return out_dict;
}

std::shared_ptr<SpkSegmentEvaluator> load_spk_segment_evaluator_py(
    const std::string& path, int32_t start_i, int32_t end_i, bool little_endian, int32_t data_type
) {
    DafCatalog catalog;
    {
        py::gil_scoped_release release;
        catalog = read_daf_catalog(path);
    }
    const auto [coverage_start_jd, coverage_end_jd] =
        spk_descriptor_coverage_jd(catalog, start_i, end_i, data_type);
    if (data_type == 13) {
        SpkType13SegmentPayload payload;
        {
            py::gil_scoped_release release;
            payload = read_spk_type13_segment_payload(path, start_i, end_i, little_endian);
        }
        return std::make_shared<SpkSegmentEvaluator>(
            coverage_start_jd,
            coverage_end_jd,
            std::move(payload.epochs_jd),
            std::move(payload.states),
            static_cast<size_t>(payload.window_size)
        );
    }

    SpkChebyshevSegmentPayload payload;
    {
        py::gil_scoped_release release;
        payload = read_spk_chebyshev_segment_payload(
            path, start_i, end_i, little_endian, data_type
        );
    }
    return std::make_shared<SpkSegmentEvaluator>(
        coverage_start_jd,
        coverage_end_jd,
        data_type,
        false,
        payload.init,
        payload.intlen,
        payload.record_count,
        payload.component_count,
        payload.coefficient_count,
        std::move(payload.coefficients)
    );
}

std::shared_ptr<NativeSpkKernelHandle> open_spk_kernel_py(const std::string& path) {
    py::gil_scoped_release release;
    return std::make_shared<NativeSpkKernelHandle>(path);
}

py::dict pck_catalog_to_py_dict(const NativePckKernelHandle& handle) {
    py::list summaries;
    for (const PckSummaryEntry& entry : handle.summaries()) {
        py::dict item;
        item["name"] = py::bytes(entry.name);
        item["descriptor"] = py::make_tuple(
            entry.start_second,
            entry.end_second,
            entry.frame_class_id,
            entry.inertial_frame_id,
            entry.data_type,
            entry.start_i,
            entry.end_i
        );
        summaries.append(std::move(item));
    }
    py::dict out;
    out["locidw"] = handle.catalog.locidw;
    out["locfmt"] = handle.catalog.locfmt;
    out["nd"] = handle.catalog.nd;
    out["ni"] = handle.catalog.ni;
    out["little_endian"] = handle.catalog.little_endian;
    out["summaries"] = std::move(summaries);
    return out;
}

std::shared_ptr<NativePckKernelHandle> open_pck_kernel_py(const std::string& path) {
    py::gil_scoped_release release;
    return std::make_shared<NativePckKernelHandle>(path);
}

py::dict read_spk_type13_segment_payload_py(const std::string& path, int32_t start_i, int32_t end_i, bool little_endian) {
    SpkType13SegmentPayload payload;
    {
        py::gil_scoped_release release;
        payload = read_spk_type13_segment_payload(path, start_i, end_i, little_endian);
    }

    py::tuple epochs(payload.state_count);
    for (size_t i = 0; i < static_cast<size_t>(payload.state_count); ++i) {
        epochs[i] = payload.epochs_jd[i];
    }

    py::list states;
    for (size_t axis = 0; axis < 6; ++axis) {
        py::tuple axis_values(payload.state_count);
        for (size_t i = 0; i < static_cast<size_t>(payload.state_count); ++i) {
            axis_values[i] = payload.states[axis * payload.state_count + i];
        }
        states.append(std::move(axis_values));
    }

    py::dict out_dict;
    out_dict["epochs_jd"] = epochs;
    out_dict["states"] = states;
    out_dict["window_size"] = payload.window_size;
    return out_dict;
}

// Eclipse-footprint bindings below accept ordinary Python sequences. Public
// doctrine and topology remain in moira/eclipse.py; this boundary has no
// NumPy or buffer-library dependency.


} // namespace

PYBIND11_MODULE(_moira_native, m) {
    m.def("_marriage_chebyshev_enclosure", &marriage_chebyshev_enclosure,
          py::arg("coefficients"), py::arg("argument"), py::arg("centered"),
          "Private exact-schedule interval Chebyshev value/derivative recurrence.");
    m.doc() = "Moira Native Backend Forge";

    // --- Evaluators ---
    py::class_<IEvaluator, std::shared_ptr<IEvaluator>>(m, "IEvaluator")
        .def("evaluate", [](IEvaluator& self, double jd) {
            double res[6];
            {
                py::gil_scoped_release release;
                self.evaluate(jd, res);
            }
            std::vector<double> out(res, res + 6);
            return vector_to_py_list(out);
        })
        .def("evaluate_batch", [](IEvaluator& self, const std::vector<double>& jds) {
            size_t count = jds.size();
            std::vector<std::vector<double>> out(count, std::vector<double>(6));
            std::vector<double> flat(count * 6);
            {
                py::gil_scoped_release release;
                self.evaluate_batch(jds.data(), count, flat.data());
            }
            for (size_t i = 0; i < count; ++i) {
                out[i].assign(flat.data() + i * 6, flat.data() + i * 6 + 6);
            }
            return out;
        });

    py::class_<TtToTdbEvaluator, IEvaluator, std::shared_ptr<TtToTdbEvaluator>>(
        m, "TtToTdbEvaluator"
    )
        .def(py::init<std::shared_ptr<IEvaluator>>(), py::arg("raw_tdb_evaluator"))
        .def_static("convert", &TtToTdbEvaluator::convert, py::arg("jd_tt"));

    py::class_<PenumbralClearanceScan>(m, "PenumbralClearanceScan")
        .def_readonly("sampled_maximum", &PenumbralClearanceScan::sampled_maximum)
        .def_readonly(
            "local_maximum_brackets",
            &PenumbralClearanceScan::local_maximum_brackets
        );

    py::class_<PenumbralClearanceScanner>(m, "PenumbralClearanceScanner")
        .def(
            py::init([](
                const std::vector<double>& epochs,
                const std::vector<std::array<double, 8>>& rows
            ) {
                std::vector<PenumbralClearanceShadow> shadows;
                shadows.reserve(rows.size());
                for (const auto& row : rows) {
                    shadows.push_back({
                        {row[0], row[1], row[2]},
                        {row[3], row[4], row[5]},
                        row[6],
                        row[7],
                    });
                }
                return PenumbralClearanceScanner(epochs, std::move(shadows));
            }),
            py::arg("epochs"),
            py::arg("shadow_rows")
        )
        .def_property_readonly("size", &PenumbralClearanceScanner::size)
        .def(
            "scan",
            [](
                const PenumbralClearanceScanner& scanner,
                const std::array<double, 3>& site_xyz_itrf_km,
                double witness_epoch,
                const std::array<double, 8>& row
            ) {
                const PenumbralClearanceShadow witness{
                    {row[0], row[1], row[2]},
                    {row[3], row[4], row[5]},
                    row[6],
                    row[7],
                };
                py::gil_scoped_release release;
                return scanner.scan(site_xyz_itrf_km, witness_epoch, witness);
            },
            py::arg("site_xyz_itrf_km"),
            py::arg("witness_epoch"),
            py::arg("witness_shadow_row")
        );

    py::class_<PenumbralEnvelopeCandidate>(m, "PenumbralEnvelopeCandidate")
        .def_readonly("azimuth_rad", &PenumbralEnvelopeCandidate::azimuth_rad)
        .def_readonly("xyz_itrf_km", &PenumbralEnvelopeCandidate::xyz_itrf_km)
        .def_readonly("latitude_deg", &PenumbralEnvelopeCandidate::latitude_deg)
        .def_readonly("longitude_deg", &PenumbralEnvelopeCandidate::longitude_deg)
        .def_readonly(
            "signed_half_chord_sq_km2",
            &PenumbralEnvelopeCandidate::signed_half_chord_sq_km2
        );

    m.def(
        "penumbral_lawful_azimuth_interval",
        [](const std::array<double, 15>& row) {
            const PenumbralEnvelopeShadow shadow{
                {
                    {row[0], row[1], row[2]},
                    {row[3], row[4], row[5]},
                    row[7],
                    row[8],
                },
                row[6],
                {row[9], row[10], row[11]},
                {row[12], row[13], row[14]},
            };
            try {
                py::gil_scoped_release release;
                eclipse_footprint_detail::validate_envelope_shadow(shadow);
                return eclipse_footprint_detail::lawful_azimuth_interval(shadow);
            } catch (const std::invalid_argument&) {
                throw;
            } catch (const std::runtime_error& error) {
                PyErr_SetString(PyExc_ArithmeticError, error.what());
                throw py::error_already_set();
            }
        },
        py::arg("shadow_row")
    );

    m.def(
        "penumbral_envelope_candidates",
        [](
            const std::array<double, 15>& row,
            const std::array<double, 8>& before_row,
            const std::array<double, 8>& after_row
        ) {
            const PenumbralEnvelopeShadow shadow{
                {
                    {row[0], row[1], row[2]},
                    {row[3], row[4], row[5]},
                    row[7],
                    row[8],
                },
                row[6],
                {row[9], row[10], row[11]},
                {row[12], row[13], row[14]},
            };
            const PenumbralClearanceShadow before{
                {before_row[0], before_row[1], before_row[2]},
                {before_row[3], before_row[4], before_row[5]},
                before_row[6],
                before_row[7],
            };
            const PenumbralClearanceShadow after{
                {after_row[0], after_row[1], after_row[2]},
                {after_row[3], after_row[4], after_row[5]},
                after_row[6],
                after_row[7],
            };
            try {
                py::gil_scoped_release release;
                return penumbral_envelope_candidates(shadow, before, after);
            } catch (const std::invalid_argument&) {
                throw;
            } catch (const std::runtime_error& error) {
                PyErr_SetString(PyExc_ArithmeticError, error.what());
                throw py::error_already_set();
            }
        },
        py::arg("shadow_row"),
        py::arg("before_shadow_row"),
        py::arg("after_shadow_row")
    );

    py::class_<ChebyshevEvaluator, IEvaluator, std::shared_ptr<ChebyshevEvaluator>>(m, "ChebyshevEvaluator")
        .def(py::init<double, double, size_t, size_t, size_t, std::vector<double>>());

    py::class_<SpkSegmentEvaluator, std::shared_ptr<SpkSegmentEvaluator>, IEvaluator>(m, "SpkSegmentEvaluator")
        .def("chebyshev_metadata", [](const SpkSegmentEvaluator& self) {
            if (self.data_type != 2 && self.data_type != 3) {
                throw std::invalid_argument("Chebyshev metadata requires an SPK type 2 or 3 segment");
            }
            return py::make_tuple(self.init, self.intlen, self.record_count,
                                  self.component_count, self.coefficient_count);
        })
        .def("chebyshev_records", [](const SpkSegmentEvaluator& self, int64_t first, int64_t count) {
            // Copy a bounded selection from this exact serving evaluator. Coefficients
            // are always returned in ascending degree order, independently of its
            // private storage order. No file is reopened and no mutable view escapes.
            if ((self.data_type != 2 && self.data_type != 3) || first < 0 ||
                count < 1 || count > 256 || static_cast<uint64_t>(first) >= self.record_count ||
                static_cast<uint64_t>(count) > self.record_count - static_cast<size_t>(first)) {
                throw std::invalid_argument("invalid bounded Chebyshev record selection");
            }
            py::tuple records(static_cast<size_t>(count));
            for (int64_t r = 0; r < count; ++r) {
                py::tuple components(self.component_count);
                for (size_t c = 0; c < self.component_count; ++c) {
                    py::tuple values(self.coefficient_count);
                    const size_t base = ((static_cast<size_t>(first + r) * self.component_count) + c) * self.coefficient_count;
                    for (size_t k = 0; k < self.coefficient_count; ++k) {
                        const size_t stored = self.coefficients_in_file_order ? k : self.coefficient_count - 1 - k;
                        values[k] = self.coefficients[base + stored];
                    }
                    components[c] = values;
                }
                records[static_cast<size_t>(r)] = components;
            }
            return records;
        }, py::arg("first"), py::arg("count") = 1)
        .def("position", [](const SpkSegmentEvaluator& self, double jd, double jd2) {
            double result[3];
            {
                py::gil_scoped_release release;
                self.position(jd, jd2, result);
            }
            py::list out;
            for (double value : result) out.append(value);
            return out;
        }, py::arg("jd"), py::arg("jd2") = 0.0)
        .def("position_and_velocity", [](const SpkSegmentEvaluator& self, double jd, double jd2) {
            double position[3];
            double velocity[3];
            {
                py::gil_scoped_release release;
                self.position_and_velocity(jd, jd2, position, velocity);
            }
            py::list pos_out;
            py::list vel_out;
            for (double value : position) pos_out.append(value);
            for (double value : velocity) vel_out.append(value);
            return py::make_tuple(pos_out, vel_out);
        }, py::arg("jd"), py::arg("jd2") = 0.0);

    py::class_<NativeSpkKernelHandle, std::shared_ptr<NativeSpkKernelHandle>>(m, "NativeSpkKernelHandle")
        .def("catalog", [](const NativeSpkKernelHandle& self) {
            return catalog_to_py_dict(self.catalog);
        })
        .def("load_segment_evaluator", [](NativeSpkKernelHandle& self, int32_t start_i, int32_t end_i, int32_t data_type) {
            std::shared_ptr<SpkSegmentEvaluator> evaluator;
            {
                py::gil_scoped_release release;
                evaluator = self.get_segment_evaluator(start_i, end_i, data_type);
            }
            return evaluator;
        })
        .def("batch_segment_position_and_velocity", [](NativeSpkKernelHandle& self, py::iterable specs, double jd, double jd2) {
            struct SegmentRequestSpec {
                int32_t start_i;
                int32_t end_i;
                int32_t data_type;
            };

            std::vector<SegmentRequestSpec> parsed_specs;
            for (py::handle spec_handle : specs) {
                auto spec = py::cast<py::tuple>(spec_handle);
                if (py::len(spec) != 3) {
                    throw std::runtime_error("segment spec must be a 3-tuple of (start_i, end_i, data_type)");
                }
                parsed_specs.push_back(SegmentRequestSpec{
                    py::cast<int32_t>(spec[0]),
                    py::cast<int32_t>(spec[1]),
                    py::cast<int32_t>(spec[2]),
                });
            }

            std::vector<std::array<double, 3>> positions(parsed_specs.size());
            std::vector<std::array<double, 3>> velocities(parsed_specs.size());
            {
                py::gil_scoped_release release;
                for (size_t index = 0; index < parsed_specs.size(); ++index) {
                    const SegmentRequestSpec& spec = parsed_specs[index];
                    self.segment_position_and_velocity(
                        spec.start_i,
                        spec.end_i,
                        spec.data_type,
                        jd,
                        jd2,
                        positions[index].data(),
                        velocities[index].data()
                    );
                }
            }

            py::list out;
            for (size_t index = 0; index < parsed_specs.size(); ++index) {
                py::list pos_out;
                py::list vel_out;
                for (double value : positions[index]) pos_out.append(value);
                for (double value : velocities[index]) vel_out.append(value);
                out.append(py::make_tuple(pos_out, vel_out));
            }
            return out;
        }, py::arg("specs"), py::arg("jd"), py::arg("jd2") = 0.0)
        .def("batch_segment_position_requests", [](NativeSpkKernelHandle& self, py::iterable requests) {
            struct SegmentPositionRequest {
                int32_t start_i;
                int32_t end_i;
                int32_t data_type;
                double jd;
                double jd2;
            };

            std::vector<SegmentPositionRequest> parsed_requests;
            for (py::handle request_handle : requests) {
                auto request = py::cast<py::tuple>(request_handle);
                if (py::len(request) != 4 && py::len(request) != 5) {
                    throw std::runtime_error(
                        "segment request must be a 4-tuple of (start_i, end_i, data_type, jd) "
                        "or 5-tuple of (start_i, end_i, data_type, jd, jd2)"
                    );
                }
                parsed_requests.push_back(SegmentPositionRequest{
                    py::cast<int32_t>(request[0]),
                    py::cast<int32_t>(request[1]),
                    py::cast<int32_t>(request[2]),
                    py::cast<double>(request[3]),
                    py::len(request) == 5 ? py::cast<double>(request[4]) : 0.0,
                });
            }

            std::vector<std::array<double, 3>> positions(parsed_requests.size());
            {
                py::gil_scoped_release release;
                for (size_t index = 0; index < parsed_requests.size(); ++index) {
                    const SegmentPositionRequest& request = parsed_requests[index];
                    self.segment_position(
                        request.start_i,
                        request.end_i,
                        request.data_type,
                        request.jd,
                        request.jd2,
                        positions[index].data()
                    );
                }
            }

            py::list out;
            for (const auto& position : positions) {
                py::list pos_out;
                for (double value : position) pos_out.append(value);
                out.append(pos_out);
            }
            return out;
        })
        .def("segment_position", [](NativeSpkKernelHandle& self, int32_t start_i, int32_t end_i, int32_t data_type, double jd, double jd2) {
            double result[3];
            {
                py::gil_scoped_release release;
                self.segment_position(start_i, end_i, data_type, jd, jd2, result);
            }
            py::list out;
            for (double value : result) out.append(value);
            return out;
        }, py::arg("start_i"), py::arg("end_i"), py::arg("data_type"), py::arg("jd"), py::arg("jd2") = 0.0)
        .def("segment_position_and_velocity", [](NativeSpkKernelHandle& self, int32_t start_i, int32_t end_i, int32_t data_type, double jd, double jd2) {
            double position[3];
            double velocity[3];
            {
                py::gil_scoped_release release;
                self.segment_position_and_velocity(start_i, end_i, data_type, jd, jd2, position, velocity);
            }
            py::list pos_out;
            py::list vel_out;
            for (double value : position) pos_out.append(value);
            for (double value : velocity) vel_out.append(value);
            return py::make_tuple(pos_out, vel_out);
        }, py::arg("start_i"), py::arg("end_i"), py::arg("data_type"), py::arg("jd"), py::arg("jd2") = 0.0)
        .def("close", [](NativeSpkKernelHandle& self) {
            py::gil_scoped_release release;
            self.close();
        })
        .def("segment_cache_size", [](NativeSpkKernelHandle& self) {
            py::gil_scoped_release release;
            return self.segment_cache_size();
        })
        .def("segment_cache_limit", [](NativeSpkKernelHandle& self) {
            return self.segment_cache_limit();
        });

    py::class_<NativePckKernelHandle, std::shared_ptr<NativePckKernelHandle>>(m, "NativePckKernelHandle")
        .def("catalog", [](const NativePckKernelHandle& self) {
            return pck_catalog_to_py_dict(self);
        })
        .def("coverage", [](const NativePckKernelHandle& self, int32_t frame_class_id) {
            const auto coverage = self.coverage_jd(frame_class_id);
            return py::make_tuple(coverage.first, coverage.second);
        }, py::arg("frame_class_id") = 31008)
        .def("euler_angles", [](NativePckKernelHandle& self, double jd_tdb, int32_t frame_class_id) {
            std::array<double, 3> angles;
            {
                py::gil_scoped_release release;
                angles = self.euler_angles(jd_tdb, frame_class_id);
            }
            return py::make_tuple(angles[0], angles[1], angles[2]);
        }, py::arg("jd_tdb"), py::arg("frame_class_id") = 31008)
        .def("rotation_matrix", [](NativePckKernelHandle& self, double jd_tdb, int32_t frame_class_id) {
            Mat3 matrix;
            {
                py::gil_scoped_release release;
                matrix = self.rotation_matrix(jd_tdb, frame_class_id);
            }
            return mat3_to_py_tuple(matrix);
        }, py::arg("jd_tdb"), py::arg("frame_class_id") = 31008)
        .def("close", [](NativePckKernelHandle& self) {
            py::gil_scoped_release release;
            self.close();
        })
        .def("segment_cache_size", [](const NativePckKernelHandle& self) {
            py::gil_scoped_release release;
            return self.segment_cache_size();
        });

    py::class_<NativePlanetaryEvaluator, std::shared_ptr<NativePlanetaryEvaluator>>(m, "NativePlanetaryEvaluator")
        .def(py::init<std::shared_ptr<NativeSpkKernelHandle>>())
        .def(
            "evaluate_all_planets_apparent_geocentric_ecliptic",
            [](const NativePlanetaryEvaluator& self,
               const std::vector<std::string>& bodies,
               py::iterable public_specs_src,
               py::dict body_specs_src,
               double epoch_tdb,
               double obliquity_deg,
               const py::sequence& rotation_matrix_src) {
                std::vector<NativePlanetaryEvaluator::SegmentSpec> public_specs;
                for (py::handle spec_handle : public_specs_src) {
                    py::tuple spec = py::cast<py::tuple>(spec_handle);
                    if (py::len(spec) != 3) {
                        throw std::runtime_error("public segment spec must be a 3-tuple");
                    }
                    public_specs.emplace_back(
                        py::cast<int32_t>(spec[0]),
                        py::cast<int32_t>(spec[1]),
                        py::cast<int32_t>(spec[2])
                    );
                }

                NativePlanetaryEvaluator::SegmentSpecMap body_specs;
                for (auto item : body_specs_src) {
                    const std::string body = py::cast<std::string>(item.first);
                    py::iterable specs_src = py::reinterpret_borrow<py::iterable>(item.second);
                    std::vector<NativePlanetaryEvaluator::SegmentSpec> specs;
                    for (py::handle spec_handle : specs_src) {
                        py::tuple spec = py::cast<py::tuple>(spec_handle);
                        if (py::len(spec) != 3) {
                            throw std::runtime_error("body segment spec must be a 3-tuple");
                        }
                        specs.emplace_back(
                            py::cast<int32_t>(spec[0]),
                            py::cast<int32_t>(spec[1]),
                            py::cast<int32_t>(spec[2])
                        );
                    }
                    body_specs.emplace(body, std::move(specs));
                }

                double matrix_vals[3][3];
                load_mat3(rotation_matrix_src, matrix_vals);
                Mat3 rotation_matrix = Mat3::identity();
                for (size_t i = 0; i < 3; ++i) {
                    for (size_t j = 0; j < 3; ++j) {
                        rotation_matrix[i][j] = matrix_vals[i][j];
                    }
                }

                std::vector<NativePlanetaryPayload> payloads;
                {
                    py::gil_scoped_release release;
                    payloads = self.evaluate_all_planets_apparent_geocentric_ecliptic(
                        bodies,
                        public_specs,
                        body_specs,
                        epoch_tdb,
                        obliquity_deg,
                        rotation_matrix
                    );
                }

                py::list out;
                for (const NativePlanetaryPayload& payload : payloads) {
                    out.append(py::make_tuple(
                        payload.name,
                        payload.longitude,
                        payload.latitude,
                        payload.distance,
                        payload.speed,
                        payload.retrograde
                    ));
                }
                return out;
            },
            py::arg("bodies"),
            py::arg("public_specs"),
            py::arg("body_specs"),
            py::arg("epoch_tdb"),
            py::arg("obliquity_deg"),
            py::arg("rotation_matrix")
        )
        .def(
            "evaluate_all_planets_apparent_with_speed",
            [](const NativePlanetaryEvaluator& self,
               const std::vector<std::string>& bodies,
               py::iterable public_specs_src,
               py::dict body_specs_src,
               double epoch_tdb,
               double obliquity_deg,
               const py::sequence& rotation_matrix_src,
               double rate_step_days,
               double rate_minus_epoch_tdb,
               double rate_minus_obliquity_deg,
               const py::sequence& rate_minus_rotation_matrix_src,
               double rate_plus_epoch_tdb,
               double rate_plus_obliquity_deg,
               const py::sequence& rate_plus_rotation_matrix_src) {
                std::vector<NativePlanetaryEvaluator::SegmentSpec> public_specs;
                for (py::handle spec_handle : public_specs_src) {
                    py::tuple spec = py::cast<py::tuple>(spec_handle);
                    if (py::len(spec) != 3) {
                        throw std::runtime_error("public segment spec must be a 3-tuple");
                    }
                    public_specs.emplace_back(
                        py::cast<int32_t>(spec[0]),
                        py::cast<int32_t>(spec[1]),
                        py::cast<int32_t>(spec[2])
                    );
                }

                NativePlanetaryEvaluator::SegmentSpecMap body_specs;
                for (auto item : body_specs_src) {
                    const std::string body = py::cast<std::string>(item.first);
                    py::iterable specs_src = py::reinterpret_borrow<py::iterable>(item.second);
                    std::vector<NativePlanetaryEvaluator::SegmentSpec> specs;
                    for (py::handle spec_handle : specs_src) {
                        py::tuple spec = py::cast<py::tuple>(spec_handle);
                        if (py::len(spec) != 3) {
                            throw std::runtime_error("body segment spec must be a 3-tuple");
                        }
                        specs.emplace_back(
                            py::cast<int32_t>(spec[0]),
                            py::cast<int32_t>(spec[1]),
                            py::cast<int32_t>(spec[2])
                        );
                    }
                    body_specs.emplace(body, std::move(specs));
                }

                auto parse_mat3 = [](const py::sequence& src) -> Mat3 {
                    double vals[3][3];
                    load_mat3(src, vals);
                    Mat3 m = Mat3::identity();
                    for (size_t i = 0; i < 3; ++i) {
                        for (size_t j = 0; j < 3; ++j) {
                            m[i][j] = vals[i][j];
                        }
                    }
                    return m;
                };

                Mat3 rot_mat = parse_mat3(rotation_matrix_src);
                Mat3 rate_minus_rot_mat = parse_mat3(rate_minus_rotation_matrix_src);
                Mat3 rate_plus_rot_mat = parse_mat3(rate_plus_rotation_matrix_src);

                std::vector<NativePlanetaryPayload> payloads;
                {
                    py::gil_scoped_release release;
                    payloads = self.evaluate_all_planets_apparent_with_speed(
                        bodies,
                        public_specs,
                        body_specs,
                        epoch_tdb,
                        obliquity_deg,
                        rot_mat,
                        rate_step_days,
                        rate_minus_epoch_tdb,
                        rate_minus_obliquity_deg,
                        rate_minus_rot_mat,
                        rate_plus_epoch_tdb,
                        rate_plus_obliquity_deg,
                        rate_plus_rot_mat
                    );
                }

                py::list out;
                for (const NativePlanetaryPayload& payload : payloads) {
                    out.append(py::make_tuple(
                        payload.name,
                        payload.longitude,
                        payload.latitude,
                        payload.distance,
                        payload.speed,
                        payload.retrograde
                    ));
                }
                return out;
            },
            py::arg("bodies"),
            py::arg("public_specs"),
            py::arg("body_specs"),
            py::arg("epoch_tdb"),
            py::arg("obliquity_deg"),
            py::arg("rotation_matrix"),
            py::arg("rate_step_days"),
            py::arg("rate_minus_epoch_tdb"),
            py::arg("rate_minus_obliquity_deg"),
            py::arg("rate_minus_rotation_matrix"),
            py::arg("rate_plus_epoch_tdb"),
            py::arg("rate_plus_obliquity_deg"),
            py::arg("rate_plus_rotation_matrix")
        );

    py::class_<Type13Evaluator, IEvaluator, std::shared_ptr<Type13Evaluator>>(m, "Type13Evaluator")
        .def(py::init<std::vector<double>, std::vector<double>, size_t>());
    m.attr("LagrangeEvaluator") = m.attr("Type13Evaluator");

    py::class_<RelativeEvaluator, IEvaluator, std::shared_ptr<RelativeEvaluator>>(m, "RelativeEvaluator")
        .def(py::init<std::shared_ptr<IEvaluator>, std::shared_ptr<IEvaluator>>());
    
    py::class_<SumEvaluator, IEvaluator, std::shared_ptr<SumEvaluator>>(m, "SumEvaluator")
        .def(py::init<std::shared_ptr<IEvaluator>, std::shared_ptr<IEvaluator>>());

    py::class_<NegateEvaluator, IEvaluator, std::shared_ptr<NegateEvaluator>>(m, "NegateEvaluator")
        .def(py::init<std::shared_ptr<IEvaluator>>());

    m.def("longitude_difference", [](std::shared_ptr<IEvaluator> t1, std::shared_ptr<IEvaluator> t2, std::shared_ptr<IEvaluator> obs, double jd) {
        return longitude_difference(*t1, *t2, *obs, jd);
    });

    m.def("longitude_difference_batch", [](std::shared_ptr<IEvaluator> t1, std::shared_ptr<IEvaluator> t2, std::shared_ptr<IEvaluator> obs, const std::vector<double>& jds) {
        std::vector<double> out(jds.size());
        {
            py::gil_scoped_release release;
            for (size_t i = 0; i < jds.size(); ++i) {
                out[i] = longitude_difference(*t1, *t2, *obs, jds[i]);
            }
        }
        return out;
    });

    m.def("ecliptic_longitude_batch", [](std::shared_ptr<IEvaluator> target, std::shared_ptr<IEvaluator> obs, const std::vector<double>& jds) {
        std::vector<double> out(jds.size());
        {
            py::gil_scoped_release release;
            NutationEpochCache cache;
            for (size_t i = 0; i < jds.size(); ++i) {
                double rt[6], ro[6];
                target->evaluate(jds[i], rt);
                obs->evaluate(jds[i], ro);
                Vec3 geo(rt[0] - ro[0], rt[1] - ro[1], rt[2] - ro[2]);
                const NutationResult nut = cache.evaluate(jds[i]);
                double lon = _true_ecliptic_longitude(geo, jds[i], nut);
                if (lon < 0.0) lon += 360.0;
                out[i] = lon;
            }
        }
        return out;
    });

    m.def("find_aspects_to_longitude", [](std::shared_ptr<IEvaluator> target, std::shared_ptr<IEvaluator> obs, double frozen_lon_deg, double aspect_deg, double a, double b, double dt) {
        auto f = [&](double jd) {
            double rt[6], ro[6];
            target->evaluate(jd, rt);
            obs->evaluate(jd, ro);
            Vec3 geo(rt[0] - ro[0], rt[1] - ro[1], rt[2] - ro[2]);
            NutationResult nut = nutation_2000r06(jd);
            double lon = _true_ecliptic_longitude(geo, jd, nut);
            double val = lon - frozen_lon_deg - aspect_deg;
            while (val > 180.0) val -= 360.0;
            while (val <= -180.0) val += 360.0;
            return val;
        };
        py::gil_scoped_release release;
        return find_roots(f, a, b, dt);
    });

    m.def("declination_batch", [](std::shared_ptr<IEvaluator> t, std::shared_ptr<IEvaluator> obs, const std::vector<double>& jds) {
        std::vector<double> out(jds.size());
        {
            py::gil_scoped_release release;
            for (size_t i = 0; i < jds.size(); ++i) {
                double r_t[6], r_o[6];
                t->evaluate(jds[i], r_t);
                obs->evaluate(jds[i], r_o);
                double x = r_t[0] - r_o[0];
                double y = r_t[1] - r_o[1];
                double z = r_t[2] - r_o[2];
                double dist = std::sqrt(x*x + y*y + z*z);
                out[i] = std::asin(z / dist) * 180.0 / 3.14159265358979323846;
            }
        }
        return out;
    });

    m.def("angular_separation", [](std::shared_ptr<IEvaluator> t1, std::shared_ptr<IEvaluator> t2, std::shared_ptr<IEvaluator> obs, double jd) {
        return angular_separation(*t1, *t2, *obs, jd);
    });

    m.def("find_conjunctions", [](std::shared_ptr<IEvaluator> t1, std::shared_ptr<IEvaluator> t2, std::shared_ptr<IEvaluator> obs, double a, double b, double dt) {
        auto f = [&](double jd) { return longitude_difference(*t1, *t2, *obs, jd); };
        py::gil_scoped_release release;
        return find_roots(f, a, b, dt);
    });

    m.def("find_aspects", [](std::shared_ptr<IEvaluator> t1, std::shared_ptr<IEvaluator> t2, std::shared_ptr<IEvaluator> obs, double aspect_deg, double a, double b, double dt) {
        auto f = [&](double jd) { 
            double diff = longitude_difference(*t1, *t2, *obs, jd);
            double val = diff - aspect_deg;
            while (val > 180.0) val -= 360.0;
            while (val <= -180.0) val += 360.0;
            return val;
        };
        py::gil_scoped_release release;
        return find_roots(f, a, b, dt);
    });

    m.def("find_stations", [](std::shared_ptr<IEvaluator> target, std::shared_ptr<IEvaluator> obs, double a, double b, double dt) {
        py::gil_scoped_release release;
        return find_stations(*target, *obs, a, b, dt);
    });

    m.def("find_ingresses", [](std::shared_ptr<IEvaluator> target, std::shared_ptr<IEvaluator> obs, double a, double b, double dt) {
        py::gil_scoped_release release;
        return find_ingresses(*target, *obs, a, b, dt);
    });

    py::class_<OccultationEvent>(m, "OccultationEvent")
        .def_readonly("t_mid", &OccultationEvent::t_mid)
        .def_readonly("separation_min", &OccultationEvent::separation_min)
        .def_readonly("t_start", &OccultationEvent::t_start)
        .def_readonly("t_end", &OccultationEvent::t_end)
        .def_readonly("is_total", &OccultationEvent::is_total);

    m.def("find_occultations",
        [](const IEvaluator& target1, double r1_km,
           const IEvaluator& target2, double r2_km,
           const IEvaluator& observer, double a, double b, double dt) {
            py::gil_scoped_release release;
            return find_occultations(target1, r1_km, target2, r2_km, observer, a, b, dt);
        },
        py::arg("target1"), py::arg("r1_km"),
        py::arg("target2"), py::arg("r2_km"),
        py::arg("observer"), py::arg("a"), py::arg("b"), py::arg("dt"));

    py::enum_<EventType>(m, "EventType")
        .value("STATION", EventType::STATION)
        .value("INGRESS", EventType::INGRESS)
        .value("OCCULTATION", EventType::OCCULTATION)
        .export_values();

    py::class_<SearchResult>(m, "SearchResult")
        .def_readonly("type", &SearchResult::type)
        .def_readonly("jd", &SearchResult::jd)
        .def_readonly("description", &SearchResult::description)
        .def_readonly("value", &SearchResult::value);

    py::class_<SearchPool>(m, "SearchPool")
        .def(py::init<>())
        .def("add_station_task", &SearchPool::add_station_task)
        .def("add_ingress_task", &SearchPool::add_ingress_task)
        .def("add_occultation_task", &SearchPool::add_occultation_task)
        .def("run", [](SearchPool& self, double a, double b, double dt) {
            py::gil_scoped_release release;
            return self.run(a, b, dt);
        }, py::arg("a"), py::arg("b"), py::arg("dt") = 0.5);


    // --- Geometry ---
    py::class_<Vec3>(m, "Vec3")
        .def(py::init<double, double, double>())
        .def_property("x", [](Vec3& v) { return v[0]; }, [](Vec3& v, double val) { v[0] = val; })
        .def_property("y", [](Vec3& v) { return v[1]; }, [](Vec3& v, double val) { v[1] = val; })
        .def_property("z", [](Vec3& v) { return v[2]; }, [](Vec3& v, double val) { v[2] = val; })
        .def("norm", &Vec3::norm)
        .def("unit", &Vec3::unit);

    py::class_<Mat3>(m, "Mat3")
        .def(py::init<>())
        .def(py::init<std::array<std::array<double, 3>, 3>>())
        .def_readwrite("data", &Mat3::data)
        .def_static("identity", &Mat3::identity)
        .def_static("mul_vec", py::overload_cast<const Mat3&, const Vec3&>(&Mat3::mul))
        .def_static("mul_mat", py::overload_cast<const Mat3&, const Mat3&>(&Mat3::mul))
        .def_static("compose", &Mat3::compose)
        .def("transpose", &Mat3::transpose)
        .def("determinant", &Mat3::determinant)
        .def("inverse", &Mat3::inverse)
        .def("is_orthonormal", &Mat3::is_orthonormal, py::arg("epsilon") = 1e-12)
        .def("orthonormalize", &Mat3::orthonormalize)
        .def_static("rot_x", &Mat3::rot_x)
        .def_static("rot_y", &Mat3::rot_y)
        .def_static("rot_z", &Mat3::rot_z);

    // --- Julian ---
    m.def("julian_day", &julian_day, py::arg("year"), py::arg("month"), py::arg("day"), py::arg("hour") = 0.0);
    m.def("calendar_from_jd", &calendar_from_jd, py::arg("jd"));
    m.def("load_naif_lsk", &load_naif_lsk_py, py::arg("path"));
    m.def("jd_utc_to_et_seconds_past_j2000", &jd_utc_to_et_seconds_past_j2000_py, py::arg("jd_utc"));
    m.def("earth_rotation_angle", &earth_rotation_angle, py::arg("jd_ut"));
    m.def("greenwich_mean_sidereal_time", &greenwich_mean_sidereal_time, py::arg("jd_ut"));
    m.def("apparent_sidereal_time", &apparent_sidereal_time, py::arg("jd_ut"), py::arg("nutation_longitude"), py::arg("obliquity"));
    m.def("set_nutation_2000a_tables", &set_nutation_2000a_tables_py,
        py::arg("ls_terms"), py::arg("pl_terms"), py::arg("ls_j0_count"), py::arg("pl_j0_count"));
    m.def("nutation_2000a", &nutation_2000a_py, py::arg("jd_tt"));
    m.def("_runtime_avx2_available", &runtime_avx2_available);
    m.def("nutation_2000r06_series_ready", &nutation_2000r06_series_ready);
    py::class_<NutationEpochCache>(m, "NutationEpochCache")
        .def(py::init<>())
        .def_property_readonly("size", &NutationEpochCache::size);
    m.def("rotation_matrix_multiply", &rotation_matrix_multiply_py, py::arg("a"), py::arg("b"));
    m.def("rotation_matrix_apply", &rotation_matrix_apply_py, py::arg("matrix"), py::arg("vector"));
    m.def("apply_aberration_velocity", &apply_aberration_velocity_py, py::arg("xyz"), py::arg("velocity"));
    m.def("apply_frame_bias", &apply_frame_bias_py, py::arg("xyz"));

    // --- Houses & Local Angles ---
    m.def("mc_from_armc", &mc_from_armc, py::arg("armc"), py::arg("obliquity"));
    m.def("asc_from_armc", &asc_from_armc, py::arg("armc"), py::arg("obliquity"), py::arg("lat"));
    m.def("vertex_from_armc", &vertex_from_armc, py::arg("armc"), py::arg("obliquity"), py::arg("lat"));
    m.def("east_point_from_armc", &east_point_from_armc, py::arg("armc"), py::arg("obliquity"));
    m.def("reduce_local_angles", [](double jd_ut, double jd_tt, double lat, double lon) {
        py::gil_scoped_release release;
        return reduce_local_angles(jd_ut, jd_tt, lat, lon);
    }, py::arg("jd_ut"), py::arg("jd_tt"), py::arg("lat"), py::arg("lon"));
    m.def("calculate_houses_cusps", [](double armc, double obliquity, double lat, double asc, double mc, const std::string& system) {
        if (system.empty()) {
            throw std::invalid_argument("House system code cannot be empty");
        }
        std::array<double, 12> cusps;
        {
            py::gil_scoped_release release;
            cusps = calculate_houses_cusps(armc, obliquity, lat, asc, mc, system[0]);
        }
        py::list out;
        for (double c : cusps) {
            out.append(c);
        }
        return out;
    }, py::arg("armc"), py::arg("obliquity"), py::arg("lat"), py::arg("asc"), py::arg("mc"), py::arg("system"));

    // --- Primary Directions (Under-Pole Arcs) ---
    py::class_<NativeSpeculumPoint>(m, "NativeSpeculumPoint")
        .def(py::init<std::string, double, double, double, double, double, double, double, bool, double, bool>(),
             py::arg("name") = "", py::arg("lon") = 0.0, py::arg("lat") = 0.0,
             py::arg("ra") = 0.0, py::arg("dec") = 0.0, py::arg("ha") = 0.0,
             py::arg("dsa") = 0.0, py::arg("nsa") = 0.0, py::arg("upper") = true,
             py::arg("f") = 0.0, py::arg("is_eastern") = false)
        .def_readwrite("name", &NativeSpeculumPoint::name)
        .def_readwrite("lon", &NativeSpeculumPoint::lon)
        .def_readwrite("lat", &NativeSpeculumPoint::lat)
        .def_readwrite("ra", &NativeSpeculumPoint::ra)
        .def_readwrite("dec", &NativeSpeculumPoint::dec)
        .def_readwrite("ha", &NativeSpeculumPoint::ha)
        .def_readwrite("dsa", &NativeSpeculumPoint::dsa)
        .def_readwrite("nsa", &NativeSpeculumPoint::nsa)
        .def_readwrite("upper", &NativeSpeculumPoint::upper)
        .def_readwrite("f", &NativeSpeculumPoint::f)
        .def_readwrite("is_eastern", &NativeSpeculumPoint::is_eastern);

    m.def("campanus_regio_sin_zenith_distance", &campanus_regio_sin_zenith_distance,
          py::arg("dec"), py::arg("ha"), py::arg("geo_lat"),
          "Branch-independent sine of zenith distance for Campanus/Regiomontanus.");
    m.def("regiomontanus_pole_height", &regiomontanus_pole_height,
          py::arg("dec"), py::arg("ha"), py::arg("geo_lat"),
          "Pole height in degrees under the Regiomontanus circle of position.");
    m.def("topocentric_pole_height", &topocentric_pole_height,
          py::arg("ha"), py::arg("dsa"), py::arg("nsa"), py::arg("upper"), py::arg("geo_lat"),
          "Pole height in degrees under the Topocentric proportional semi-arc law.");
    m.def("under_pole_w", &under_pole_w,
          py::arg("ra"), py::arg("dec"), py::arg("pole_deg"), py::arg("is_eastern"),
          "Oblique ascension/descension under pole height in degrees [0, 360).");
    m.def("under_pole_arc_native", &under_pole_arc_native,
          py::arg("sig_ra"), py::arg("sig_dec"), py::arg("sig_is_eastern"),
          py::arg("prom_ra"), py::arg("prom_dec"), py::arg("pole_deg"),
          "Direct arc under pole in degrees [0, 360).");
    m.def("regiomontanus_under_pole_arc", &regiomontanus_under_pole_arc,
          py::arg("sig_ra"), py::arg("sig_dec"), py::arg("sig_ha"), py::arg("sig_is_eastern"),
          py::arg("prom_ra"), py::arg("prom_dec"), py::arg("geo_lat"),
          "Regiomontanus under-pole direct arc in degrees [0, 360).");
    m.def("topocentric_under_pole_arc", &topocentric_under_pole_arc,
          py::arg("sig_ra"), py::arg("sig_dec"), py::arg("sig_ha"), py::arg("sig_dsa"), py::arg("sig_nsa"),
          py::arg("sig_upper"), py::arg("sig_is_eastern"),
          py::arg("prom_ra"), py::arg("prom_dec"), py::arg("geo_lat"),
          "Topocentric under-pole direct arc in degrees [0, 360).");
    m.def("compute_under_pole_pair_arcs", [](py::handle sig_obj, py::handle prom_obj, double geo_lat, const std::string& method_str) {
        if (method_str.empty()) throw std::invalid_argument("Primary direction method code cannot be empty");
        char method_code = method_str[0];
        auto extract_pt = [](py::handle item) {
            if (py::isinstance<NativeSpeculumPoint>(item)) {
                return item.cast<NativeSpeculumPoint>();
            }
            NativeSpeculumPoint pt;
            if (py::hasattr(item, "name")) pt.name = item.attr("name").cast<std::string>();
            pt.ra = item.attr("ra").cast<double>();
            pt.dec = item.attr("dec").cast<double>();
            pt.ha = item.attr("ha").cast<double>();
            if (py::hasattr(item, "dsa") && !item.attr("dsa").is_none()) pt.dsa = item.attr("dsa").cast<double>();
            if (py::hasattr(item, "nsa") && !item.attr("nsa").is_none()) pt.nsa = item.attr("nsa").cast<double>();
            if (py::hasattr(item, "upper") && !item.attr("upper").is_none()) pt.upper = item.attr("upper").cast<bool>();
            if (py::hasattr(item, "f") && !item.attr("f").is_none()) pt.f = item.attr("f").cast<double>();
            if (py::hasattr(item, "is_eastern")) {
                pt.is_eastern = item.attr("is_eastern").cast<bool>();
            } else {
                pt.is_eastern = (pt.ha < 0.0);
            }
            return pt;
        };
        NativeSpeculumPoint sig = extract_pt(sig_obj);
        NativeSpeculumPoint prom = extract_pt(prom_obj);
        return compute_under_pole_pair_arcs(sig, prom, geo_lat, method_code);
    }, py::arg("sig"), py::arg("prom"), py::arg("geo_lat"), py::arg("method"),
    "Compute traditional (direct, converse) primary direction under-pole arcs for a single pair.");
    m.def("compute_under_pole_arcs_matrix", [](const py::sequence& points_seq, double geo_lat, const std::string& method_str) {
        if (method_str.empty()) throw std::invalid_argument("Primary direction method code cannot be empty");
        char method_code = method_str[0];
        std::vector<NativeSpeculumPoint> points;
        points.reserve(py::len(points_seq));
        for (auto item_handle : points_seq) {
            py::object item = py::reinterpret_borrow<py::object>(item_handle);
            if (py::isinstance<NativeSpeculumPoint>(item)) {
                points.push_back(item.cast<NativeSpeculumPoint>());
            } else {
                NativeSpeculumPoint pt;
                if (py::hasattr(item, "name")) pt.name = item.attr("name").cast<std::string>();
                pt.ra = item.attr("ra").cast<double>();
                pt.dec = item.attr("dec").cast<double>();
                pt.ha = item.attr("ha").cast<double>();
                if (py::hasattr(item, "dsa") && !item.attr("dsa").is_none()) pt.dsa = item.attr("dsa").cast<double>();
                if (py::hasattr(item, "nsa") && !item.attr("nsa").is_none()) pt.nsa = item.attr("nsa").cast<double>();
                if (py::hasattr(item, "upper") && !item.attr("upper").is_none()) pt.upper = item.attr("upper").cast<bool>();
                if (py::hasattr(item, "f") && !item.attr("f").is_none()) pt.f = item.attr("f").cast<double>();
                if (py::hasattr(item, "is_eastern")) {
                    pt.is_eastern = item.attr("is_eastern").cast<bool>();
                } else {
                    pt.is_eastern = (pt.ha < 0.0);
                }
                points.push_back(pt);
            }
        }
        std::vector<std::vector<std::pair<double, double>>> matrix;
        {
            py::gil_scoped_release release;
            matrix = compute_under_pole_arcs_matrix(points, geo_lat, method_code);
        }
        return matrix;
    }, py::arg("points"), py::arg("geo_lat"), py::arg("method"),
    "Compute NxN matrix of (direct, converse) primary direction under-pole arcs with GIL released.");

    m.def("placidian_required_ha", &placidian_required_ha,
          py::arg("f"), py::arg("dsa"), py::arg("nsa"),
          "Required hour angle for temporal fraction under Placidian semi-arcs.");
    m.def("placidian_mundane_arc", &placidian_mundane_arc,
          py::arg("sig_f"), py::arg("prom_ha"), py::arg("prom_dsa"), py::arg("prom_nsa"),
          "Direct Placidian mundane primary direction arc in degrees [0, 360).");
    m.def("compute_placidian_pair_arcs", [](py::handle sig_obj, py::handle prom_obj, const std::string& converse_mode_str) {
        char conv_mode = converse_mode_str.empty() ? 'T' : converse_mode_str[0];
        auto extract_pt = [](py::handle item) {
            if (py::isinstance<NativeSpeculumPoint>(item)) {
                return item.cast<NativeSpeculumPoint>();
            }
            NativeSpeculumPoint pt;
            if (py::hasattr(item, "name")) pt.name = item.attr("name").cast<std::string>();
            pt.ra = item.attr("ra").cast<double>();
            pt.dec = item.attr("dec").cast<double>();
            pt.ha = item.attr("ha").cast<double>();
            if (py::hasattr(item, "dsa") && !item.attr("dsa").is_none()) pt.dsa = item.attr("dsa").cast<double>();
            if (py::hasattr(item, "nsa") && !item.attr("nsa").is_none()) pt.nsa = item.attr("nsa").cast<double>();
            if (py::hasattr(item, "upper") && !item.attr("upper").is_none()) pt.upper = item.attr("upper").cast<bool>();
            if (py::hasattr(item, "f") && !item.attr("f").is_none()) pt.f = item.attr("f").cast<double>();
            if (py::hasattr(item, "is_eastern")) {
                pt.is_eastern = item.attr("is_eastern").cast<bool>();
            } else {
                pt.is_eastern = (pt.ha < 0.0);
            }
            return pt;
        };
        NativeSpeculumPoint sig = extract_pt(sig_obj);
        NativeSpeculumPoint prom = extract_pt(prom_obj);
        return compute_placidian_pair_arcs(sig, prom, conv_mode);
    }, py::arg("sig"), py::arg("prom"), py::arg("converse_mode") = "T",
    "Compute Placidian (direct, converse) primary direction arcs for a single pair.");

    m.def("compute_placidian_arcs_matrix", [](const py::sequence& points_seq, const std::string& converse_mode_str) {
        char conv_mode = converse_mode_str.empty() ? 'T' : converse_mode_str[0];
        std::vector<NativeSpeculumPoint> points;
        points.reserve(py::len(points_seq));
        for (auto item_handle : points_seq) {
            py::object item = py::reinterpret_borrow<py::object>(item_handle);
            if (py::isinstance<NativeSpeculumPoint>(item)) {
                points.push_back(item.cast<NativeSpeculumPoint>());
            } else {
                NativeSpeculumPoint pt;
                if (py::hasattr(item, "name")) pt.name = item.attr("name").cast<std::string>();
                pt.ra = item.attr("ra").cast<double>();
                pt.dec = item.attr("dec").cast<double>();
                pt.ha = item.attr("ha").cast<double>();
                if (py::hasattr(item, "dsa") && !item.attr("dsa").is_none()) pt.dsa = item.attr("dsa").cast<double>();
                if (py::hasattr(item, "nsa") && !item.attr("nsa").is_none()) pt.nsa = item.attr("nsa").cast<double>();
                if (py::hasattr(item, "upper") && !item.attr("upper").is_none()) pt.upper = item.attr("upper").cast<bool>();
                if (py::hasattr(item, "f") && !item.attr("f").is_none()) pt.f = item.attr("f").cast<double>();
                if (py::hasattr(item, "is_eastern")) {
                    pt.is_eastern = item.attr("is_eastern").cast<bool>();
                } else {
                    pt.is_eastern = (pt.ha < 0.0);
                }
                points.push_back(pt);
            }
        }
        std::vector<std::vector<std::pair<double, double>>> matrix;
        {
            py::gil_scoped_release release;
            matrix = compute_placidian_arcs_matrix(points, conv_mode);
        }
        return matrix;
    }, py::arg("points"), py::arg("converse_mode") = "T",
    "Compute NxN matrix of Placidian (direct, converse) primary direction arcs with GIL released.");


    // --- Interpolation ---
    m.def("horner", [](const py::sequence& coeffs, double x) {
        std::vector<double> coeff_vec = load_double_vector(coeffs, "Horner coefficients");
        double res = 0.0;
        for (int i = static_cast<int>(coeff_vec.size()) - 1; i >= 0; --i) {
            res = res * x + coeff_vec[static_cast<size_t>(i)];
        }
        return res;
    }, py::arg("coeffs"), py::arg("x"));
    m.def("lagrange_interpolate", [](const py::sequence& x_pts, const py::sequence& y_pts, double x) {
        std::vector<double> x_vec = load_double_vector(x_pts, "Lagrange x points");
        std::vector<double> y_vec = load_double_vector(y_pts, "Lagrange y points");
        if (x_vec.size() != y_vec.size()) throw std::runtime_error("Lagrange points must have same size");
        return lagrange_interpolate(x_vec.data(), y_vec.data(), x_vec.size(), x);
    }, py::arg("x_pts"), py::arg("y_pts"), py::arg("x"));
    
    m.def("spk_chebyshev_record", &spk_chebyshev_record_py, py::arg("coeff_record"), py::arg("s"));
    m.def("spk_chebyshev_record_with_derivative", &spk_chebyshev_record_with_derivative_py, py::arg("coeff_record"), py::arg("s"), py::arg("derivative_scale"));
    
    m.def("spk_chebyshev_series_record", &spk_chebyshev_series_record_py, py::arg("coefficients"), py::arg("record_index"), py::arg("s"));
    m.def("spk_chebyshev_series_record_with_derivative", &spk_chebyshev_series_record_with_derivative_py, py::arg("coefficients"), py::arg("record_index"), py::arg("s"), py::arg("derivative_scale"));
    
    m.def("spk_chebyshev_series_bulk_evaluate", &spk_chebyshev_series_bulk_evaluate_py, 
          py::arg("coeff_series"), py::arg("record_indices"), py::arg("s_values"), py::arg("derivative_scales"), py::arg("need_rates"));

    m.def("spk_type13_record", [](const py::sequence& epochs, const py::sequence& states, size_t window_size, double jd) {
        std::vector<double> epochs_vec = load_double_vector(epochs, "Type13 epochs");
        DenseMatrixView states_view = load_double_matrix(states, "Type13 states");
        if (states_view.rows != 6 || states_view.cols != epochs_vec.size()) {
            throw std::invalid_argument(
                "Type13 states must have shape (6, epoch_count)"
            );
        }
        validate_type13_payload(
            epochs_vec, states_view.values, window_size, "spk_type13_record"
        );
        std::vector<double> result(6);
        spk_type13_record_inplace(
            epochs_vec.data(),
            states_view.values.data(),
            epochs_vec.size(),
            window_size,
            jd,
            result.data()
        );
        return vector_to_py_list(result);
    }, py::arg("epochs"), py::arg("states"), py::arg("window_size"), py::arg("jd"));

    m.def("read_daf_catalog", &read_daf_catalog_py, py::arg("path"));
    m.def("open_spk_kernel", &open_spk_kernel_py, py::arg("path"));
    m.def("open_pck_kernel", &open_pck_kernel_py, py::arg("path"));
    m.def("read_spk_chebyshev_segment_payload", &read_spk_chebyshev_segment_payload_py, py::arg("path"), py::arg("start_i"), py::arg("end_i"), py::arg("little_endian"), py::arg("data_type"));
    m.def("load_spk_segment_evaluator", &load_spk_segment_evaluator_py, py::arg("path"), py::arg("start_i"), py::arg("end_i"), py::arg("little_endian"), py::arg("data_type"));
    m.def("read_spk_type13_segment_payload", &read_spk_type13_segment_payload_py, py::arg("path"), py::arg("start_i"), py::arg("end_i"), py::arg("little_endian"));

    // --- Solvers ---
    m.def("brent_root", &brent_root, py::arg("f"), py::arg("a"), py::arg("b"), py::arg("tol") = 1e-12, py::arg("max_iter") = 100);
    m.def("brent_minimize", &brent_minimize, py::arg("f"), py::arg("a"), py::arg("b"), py::arg("tol") = 1e-12, py::arg("max_iter") = 100);
    m.def("newton_safe", &newton_safe, py::arg("f"), py::arg("df"), py::arg("a"), py::arg("b"), py::arg("tol") = 1e-12, py::arg("max_iter") = 100);
    m.def("find_roots", &find_roots, py::arg("f"), py::arg("a"), py::arg("b"), py::arg("dt"), py::arg("tol") = 1e-12);
    m.def("find_extrema", &find_extrema, py::arg("f"), py::arg("a"), py::arg("b"), py::arg("dt"), py::arg("tol") = 1e-12);
    m.def("solve_light_time", &solve_light_time, py::arg("target_ephemeris"), py::arg("observer_pos"), py::arg("t_obs"), py::arg("initial_tau") = 0.0, py::arg("tol") = 1e-12, py::arg("max_iter") = 10);

    // --- Events & Separation ---
    m.def("angular_separation", py::overload_cast<const IEvaluator&, const IEvaluator&, const IEvaluator&, double>(&angular_separation));
    m.def("angular_separation", py::overload_cast<const Vec3&, const Vec3&>(&angular_separation));
    
    m.def("find_solar_eclipses", [](std::shared_ptr<IEvaluator> sun, std::shared_ptr<IEvaluator> moon, double jd_start, double jd_end, double r_sun_km, double r_moon_km, double dt_days) {
        py::gil_scoped_release release;
        return find_solar_eclipses(sun, moon, jd_start, jd_end, r_sun_km, r_moon_km, dt_days);
    }, py::arg("sun"), py::arg("moon"), py::arg("jd_start"), py::arg("jd_end"), py::arg("r_sun_km"), py::arg("r_moon_km"), py::arg("dt_days"));
    m.def("find_lunar_eclipses", [](std::shared_ptr<IEvaluator> sun, std::shared_ptr<IEvaluator> moon, double jd_start, double jd_end, double r_sun_km, double r_moon_km, double r_earth_km, double dt_days) {
        py::gil_scoped_release release;
        return find_lunar_eclipses(sun, moon, jd_start, jd_end, r_sun_km, r_moon_km, r_earth_km, dt_days);
    }, py::arg("sun"), py::arg("moon"), py::arg("jd_start"), py::arg("jd_end"), py::arg("r_sun_km"), py::arg("r_moon_km"), py::arg("r_earth_km"), py::arg("dt_days"));
    m.def("find_solar_syzygy_candidates", [](std::shared_ptr<IEvaluator> sun, std::shared_ptr<IEvaluator> moon, double jd_start_tt, double jd_end_tt, double max_separation_deg, double step_days) {
        py::gil_scoped_release release;
        return find_solar_syzygy_candidates(sun, moon, jd_start_tt, jd_end_tt, max_separation_deg, step_days);
    }, py::arg("sun"), py::arg("moon"), py::arg("jd_start_tt"), py::arg("jd_end_tt"), py::arg("max_separation_deg") = 3.0, py::arg("step_days") = 2.0);
    m.def("find_lunar_syzygy_candidates", [](std::shared_ptr<IEvaluator> sun, std::shared_ptr<IEvaluator> moon, double jd_start_tt, double jd_end_tt, double max_separation_deg, double step_days) {
        py::gil_scoped_release release;
        return find_lunar_syzygy_candidates(sun, moon, jd_start_tt, jd_end_tt, max_separation_deg, step_days);
    }, py::arg("sun"), py::arg("moon"), py::arg("jd_start_tt"), py::arg("jd_end_tt"), py::arg("max_separation_deg") = 3.0, py::arg("step_days") = 2.0);
    
    py::class_<Event>(m, "Event")
        .def_readwrite("type", &Event::type)
        .def_readwrite("t_mid", &Event::t_mid)
        .def_readwrite("t_start", &Event::t_start)
        .def_readwrite("t_end", &Event::t_end)
        .def_readwrite("value", &Event::value)
        .def_readwrite("description", &Event::description);

    // --- Coordinates ---
    m.def("radec_to_vec3", &radec_to_vec3, py::arg("ra_deg"), py::arg("dec_deg"), py::arg("dist") = 1.0);
    m.def("vec3_to_radec", &vec3_to_radec, py::arg("v"));
    m.def("lonlat_to_vec3", &lonlat_to_vec3, py::arg("lon_deg"), py::arg("lat_deg"), py::arg("dist") = 1.0);
    m.def("vec3_to_lonlat", &vec3_to_lonlat, py::arg("v"));
    m.def("vec3_to_lonlat_signed", &vec3_to_lonlat_signed, py::arg("v"));
    m.def("geodetic_to_cartesian_wgs84", &geodetic_to_cartesian_wgs84,
        py::arg("lon_deg"), py::arg("lat_deg"), py::arg("elevation_m") = 0.0);
    m.def("ecliptic_to_equatorial", &ecliptic_to_equatorial, py::arg("lon_deg"), py::arg("lat_deg"), py::arg("obliquity_deg"));
    m.def("equatorial_to_ecliptic", &equatorial_to_ecliptic, py::arg("ra_deg"), py::arg("dec_deg"), py::arg("obliquity_deg"));

    // --- Topocentric ---
    py::class_<TopocentricEvaluator, IEvaluator, std::shared_ptr<TopocentricEvaluator>>(m, "TopocentricEvaluator")
        .def(py::init<std::shared_ptr<IEvaluator>, double, double, double>());

    py::class_<FixedStarEvaluator, IEvaluator, std::shared_ptr<FixedStarEvaluator>>(m, "FixedStarEvaluator")
        .def(py::init<double, double, double, double, double, double>(),
             py::arg("ra_deg"), py::arg("dec_deg"), py::arg("pmra_mas"), py::arg("pmdec_mas"), py::arg("parallax"), py::arg("rv"));

    // Legacy observer-grid cartography is intentionally not admitted here.
    // Its magnitude/raster products are distinct from the public footprint.

    // --- LOLA (Lunar Orbiter Laser Altimeter) ---
    
    py::class_<lola::SphericalCoords>(m, "SphericalCoords")
        .def_readwrite("lon_deg", &lola::SphericalCoords::lon_deg)
        .def_readwrite("lat_deg", &lola::SphericalCoords::lat_deg)
        .def_readwrite("radius_km", &lola::SphericalCoords::radius_km);

    py::class_<CubicSmoothingSpline>(m, "CubicSmoothingSpline")
        .def(py::init<const std::vector<double>&, const std::vector<double>&, double>(),
             py::arg("x"), py::arg("y"), py::arg("p"))
        .def("evaluate", &CubicSmoothingSpline::evaluate, py::arg("target"))
        .def("get_knots", &CubicSmoothingSpline::get_knots)
        .def("get_y_hat", &CubicSmoothingSpline::get_y_hat)
        .def("get_M", &CubicSmoothingSpline::get_M);

    py::class_<lola::SkyPlaneProjection>(m, "SkyPlaneProjection")
        .def_readwrite("east_km", &lola::SkyPlaneProjection::east_km)
        .def_readwrite("north_km", &lola::SkyPlaneProjection::north_km)
        .def_readwrite("radius_km", &lola::SkyPlaneProjection::radius_km)
        .def_readwrite("pa_deg", &lola::SkyPlaneProjection::pa_deg);

    py::class_<lola::SkyPlanePaBinMaxima>(m, "SkyPlanePaBinMaxima")
        .def_readwrite("bin_indices", &lola::SkyPlanePaBinMaxima::bin_indices)
        .def_readwrite(
            "bin_centers_unwrapped_deg",
            &lola::SkyPlanePaBinMaxima::bin_centers_unwrapped_deg)
        .def_readwrite("radii_km", &lola::SkyPlanePaBinMaxima::radii_km)
        .def_readwrite("point_indices", &lola::SkyPlanePaBinMaxima::point_indices)
        .def_readwrite("bin_count", &lola::SkyPlanePaBinMaxima::bin_count)
        .def_readwrite(
            "admitted_source_point_count",
            &lola::SkyPlanePaBinMaxima::admitted_source_point_count);

    py::class_<lola::MaxPerBin>(m, "MaxPerBin")
        .def_readwrite("bins", &lola::MaxPerBin::bins)
        .def_readwrite("radii_km", &lola::MaxPerBin::radii_km)
        .def_readwrite("point_indices", &lola::MaxPerBin::point_indices);

    py::class_<lola::Point2D>(m, "Point2D")
        .def(py::init<double, double>())
        .def_readwrite("x", &lola::Point2D::x)
        .def_readwrite("y", &lola::Point2D::y);

    py::class_<lola::LolaPointCloud>(m, "LolaPointCloud")
        .def(py::init<const std::vector<double>&, const std::vector<double>&, const std::vector<double>&>())
        .def(py::init<>())
        .def("size", &lola::LolaPointCloud::size)
        .def("x_data", [](const lola::LolaPointCloud& self) { return reinterpret_cast<uintptr_t>(self.x_data()); })
        .def("y_data", [](const lola::LolaPointCloud& self) { return reinterpret_cast<uintptr_t>(self.y_data()); })
        .def("z_data", [](const lola::LolaPointCloud& self) { return reinterpret_cast<uintptr_t>(self.z_data()); })
        .def("get_x", &lola::LolaPointCloud::x_list)
        .def("get_y", &lola::LolaPointCloud::y_list)
        .def("get_z", &lola::LolaPointCloud::z_list)
        .def("filter_by_visibility", &lola::LolaPointCloud::filter_by_visibility)
        .def("filter_by_position_angle", &lola::LolaPointCloud::filter_by_position_angle)
        .def("filter_by_radius", &lola::LolaPointCloud::filter_by_radius)
        .def("filter_combined", &lola::LolaPointCloud::filter_combined)
        .def("to_spherical", &lola::LolaPointCloud::to_spherical)
        .def(
            "project_to_sky_plane",
            &lola::LolaPointCloud::project_to_sky_plane,
            py::arg("observer_dir"),
            py::arg("sky_east"),
            py::arg("sky_north"),
            py::arg("observer_distance_km") =
                std::numeric_limits<double>::infinity())
        .def(
            "project_max_radius_per_pa_bin",
            &lola::LolaPointCloud::project_max_radius_per_pa_bin,
            py::arg("observer_dir"),
            py::arg("sky_east"),
            py::arg("sky_north"),
            py::arg("pa_lower_unwrapped_deg"),
            py::arg("pa_upper_unwrapped_deg"),
            py::arg("bin_width_deg"),
            py::arg("observer_distance_km") =
                std::numeric_limits<double>::infinity(),
            py::arg("raw_radius_min_km") = py::none(),
            py::arg("raw_radius_max_km") = py::none(),
            py::call_guard<py::gil_scoped_release>());

    m.def("normalize_vectors_bulk", [](const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z) {
        size_t count = x.size();
        if (y.size() != count || z.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> ox(count), oy(count), oz(count);
        {
            py::gil_scoped_release release;
            lola::normalize_vectors_bulk(x.data(), y.data(), z.data(), ox.data(), oy.data(), oz.data(), count);
        }
        return py::make_tuple(ox, oy, oz);
    }, py::arg("x"), py::arg("y"), py::arg("z"));

    m.def("dot_product_bulk", [](const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const Vec3& reference) {
        size_t count = x.size();
        if (y.size() != count || z.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> results(count);
        {
            py::gil_scoped_release release;
            lola::dot_product_bulk(x.data(), y.data(), z.data(), reference, results.data(), count);
        }
        return results;
    }, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("reference"));

    m.def("cross_product_bulk", [](const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const Vec3& reference) {
        size_t count = x.size();
        if (y.size() != count || z.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> ox(count), oy(count), oz(count);
        {
            py::gil_scoped_release release;
            lola::cross_product_bulk(x.data(), y.data(), z.data(), reference, ox.data(), oy.data(), oz.data(), count);
        }
        return py::make_tuple(ox, oy, oz);
    }, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("reference"));

    m.def("project_onto_plane_bulk", [](const std::vector<double>& x_in, const std::vector<double>& y_in, const std::vector<double>& z_in, const Vec3& plane_normal) {
        size_t count = x_in.size();
        if (y_in.size() != count || z_in.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> ox(count), oy(count), oz(count);
        {
            py::gil_scoped_release release;
            lola::project_onto_plane_bulk(x_in.data(), y_in.data(), z_in.data(), plane_normal, ox.data(), oy.data(), oz.data(), count);
        }
        return py::make_tuple(ox, oy, oz);
    }, py::arg("x_in"), py::arg("y_in"), py::arg("z_in"), py::arg("plane_normal"));

    m.def("cartesian_to_spherical_bulk", [](const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z) {
        size_t count = x.size();
        if (y.size() != count || z.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> lon(count), lat(count), rad(count);
        {
            py::gil_scoped_release release;
            lola::cartesian_to_spherical_bulk(x.data(), y.data(), z.data(), lon.data(), lat.data(), rad.data(), count);
        }
        return py::make_tuple(lon, lat, rad);
    }, py::arg("x"), py::arg("y"), py::arg("z"));

    m.def("spherical_to_cartesian_bulk", [](const std::vector<double>& lon_deg, const std::vector<double>& lat_deg, const std::vector<double>& radius_km) {
        size_t count = lon_deg.size();
        if (lat_deg.size() != count || radius_km.size() != count) throw std::invalid_argument("LOLA: Coordinate vectors must have the same size");
        std::vector<double> x(count), y(count), z(count);
        {
            py::gil_scoped_release release;
            lola::spherical_to_cartesian_bulk(lon_deg.data(), lat_deg.data(), radius_km.data(), x.data(), y.data(), z.data(), count);
        }
        return py::make_tuple(x, y, z);
    }, py::arg("lon_deg"), py::arg("lat_deg"), py::arg("radius_km"));

    m.def("normalize_longitude_bulk", [](const std::vector<double>& lon_deg) {
        std::vector<double> out = lon_deg;
        {
            py::gil_scoped_release release;
            lola::normalize_longitude_bulk(out.data(), out.size());
        }
        return out;
    }, py::arg("lon_deg"));

    m.def("bin_by_position_angle", [](const std::vector<double>& pa_deg, double target_pa_deg, double bin_width_deg) {
        py::gil_scoped_release release;
        return lola::bin_by_position_angle(pa_deg.data(), target_pa_deg, bin_width_deg, pa_deg.size());
    }, py::arg("pa_deg"), py::arg("target_pa_deg"), py::arg("bin_width_deg"));

    m.def("select_max_radius_per_bin", [](const std::vector<int>& bin_indices, const std::vector<double>& radius_km) {
        if (bin_indices.size() != radius_km.size()) throw std::invalid_argument("LOLA: bin_indices and radius_km must have the same size");
        py::gil_scoped_release release;
        return lola::select_max_radius_per_bin(bin_indices.data(), radius_km.data(), bin_indices.size());
    }, py::arg("bin_indices"), py::arg("radius_km"));

    m.def("lexsort_by_bin_and_radius", [](const std::vector<int>& bin_indices, const std::vector<double>& radius_km) {
        if (bin_indices.size() != radius_km.size()) throw std::invalid_argument("LOLA: bin_indices and radius_km must have the same size");
        py::gil_scoped_release release;
        return lola::lexsort_by_bin_and_radius(bin_indices.data(), radius_km.data(), bin_indices.size());
    }, py::arg("bin_indices"), py::arg("radius_km"));

    m.def("convex_hull_2d", &lola::convex_hull_2d, py::arg("points"));
    
    m.def("ray_hull_intersection", &lola::ray_hull_intersection, 
        py::arg("hull"), py::arg("position_angle_deg"), py::arg("fallback_radius_km"));

    // --- Visibility ---
    py::class_<PhysicalVisibilityDirectExtinctionKernel>(
        m, "_PhysicalVisibilityDirectExtinctionKernel"
    )
        .def(
            py::init<std::vector<double>, std::size_t>(),
            py::arg("extinction_magnitude"),
            py::arg("spectral_bin_count")
        )
        .def(
            "interpolate",
            [](const PhysicalVisibilityDirectExtinctionKernel& kernel,
               std::size_t low,
               std::size_t high,
               double fraction) {
                PhysicalVisibilityDirectExtinction result;
                {
                    py::gil_scoped_release release;
                    result = kernel.interpolate(low, high, fraction);
                }
                return py::make_tuple(
                    vector_to_py_tuple(result.extinction_magnitude),
                    vector_to_py_tuple(result.transmission)
                );
            },
            py::arg("low"),
            py::arg("high"),
            py::arg("fraction")
        );

    m.def(
        "_physical_visibility_resolve_response_weights",
        &physical_visibility_response_weights_py,
        py::arg("band_wavelength_nm"),
        py::arg("band_differential_magnitude"),
        py::arg("spectral_bin_start_nm"),
        py::arg("base_scotopic_to_photopic_ratio"),
        py::arg("base_photopic"),
        py::arg("base_scotopic"),
        "Resolve one doctrine-free physical-visibility response-weight kernel."
    );

    py::class_<HeliacalEvent>(m, "HeliacalEvent")
        .def_readonly("event_kind", &HeliacalEvent::event_kind)
        .def_readonly("jd_ut", &HeliacalEvent::jd_ut)
        .def_readonly("is_found", &HeliacalEvent::is_found)
        .def_readonly("arcus_visionis", &HeliacalEvent::arcus_visionis)
        .def_readonly("elongation", &HeliacalEvent::elongation)
        .def_readonly("star_altitude", &HeliacalEvent::star_altitude)
        .def_readonly("day_offset", &HeliacalEvent::day_offset);

    m.def("arcus_visionis", &arcus_visionis, 
        py::arg("mag"), py::arg("limiting_mag"), py::arg("extinction_k"),
        "Compute required Arcus Visionis for visibility via Schoch/Ptolemy logic.");

    m.def("target_topocentric_altitude",
        [](const IEvaluator& target_eval, double jd_ut, double lat_deg, double lon_deg,
           double pressure_mbar, double temperature_c, bool use_refraction, double delta_t,
           const IEvaluator* earth_eval) {
            py::gil_scoped_release release;
            return target_topocentric_altitude(target_eval, jd_ut, lat_deg, lon_deg,
                pressure_mbar, temperature_c, use_refraction, delta_t, earth_eval);
        },
        py::arg("target_eval"), py::arg("jd_ut"), py::arg("lat_deg"), py::arg("lon_deg"),
        py::arg("pressure_mbar") = 1013.25, py::arg("temperature_c") = 10.0, py::arg("use_refraction") = true,
        py::arg("delta_t") = 64.184, py::arg("earth_eval") = nullptr,
        "Compute topocentric altitude with optional annual aberration via Earth evaluator.");

    m.def("find_sun_at_alt",
        [](const IEvaluator& sun_eval, double jd_midnight, double lat_deg, double lon_deg,
           double target_alt, bool morning, double delta_t, const IEvaluator* earth_eval) {
            py::gil_scoped_release release;
            return find_sun_at_alt(sun_eval, jd_midnight, lat_deg, lon_deg,
                target_alt, morning, delta_t, earth_eval);
        },
        py::arg("sun_eval"), py::arg("jd_midnight"), py::arg("lat_deg"), py::arg("lon_deg"),
        py::arg("target_alt"), py::arg("morning"), py::arg("delta_t") = 64.184,
        py::arg("earth_eval") = nullptr,
        "Solve for JD within a half-day window where the Sun reaches a target altitude.");

    m.def("search_heliacal_rising",
        [](const IEvaluator& star_eval, const IEvaluator& sun_eval, double jd_start,
           double lat, double lon, double arcus_visionis_val, int search_days,
           double delta_t, const IEvaluator* earth_eval,
           double delta_t_rate_seconds_per_day, NutationEpochCache* nutation_cache) {
            py::gil_scoped_release release;
            return search_heliacal_rising(star_eval, sun_eval, jd_start, lat, lon,
                arcus_visionis_val, search_days, delta_t, earth_eval,
                delta_t_rate_seconds_per_day, nutation_cache);
        },
        py::arg("star_eval"), py::arg("sun_eval"), py::arg("jd_start"),
        py::arg("lat"), py::arg("lon"), py::arg("arcus_visionis_val"),
        py::arg("search_days"), py::arg("delta_t") = 64.184, py::arg("earth_eval") = nullptr,
        py::arg("delta_t_rate_seconds_per_day") = 0.0,
        py::arg("nutation_cache") = nullptr,
        "Native fast-path search for star heliacal rising (morning appearance).");

    m.def("search_heliacal_setting",
        [](const IEvaluator& star_eval, const IEvaluator& sun_eval, double jd_start,
           double lat, double lon, double arcus_visionis_val, int search_days,
           double delta_t, const IEvaluator* earth_eval,
           double setting_elongation_threshold, double setting_visibility_factor,
           double delta_t_rate_seconds_per_day, NutationEpochCache* nutation_cache) {
            py::gil_scoped_release release;
            return search_heliacal_setting(star_eval, sun_eval, jd_start, lat, lon,
                arcus_visionis_val, search_days, delta_t, earth_eval,
                setting_elongation_threshold, setting_visibility_factor,
                delta_t_rate_seconds_per_day, nutation_cache);
        },
        py::arg("star_eval"), py::arg("sun_eval"), py::arg("jd_start"),
        py::arg("lat"), py::arg("lon"), py::arg("arcus_visionis_val"),
        py::arg("search_days"), py::arg("delta_t") = 64.184, py::arg("earth_eval") = nullptr,
        py::arg("setting_elongation_threshold") = 12.0,
        py::arg("setting_visibility_factor") = 0.5,
        py::arg("delta_t_rate_seconds_per_day") = 0.0,
        py::arg("nutation_cache") = nullptr,
        "Native fast-path search for star heliacal setting (morning disappearance).");

    m.def("heliacal_signed_elongation", &heliacal_signed_elongation,
        py::arg("star_eval"), py::arg("sun_eval"), py::arg("jd_ut"),
        py::arg("delta_t") = 64.184, py::arg("earth_eval") = nullptr,
        py::arg("nutation_cache") = nullptr,
        "Compute signed ecliptic elongation between star and Sun.");

    m.def("mean_obliquity_p03", &mean_obliquity_p03, py::arg("jd_tt"),
        "IAU 2006 Mean Obliquity (P03 model).");

    // ── Harmogram acceleration kernels ────────────────────────────────────────

    m.def("harmogram_compute_components",
        [](const std::vector<double>& longitudes_deg,
           const std::vector<int>&    harmonics,
           bool                       raw_sum)
        {
            py::gil_scoped_release release;
            return harmogram_compute_components(longitudes_deg, harmonics, raw_sum);
        },
        py::arg("longitudes_deg"), py::arg("harmonics"), py::arg("raw_sum"),
        "Compute harmonic vector components (amplitude, phase) for a set of "
        "ecliptic longitudes. Returns list of (harmonic, amplitude, phase_deg).");

    m.def("harmogram_trace_batch",
        [](const std::vector<std::vector<double>>& samples_source_lons,
           const std::vector<std::vector<double>>& samples_target_lons,
           bool                                    same_source_target,
           bool                                    ordered,
           bool                                    include_self,
           bool                                    raw_sum,
           double                                  h0_contribution,
           const std::vector<int>&                 harmonics,
           const std::vector<std::tuple<int, double, double>>& intensity_components)
        {
            HarmogramBatchResult res;
            {
                py::gil_scoped_release release;
                res = harmogram_trace_batch(
                    samples_source_lons, samples_target_lons,
                    same_source_target, ordered, include_self, raw_sum,
                    h0_contribution, harmonics, intensity_components
                );
            }
            return py::make_tuple(res.total_strengths, res.sample_components);
        },
        py::arg("samples_source_lons"),
        py::arg("samples_target_lons"),
        py::arg("same_source_target"),
        py::arg("ordered"),
        py::arg("include_self"),
        py::arg("raw_sum"),
        py::arg("h0_contribution"),
        py::arg("harmonics"),
        py::arg("intensity_components"),
        "Batch ZA-parts Fourier + intensity projection for N time samples, "
        "parallelised via OpenMP. Returns (total_strengths[N], "
        "sample_components[N][H]) where each component is (harmonic, amplitude, phase_deg).");

    m.def("harmogram_intensity_components",
        [](int                       harmonic_number,
           int                       harmonic_start,
           int                       harmonic_stop,
           int                       sample_count,
           const std::string&        orb_mode,
           double                    orb_width_deg,
           bool                      include_conjunction,
           const std::string&        orb_scaling_mode,
           double                    gaussian_width_deg,
           bool                      gaussian_fwhm_mode)
        {
            std::tuple<double, std::vector<std::tuple<int, double, double>>> result;
            {
                py::gil_scoped_release release;
                result = harmogram_intensity_components(
                    harmonic_number, harmonic_start, harmonic_stop,
                    sample_count, orb_mode, orb_width_deg, include_conjunction,
                    orb_scaling_mode, gaussian_width_deg, gaussian_fwhm_mode
                );
            }
            auto [h0_amp, comps] = result;
            return py::make_tuple(h0_amp, comps);
        },
        py::arg("harmonic_number"),
        py::arg("harmonic_start"),
        py::arg("harmonic_stop"),
        py::arg("sample_count"),
        py::arg("orb_mode"),
        py::arg("orb_width_deg"),
        py::arg("include_conjunction"),
        py::arg("orb_scaling_mode"),
        py::arg("gaussian_width_deg") = 0.0,
        py::arg("gaussian_fwhm_mode") = false,
        "Sample the intensity orb function and compute its DFT spectrum. "
        "Returns (h0_amplitude, [(harmonic, amplitude, phase_deg), ...]).");
}