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

// =============================================================================
//  SPHERE PATTERN WRAP (v21 - Fibonacci, Adaptive Symmetry & Polar Injection)
// =============================================================================

export enum OutputMode
{
    annotation {"Name" : "Solid" } SOLID,
    annotation {"Name" : "Surface" } SURFACE
}

export enum BooleanMode
{
    annotation {"Name" : "Add" } EMBOSS,
    annotation {"Name" : "Subtract" } DEBOSS
}

export enum SpherePatternArrangement
{
    annotation { "Name" : "Golden Spiral (Fibonacci)" } FIBONACCI,
    annotation { "Name" : "Circular Rings" } RINGS,
    annotation { "Name" : "Honeycomb" } HONEYCOMB,
    annotation { "Name" : "UV Grid" } UV_GRID
}

export enum PatternProfileShape
{
    annotation { "Name" : "Circle" } CIRCLE,
    annotation { "Name" : "Square" } SQUARE,
    annotation { "Name" : "Diamond" } DIAMOND,
    annotation { "Name" : "Hexagon" } HEXAGON
}

annotation { "Feature Type Name" : "Sphere Pattern Wrap" }
export const spherePatternWrap = defineFeature(function(context is Context, id is Id, definition is map)
    precondition
    {
        annotation { "Name" : "Output", "UIHint" : [UIHint.HORIZONTAL_ENUM, UIHint.REMEMBER_PREVIOUS_VALUE] }
        definition.outputMode is OutputMode;

        if (definition.outputMode == OutputMode.SOLID)
        {
            annotation { "Name" : "Operation", "UIHint" : [UIHint.HORIZONTAL_ENUM, UIHint.REMEMBER_PREVIOUS_VALUE] }
            definition.booleanMode is BooleanMode;
        }

        annotation { "Name" : "Target spherical face", "Filter" : EntityType.FACE, "MaxNumberOfPicks" : 1 }
        definition.targetFace is Query;

        annotation { "Name" : "Profile shape", "UIHint" : UIHint.REMEMBER_PREVIOUS_VALUE }
        definition.profileShape is PatternProfileShape;

        annotation { "Name" : "Profile size (width)" }
        isLength(definition.profileSize, { (millimeter) : [0.01, 5.0, 500.0], (inch) : 0.2 } as LengthBoundSpec);

        if (definition.profileShape == PatternProfileShape.DIAMOND)
        {
            annotation { "Name" : "Diamond ratio", "Description" : "Width : height ratio. 1.0 = symmetric, >1 = wider, <1 = taller" }
            isReal(definition.diamondRatio, { (unitless) : [0.2, 1.0, 5.0] } as RealBoundSpec);
        }

        annotation { "Name" : "Pattern arrangement", "UIHint" : UIHint.REMEMBER_PREVIOUS_VALUE }
        definition.patternArrangement is SpherePatternArrangement;

        annotation { "Name" : "Equatorial count", "Description": "Number of shapes around the equator" }
        isInteger(definition.equatorialCount, { (unitless) : [3, 9, 500] } as IntegerBoundSpec);

        annotation { "Name" : "Latitude limit (degrees)", "Description": "Stop patterning this close to the poles (90 = full sphere)", "Default" : 90 }
        isReal(definition.latLimit, { (unitless) : [0, 90, 90] } as RealBoundSpec);

        if (definition.outputMode == OutputMode.SOLID)
        {
            annotation { "Name" : "Depth" }
            isLength(definition.solidDepth, { (millimeter) : [0.01, 0.5, 50.0], (inch) : 0.08 } as LengthBoundSpec);
        }

        if (definition.outputMode == OutputMode.SURFACE)
        {
            annotation { "Name" : "Split target face", "Default" : true }
            definition.splitFace is boolean;
        }

        annotation { "Name" : "Keep wrap curves", "Default" : false }
        definition.keepCurves is boolean;

        annotation { "Name" : "Skip locations", "Filter" : EntityType.VERTEX || BodyType.MATE_CONNECTOR || BodyType.POINT, "MaxNumberOfPicks" : 200 }
        definition.skipLocations is Query;
    }
    {
        if (isQueryEmpty(context, definition.targetFace))
            throw regenError("Select a spherical face.", ["targetFace"]);

        wrapGeneratedPattern(context, id, definition);
    });

