FeatureScript 3008;
import(path : "onshape/std/geometry.fs", version : "3008.0");

// =====================================================================
// GOLF BALL CUSTOM FEATURE
// =====================================================================

export enum GolfBallPatternType
{
    annotation { "Name" : "Icosahedral (geodesic mosaic)" }
    ICOSAHEDRAL,
    annotation { "Name" : "Hexagonal (latitude rows)" }
    HEXAGONAL,
    annotation { "Name" : "Tetrahedral (geodesic mosaic)" }
    TETRAHEDRAL
}

annotation { "Feature Type Name" : "Golf Ball Pattern" }
export const golfBall = defineFeature(function(context is Context, id is Id, definition is map)
    precondition
    {
        annotation { "Name" : "Center point (optional - defaults to origin)",
                     "Filter" : EntityType.VERTEX,
                     "MaxNumberOfPicks" : 1 }
        definition.centerPoint is Query;

        annotation { "Name" : "Ball diameter" }
        isLength(definition.ballDiameter, { (millimeter) : [10, 42.67, 80] } as LengthBoundSpec);

        annotation { "Name" : "Dimple diameter" }
        isLength(definition.dimpleDiameter, { (millimeter) : [3.3, 4.1, 4.1] } as LengthBoundSpec);

        annotation { "Name" : "Dimple depth" }
        isLength(definition.dimpleDepth, { (millimeter) : [0.2, 0.25, 0.3] } as LengthBoundSpec);

        annotation { "Name" : "Target dimple count" }
        isInteger(definition.dimpleCount, { (unitless) : [300, 322, 500] } as IntegerBoundSpec);

        annotation { "Name" : "Pattern type" }
        definition.patternType is GolfBallPatternType;
    }
    {
        const R = definition.ballDiameter / 2;
        var centerPos = vector(0, 0, 0) * meter;
        if (size(evaluateQuery(context, definition.centerPoint)) > 0)
        {
            centerPos = evVertexPoint(context, { "vertex" : definition.centerPoint });
        }

        // ---- 1. Base sphere ----------------------------------------
        const ballId = id + "ball";
        opSphere(context, ballId, {
            "center" : centerPos,
            "radius" : R
        });
        const ballBody = qCreatedBy(ballId, EntityType.BODY);

        // ---- 2. Dimple direction lattice ----------------------------
        var directions is array = [];
        if (definition.patternType == GolfBallPatternType.ICOSAHEDRAL)
        {
            directions = icosahedralDirections(definition.dimpleCount);
            directions = filterMinSpacing(directions, minSpacingDotFor(size(directions)));
        }
        else if (definition.patternType == GolfBallPatternType.TETRAHEDRAL)
        {
            directions = tetrahedralDirections(definition.dimpleCount);
            directions = filterMinSpacing(directions, minSpacingDotFor(size(directions)));
        }
        else
        {
            directions = hexagonalDirections(definition.dimpleCount);
        }

        // ---- 3. Dimple cutting-sphere radius (spherical cap math) ---
        // For a spherical cap of chord radius a = dimpleDiameter/2 and
        // depth h, the generating sphere radius is:
        //     r_d = (a^2 + h^2) / (2h)
        const a = definition.dimpleDiameter / 2;
        const h = definition.dimpleDepth;
        const dimpleR = (a * a + h * h) / (2 * h);

        // The cutting sphere's center sits along the radial direction
        // at distance (R - h + r_d) from the ball center, so its near
        // surface dents the ball to exactly depth h at the dimple
        // center, and its far surface meets the nominal ball surface
        // at radius a (chord radius) - i.e. dimple diameter is exact.
        const centerOffset = R - h + dimpleR;

        // ---- 4. Create one cutting sphere per lattice point ---------
        var toolIds is array = [];
        for (var i = 0; i < size(directions); i += 1)
        {
            const dir = directions[i];
            const dimpleCenter = centerPos + dir * centerOffset;
            const toolId = id + ("dimple" ~ i);
            opSphere(context, toolId, {
                "center" : dimpleCenter,
                "radius" : dimpleR
            });
            toolIds = append(toolIds, toolId);
        }

        var toolQueries is array = mapArray(toolIds, function(tId)
        {
            return qCreatedBy(tId, EntityType.BODY);
        });
        const toolQuery = qUnion(toolQueries);

        // ---- 5. Subtract all dimples from the ball in one boolean ---
        opBoolean(context, id + "cut", {
            "tools" : toolQuery,
            "targets" : ballBody,
            "operationType" : BooleanOperationType.SUBTRACTION
        });
    });

