; $Id: cir_3pnt.pro,v 1.13 2004/01/21 15:54:49 scottm Exp $
;
; Copyright (c) 1992-2004, Research Systems, Inc.  All rights reserved.
;   Unauthorized reproduction prohibited.
;

PRO cir_3pnt, xIn, yIn, r, x0, y0

;+
; NAME:
;   CIR_3PNT
;
; PURPOSE:
;   This procedure returns the radius and center of a circle,
;   given 3 points on the circle. This is analogous to finding
;   the circumradius and circumcircle of a triangle; the center
;   of the circumcircle is the point at which the three perpendicular
;   bisectors of the triangle formed by the points meet.
;
; CATEGORY:
;   Analytical geometry.
;
; CALLING SEQUENCE:
;   CIR_3PNT, X, Y, R, X0, Y0
;
; INPUTS:
;   X: A three-element vector containing the X-coordinates of the points.
;   Y: A three-element vector containing the Y-coordinates of the points.
;
; OUTPUTS:
;   R: The radius of the circle. The procedure returns 0.0 if the points
;      are co-linear.
;   X0, Y0: The coordinates of the center of the circle. The procedure
;           returns 0.0 if the points are co-linear.
;
; PROCEDURE:
;   Derived from Glasser, ed.  Graphics Gems, Volume 1, Page 22.
;
; EXAMPLE:
;   X = [1.0, 2.0, 3.0]
;   Y = [1.0, 2.0, 1.0]
;   CIR_3PNT, X, Y, R, X0, Y0
;   Print, 'The radius is: ', R
;   Print, 'The center of the circle is at: ', X0, Y0

; MODIFICATION HISTORY:
;   Written by: DMS, RSI, Nov, 1992.
;                       SJL, Nov, 1997.  - Formatting.
;   Modified: CT, RSI, August 2002. Corrected abs(c) test.
;             Also don't stomp on input vars.
;-

    x = double(xIn)
    y = double(yIn)

    d1 = (x[2]-x[0]) * (x[1]-x[0]) + (y[2]-y[0]) * (y[1]-y[0])
    d2 = (x[2]-x[1]) * (x[0]-x[1]) + (y[2]-y[1]) * (y[0]-y[1])
    d3 = (x[0]-x[2]) * (x[1]-x[2]) + (y[0]-y[2]) * (y[1]-y[2])

    c1 = d2 * d3
    c2 = d3 * d1
    c3 = d1 * d2
    c = c1 + c2 + c3

    if (abs(c) lt 1e-14) then begin ;Colinear?
        r =  0d
        x0 = 0d
        y0 = 0d
        return
    endif

    r = sqrt((d1 + d2)*(d2 + d3) * (d3 + d1)/c)/2.
    v = [ c2 + c3, c3 + c1, c1 + c2 ]/(2.*c)
    x0 = total(v * x)
    y0 = total(v * y)

end