// =============================================================================
//  GEOMETRY HELPERS
// =============================================================================

function getSphereInfo(context is Context, face is Query) returns map
{
    try silent
    {
        var surf = evSurfaceDefinition(context, {"face": face});
        if (surf.surfaceType == SurfaceType.SPHERE)
        {
            return {
                "center" : surf.coordSystem.origin,
                "R"      : surf.radius,
                "zAxis"  : surf.coordSystem.zAxis,
                "xAxis"  : surf.coordSystem.xAxis
            };
        }
    }
    var sg = fitSphereFromFace(context, face);
    return {"center": sg.centre, "R": sg.R, "zAxis": vector(0,0,1), "xAxis": vector(1,0,0)};
}

function mapToSphere(centerData is map, x is ValueWithUnits, y is ValueWithUnits) returns Vector
{
    var raw    = centerData.point + x * centerData.U + y * centerData.V;
    var dir    = raw - centerData.sphereCenter;
    var dirLen = norm(dir);
    if (dirLen > 1e-9 * meter)
        return centerData.sphereCenter + centerData.sphereR * (dir / dirLen);
    return centerData.point;
}

function getShapeRatio(definition is map) returns number
{
    var ratio = 1.0;
    if (definition.profileShape == PatternProfileShape.DIAMOND)
        try silent { ratio = definition.diamondRatio; }
    return ratio;
}

function cornerPointsForShape(r is ValueWithUnits, shape is PatternProfileShape, ratio is number) returns array
{
    if (shape == PatternProfileShape.SQUARE)
    {
        var a = r * 0.707106781;
        return [ vector(a,a), vector(-a,a), vector(-a,-a), vector(a,-a) ];
    }
    else if (shape == PatternProfileShape.DIAMOND)
    {
        var dx = r * sqrt(ratio);
        var dy = r / sqrt(ratio);
        return [ vector(dx, 0*meter), vector(0*meter, dy), vector(-dx, 0*meter), vector(0*meter, -dy) ];
    }
    else if (shape == PatternProfileShape.HEXAGON)
    {
        var corners = [];
        for (var i = 0; i < 6; i += 1)
        {
            var t = (i * PI / 3) * radian;
            corners = append(corners, vector(r*cos(t), r*sin(t)));
        }
        return corners;
    }
    return [];
}

function generateProfilePoints(centerData is map, r is ValueWithUnits, shape is PatternProfileShape, ratio is number, onPlane is boolean) returns array
{
    var pts = [];
    if (shape == PatternProfileShape.CIRCLE)
    {
        var nCirc = 48; // Efficient spline extrusion sample density
        for (var i = 0; i < nCirc; i += 1)
        {
            var t = (i * 2 * PI / nCirc) * radian;
            if (onPlane)
                pts = append(pts, centerData.point + r * cos(t) * centerData.U + r * sin(t) * centerData.V);
            else
                pts = append(pts, mapToSphere(centerData, r*cos(t), r*sin(t)));
        }
    }
    else
    {
        var corners = cornerPointsForShape(r, shape, ratio);
        var nEdges  = size(corners);
        var samples = 8;
        for (var i = 0; i < nEdges; i += 1)
        {
            var c1 = corners[i];
            var c2 = corners[(i+1) % nEdges];
            for (var j = 0; j < samples; j += 1)
            {
                var t = j / samples;
                var x = c1[0]*(1-t) + c2[0]*t;
                var y = c1[1]*(1-t) + c2[1]*t;
                if (onPlane)
                    pts = append(pts, centerData.point + x * centerData.U + y * centerData.V);
                else
                    pts = append(pts, mapToSphere(centerData, x, y));
            }
        }
    }
    return pts;
}

