use bone_kernel::{
    Curve3, Degree, KnotVector, NurbsCurve3, NurbsSurface3, RationalControlPoint, Surface3, TriMesh,
};
use bone_types::{AngleTolerance, ChordHeightTolerance, Parameter, Point3, Tolerance, Vec3};
use core::f64::consts::FRAC_1_SQRT_2;
use uom::si::length::millimeter;

const TOL: Tolerance = Tolerance::new(1.0e-9);

fn p(value: f64) -> Parameter {
    Parameter::new(value)
}

fn control(x: f64, y: f64, z: f64, weight: f64) -> RationalControlPoint {
    let Ok(point) = RationalControlPoint::new(Point3::from_mm(x, y, z), weight, TOL) else {
        panic!("control point weight is positive");
    };
    point
}

fn degree(value: usize) -> Degree {
    let Ok(degree) = Degree::new(value) else {
        panic!("degree is positive");
    };
    degree
}

fn knots(values: &[f64]) -> KnotVector {
    let Ok(vector) = KnotVector::new(values.to_vec()) else {
        panic!("knot vector is finite, monotone, and spread");
    };
    vector
}

fn line() -> NurbsCurve3 {
    let Ok(curve) = NurbsCurve3::new(
        degree(1),
        knots(&[0.0, 0.0, 1.0, 1.0]),
        vec![control(1.0, 2.0, 3.0, 1.0), control(4.0, 6.0, 8.0, 1.0)],
    ) else {
        panic!("a degree-one curve with two control points is valid");
    };
    curve
}

fn quarter_circle() -> NurbsCurve3 {
    let Ok(curve) = NurbsCurve3::new(
        degree(2),
        knots(&[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]),
        vec![
            control(1.0, 0.0, 0.0, 1.0),
            control(1.0, 1.0, 0.0, FRAC_1_SQRT_2),
            control(0.0, 1.0, 0.0, 1.0),
        ],
    ) else {
        panic!("the rational quadratic quarter circle is valid");
    };
    curve
}

#[test]
fn degree_one_curve_is_a_straight_line() {
    let curve = line();
    let start = curve.evaluate(p(0.0));
    let mid = curve.evaluate(p(0.5));
    let end = curve.evaluate(p(1.0));
    assert!((start - Point3::from_mm(1.0, 2.0, 3.0)).norm_mm() < 1.0e-9);
    assert!((end - Point3::from_mm(4.0, 6.0, 8.0)).norm_mm() < 1.0e-9);
    assert!((mid - Point3::from_mm(2.5, 4.0, 5.5)).norm_mm() < 1.0e-9);
    let tangent = curve.derivative(p(0.5));
    assert!((tangent - Vec3::from_mm(3.0, 4.0, 5.0)).norm_mm() < 1.0e-6);
}

#[test]
fn rational_quadratic_traces_the_unit_circle() {
    let curve = quarter_circle();
    (0..=16).for_each(|i| {
        let t = f64::from(i) / 16.0;
        let point = curve.evaluate(p(t));
        let (x, y, z) = point.coords_mm();
        assert!(
            (x.hypot(y) - 1.0).abs() < 1.0e-9,
            "every point sits on the unit circle: ({x}, {y}) radius {}",
            x.hypot(y)
        );
        assert!(z.abs() < 1.0e-12);
    });
    let mid = curve.evaluate(p(0.5));
    assert!((mid - Point3::from_mm(FRAC_1_SQRT_2, FRAC_1_SQRT_2, 0.0)).norm_mm() < 1.0e-9);
}

#[test]
fn circle_tangent_is_perpendicular_to_the_radius() {
    let curve = quarter_circle();
    [0.0, 0.25, 0.5, 0.75, 1.0].into_iter().for_each(|t| {
        let radius = curve.evaluate(p(t)) - Point3::origin();
        let tangent = curve.derivative(p(t));
        assert!(
            radius.dot_mm2(tangent).abs() < 1.0e-6,
            "tangent is perpendicular to the radius at t={t}"
        );
    });
}

#[test]
fn curve_bounding_box_covers_the_control_hull() {
    let curve = quarter_circle();
    let bbox = curve.bounding_box();
    (0..=32).for_each(|i| {
        let t = f64::from(i) / 32.0;
        assert!(bbox.contains(curve.evaluate(p(t))));
    });
}

#[test]
fn curve_tessellation_stays_within_chord_tolerance() {
    let curve = quarter_circle();
    let chord = ChordHeightTolerance::from_mm(0.01);
    let points = curve.tessellate(chord);
    assert!(points.len() >= 2);
    points.iter().for_each(|point| {
        let (x, y, _) = point.coords_mm();
        assert!(
            (x.hypot(y) - 1.0).abs() < 1.0e-9,
            "tessellated points stay on the circle"
        );
    });
    assert_eq!(curve.tessellate(chord), curve.tessellate(chord));
}

