#![allow(clippy::excessive_precision)]

use sidereon_core::astro::bodies::{
    moon_az_el, observe, observe_spk_body, sun_az_el, ObserveError, ObserveOptions, Refraction,
    Target,
};
use sidereon_core::astro::constants::earth::OMEGA_E_DOT_RAD_S;
use sidereon_core::astro::constants::physics::SPEED_OF_LIGHT_M_S;
use sidereon_core::astro::constants::time::{J2000_JD, SECONDS_PER_DAY};
use sidereon_core::astro::constants::units::M_PER_KM;
use sidereon_core::astro::frames::transforms::{
    gcrs_to_itrs_matrix, gcrs_to_true_of_date_matrix, geodetic_to_itrs,
    greenwich_apparent_sidereal_time_radians, itrs_to_gcrs_matrix, mat3_vec3_mul,
    GeodeticStationKm,
};
use sidereon_core::astro::math::mat3::{inline_rxr, Mat3};
use sidereon_core::astro::math::vec3::{add3, dot3, norm3, scale3, sub3};
use sidereon_core::astro::passes::UtcInstant;
use sidereon_core::astro::spk::Spk;
use sidereon_core::astro::time::TimeScales;

const ARCSEC_PER_RADIAN: f64 = 206_264.806_247_096_36;
const C_KM_S: f64 = SPEED_OF_LIGHT_M_S / M_PER_KM;

#[derive(Clone, Copy)]
struct GoldenEquatorial {
    ra_deg: f64,
    dec_deg: f64,
    distance_km: f64,
}

#[derive(Clone, Copy)]
struct GoldenHorizontal {
    az_deg: f64,
    el_deg: f64,
    range_km: f64,
}

#[derive(Clone, Copy)]
struct GoldenEcliptic {
    lon_deg: f64,
    lat_deg: f64,
    distance_km: f64,
}

#[derive(Clone, Copy)]
struct ObserveGolden {
    name: &'static str,
    target_naif: i32,
    station: fn() -> GeodeticStationKm,
    utc: (i32, i32, i32, i32, i32, i32),
    astrometric: GoldenEquatorial,
    apparent_icrs: GoldenEquatorial,
    apparent: GoldenEquatorial,
    horizontal: GoldenHorizontal,
    ecliptic: GoldenEcliptic,
}

type FrameGoldenCase = ((i32, i32, i32, i32, i32, f64), Mat3);

