Fit a circle to 2D coords (least squares method)

Function FitCircleTo2DCoords(w)
    Wave w
    // requires 2D numeric wave with two columnns corresponding to xy coords
    if(Dimsize(w,1) != 2)
        return -1
    endif
    // make two 1D waves for x and y coords
    SplitWave/O/NAME="xW;yW;" w
    Wave xW,yW
    // solve in terms of u and v coordinates
    Duplicate/O xw, uW
    Duplicate/O yw, vW
    uW[] = xw[p] - mean(xw)
    vW[] = yw[p] - mean(yw)
    // sigma calcs - store as variables for clean code
    Variable Su = sum(uW)
    Variable Sv = sum(vW)
    MatrixOp/O/FREE uu = uW * uW
    Variable Suu = sum(uu)
    MatrixOp/O/FREE vv = vW * vW
    Variable Svv = sum(vv)
    MatrixOp/O/FREE uv = uW * vW
    Variable Suv = sum(uv)
    MatrixOp/O/FREE uuu = uW * uW * uW
    Variable Suuu = sum(uuu)
    MatrixOp/O/FREE uvv = uW * vW * vW
    Variable Suvv = sum(uvv)
    MatrixOp/O/FREE vvv = vW * vW * vW
    Variable Svvv = sum(vvv)
    MatrixOp/O/FREE vuu = vW * uW * uW
    Variable Svuu = sum(vuu)
    // linear system
    Make/O/N=(2,2) matA = {{Suu,Suv},{Suv,Svv}}
    Make/O/N=(2,1) matB = {0.5 * (Suuu + Suvv),0.5 * (Svvv + Svuu)}
    // Solve it. matB is overwritten with the result uc,vc
    MatrixLinearSolve/O matA matB
    // transform back to x y coordinate system to get xc,yc
    Duplicate/O matB, matC
    // matC is saved and contains the centre of the circle
    matC[0] = matB[0] + mean(xW)
    matC[1] = matB[1] + mean(yW)
    // alpha is the radius squared
    Variable alpha = matB[0]^2 + matB[1]^2 + ( (Suu + Svv) / numpnts(xW) )
    // clean up 1D waves
    KillWaves/Z xW,yW
   
    // return radius
    return sqrt(alpha)
End