function buildPolylineEdges(context is Context, idBase is Id, pts is array, closed is boolean) returns Query
{
    var n        = size(pts);
    var segCount = closed ? n : (n - 1);
    var eqs      = [];
    for (var i = 0; i < segCount; i += 1)
    {
        var segId = idBase + ("seg" ~ toString(i));
        opFitSpline(context, segId, {
            "points"    : [pts[i], pts[(i+1) % n]],
            "closed"    : false,
            "tolerance" : 1e-7 * meter
        });
        eqs = append(eqs, qCreatedBy(segId, EntityType.EDGE));
    }
    return qUnion(eqs);
}

function buildProfileEdges(context is Context, idBase is Id, pts is array, shape is PatternProfileShape) returns Query
{
    if (shape == PatternProfileShape.CIRCLE)
    {
        var n           = size(pts);
        var numSegments = 8;
        var ptsPerSeg   = n / numSegments;
        var eqs         = [];
        for (var s = 0; s < numSegments; s += 1)
        {
            var segPts = [];
            for (var k = 0; k <= ptsPerSeg; k += 1)
                segPts = append(segPts, pts[(s * ptsPerSeg + k) % n]);
            var segId = idBase + ("circSeg" ~ toString(s));
            opFitSpline(context, segId, {
                "points"    : segPts,
                "closed"    : false,
                "tolerance" : 1e-7 * meter
            });
            eqs = append(eqs, qCreatedBy(segId, EntityType.EDGE));
        }
        return qUnion(eqs);
    }
    return buildPolylineEdges(context, idBase, pts, true);
}

function getQueryPoint(context is Context, query is Query)
{
    try silent
    {
        return evVertexPoint(context, { "vertex" : query });
    }
    try silent
    {
        var cs = evMateConnector(context, { "mateConnector" : query });
        return cs.origin;
    }
    try silent
    {
        return evApproximateCentroid(context, { "entities" : query });
    }
    return undefined;
}

// =============================================================================
//  PATTERN CENTER GENERATORS
// =============================================================================

function generateFibonacciCenters(context is Context, id is Id, definition is map, sInfo is map, minSep is ValueWithUnits, skipPoints is array) returns array
{
    var C = sInfo.center;
    var R = sInfo.R;
    var Z = sInfo.zAxis;
    var X = sInfo.xAxis;
    var Y = cross(Z, X);

    var count = definition.equatorialCount;
    var latMax = definition.latLimit * degree;

    // Estimate the total number of points on a full sphere to match the equatorial density.
    var d = (2 * PI * R) / count;
    var areaPerPoint = (sqrt(3) / 2) * d * d; // Hexagonal grid density estimate
    var sphereArea = 4 * PI * R * R;
    var totalPointsEstimate = max(3, round(sphereArea / areaPerPoint));

    var centers = [];
    var goldenRatio = (1 + sqrt(5)) / 2;
    var goldenAngle = 2 * PI * (1 - 1 / goldenRatio) * radian; // Safe golden angle with Angle units!

    var skipTolerance = definition.profileSize * 1.5;

    for (var i = 0; i < totalPointsEstimate; i += 1)
    {
        // Even z distribution spanning the full sphere height
        var t = i / (totalPointsEstimate - 1);
        var z = 1.0 - 2.0 * t;
        var radiusAtZ = sqrt(max(0.0, 1.0 - z * z));
        var theta = i * goldenAngle; // goldenAngle has units, so theta is an Angle!

        // Robust, unit-safe latitude restriction comparison
        if (abs(z) > sin(latMax))
            continue;

        var dir = radiusAtZ * cos(theta) * X + radiusAtZ * sin(theta) * Y + z * Z; // cos/sin now execute flawlessly!
        var pt = C + R * dir;

        // Check if user requested to skip this location
        var shouldSkip = false;
        for (var sp in skipPoints)
        {
            if (norm(pt - sp) < skipTolerance)
            {
                shouldSkip = true;
                break;
            }
        }
        if (shouldSkip)
            continue;

        // Global overlap validation check
        var tooClose = false;
        for (var c in centers)
        {
            if (norm(pt - c.point) < minSep)
            {
                tooClose = true;
                break;
            }
        }

        if (!tooClose)
        {
            // Tangent vectors for coordinates
            var vDir = -z * cos(theta) * X - z * sin(theta) * Y + radiusAtZ * Z;
            var vLen = norm(vDir);
            var V = (vLen > 1e-9) ? vDir / vLen : X;
            var U = cross(V, dir);
            var uLen = norm(U);
            U = (uLen > 1e-9) ? U / uLen : Y;

            centers = append(centers, {
                "point": pt,
                "dir": dir,
                "U": U,
                "V": V,
                "sphereCenter": C,
                "sphereR": R
            });
        }
    }
    return centers;
}

