Primitive sphere creation from points

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


annotation { "Feature Type Name" : "Sphere 3 Point" }
export const sphereThreePoint = defineFeature(function(context is Context, id is Id, definition is map)
    precondition
    {
         annotation { "Name" : "Points", "Filter" : QueryFilterCompound.ALLOWS_VERTEX, "MaxNumberOfPicks" : 3 }
         definition.points is Query;
         
    }
    {
        var pts = [undefined,undefined,undefined];
        for(var n=0;n<3;n+=1)
        {
            pts[n] =evVertexPoint(context, {
                    "vertex" : qNthElement(definition.points, n)
            });
        }
      
       var cir =getCircleFromPoints(pts[0],pts[1],pts[2]);
       
       opSphere(context, id, {
               "radius" : cir.radius,
               "center" : cir.coordSystem.origin
       });
       
    }); 
function getCircleFromPoints(a is Vector, b is Vector, c is Vector)
{
    try silent
    { 
        var Xa = a[0] / millimeter;
        var Ya = a[1] / millimeter;
        var Za = a[2] / millimeter;
        var Xb = b[0] / millimeter;
        var Yb = b[1] / millimeter;
        var Zb = b[2] / millimeter;
        var Xc = c[0] / millimeter;
        var Yc = c[1] / millimeter;
        var Zc = c[2] / millimeter;
        // Lengths of AB, AC, AC
        var AB = (((Xa - Xb) ^ 2) + ((Ya - Yb) ^ 2) + ((Za - Zb) ^ 2)) ^ 0.5;
        var BC = (((Xb - Xc) ^ 2) + ((Yb - Yc) ^ 2) + ((Zb - Zc) ^ 2)) ^ 0.5;
        var AC = (((Xa - Xc) ^ 2) + ((Ya - Yc) ^ 2) + ((Za - Zc) ^ 2)) ^ 0.5;
        //Direction cosines of AB(ABi,ABj,ABk)
        var ABi = (Xb - Xa) / AB;
        var ABj = (Yb - Ya) / AB;
        var ABk = (Zb - Za) / AB;
        //'   Direction cosines of AC(ACi,ACj,ACk)
        //var ACi = (Xc - Xa) / AC;
        //var ACj = (Yc - Ya) / AC;
        //var ACk = (Zc - Za) / AC;
        //'   Cosine of angle BAC
        var cosBAC = (AB ^ 2 + AC ^ 2 - BC ^ 2) / (2 * AB * AC);
        var AD = cosBAC * AC;
        var CD = (AC ^ 2 - AD ^ 2) ^ 0.5;
        //'   Position of point D, which is C projected normally onto AB
        var Xd = Xa + (AD * ABi);
        var Yd = Ya + (AD * ABj);
        var Zd = Za + (AD * ABk);
        //'   Direction cosines of CD(Cdi,CDj,CDk)
        var CDi = (Xc - Xd) / CD;
        var CDj = (Yc - Yd) / CD;
        var CDk = (Zc - Zd) / CD;
        //'   Direction cosines of normal to AB and CD
        //'   — to be used for rotations of circle centre
        var Ni = (ABk * CDj) - (ABj * CDk);
        var Nj = (ABi * CDk) - (ABk * CDi);
        var Nk = (ABj * CDi) - (ABi * CDj);
        //'   # Diameter of circumscribed circle of a triangle is equal to the
        //'   the length of any side divided by sine of the opposite angle.
        //'   This is done in a coordinate system where X is colinear with AB, Y is // to CD,
        //'   and Z is the normal (N) to X and Y, and the origin is point A
        //'         R = D / 2
        var sinBAC = (1 - cosBAC ^ 2) ^ 0.5;
        var R = (BC / sinBAC) / 2;
        //'   Centre of circumscribed circle is point E
        var X2e = AB / 2;
        var Y2e = (R ^ 2 - X2e ^ 2) ^ 0.5;
        var Z2e = 0;
        //'   Transform matrix
        //'                    Rotations                 Translations
        //'             ——————————————————————————————————————————————
        //'               ABi  ,   ABj  ,  ABk                 Xa
        //'               CDi  ,   CDj  ,  CDk                 Ya
        //'                Ni  ,    Nj  ,   Nk                 Za
        //'             ——————————————————————————————————————————————
        //'   Position of circle centre in absolute axis system
        var X_centre = Xa + (X2e * ABi) + (Y2e * CDi) + (Z2e * Ni);
        var Y_centre = Ya + (X2e * ABj) + (Y2e * CDj) + (Z2e * Nj);
        var Z_centre = Za + (X2e * ABk) + (Y2e * CDk) + (Z2e * Nk);
        //'   ———————————————————————————————————————————————————————————————————

        var circleCenter = vector(X_centre * millimeter, Y_centre * millimeter, Z_centre * millimeter);
        var radius = pointDistance(circleCenter, a);
        var plane = try(planeFrom3Points(a, b, c));
        if (plane != undefined)
        {
            plane.origin = circleCenter;
            var coordsys = coordSystem(plane);
            return circle(coordsys, radius);
        }
    }
    return undefined;
}
export function pointDistance(vec0 is Vector, vec1 is Vector) returns ValueWithUnits
{ 
    var vec = [vec0[0].value - vec1[0].value, vec0[1].value - vec1[1].value, vec0[2].value - vec1[2].value]; //makeArray(size(vec0), undefined);
    return   sqrt(vec[0]*vec[0] + vec[1] *vec[1] + vec[2] *vec[2]) * meter;  
     
}
export function planeFrom3Points(pt0 is Vector, pt1 is Vector, pt2 is Vector)
{
    
    var points = [pt0, pt1, pt2];
    var normal = cross(points[2] - points[0], points[1] - points[0]);
    try silent
    {
        //check if the 3 points are on a line;
        var vec1 = points[0] - points[1];
        var vec2 = points[2] - points[1];
        if (parallelVectors(normalize(vec1), normalize(vec2)))
        {
            return undefined;
        }
    }
    try silent
    {
        var ret = plane(points[0], normalize(normal), normalize(points[1] - points[0]));
        return ret;
    }
    return undefined;
}