// Generated by tests/fixtures/bodies/gen_observe_golden.py
// Source: observe_de.bsp with Skyfield, Sun-only deflection.
// Horizons spot check for Mars barycenter from Greenwich at
// 2024-01-01T00:00:00 UTC: apparent of-date RA/Dec agrees within
// the sub-arcsecond tolerance used below.
const OBSERVE_GOLDEN: &[ObserveGolden] = &[
    ObserveGolden {
        name: "mars_greenwich",
        target_naif: 4,
        station: station_greenwich,
        utc: (2024, 1, 1, 0, 0, 0),
        astrometric: GoldenEquatorial {
            ra_deg: 2.66695691769370058e+02,
            dec_deg: -2.39525070436205283e+01,
            distance_km: 3.62601889498879015e+08,
        },
        apparent_icrs: GoldenEquatorial {
            ra_deg: 2.66689454311101997e+02,
            dec_deg: -2.39523599806567375e+01,
            distance_km: 3.62601889498881221e+08,
        },
        apparent: GoldenEquatorial {
            ra_deg: 2.67054529373346782e+02,
            dec_deg: -2.39618517989077198e+01,
            distance_km: 3.62601889498880446e+08,
        },
        horizontal: GoldenHorizontal {
            az_deg: 2.50298037597854872e+01,
            el_deg: -6.07013533515289936e+01,
            range_km: 3.62601889498880506e+08,
        },
        ecliptic: GoldenEcliptic {
            lon_deg: 2.67308452801835415e+02,
            lat_deg: -5.50968260129116927e-01,
            distance_km: 3.62601889498880446e+08,
        },
    },
    ObserveGolden {
        name: "jupiter_greenwich",
        target_naif: 5,
        station: station_greenwich,
        utc: (2024, 1, 1, 0, 0, 0),
        astrometric: GoldenEquatorial {
            ra_deg: 3.33611390561985246e+01,
            dec_deg: 1.21507431073479779e+01,
            distance_km: 6.70420690887020946e+08,
        },
        apparent_icrs: GoldenEquatorial {
            ra_deg: 3.33635997421848529e+01,
            dec_deg: 1.21514791331984551e+01,
            distance_km: 6.70420690887020946e+08,
        },
        apparent: GoldenEquatorial {
            ra_deg: 3.36852066073541820e+01,
            dec_deg: 1.22635965691107600e+01,
            distance_km: 6.70420690887021422e+08,
        },
        horizontal: GoldenHorizontal {
            az_deg: 2.59071073620473896e+02,
            el_deg: 2.41550799692646976e+01,
            range_km: 6.70420690887021422e+08,
        },
        ecliptic: GoldenEcliptic {
            lon_deg: 3.55819952432663129e+01,
            lat_deg: -1.18573982460513649e+00,
            distance_km: 6.70420690887021422e+08,
        },
    },
    ObserveGolden {
        name: "mars_sydney",
        target_naif: 4,
        station: station_sydney,
        utc: (2024, 1, 1, 0, 0, 0),
        astrometric: GoldenEquatorial {
            ra_deg: 2.66695785862900038e+02,
            dec_deg: -2.39518562573648488e+01,
            distance_km: 3.62590229819242179e+08,
        },
        apparent_icrs: GoldenEquatorial {
            ra_deg: 2.66689685398737652e+02,
            dec_deg: -2.39517065278724672e+01,
            distance_km: 3.62590229819244266e+08,
        },
        apparent: GoldenEquatorial {
            ra_deg: 2.67054758665140866e+02,
            dec_deg: -2.39611978130430501e+01,
            distance_km: 3.62590229819243550e+08,
        },
        horizontal: GoldenHorizontal {
            az_deg: 5.82535677038007975e+01,
            el_deg: 7.31077532315459422e+01,
            range_km: 3.62590229819243550e+08,
        },
        ecliptic: GoldenEcliptic {
            lon_deg: 2.67308648931098048e+02,
            lat_deg: -5.50310128197137427e-01,
            distance_km: 3.62590229819243491e+08,
        },
    },
    ObserveGolden {
        name: "jupiter_sydney",
        target_naif: 5,
        station: station_sydney,
        utc: (2024, 1, 1, 0, 0, 0),
        astrometric: GoldenEquatorial {
            ra_deg: 3.33617406289273077e+01,
            dec_deg: 1.21513486940236071e+01,
            distance_km: 6.70428149882314205e+08,
        },
        apparent_icrs: GoldenEquatorial {
            ra_deg: 3.33641188866043379e+01,
            dec_deg: 1.21520642281249742e+01,
            distance_km: 6.70428149882314205e+08,
        },
        apparent: GoldenEquatorial {
            ra_deg: 3.36857267366724642e+01,
            dec_deg: 1.22641810147220589e+01,
            distance_km: 6.70428149882314920e+08,
        },
        horizontal: GoldenHorizontal {
            az_deg: 1.13074149464185325e+02,
            el_deg: -4.95146078485224379e+01,
            range_km: 6.70428149882314920e+08,
        },
        ecliptic: GoldenEcliptic {
            lon_deg: 3.55826684676442540e+01,
            lat_deg: -1.18535659146729366e+00,
            distance_km: 6.70428149882314920e+08,
        },
    },
    ObserveGolden {
        name: "mars_near_conjunction",
        target_naif: 4,
        station: station_greenwich,
        utc: (2023, 11, 18, 0, 0, 0),
        astrometric: GoldenEquatorial {
            ra_deg: 2.32743352531053091e+02,
            dec_deg: -1.91525143806728266e+01,
            distance_km: 3.77952647406406105e+08,
        },
        apparent_icrs: GoldenEquatorial {
            ra_deg: 2.32737582732398380e+02,
            dec_deg: -1.91516960830259819e+01,
            distance_km: 3.77952647427149713e+08,
        },
        apparent: GoldenEquatorial {
            ra_deg: 2.33077725278055482e+02,
            dec_deg: -1.92331369740263050e+01,
            distance_km: 3.77952647427149773e+08,
        },
        horizontal: GoldenHorizontal {
            az_deg: 6.53890721526704244e+00,
            el_deg: -5.76236365760795550e+01,
            range_km: 3.77952647427149713e+08,
        },
        ecliptic: GoldenEcliptic {
            lon_deg: 2.35444588446476416e+02,
            lat_deg: -1.13926919166224633e-01,
            distance_km: 3.77952647427149832e+08,
        },
    },
];