#[test]
fn curve_closest_point_lands_on_the_curve() {
    let curve = quarter_circle();
    let near = curve.closest_point(Point3::from_mm(2.0, 2.0, 0.0), TOL);
    let (x, y, _) = near.point().coords_mm();
    assert!((x.hypot(y) - 1.0).abs() < 1.0e-6);
    assert!((near.point() - Point3::from_mm(FRAC_1_SQRT_2, FRAC_1_SQRT_2, 0.0)).norm_mm() < 0.05);
}

fn bilinear_plane() -> NurbsSurface3 {
    let Ok(surface) = NurbsSurface3::new(
        degree(1),
        degree(1),
        knots(&[0.0, 0.0, 1.0, 1.0]),
        knots(&[0.0, 0.0, 1.0, 1.0]),
        vec![
            vec![control(0.0, 0.0, 0.0, 1.0), control(0.0, 4.0, 0.0, 1.0)],
            vec![control(3.0, 0.0, 0.0, 1.0), control(3.0, 4.0, 0.0, 1.0)],
        ],
    ) else {
        panic!("a 2x2 bilinear patch is valid");
    };
    surface
}

fn quarter_cylinder() -> NurbsSurface3 {
    let Ok(surface) = NurbsSurface3::new(
        degree(2),
        degree(1),
        knots(&[0.0, 0.0, 0.0, 1.0, 1.0, 1.0]),
        knots(&[0.0, 0.0, 1.0, 1.0]),
        vec![
            vec![control(1.0, 0.0, 0.0, 1.0), control(1.0, 0.0, 5.0, 1.0)],
            vec![
                control(1.0, 1.0, 0.0, FRAC_1_SQRT_2),
                control(1.0, 1.0, 5.0, FRAC_1_SQRT_2),
            ],
            vec![control(0.0, 1.0, 0.0, 1.0), control(0.0, 1.0, 5.0, 1.0)],
        ],
    ) else {
        panic!("the rational quarter cylinder is valid");
    };
    surface
}

#[test]
fn bilinear_patch_is_planar() {
    let surface = bilinear_plane();
    let center = surface.evaluate(p(0.5), p(0.5));
    assert!((center - Point3::from_mm(1.5, 2.0, 0.0)).norm_mm() < 1.0e-9);
    let normal = surface.normal(p(0.5), p(0.5));
    let (nx, ny, nz) = normal.components();
    assert!(nx.abs() < 1.0e-9 && ny.abs() < 1.0e-9 && nz.abs() > 0.999);
    let (du, dv) = surface.partials(p(0.5), p(0.5));
    assert!((du - Vec3::from_mm(3.0, 0.0, 0.0)).norm_mm() < 1.0e-6);
    assert!((dv - Vec3::from_mm(0.0, 4.0, 0.0)).norm_mm() < 1.0e-6);
}

#[test]
fn patch_membership_respects_the_surface() {
    let surface = bilinear_plane();
    let on = Tolerance::new(1.0e-6);
    assert!(surface.contains_point(surface.evaluate(p(0.5), p(0.5)), on));
    assert!(!surface.contains_point(Point3::from_mm(1.5, 2.0, 5.0), on));
}

#[test]
fn patch_membership_holds_between_grid_nodes() {
    let surface = quarter_cylinder();
    let on = Tolerance::new(1.0e-6);
    [(0.37, 0.62), (0.08, 0.5), (0.93, 0.21)]
        .into_iter()
        .for_each(|(u, v)| {
            let point = surface.evaluate(p(u), p(v));
            assert!(
                surface.contains_point(point, on),
                "a point at ({u}, {v}) lies on the surface"
            );
            let (nx, ny, nz) = surface.normal(p(u), p(v)).components();
            let lifted = point + Vec3::from_mm(nx, ny, nz) * 0.5;
            assert!(
                !surface.contains_point(lifted, on),
                "a half-millimeter normal offset at ({u}, {v}) is off the surface"
            );
        });
}

#[test]
fn rational_surface_traces_a_cylinder() {
    let surface = quarter_cylinder();
    (0..=8).for_each(|around| {
        (0..=4).for_each(|along| {
            let su = f64::from(around) / 8.0;
            let sv = f64::from(along) / 4.0;
            let point = surface.evaluate(p(su), p(sv));
            let (x, y, z) = point.coords_mm();
            assert!(
                (x.hypot(y) - 1.0).abs() < 1.0e-9,
                "every point sits on the unit cylinder: radius {}",
                x.hypot(y)
            );
            assert!((z - sv * 5.0).abs() < 1.0e-9);
        });
    });
}

#[test]
fn cylinder_normal_points_radially_outward() {
    let surface = quarter_cylinder();
    let point = surface.evaluate(p(0.5), p(0.5));
    let normal = surface.normal(p(0.5), p(0.5));
    let (px, py, _) = point.coords_mm();
    let (nx, ny, nz) = normal.components();
    assert!(nz.abs() < 1.0e-6, "the normal has no axial component");
    let radial = px * nx + py * ny;
    assert!(radial.abs() > 0.999, "the normal is radial");
}

