#| Doug Hoyte, Jan 2007
 |
 | This program takes shape data specified as a lisp
 | "vertex spec" and compiles C code for use in a
 | plotting program.
 |
 | Usage:
 |
 | (compile-shape icosahedron)
 | (compile-shape (prism N))
 | (compile-shape (antiprism N))
 | (compile-shape (pyramid N))
 | (compile-shape (bipyramid N))
 | (compile-shape icosahedron dodecahedron) ; Combine polyhedra into 1 plot
 |
 | It can also determine the geometric properties of many polyhedra:
 |
 | (vspec-info cube)
 |   Vertices = 8  Edges = 12  Faces = 6
 | (vspec-info icosahedron)
 |   Vertices = 12  Edges = 30  Faces = 20
 | (vspec-info dodecahedron)
 |   Vertices = 20  Edges = 30  Faces = 12
 | (vspec-info (prism 5))
 |   Vertices = 10  Edges = 15  Faces = 7
 | (vspec-info (antiprism 5))
 |   Vertices = 10  Edges = 20  Faces = 12
 | (vspec-info dodecahedron)
 |   Vertices = 20  Edges = 30  Faces = 12
 | (vspec-info (pyramid 9))
 |   Vertices = 10  Edges = 18  Faces = 10
 | (vspec-info (bipyramid 123))
 |   Vertices = 125  Edges = 369  Faces = 246
 |#


;; Hardcoded vertex specs:

(defconstant cube
  '(((0 1) (0 1) (0 1))))

(defconstant icosahedron
  '(((1 -1) (phi (- phi)) (0))
    ((0) (1 -1) (phi (- phi)))
    ((phi (- phi)) (0) (1 -1))))

(defconstant dodecahedron
  '(((1 -1) (1 -1) (1 -1))
    ((0) (phi^-1 (- phi^-1)) (phi (- phi)))
    ((phi^-1 (- phi^-1)) (phi (- phi)) (0))
    ((phi (- phi)) (0) (phi^-1 (- phi^-1)))))

(defconstant cuboctahedron
  '(((1 -1) (1 -1) (0))
    ((1 -1) (0) (1 -1))
    ((0) (1 -1) (1 -1))))

(defconstant truncated-hexahedron
  '(((gam (- gam)) (1 -1) (1 -1))
    ((1 -1) (gam (- gam)) (1 -1))
    ((1 -1) (1 -1) (gam (- gam)))))

(defconstant truncated-octahedron
  '(((0) (1 -1) (2 -2))
    ((0) (2 -2) (1 -1))
    ((1 -1) (0) (2 -2))
    ((1 -1) (2 -2) (0))
    ((2 -2) (0) (1 -1))
    ((2 -2) (1 -1) (0))))

(defconstant gam (- (sqrt 2) 1))
(defconstant phi (/ (+ (sqrt 5) 1) 2))
(defconstant phi^-1 (/ 1 phi))



;; Dynamic vertex specs:

(defun prism (n)
  (let ((dist (cart-dist (list (cos (/ (* 2 pi) n)) (sin (/ (* 2 pi) n)))
			 (list (cos (* 2 (/ (* 2 pi) n))) (sin (* 2 (/ (* 2 pi) n)))))))
    (loop for i from 1 to n collect
      (list
        (list (cos (* i (/ (* 2 pi) n))))
        (list (sin (* i (/ (* 2 pi) n))))
        (list 0 dist)))))

(defun antiprism (n)
  (loop for k from 0 to (- (* 2 n) 1) collect
    (list (list (cos (/ (* k pi) n)))
	  (list (sin (/ (* k pi) n)))
	  (list (* (expt -1 k) (sqrt (/ (- (cos (/ pi n)) (cos (/ (* 2 pi) n))) 2)))))))


(defun pyramid (n)
  (let ((dist (cart-dist (list (cos (/ (* 2 pi) n)) (sin (/ (* 2 pi) n)))
			 (list (cos (* 2 (/ (* 2 pi) n))) (sin (* 2 (/ (* 2 pi) n)))))))
    (cons
      (list '(0) '(0) (list (sqrt (abs (- (* dist dist) 1)))))
      (loop for i from 1 to n collect
        (list
          (list (cos (* i (/ (* 2 pi) n))))
          (list (sin (* i (/ (* 2 pi) n))))
          (list 0))))))

;; Relies on the first vertex in a pyramid being the peak (as returned in #'pyramid above)
(defun bipyramid (n)
  (let ((p (pyramid n)))
    (cons (list '(0) '(0) (list (- (caaddr (car p))))) p)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;; n-dimensional cartesian distance
(defun cart-dist (v1 v2)
  (sqrt (apply #'+ (mapcar (lambda (a b) (expt (- a b) 2)) v1 v2))))

(defun compile-shape (&rest vspecs)
  (let ((verts nil)
	(edges nil))
    (loop for vspec in vspecs do
      (setf verts (append verts (expand-vspec vspec)))
      (setf edges (append edges (calc-edges (expand-vspec vspec)))))
    (compile-to-c verts edges)))

(defun vspec-info (vspec)
  (let* ((vertlist (expand-vspec vspec))
	 (verts (length vertlist))
	 (edges (length (calc-edges vertlist)))
	 (faces (- (+ 2 edges) verts))) ; All our shapes have euler characteristic = 2
    (format t "Vertices = ~a  Edges = ~a  Faces = ~a" verts edges faces)))
  

(defun expand-vspec (vspec)
  (let ((o nil))
    (loop for v in vspec do
      (dolist (x (car v))
        (dolist (y (cadr v))
          (dolist (z (caddr v))
	    (setf o (cons (list (eval x) (eval y) (eval z)) o))))))
    o))


(defconstant *tolerance* 0.01)

(defun calc-edges (vs)
  (let* ((o nil)
	 (numvs (length vs))
	 (vhash (make-hash-table :test 'equal))
         (a (make-array numvs :initial-contents vs)))
    (loop for i from 0 to (- numvs 1) do
      (let ((min 1000000)
	    (cur nil))
        (loop for j from 0 to (- numvs 1) unless (= i j) do
          (let ((dist (cart-dist (aref a i) (aref a j))))
	    (if (< (+ dist *tolerance*) min)
	      (progn
		(setf min dist)
		(setf cur (list (list i j))))
	      (if (< (abs (- min dist)) *tolerance*)
		(setf cur (cons (list i j) cur))))))

	(loop for (i j) in cur do
	  (if (not (gethash (sort (list i j) #'<) vhash))
	    (progn
	      (setf o (cons (list i j) o))
	      (setf (gethash (sort (list i j) #'<) vhash) t))))))
    o))


(defun compile-to-c (vs es)
  (format t "float vX[] = { ~{~,4f~^, ~} };~%" (mapcar #'car vs))
  (format t "float vY[] = { ~{~,4f~^, ~} };~%" (mapcar #'cadr vs))
  (format t "float vZ[] = { ~{~,4f~^, ~} };~%~%" (mapcar #'caddr vs))
  (format t "int eS[] = { ~{~a~^, ~} };~%" (mapcar #'car es))
  (format t "int eE[] = { ~{~a~^, ~} };~%" (mapcar #'cadr es)))