// Generated from pre-refactor commit 619f449 with:
// cargo build -p sidereon-core; rustc /tmp/dump_sid_matrix.rs ...
const FRAME_GOLDEN: &[FrameGoldenCase] = &[
    (
        (2000, 1, 1, 12, 0, 0.0),
        [
            [
                1.81585122847680802e-1,
                -9.83375229832204600e-1,
                -2.26462012731456990e-5,
            ],
            [
                9.83375229723779998e-1,
                1.81585122100325458e-1,
                3.15834393603663929e-5,
            ],
            [
                -2.69461587142590570e-5,
                -2.80047961012289743e-5,
                9.99999999244818083e-1,
            ],
        ],
    ),
    (
        (2020, 6, 24, 12, 34, 56.0),
        [
            [
                -2.01022163818197153e-1,
                9.79586611856021472e-1,
                3.99407422332429831e-4,
            ],
            [
                -9.79584737922584714e-1,
                -2.01022560513757581e-1,
                1.91608810716750451e-3,
            ],
            [
                1.95726415964563851e-3,
                -6.07723776638428277e-6,
                9.99998084538203713e-1,
            ],
        ],
    ),
    (
        (2024, 1, 1, 0, 0, 0.0),
        [
            [
                -1.70985862949848105e-1,
                9.85273414718201956e-1,
                3.64583089080899262e-4,
            ],
            [
                -9.85270747182341422e-1,
                -1.70986248483862069e-1,
                2.29294050648503125e-3,
            ],
            [
                2.32151201723513039e-3,
                3.28473585994561317e-5,
                9.99997304747870186e-1,
            ],
        ],
    ),
    (
        (2026, 7, 1, 0, 0, 0.0),
        [
            [
                1.51667472875174736e-1,
                -9.88431507364607498e-1,
                -3.64582757404158340e-4,
            ],
            [
                9.88428180113040566e-1,
                1.51667908368473991e-1,
                -2.56482544345893642e-3,
            ],
            [
                2.59044978345025550e-3,
                2.86367219681365483e-5,
                9.99996644369298804e-1,
            ],
        ],
    ),
];

