float dist = chf("dist"); int divisions = chi("Divisions");
int corners[] = expandpointgroup(0, chs("group"));


int PrimPts[] = primpoints(0, @primnum), PrimPts_new[]=PrimPts;
int primCorners[]; //intersection of  current prim pts  &  selected pts group
foreach (int pt; PrimPts)
    if (find(corners, pt) >= 0 && chs("group") != "") 
        append(primCorners, pt); 

if(len(primCorners) > 0){
    foreach(int i; primCorners){
        vector corner = point(0, "P", i);

        int p1 = neighbour(0, i, 0);
        int p2 = neighbour(0, i, 1);
        vector P1 = point(0, 'P', p1); //2 neighbouring Pts
        vector P2 = point(0, 'P', p2);
        vector dir1 = normalize(corner - P1); //2 tangents' dirs
        vector dir2 = normalize(corner - P2);
        vector normal = normalize(cross(dir1, dir2)); //normal of the plane

        vector r1 = -normalize(cross(dir1, normal)); //2 radius dirs
        vector r2 = normalize(cross(dir2, normal));
        vector T1 = corner - dir1 * dist; //2 tangency Pts
        vector T2 = corner - dir2 * dist;
        vector C = T1 - r1 * length(cross(T2 - T1, r2)) / length(cross(r1, r2)); //centre of the circle
        float radius = distance(C, T1); 

        vector x = r1; vector z = normal; 
        vector y = normalize(cross(z, x)); 
        float angle = atan2(dot(r2, y), dot(r2, x)); //azimuth angle of r2
        float stepAngle = angle / float(divisions);

        int newPts[] = {};
        for (int j = 0; j <= divisions; ++j) {
            float jAngle = j * stepAngle; //Polar coord: set(jAngle, r, 0)
            vector j_on_plane = set(cos(jAngle) * radius, sin(jAngle) * radius, 0); //Cartesian coord 2D
            //maps 2D coord on the plane into world 3D coord
            matrix M = ident();
            M = set(
                set(x.x, x.y, x.z),
                set(y.x, y.y, y.z),
                set(z.x, z.y, z.z)
            );
            vector j_in_3D = j_on_plane * M + C;

            int newPt = addpoint(0, j_in_3D);
            push(newPts, newPt);
        }
        int insert = find(PrimPts, i);
        int insert_1 = find(PrimPts_new, i);

        if(PrimPts[(insert+1) % len(PrimPts)] == p1) //match new pts' order
            newPts = reverse(newPts);
        removeindex(PrimPts_new, insert_1);
        insert(PrimPts_new, insert_1, newPts);
    }
    removeprim(0, @primnum, 0);
    foreach(int i; primCorners) removepoint(0, i);
    addprim(0, "poly", PrimPts_new);
}