! For a given set of 3 or more points in 2 dimensions,
measure the circularity of the points by 
finding two concentric circles,
 a) a smallest enclosing circle, and
 b) a largest enclosed circle.
The difference in squared radius of the two circles is a measure of circularity.
If the points are precisely on a circle,
then the difference in radius will be zero;

! Keywords: Chart, Chebychev, Circularity, Gass, Harary, LINGO, MinMax, 
  Quality control, Roundness, Scatter plot, Sphere, Witzgall;

! Ref: Gass, S., C. Witzgall, and H.H. Harary (1998), 
"Fitting Circles and Spheres to Coordinate Measuring Machine Data,"
 International J. of Flexible Manufacturing Systems, v. 10, pp 5-25.;
SETS:
  POINT: X, Y, R;  ! The data points;
  POINTG: XG, YG;  ! Some reference points to graph;
ENDSETS
DATA:
!Case01; 
 X, Y =
 0  0
 4  1
 2  3
 0  1
;
!Case02 
  X,   Y =
 100  100
 104  101
 102  103
 100  101
;
!Case03 
   X,   Y =
 -100  100
  -96  101
  -98  103
 -100  101
;
!Case04
     X,       Y =
  0.8660254   0.5
    0         1
 -0.8660254   0.5
 -0.8660254  -0.5
    0        -1
  0.8660254   -0.5
;
!Case05
     X,       Y =
  7.7949     4.5
    0        8.996
 -7.7953     4.4987
 -7.7963    -4.4998
    0       -9.0031
  7.7939    -4.5001
;
!Case06
    X,       Y =
    2         9
    6         7
    7         6
;

ENDDATA
! Given a set of points ( X(k), Y( k), we want to find a reference circle,
 with center ( X0, Y0) and  radius R0 such that the given points are close
 to being on this reference circle.
 The difference between 
     the squared radius of point k from (X0, Y0) and 
     the squared radius of the reference circle is:
    ( X(k) - X0)^2 + (Y( k) - Y0)^2 - R0^2
 =   X( k)^2 + Y( k)^2 - 2* X( k)* X0 - 2* Y( k)* Y0  + X0^2 + Y0^2 - R0^2
 =   X( k)^2 + Y( k)^2 - 2* X( k)* X0 - 2* Y( k)* Y0  - RHO,
 where
   RHO = R0^2 - X0^2 - Y0^2;
! We want to find the X0, Y0, and R0 that minimize the maximum of the difference
between the squared nominal radius, the squared radius of enclosing circle,
and squared radius of the enclosed circle;

SUBMODEL FitCirc:
MIN = T;
@FOR( POINT( k):
  R( k) = X( k)^2 + Y( k)^2 - 2* X( k)* X0 - 2* Y( k)* Y0 - RHO;
 [OUTR] T >=   R( k); ! At least 1 of these must be binding;
 [INNR] T >= - R( k); ! At least 1 of these must be binding, but not same k;
  @FREE( R( k)); ! Unconstrained in sign;
    );
  @FREE( X0);
  @FREE( Y0);
  @FREE( RHO);
! Where we should have additionally the nonlinear constraint;
! [NONL]   RHO = R0^2 - X0^2 - Y0^2;
! We want to show that if the nonlinear constraint [NONL] is dropped,
  then at an optimal solution there is a feasible solution to [NONL],
  so dropping it is justified.
The argument is as follows.
 At an optimum, at least one of the [INNR] constraints must hold as an =, so
   T = - X( k)^2 -Y( k)^2 + 2* X( k)* X0 + 2* Y( k)* Y0 + RHO
 Now T >= 0, so
   0 <= - X( k)^2 -Y( k)^2 + 2* X( k)* X0 + 2* Y( k)* Y0 + RHO
 or  
   X( k)^2  +  Y( k)^2 - 2* X( k)* X0 - 2* Y( k)* Y0 <= RHO
 or  
   X( k)^2  +  Y( k)^2 - 2* X( k)* X0 - 2* Y( k)* Y0  + X0^2 + Y0^2 <= RHO + X0^2 + Y0^2, 
 or  
   ( X( k) - X0)^2  + ( Y( k) - Y0) ^2 <= RHO + X0^2 + Y0^2,
 or 
     0 <= RHO + X0^2 + Y0^2,
 thus we can solve:
   R0^2 = RHO + X0^2 + Y0^2, 
 Because the RHS >= 0;
ENDSUBMODEL

CALC:
  @SOLVE( FitCIrc);
  
  R0 = ( RHO + X0^2 + Y0^2)^0.5; ! Get fitted radius;
 ! Display a graph;
 ! Generate a set of points on a circle;
 NPTG = 64;
 POINTG = 1..NPTG;
 
! Generate some points on the reference circle;
 @FOR( POINTG( j):
    XG( j) = X0 + R0* @COS( 2*@PI()*(j-1)/NPTG);
    YG( j) = Y0 + R0* @SIN( 2*@PI()*(j-1)/NPTG);
     );
! Generate a scatter plot of the two sets of points;
 @CHARTSCATTER( 'Fit a Circle to a Set of Points, Tau= '+@FORMAT( T,"5.2f"), 'x-axis', 'y-axis', 
     'Data', X, Y,
     'Fitted', XG, YG);
ENDCALC