#[test]
fn spk_observe_matches_skyfield_golden() {
    let kernel = observe_kernel();
    for case in OBSERVE_GOLDEN {
        let station = (case.station)();
        let time = utc(case.utc);
        let got = observe_spk_body(&station, time, &kernel, case.target_naif)
            .unwrap_or_else(|err| panic!("{}: {err}", case.name));

        assert_equatorial(
            case.name,
            "astrometric",
            got.astrometric,
            case.astrometric,
            1.0,
        );
        assert_equatorial(
            case.name,
            "apparent_icrs",
            got.apparent_icrs,
            case.apparent_icrs,
            1.0,
        );
        assert_equatorial(case.name, "apparent", got.apparent, case.apparent, 1.0);
        assert_angle_arcsec(
            case.name,
            "az",
            got.horizontal.azimuth_deg,
            case.horizontal.az_deg,
            2.0,
        );
        assert_angle_arcsec(
            case.name,
            "el",
            got.horizontal.elevation_deg,
            case.horizontal.el_deg,
            2.0,
        );
        assert_relative(
            case.name,
            "horizontal range",
            got.horizontal.range_km,
            case.horizontal.range_km,
            1.0e-8,
        );
        assert_angle_arcsec(
            case.name,
            "ecliptic lon",
            got.ecliptic.longitude_deg,
            case.ecliptic.lon_deg,
            1.0,
        );
        assert_linear_arcsec(
            case.name,
            "ecliptic lat",
            got.ecliptic.latitude_deg,
            case.ecliptic.lat_deg,
            1.0,
        );
        assert_relative(
            case.name,
            "ecliptic distance",
            got.ecliptic.distance_km,
            case.ecliptic.distance_km,
            1.0e-8,
        );
        assert!(!got.reduced, "{} reduced flag", case.name);
    }
}

#[test]
fn hour_angle_reconstructs_altaz_without_polar_motion() {
    let kernel = observe_kernel();
    let case = OBSERVE_GOLDEN
        .iter()
        .find(|case| case.name == "mars_sydney")
        .expect("golden case");
    let station = (case.station)();
    let got = observe_spk_body(&station, utc(case.utc), &kernel, case.target_naif)
        .expect("valid observation");

    let lat = station.latitude_deg.to_radians();
    let dec = got.apparent.declination_deg.to_radians();
    let h = got.hour_angle_deg.to_radians();
    let sin_alt = lat.sin() * dec.sin() + lat.cos() * dec.cos() * h.cos();
    let alt = sin_alt.clamp(-1.0, 1.0).asin();
    let east = -dec.cos() * h.sin();
    let north = dec.sin() * lat.cos() - dec.cos() * lat.sin() * h.cos();
    let az = wrap_360(east.atan2(north).to_degrees());

    assert!(
        (alt.to_degrees() - got.horizontal.elevation_deg).abs() < 1.0e-6,
        "alt {} vs {}",
        alt.to_degrees(),
        got.horizontal.elevation_deg
    );
    assert!(
        angle_diff_deg(az, got.horizontal.azimuth_deg).abs() < 1.0e-6,
        "az {az} vs {}",
        got.horizontal.azimuth_deg
    );
}

#[test]
fn gcrs_to_true_of_date_refactor_preserves_gcrs_to_itrs() {
    for (utc_fields, expected) in FRAME_GOLDEN {
        let ts = TimeScales::from_utc(
            utc_fields.0,
            utc_fields.1,
            utc_fields.2,
            utc_fields.3,
            utc_fields.4,
            utc_fields.5,
        )
        .expect("valid UTC");
        let got = gcrs_to_itrs_matrix(&ts).expect("valid matrix");
        assert_matrix_near("direct", &got, expected, 1.0e-12);

        let tod = gcrs_to_true_of_date_matrix(&ts).expect("valid true-of-date matrix");
        let gast = greenwich_apparent_sidereal_time_radians(&ts).expect("valid GAST");
        let rebuilt = inline_rxr(&rot_z(-gast), &tod);
        assert_matrix_near("rebuilt", &rebuilt, expected, 1.0e-12);
    }

    let ts = TimeScales::from_utc(2024, 1, 1, 0, 0, 0.0).expect("valid UTC");
    let unit = unit_from_radec(267.0545293733468, -23.96185179890772);
    let itrs = mat3_vec3_mul(&gcrs_to_itrs_matrix(&ts).expect("valid matrix"), &unit)
        .expect("GCRS to ITRS");
    let roundtrip = mat3_vec3_mul(&itrs_to_gcrs_matrix(&ts).expect("valid matrix"), &itrs)
        .expect("ITRS to GCRS");
    for i in 0..3 {
        assert!((roundtrip[i] - unit[i]).abs() < 1.0e-12);
    }
}