// =====================================================================
// Hexagonal latitude-ring lattice
//
// Each ring's rotation (offsetPhi) is chosen by bestRingOffset() below,
// which searches candidate rotations and picks whichever maximizes the
// minimum angular distance to the PREVIOUS ring's points. A fixed
// golden-angle rotation avoids long-range periodic repeats across many
// rings, but says nothing about any single adjacent pair - two
// neighboring rings can still land in a close phase by coincidence
// (worse when their point counts differ by rounding), producing a
// tight local cluster like the one circled in testing. Explicitly
// optimizing each ring against its neighbor removes that coincidence
// instead of hoping around it.
// =====================================================================
function hexagonalDirections(targetCount is number) returns array
{
    var numRings = round(sqrt(targetCount * PI / 4));
    if (numRings < 1)
    {
        numRings = 1;
    }

    // Ring weight ~ sin(theta), i.e. proportional to the ring's
    // circumference, so equatorial rings get more points than rings
    // near the top/bottom.
    var thetas is array = [];
    var weights is array = [];
    var weightSum = 0;
    for (var i = 0; i < numRings; i += 1)
    {
        const theta = (PI * (i + 0.5) / numRings) * radian;
        thetas = append(thetas, theta);
        const w = sin(theta);
        weights = append(weights, w);
        weightSum = weightSum + w;
    }

    var points is array = [];
    var prevRingPoints is array = [];
    for (var i = 0; i < numRings; i += 1)
    {
        const theta = thetas[i];
        var n = round(targetCount * weights[i] / weightSum);
        if (n < 3)
        {
            n = 3;
        }

        var offsetPhi = 0 * radian;
        if (i > 0)
        {
            offsetPhi = bestRingOffset(prevRingPoints, n, theta);
        }

        var ringPoints is array = [];
        for (var j = 0; j < n; j += 1)
        {
            const phi = (2 * PI * j / n) * radian + offsetPhi;
            const x = sin(theta) * cos(phi);
            const y = sin(theta) * sin(phi);
            const z = cos(theta);
            const p = vector(x, y, z);
            ringPoints = append(ringPoints, p);
            points = append(points, p);
        }
        prevRingPoints = ringPoints;
    }
    return points;
}

// =====================================================================
// bestRingOffset
// Searches candidate rotations of a ring (n points at fixed latitude
// theta) over one full point-to-point period (2*PI/n - beyond that the
// pattern just repeats), and returns whichever rotation maximizes the
// SMALLEST angular distance from any point in this ring to any point
// in the previous ring. This directly prevents the "two rings happen
// to phase-align" clustering that a fixed golden-angle offset can't
// guard against for a specific adjacent pair.
// =====================================================================
function bestRingOffset(prevPoints is array, n is number, theta is ValueWithUnits) returns ValueWithUnits
{
    if (size(prevPoints) == 0)
    {
        return 0 * radian;
    }

    const period = 2 * PI / n;
    const numCandidates = 48;

    var bestOffset = 0 * radian;
    var bestMinDist = -1 * radian;

    for (var c = 0; c < numCandidates; c += 1)
    {
        const candidateOffset = (period * c / numCandidates) * radian;

        var minDist = 1000 * radian;
        for (var j = 0; j < n; j += 1)
        {
            const phi = (2 * PI * j / n) * radian + candidateOffset;
            const x = sin(theta) * cos(phi);
            const y = sin(theta) * sin(phi);
            const z = cos(theta);
            const p = vector(x, y, z);

            for (var k = 0; k < size(prevPoints); k += 1)
            {
                var d = dotProduct(p, prevPoints[k]);
                if (d > 1) { d = 1; }
                if (d < -1) { d = -1; }
                const ang = acos(d);
                if (ang < minDist)
                {
                    minDist = ang;
                }
            }
        }

        if (minDist > bestMinDist)
        {
            bestMinDist = minDist;
            bestOffset = candidateOffset;
        }
    }

    return bestOffset;
}

// =====================================================================
// Spherical linear interpolation (slerp)
// =====================================================================
function slerp(A is Vector, B is Vector, t is number) returns Vector
{
    var d = dotProduct(A, B);
    if (d > 1) { d = 1; }
    if (d < -1) { d = -1; }
    const omega = acos(d);
    const sinOmega = sin(omega);

    if (sinOmega < 0.000001)
    {
        const blend = A + (B - A) * t;
        return blend / norm(blend);
    }

    const coeffA = sin((1 - t) * omega) / sinOmega;
    const coeffB = sin(t * omega) / sinOmega;
    return A * coeffA + B * coeffB;
}