function wrapGeneratedPattern(context is Context, id is Id, definition is map)
{
    var sInfo = getSphereInfo(context, definition.targetFace);
    var C = sInfo.center;
    var R = sInfo.R;
    var Z = sInfo.zAxis;
    var X = sInfo.xAxis;
    var Y = cross(Z, X);

    var count = definition.equatorialCount;
    if (definition.profileSize > (2 * PI * R / count))
        reportFeatureInfo(context, id, "Warning: Profile size exceeds equator spacing.");
    if (definition.profileSize > R * 0.6)
        reportFeatureInfo(context, id, "Warning: Profile size is large relative to the sphere radius.");

    // Parse skip locations
    var skipPoints = [];
    if (!isQueryEmpty(context, definition.skipLocations))
    {
        var skippers = evaluateQuery(context, definition.skipLocations);
        for (var sk in skippers)
        {
            var p = getQueryPoint(context, sk);
            if (p != undefined)
            {
                skipPoints = append(skipPoints, p);
            }
        }
    }

    var centers = [];
    var minSep = definition.profileSize * 1.05; // 5% safety margin to guarantee no touches or overlaps
    var skipTolerance = definition.profileSize * 1.5;

    if (definition.patternArrangement == SpherePatternArrangement.FIBONACCI)
    {
        centers = generateFibonacciCenters(context, id, definition, sInfo, minSep, skipPoints);
    }
    else
    {
        var latMax = definition.latLimit * degree;
        var dTheta = (2 * PI / count) * radian;
        var dPhi   = dTheta;
        if (definition.patternArrangement == SpherePatternArrangement.HONEYCOMB)
            dPhi = dTheta * (sqrt(3) / 2);

        var jMax          = floor(latMax / dPhi);

        for (var j = -jMax; j <= jMax; j += 1)
        {
            var phi    = j * dPhi;
            var K      = count;
            var offset = 0 * radian;

            if (definition.patternArrangement == SpherePatternArrangement.RINGS)
                K = max(1, round(count * cos(phi)));
            else if (definition.patternArrangement == SpherePatternArrangement.HONEYCOMB)
            {
                K = max(1, round(count * cos(phi)));
                if (abs(j) % 2 == 1) offset = (PI / K) * radian;
            }

            // Reject sparse polar rings (K < 3) to enforce strict rotational symmetry.
            if (K < 3)
                continue;

            var dThetaRow  = (2 * PI / K) * radian;

            for (var i = 0; i < K; i += 1)
            {
                var theta = i * dThetaRow + offset;
                var dir   = cos(phi)*cos(theta)*X + cos(phi)*sin(theta)*Y + sin(phi)*Z;
                var pt    = C + R * dir;

                // Check if user requested to skip this location manually
                var shouldSkip = false;
                for (var sp in skipPoints)
                {
                    if (norm(pt - sp) < skipTolerance)
                    {
                        shouldSkip = true;
                        break;
                    }
                }
                if (shouldSkip)
                    continue;

                // We bypass adjacent-row culling to naturally distribute shapes.
                var vDir = -sin(phi)*cos(theta)*X - sin(phi)*sin(theta)*Y + cos(phi)*Z;
                var vLen = norm(vDir);
                var V    = (vLen > 1e-9) ? vDir/vLen : X;
                var U    = cross(V, dir);
                var uLen = norm(U);
                U        = (uLen > 1e-9) ? U/uLen : Y;

                centers = append(centers, { "point": pt, "dir": dir, "U": U, "V": V, "sphereCenter": C, "sphereR": R });
            }
        }

        // ---- Step 4.5: Explicit Dynamic Polar Injection ----
        // If the latitude limit spans the full sphere, explicitly test and place single 
        // entities exactly at the North Pole and South Pole.
        if (definition.latLimit >= 89.9)
        {
            var poles = [ { "pt": C + R * Z, "dir": Z }, { "pt": C - R * Z, "dir": -Z } ];
            for (var pole in poles)
            {
                var pt = pole.pt;
                var dir = pole.dir;
                
                // Check if poles are marked for skip
                var shouldSkip = false;
                for (var sp in skipPoints)
                {
                    if (norm(pt - sp) < skipTolerance)
                    {
                        shouldSkip = true;
                        break;
                    }
                }
                if (shouldSkip)
                    continue;

                // Spacing safety buffer check against nearest active rings
                var tooClose = false;
                for (var c in centers)
                {
                    if (norm(pt - c.point) < minSep)
                    {
                        tooClose = true;
                        break;
                    }
                }

                if (!tooClose)
                {
                    var U = X;
                    var V = Y;
                    if (dir[2] < 0)
                    {
                        U = -X; // Ensure consistent right-handed coordinate system
                    }
                    centers = append(centers, { "point": pt, "dir": dir, "U": U, "V": V, "sphereCenter": C, "sphereR": R });
                }
            }
        }
    }

    if (size(centers) == 0)
    {
        reportFeatureWarning(context, id, "No pattern locations were generated - check parameters.");
        return;
    }

    if (definition.outputMode == OutputMode.SOLID)
        applySplitAndOffset(context, id, definition, centers);
    else
        applySurfaceSplit(context, id, definition, centers);
}