#[test]
fn reduced_sun_observation_reproduces_existing_sun_az_el() {
    let station = station_greenwich();
    let time = UtcInstant::from_utc(2024, 6, 20, 12, 1, 42, 0).expect("valid UTC");
    let legacy = sun_az_el(&station, time).expect("sun look angle");
    let full = observe(&station, time, Target::Sun, ObserveOptions::default())
        .expect("reduced sun observation");

    assert_eq!(
        full.horizontal.azimuth_deg.to_bits(),
        legacy.azimuth_deg.to_bits()
    );
    assert_eq!(
        full.horizontal.elevation_deg.to_bits(),
        legacy.elevation_deg.to_bits()
    );
    assert_eq!(
        full.horizontal.range_km.to_bits(),
        legacy.range_km.to_bits()
    );
    assert!(full.reduced);

    let moon = observe(&station, time, Target::Moon, ObserveOptions::default())
        .expect("reduced moon observation");
    let legacy_moon = moon_az_el(&station, time).expect("moon look angle");
    assert_eq!(
        moon.horizontal.azimuth_deg.to_bits(),
        legacy_moon.azimuth_deg.to_bits()
    );
    assert_eq!(
        moon.horizontal.elevation_deg.to_bits(),
        legacy_moon.elevation_deg.to_bits()
    );
    assert_eq!(
        moon.horizontal.range_km.to_bits(),
        legacy_moon.range_km.to_bits()
    );
    assert!(moon.reduced);
}

#[test]
fn correction_toggles_have_expected_scales() {
    let kernel = observe_kernel();
    let station = station_greenwich();
    let time = UtcInstant::from_utc(2024, 1, 1, 0, 0, 0, 0).expect("valid UTC");

    let full = observe_spk_body(&station, time, &kernel, 4).expect("full observation");
    let no_aberration = observe(
        &station,
        time,
        Target::Spk {
            kernel: &kernel,
            naif_id: 4,
        },
        ObserveOptions {
            aberration: false,
            ..ObserveOptions::default()
        },
    )
    .expect("observation without aberration");
    let shift = equatorial_sep_arcsec_between(full.apparent_icrs, no_aberration.apparent_icrs);
    let estimate = aberration_estimate_arcsec(&station, time, &kernel, no_aberration.apparent_icrs);
    assert!(shift > 1.0 && shift < 30.0, "aberration shift {shift}");
    assert!(
        (shift - estimate).abs() < 0.5,
        "aberration shift {shift} estimate {estimate}"
    );

    let conjunction_time = UtcInstant::from_utc(2023, 11, 18, 0, 0, 0, 0).expect("valid UTC");
    let with_deflection =
        observe_spk_body(&station, conjunction_time, &kernel, 4).expect("full conjunction");
    let without_deflection = observe(
        &station,
        conjunction_time,
        Target::Spk {
            kernel: &kernel,
            naif_id: 4,
        },
        ObserveOptions {
            deflection: false,
            ..ObserveOptions::default()
        },
    )
    .expect("conjunction without deflection");
    let deflection_shift = equatorial_sep_arcsec_between(
        with_deflection.apparent_icrs,
        without_deflection.apparent_icrs,
    );
    let elongation = equatorial_sep_arcsec_between(
        without_deflection.apparent_icrs,
        observe_spk_body(&station, conjunction_time, &kernel, 10)
            .expect("sun SPK observation")
            .apparent_icrs,
    ) / ARCSEC_PER_RADIAN;
    let bound_arcsec = 4.0 * 132_712_440_041.939_38
        / (C_KM_S * C_KM_S * 149_597_870.700 * elongation)
        * ARCSEC_PER_RADIAN;
    assert!(
        deflection_shift > 0.01 && deflection_shift < bound_arcsec * 2.0,
        "deflection shift {deflection_shift}, bound {bound_arcsec}"
    );
}

