;;; -*- Mode:LISP; Package:VIF; Base:10; Readtable:CL -*-

(defstruct (line :named)
  slope y-intercept)

(defstruct (singular-line :named (:constructor make-singular-line (x-intercept)))
  x-intercept)

(defun y-from-x-and-line (x line)
  (if (typep line 'singular-line)
      nil
    (+ (* x (line-slope line)) (line-y-intercept line))))

(defun line-from-point-and-slope (x1 y1 slope)
  (make-line :slope slope :y-intercept (- y1 (* slope x1))))

(defun line-from-two-points (x1 y1 x2 y2)
  (if (= x1 x2)
      (make-singular-line x1)
    (let ((slope (/ (float (- y2 y1)) (float (- x2 x1)))))
      (line-from-point-and-slope x1 y1 slope))))

(defun line-through-point-perpendicular-to-line (x1 y1 line)
  (cond ((= (line-slope line) 0)		;horizontal
	 (make-singular-line x1))
	((typep line 'singular-line)
	 (make-line :slope 0 :y-intercept y1))
	(t (line-from-point-and-slope x1 y1 (- (/ (line-slope line)))))))

(defun line-parallel-p (line1 line2)
  (cond ((and (typep line1 'singular-line)
	      (typep line2 'singular-line)) t)
	((= (line-slope line1)
	    (line-slope line2)) t)
	(t nil)))

(defun lines-intersect-at (line1 line2)
  (cond ((line-parallel-p line1 line2)
	 nil)
	((typep line1 'singular-line)
	 (let* ((x (singular-line-x-intercept line1))
		(y (y-from-x-and-line x line2)))
	   (list x y)))
	((typep line2 'singular-line)
	 (let* ((x (singular-line-x-intercept line2))
		(y (y-from-x-and-line x line1)))
	   (list x y)))
	(t (let* ((slope1 (line-slope line1))
		  (slope2 (line-slope line2))
		  (intercept1 (line-y-intercept line1))
		  (intercept2 (line-y-intercept line2))
		  (x (/ (- intercept2 intercept1)
			(- slope1 slope2)))
		  (y (+ (* slope1 x) intercept1)))
	     (list x y)))))

(defun line-segment-length (x1 y1 x2 y2)
  (flet ((square (n)
	   (* n n)))
    (sqrt (float (+ (square (- x2 x1)) (square (- y2 y1)))))))

(defun line-segment-middle-point (x1 y1 x2 y2)
  (list (/ (float (+ x1 x2)) 2)
	(/ (float (+ y1 y2)) 2)))

(defun circle-center-and-radius-from-3-points (x1 y1 x2 y2 x3 y3)
  (flet ((circle-radius (x1 y1 x2 y2)
	   (let ((midpoint (line-segment-middle-point x1 y1 x2 y2)))
	     (line-through-point-perpendicular-to-line (first midpoint) (second midpoint)
						       (line-from-two-points x1 y1 x2 y2)))))
    (let* ((radius1 (circle-radius x1 y1 x2 y2))
	   (radius2 (circle-radius x2 y2 x3 y3))
	   (circle-center (lines-intersect-at radius1 radius2))
	   (circle-radius (line-segment-length x1 y1
					       (first circle-center) (second circle-center))))
      (values (first circle-center)(second circle-center) circle-radius))))