#[test]
fn surface_tessellation_winding_follows_normals() {
    let surface = quarter_cylinder();
    let mesh: TriMesh = surface.tessellate(
        ChordHeightTolerance::from_mm(0.05),
        AngleTolerance::from_radians(0.3),
    );
    let verts = mesh.vertices();
    mesh.triangles().iter().for_each(|&[a, b, c]| {
        let va = verts[a as usize];
        let vb = verts[b as usize];
        let vc = verts[c as usize];
        let face = (vb.position() - va.position()).cross(vc.position() - va.position());
        [va, vb, vc].iter().for_each(|vert| {
            let (nx, ny, nz) = vert.normal().components();
            let dot = face.dot_mm2(Vec3::from_mm(nx, ny, nz));
            assert!(dot > 0.0, "triangle [{a},{b},{c}] winds with its normal");
        });
    });
}

#[test]
fn surface_normal_matches_partial_cross_product() {
    let surface = quarter_cylinder();
    [(0.0, 0.5), (0.5, 0.25), (1.0, 1.0), (0.3, 0.7)]
        .into_iter()
        .for_each(|(u, v)| {
            let (du, dv) = surface.partials(p(u), p(v));
            let Ok(cross) = du.cross(dv).try_normalize(TOL) else {
                panic!("partials are non-degenerate at ({u}, {v})");
            };
            let normal = surface.normal(p(u), p(v));
            let (cx, cy, cz) = cross.components();
            let (nx, ny, nz) = normal.components();
            assert!(
                (cx - nx).abs() < 1.0e-9 && (cy - ny).abs() < 1.0e-9 && (cz - nz).abs() < 1.0e-9,
                "normal agrees with normalize(du x dv) at ({u}, {v})"
            );
        });
}

#[test]
fn curve_closest_point_refines_past_the_seed_grid() {
    let curve = quarter_circle();
    let target = curve.evaluate(p(0.37));
    let near = curve.closest_point(target, TOL);
    assert!(
        near.distance().get::<millimeter>() < 1.0e-7,
        "a point on the curve projects back onto it within tolerance"
    );
    assert!((near.parameter().value() - 0.37).abs() < 1.0e-4);
}

#[test]
fn surface_tessellation_responds_to_angle_tolerance() {
    let surface = quarter_cylinder();
    let slack_chord = ChordHeightTolerance::from_mm(1.0e3);
    let coarse = surface
        .tessellate(slack_chord, AngleTolerance::from_radians(1.5))
        .vertices()
        .len();
    let fine = surface
        .tessellate(slack_chord, AngleTolerance::from_radians(0.05))
        .vertices()
        .len();
    assert!(
        fine > coarse,
        "tightening the angle tolerance refines the mesh: coarse={coarse}, fine={fine}"
    );
}

#[test]
fn high_multiplicity_knots_evaluate_finite() {
    let Ok(curve) = NurbsCurve3::new(
        degree(2),
        knots(&[0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0]),
        vec![
            control(0.0, 0.0, 0.0, 1.0),
            control(1.0, 1.0, 0.0, 1.0),
            control(2.0, 0.0, 0.0, 1.0),
            control(3.0, 1.0, 0.0, 1.0),
            control(4.0, 0.0, 0.0, 1.0),
            control(5.0, 1.0, 0.0, 1.0),
        ],
    ) else {
        panic!("the knot count matches six degree-two control points");
    };
    (0..=20).for_each(|i| {
        let point = curve.evaluate(p(f64::from(i) / 20.0));
        let (x, y, z) = point.coords_mm();
        assert!(
            x.is_finite() && y.is_finite() && z.is_finite(),
            "a knot of multiplicity equal to the degree never divides by zero"
        );
    });
}

#[test]
fn degenerate_nurbs_inputs_are_rejected() {
    assert!(Degree::new(0).is_err());
    assert!(KnotVector::new(vec![1.0]).is_err());
    assert!(KnotVector::new(vec![0.0, 0.0]).is_err());
    assert!(KnotVector::new(vec![1.0, 0.0]).is_err());
    assert!(RationalControlPoint::new(Point3::origin(), 0.0, TOL).is_err());
    assert!(RationalControlPoint::new(Point3::origin(), f64::NAN, TOL).is_err());
    assert!(
        NurbsCurve3::new(
            degree(2),
            knots(&[0.0, 0.0, 1.0, 1.0]),
            vec![control(0.0, 0.0, 0.0, 1.0), control(1.0, 0.0, 0.0, 1.0)],
        )
        .is_err(),
        "two control points cannot carry a degree-two curve"
    );
    assert!(
        NurbsCurve3::new(
            degree(2),
            knots(&[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]),
            vec![
                control(0.0, 0.0, 0.0, 1.0),
                control(1.0, 0.0, 0.0, 1.0),
                control(2.0, 0.0, 0.0, 1.0),
            ],
        )
        .is_err(),
        "an unclamped knot vector does not interpolate its endpoints"
    );
}