// =============================================================================
//  SHARED: build wire loops per center
// =============================================================================

function buildAllWires(context is Context, id is Id, definition is map, centers is array, idPrefix is string, onPlane is boolean) returns map
{
    var radiusVal   = definition.profileSize / 2.0;
    var ratio       = getShapeRatio(definition);
    var edgeQueries = [];
    var wireBodies  = [];
    var idx         = 0;

    for (var c in centers)
    {
        var pts = generateProfilePoints(c, radiusVal, definition.profileShape, ratio, onPlane);
        if (size(pts) >= 3)
        {
            var loopId = id + (idPrefix ~ toString(idx));
            try silent
            {
                var edges = buildProfileEdges(context, loopId, pts, definition.profileShape);
                edgeQueries = append(edgeQueries, edges);
                wireBodies  = append(wireBodies, qOwnerBody(edges));
            }
        }
        idx += 1;
    }
    return { "edgeQueries": edgeQueries, "wireBodies": wireBodies };
}

// =============================================================================
//  SOLID MODE (Extruded Sheets Split & Offset Face)
// =============================================================================

function applySplitAndOffset(context is Context, id is Id, definition is map, centers is array)
{
    var depth    = definition.solidDepth;
    var isDeboss = (definition.booleanMode == BooleanMode.DEBOSS);

    // ---- Step 1: Generate flat wire profiles on each tangent plane ----
    var wireData    = buildAllWires(context, id, definition, centers, "flat_wire", true);
    var flatEdges   = wireData.edgeQueries;
    var flatBodies  = wireData.wireBodies;

    if (size(flatEdges) == 0)
    {
        reportFeatureWarning(context, id, "No profile curves could be generated.");
        return;
    }

    // ---- Step 2: Extrude the wires into surface (sheet) bodies cutting the sphere ----
    var sheetBodies = [];
    var extrudeDistance = max(5.0 * millimeter, depth * 5); 

    for (var idx = 0; idx < size(centers); idx += 1)
    {
        var c = centers[idx];
        var extrudeId = id + ("sheet_extrude" ~ toString(idx));
        try silent
        {
            opExtrude(context, extrudeId, {
                "entities"   : flatEdges[idx],
                "direction"  : c.dir,
                "endBound"   : BoundingType.BLIND,
                "endDepth"   : extrudeDistance,
                "startBound" : BoundingType.BLIND,
                "startDepth" : extrudeDistance
            });
            sheetBodies = append(sheetBodies, qCreatedBy(extrudeId, EntityType.BODY));
        }
    }

    if (size(sheetBodies) == 0)
    {
        reportFeatureWarning(context, id, "Failed to extrude slicing sheets from profiles.");
        cleanupWires(context, id, flatBodies, false);
        return;
    }

    // ---- Step 3: Run Split Face using the extruded sheets as splitting tools ----
    var splitId = id + "solid_split";
    try
    {
        opSplitFace(context, splitId, {
            "faceTargets" : definition.targetFace,
            "bodyTools"   : qUnion(sheetBodies)
        });
    }
    catch (error)
    {
        reportFeatureWarning(context, id, "Split face operation failed: " ~ error.message);
        cleanupWires(context, id, flatBodies, false);
        try silent { opDeleteBodies(context, id + "cleanSheets", { "entities" : qUnion(sheetBodies) }); }
        return;
    }

    // ---- Step 4: Isolate pattern faces using centroid proximity mapping ----
    var bornArray = evaluateQuery(context, definition.targetFace);
    if (size(bornArray) <= 1)
    {
        reportFeatureWarning(context, id, "Split produced no new faces. Ensure the profile size and arrangement are correct.");
        cleanupWires(context, id, flatBodies, false);
        try silent { opDeleteBodies(context, id + "cleanSheets", { "entities" : qUnion(sheetBodies) }); }
        return;
    }

    var patchFaces = [];
    var maxCentroidDist = definition.profileSize * 0.45; // Centroid of a valid cap matches its center point extremely closely

    for (var f in bornArray)
    {
        try silent
        {
            var centroid = evApproximateCentroid(context, { "entities": f });
            var matched = false;
            
            // Map centroid to see if it belongs to one of our pattern centers
            for (var c in centers)
            {
                if (norm(centroid - c.point) < maxCentroidDist)
                {
                    matched = true;
                    break;
                }
            }

            if (matched)
            {
                patchFaces = append(patchFaces, f);
            }
        }
    }

    if (size(patchFaces) == 0)
    {
        reportFeatureWarning(context, id, "Could not isolate small pattern faces. Try smaller profile sizes.");
        cleanupWires(context, id, flatBodies, false);
        try silent { opDeleteBodies(context, id + "cleanSheets", { "entities" : qUnion(sheetBodies) }); }
        return;
    }

    // ---- Step 5: Apply native Move Face/Offset to achieve 3D Emboss/Deboss ----
    var offsetDist = isDeboss ? -depth : depth;
    try
    {
        opOffsetFace(context, id + "solid_offset", {
            "moveFaces"      : qUnion(patchFaces),
            "offsetDistance" : offsetDist
        });
    }
    catch (error)
    {
        reportFeatureWarning(context, id, "Move/Offset face failed: " ~ error.message ~ ". Try a smaller Depth.");
    }

    // ---- Step 6: Clean up temporary flat wires and extruded sheet bodies ----
    cleanupWires(context, id, flatBodies, definition.keepCurves);
    try silent
    {
        opDeleteBodies(context, id + "cleanSheets", {
            "entities" : qUnion(sheetBodies)
        });
    }
}