// =====================================================================
// Double-slerp triangle subdivision
// =====================================================================
function subdivideTriangleSlerp(A is Vector, B is Vector, C is Vector, frequency is number) returns array
{
    var points is array = [];
    for (var i = 0; i <= frequency; i += 1)
    {
        const S = slerp(A, B, i / frequency);
        if (i == frequency)
        {
            points = append(points, S);
        }
        else
        {
            const E = slerp(B, C, (frequency - i) / frequency);
            const rowCount = frequency - i;
            for (var j = 0; j <= rowCount; j += 1)
            {
                points = append(points, slerp(S, E, j / rowCount));
            }
        }
    }
    return points;
}

// =====================================================================
// Minimum-spacing filter
// =====================================================================
function filterMinSpacing(points is array, minDot is number) returns array
{
    var accepted is array = [];
    for (var idx = 0; idx < size(points); idx += 1)
    {
        const candidate = points[idx];
        var tooClose = false;
        for (var jdx = 0; jdx < size(accepted); jdx += 1)
        {
            if (dotProduct(candidate, accepted[jdx]) > minDot)
            {
                tooClose = true;
            }
        }
        if (!tooClose)
        {
            accepted = append(accepted, candidate);
        }
    }
    return accepted;
}

function minSpacingDotFor(count is number) returns number
{
    const avgSpacingRadians = sqrt(4 * PI / max(1, count));
    return cos(avgSpacingRadians * 0.2 * radian);
}

function dotProduct(a is Vector, b is Vector) returns number
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

// =====================================================================
// Icosahedral geodesic lattice
// =====================================================================
function icosahedralDirections(targetCount is number) returns array
{
    const rawF = sqrt(max(1, (targetCount - 2) / 10));
    var frequency = round(rawF);
    if (frequency < 1)
    {
        frequency = 1;
    }

    return subdivideAndDedup(icosahedronVertices(), icosahedronFaces(), frequency);
}

function icosahedronVertices() returns array
{
    const phi = (1 + sqrt(5)) / 2;
    const raw = [
        vector(-1, phi, 0), vector(1, phi, 0), vector(-1, -phi, 0), vector(1, -phi, 0),
        vector(0, -1, phi), vector(0, 1, phi), vector(0, -1, -phi), vector(0, 1, -phi),
        vector(phi, 0, -1), vector(phi, 0, 1), vector(-phi, 0, -1), vector(-phi, 0, 1)
    ];
    return mapArray(raw, function(v) { return v / norm(v); });
}

function icosahedronFaces() returns array
{
    return [
        [0, 11, 5], [0, 5, 1], [0, 1, 7], [0, 7, 10], [0, 10, 11],
        [1, 5, 9], [5, 11, 4], [11, 10, 2], [10, 7, 6], [7, 1, 8],
        [3, 9, 4], [3, 4, 2], [3, 2, 6], [3, 6, 8], [3, 8, 9],
        [4, 9, 5], [2, 4, 11], [6, 2, 10], [8, 6, 7], [9, 8, 1]
    ];
}

function pointKey(p is Vector) returns string
{
    const scale = 100000;
    const rx = round(p[0] * scale);
    const ry = round(p[1] * scale);
    const rz = round(p[2] * scale);
    return (rx ~ "_") ~ (ry ~ "_") ~ (rz ~ "");
}

function subdivideAndDedup(verts is array, faces is array, frequency is number) returns array
{
    var seen is map = {};
    var points is array = [];

    for (var fIdx = 0; fIdx < size(faces); fIdx += 1)
    {
        const tri = faces[fIdx];
        const A = verts[tri[0]];
        const B = verts[tri[1]];
        const C = verts[tri[2]];

        const facePoints = subdivideTriangleSlerp(A, B, C, frequency);
        for (var pIdx = 0; pIdx < size(facePoints); pIdx += 1)
        {
            const p = facePoints[pIdx];
            const key = pointKey(p);
            if (seen[key] == undefined)
            {
                seen[key] = true;
                points = append(points, p);
            }
        }
    }
    return points;
}

// =====================================================================
// Tetrahedral geodesic lattice
// =====================================================================
function tetrahedronVertices() returns array
{
    const raw = [
        vector(1, 1, 1), vector(1, -1, -1), vector(-1, 1, -1), vector(-1, -1, 1)
    ];
    return mapArray(raw, function(v) { return v / norm(v); });
}

function tetrahedronFaces() returns array
{
    return [
        [0, 1, 2], [0, 3, 1], [0, 2, 3], [1, 3, 2]
    ];
}

function tetrahedralDirections(targetCount is number) returns array
{
    const rawF = sqrt(max(1, (targetCount - 2) / 2));
    var frequency = round(rawF * 1.3);
    if (frequency < 1)
    {
        frequency = 1;
    }

    return subdivideAndDedup(tetrahedronVertices(), tetrahedronFaces(), frequency);
}