#[test]
fn refraction_matches_skyfield_bennett_reference() {
    let kernel = observe_kernel();
    let station = station_sydney();
    let time = UtcInstant::from_utc(2024, 1, 1, 0, 0, 0, 0).expect("valid UTC");
    let got = observe(
        &station,
        time,
        Target::Spk {
            kernel: &kernel,
            naif_id: 4,
        },
        ObserveOptions {
            refraction: Some(Refraction {
                pressure_mbar: 1013.25,
                temperature_c: 10.0,
            }),
            ..ObserveOptions::default()
        },
    )
    .expect("refracted observation");
    assert!(
        (got.horizontal.elevation_deg - 73.11279565959617).abs() * 60.0 < 0.1,
        "refracted elevation {}",
        got.horizontal.elevation_deg
    );
}

#[test]
fn non_finite_inputs_return_typed_errors() {
    let kernel = observe_kernel();
    let station = station_greenwich();
    let time = UtcInstant::from_utc(2024, 1, 1, 0, 0, 0, 0).expect("valid UTC");
    let err = observe(
        &station,
        time,
        Target::BarycentricState {
            kernel: &kernel,
            position_km: [f64::NAN, 0.0, 0.0],
            velocity_km_s: [0.0, 0.0, 0.0],
        },
        ObserveOptions::default(),
    )
    .unwrap_err();
    assert!(matches!(err, ObserveError::NonFinite));
}

fn observe_kernel() -> Spk {
    Spk::from_bytes(include_bytes!("fixtures/bodies/observe_de.bsp")).expect("fixture SPK")
}

fn utc(fields: (i32, i32, i32, i32, i32, i32)) -> UtcInstant {
    UtcInstant::from_utc(
        fields.0, fields.1, fields.2, fields.3, fields.4, fields.5, 0,
    )
    .expect("valid UTC")
}

fn station_greenwich() -> GeodeticStationKm {
    GeodeticStationKm {
        latitude_deg: 51.4779,
        longitude_deg: -0.0015,
        altitude_km: 0.046,
    }
}

fn station_sydney() -> GeodeticStationKm {
    GeodeticStationKm {
        latitude_deg: -33.8568,
        longitude_deg: 151.2153,
        altitude_km: 0.040,
    }
}

fn assert_equatorial(
    case: &str,
    field: &str,
    actual: sidereon_core::astro::bodies::Equatorial,
    expected: GoldenEquatorial,
    max_arcsec: f64,
) {
    let sep = equatorial_sep_arcsec(actual, expected);
    assert!(sep < max_arcsec, "{case} {field} separation {sep} arcsec");
    assert_relative(
        case,
        field,
        actual.distance_km,
        expected.distance_km,
        1.0e-8,
    );
}

fn equatorial_sep_arcsec(
    actual: sidereon_core::astro::bodies::Equatorial,
    expected: GoldenEquatorial,
) -> f64 {
    let a = unit_from_radec(actual.right_ascension_deg, actual.declination_deg);
    let b = unit_from_radec(expected.ra_deg, expected.dec_deg);
    dot3(a, b).clamp(-1.0, 1.0).acos() * ARCSEC_PER_RADIAN
}

fn equatorial_sep_arcsec_between(
    a: sidereon_core::astro::bodies::Equatorial,
    b: sidereon_core::astro::bodies::Equatorial,
) -> f64 {
    let a = unit_from_radec(a.right_ascension_deg, a.declination_deg);
    let b = unit_from_radec(b.right_ascension_deg, b.declination_deg);
    dot3(a, b).clamp(-1.0, 1.0).acos() * ARCSEC_PER_RADIAN
}