// =============================================================================
//  SURFACE MODE
// =============================================================================

function applySurfaceSplit(context is Context, id is Id, definition is map, centers is array)
{
    var wireData    = buildAllWires(context, id, definition, centers, "surf", false);
    var edgeQueries = wireData.edgeQueries;
    var wireBodies  = wireData.wireBodies;

    if (size(edgeQueries) == 0)
    {
        reportFeatureWarning(context, id, "No projected curves were generated.");
        return;
    }

    var wireEdges = qUnion(edgeQueries);

    if (definition.splitFace)
    {
        try
        {
            opSplitFace(context, id + "split", {
                "faceTargets" : definition.targetFace,
                "edgeTools"   : wireEdges
            });
        }
        catch (error) { reportFeatureWarning(context, id, "Split failed: " ~ error.message); }
    }

    cleanupWires(context, id, wireBodies, definition.keepCurves);
}

// =============================================================================
//  SHARED: wire cleanup helper
// =============================================================================

function cleanupWires(context is Context, id is Id, wireBodies is array, keepCurves is boolean)
{
    if (!keepCurves && size(wireBodies) > 0)
        try silent { opDeleteBodies(context, id + "wClean", { "entities": qUnion(wireBodies) }); }
}

// =============================================================================
//  FALLBACK: fit sphere from sampled face normals
// =============================================================================

