/*****************************************************************************
*                                                                            *
*   CIRCFIT.SPL  Copyright (C) 2020 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Least squares fit to a circle                              *
*                                                                            *
*   Revisions:    19 Jun 2020  RRR  Creation                                 *
*                                                                            *
*****************************************************************************/

#if @HELP_CIRCFIT

    CIRCFIT

    Purpose: Computes the least squares fit of a circle.

    Syntax:  CIRCFIT(x, y)

             (circ, coef) = CIRCFIT(x, y)

              x - A series, the X values.

              y - A series, the Y values.

    Alternate Syntax:  
    
             CIRCFIT(xy)

             (circ, coef) = CIRCFIT(xy)

              xy - An XY series, the input circle.

    Returns: An XY series, the fitted circle.

             (circ, coef) = CIRCFIT(x, y) returns the fitted circle and
             the coefficients as two series. The coefficients are of the
             form:

             coef == {xc, yc, r}

             where XC, and YC are the center coordinates and R is the radius.

    Example:
             W1: linspace(0, 2*pi, 100)
             W2: 5 * cos(W1) + grand(length(w1), 1)
             W3: 5 * sin(W1) + grand(length(w1), 1)
             W4: xy(w2, w3);setaspect(1)
             W5: circfit(w2, w3);setaspect(1);overp(w4);setx(-10, 10);sety(-10, 10)

             W4 contains a noisy circle. W5 computes the least squares fit and
             overplots the original circle.
 
    Example:
             W1: linspace(0, 2*pi, 100)
             W2: 5 * cos(W1) + grand(length(w1), 1)
             W3: 5 * sin(W1) + grand(length(w1), 1)
             W4: xy(w2, w3);setaspect(1)
             W5: circfit(w4);setaspect(1);overp(w4);setx(-10, 10);sety(-10, 10)

             Same as above except the input to CIRCFIT is an XY series.

    Example:
             W1: linspace(0, 2*pi, 100);
             W2: 3 * cos(w1) + 5
             W3: 3 * sin(w1) + 7
             W4: xy(w2, w3);setaspect(1)
             W5: w2[1..3]
             W6: w3[1..3]
             W7: (data, coef) = circfit(w5, w6);coef;table
             W8: xy(w7[3] * cos(w1) + w7[1], w7[3] * sin(w1) + w7[2])

             W4 contains a circle of radius 3 with a center at {5, 7}.

             W7 performs a three point fit of the circle where the resulting
             coefficients COEF == {xc, yc, r} == {5, 7, 3}.

             W8 computes and displays the resulting fitted circle.

    Remarks:
             CIRCFIT uses the \^ matrix solve operator to perform the least
             squares fit.

    See Also:
             \^
             Linfit
             Sinfit
             Sinfit3
             Sinfit4
             XY
#endif


/* least squares fit to a circle */
circfit(x, y)
{
        local a, xc, yc, r, coef, th, circ, xe, ye;
        
        if (argc < 2)
        {
                if (argc < 1)
                {
                        x = {};
                }

                y = yvals(x);
                x = xvals(x);
        }

        if (isempty(x))
        {
                error(sprintf("%s = XY or X and Y input series required", __FUNC__));
        }

        /* least squares fit */
        a = ravel(x, y, ones(size(x))) \^ (-(x^2 + y^2));

        /* circle center */
        xc = -0.5 * a[1];
        yc = -0.5 * a[2];

        /* circle radius */
        r = sqrt((a[1]^2 + a[2]^2) / 4 - a[3]);

        /* group center and radius */
        coef = {xc, yc, r};

        /* compute resulting circle */
        th = linspace(0, 2*pi, length(x));
        xe = r * cos(th) + xc;
        ye = r * sin(th) + yc;

        circ = xy(xe, ye);

        return(circ, coef);
}