fn assert_angle_arcsec(case: &str, field: &str, actual: f64, expected: f64, max_arcsec: f64) {
    let diff = angle_diff_deg(actual, expected).abs() * 3600.0;
    assert!(diff < max_arcsec, "{case} {field} diff {diff} arcsec");
}

fn assert_linear_arcsec(case: &str, field: &str, actual: f64, expected: f64, max_arcsec: f64) {
    let diff = (actual - expected).abs() * 3600.0;
    assert!(diff < max_arcsec, "{case} {field} diff {diff} arcsec");
}

fn assert_relative(case: &str, field: &str, actual: f64, expected: f64, max_rel: f64) {
    let rel = (actual - expected).abs() / expected.abs().max(1.0);
    assert!(rel < max_rel, "{case} {field} relative diff {rel}");
}

fn unit_from_radec(ra_deg: f64, dec_deg: f64) -> [f64; 3] {
    let ra = ra_deg.to_radians();
    let dec = dec_deg.to_radians();
    [dec.cos() * ra.cos(), dec.cos() * ra.sin(), dec.sin()]
}

fn angle_diff_deg(actual: f64, expected: f64) -> f64 {
    let diff = (actual - expected + 180.0).rem_euclid(360.0) - 180.0;
    if diff <= -180.0 {
        diff + 360.0
    } else {
        diff
    }
}

fn wrap_360(deg: f64) -> f64 {
    deg.rem_euclid(360.0)
}

fn rot_z(angle: f64) -> Mat3 {
    let c = angle.cos();
    let s = angle.sin();
    [[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]
}

fn assert_matrix_near(label: &str, actual: &Mat3, expected: &Mat3, tolerance: f64) {
    for i in 0..3 {
        for j in 0..3 {
            assert!(
                (actual[i][j] - expected[i][j]).abs() < tolerance,
                "{label}[{i}][{j}] {} vs {}",
                actual[i][j],
                expected[i][j]
            );
        }
    }
}

fn aberration_estimate_arcsec(
    station: &GeodeticStationKm,
    time: UtcInstant,
    kernel: &Spk,
    direction: sidereon_core::astro::bodies::Equatorial,
) -> f64 {
    let ts = time.time_scales();
    let et = (ts.jd_tdb - J2000_JD) * SECONDS_PER_DAY;
    let (x, y, z) = geodetic_to_itrs(
        station.latitude_deg,
        station.longitude_deg,
        station.altitude_km,
    )
    .expect("station ITRS");
    let station_itrs = [x, y, z];
    let v_itrs = [
        -OMEGA_E_DOT_RAD_S * station_itrs[1],
        OMEGA_E_DOT_RAD_S * station_itrs[0],
        0.0,
    ];
    let v_geo = mat3_vec3_mul(
        &itrs_to_gcrs_matrix(&ts).expect("ITRS to GCRS matrix"),
        &v_itrs,
    )
    .expect("station velocity");
    let v_earth = spk_velocity(kernel, 399, 0, et);
    let beta = scale3(add3(v_earth, v_geo), 1.0 / C_KM_S);
    let p = unit_from_radec(direction.right_ascension_deg, direction.declination_deg);
    let beta_parallel = scale3(p, dot3(p, beta));
    norm3(sub3(beta, beta_parallel)) * ARCSEC_PER_RADIAN
}

fn spk_velocity(kernel: &Spk, target: i32, center: i32, et: f64) -> [f64; 3] {
    let state = kernel
        .spk_state(target, center, et)
        .expect("SPK state for velocity");
    if let Some(velocity) = state.velocity_km_s {
        return velocity;
    }
    let before = kernel
        .spk_state(target, center, et - 1.0)
        .expect("SPK state before")
        .position_km;
    let after = kernel
        .spk_state(target, center, et + 1.0)
        .expect("SPK state after")
        .position_km;
    scale3(sub3(after, before), 0.5)
}