function fitSphereFromFace(context is Context, face is Query) returns map
{
    var uvSamples = [
        vector(0.2,0.2), vector(0.5,0.2), vector(0.8,0.2),
        vector(0.2,0.5), vector(0.5,0.5), vector(0.8,0.5),
        vector(0.2,0.8), vector(0.5,0.8), vector(0.8,0.8)
    ];
    var origins = []; var normals = [];
    for (var uv in uvSamples)
    {
        try silent
        {
            var tp = evFaceTangentPlane(context, { "face": face, "parameter": uv });
            origins = append(origins, tp.origin);
            normals = append(normals, tp.normal);
        }
    }
    if (size(origins) < 3) throw regenError("Cannot determine sphere geometry.", ["targetFace"]);

    var A  = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0];
    var bx = 0.0*meter; var by = 0.0*meter; var bz = 0.0*meter;
    for (var i = 0; i < size(origins); i += 1)
    {
        var p  = origins[i]; var n = normals[i];
        var nx = n[0]; var ny = n[1]; var nz = n[2];
        A[0] += 1-nx*nx; A[1] += -nx*ny;  A[2] += -nx*nz;
        A[3] += -ny*nx;  A[4] += 1-ny*ny; A[5] += -ny*nz;
        A[6] += -nz*nx;  A[7] += -nz*ny;  A[8] += 1-nz*nz;
        bx += (1-nx*nx)*p[0] + (-nx*ny)*p[1] + (-nx*nz)*p[2];
        by += (-ny*nx)*p[0]  + (1-ny*ny)*p[1] + (-ny*nz)*p[2];
        bz += (-nz*nx)*p[0]  + (-nz*ny)*p[1]  + (1-nz*nz)*p[2];
    }
    var det = A[0]*(A[4]*A[8]-A[5]*A[7]) - A[1]*(A[3]*A[8]-A[5]*A[6]) - A[2]*(A[3]*A[7]-A[4]*A[6]);
    var centre;
    if (abs(det) < 1e-10)
    {
        var cx = 0.0*meter; var cy = 0.0*meter; var cz = 0.0*meter;
        for (var p in origins) { cx += p[0]; cy += p[1]; cz += p[2]; }
        centre = vector(cx/size(origins), cy/size(origins), cz/size(origins));
    }
    else
    {
        var Cx = (bx*(A[4]*A[8]-A[5]*A[7]) - A[1]*(by*A[8]-A[5]*bz) + A[2]*(by*A[7]-A[4]*bz)) / det;
        var Cy = (A[0]*(by*A[8]-A[5]*bz)   - bx*(A[3]*A[8]-A[5]*A[6]) + A[2]*(A[3]*bz-by*A[6])) / det;
        var Cz = (A[0]*(A[4]*bz-by*A[7])   - A[1]*(A[3]*bz-by*A[6])   + bx*(A[3]*A[7]-A[4]*A[6])) / det;
        centre = vector(Cx, Cy, Cz);
    }
    var Rv = 0.0*meter;
    for (var p in origins) Rv += norm(p - centre);
    return {"centre": centre, "R": Rv / size(origins)};
}