;; The cube and tetrahedron routines have been ported from the Python
;; GTS library; The ASCII figures appearing below have been taken from
;; the Python GTS library source code.

;;
;; Returns a cube of side length 2 centered at the origin.
;;    
;;              
;;        v8 +------+ v5
;;          /      /|
;;         /    v1/ |
;;     v4 +------+  |
;;        |      |  + v6
;;        |(v7)  | /
;;        |      |/
;;     v3 +------+ v2
;; 
 
(define (cube)

  (let (;; the vertices
	(v1 (vertex-new  1.0  1.0  1.0))
	(v2 (vertex-new  1.0  -1.0 1.0))
	(v3 (vertex-new  -1.0 -1.0 1.0))
	(v4 (vertex-new  -1.0 1.0  1.0))
	(v5 (vertex-new  1.0  1.0 -1.0))
	(v6 (vertex-new  1.0  -1.0 -1.0))
	(v7 (vertex-new -1.0  -1.0 -1.0))
	(v8 (vertex-new -1.0   1.0 -1.0))
	)
    (let (
	  ;; the edges
	  (e12 (edge-new v1 v2))
	  (e23 (edge-new v2 v3))
	  (e34 (edge-new v3 v4))
	  (e14 (edge-new v1 v4))

	  (e56 (edge-new v5 v6))
	  (e67 (edge-new v6 v7))
	  (e78 (edge-new v7 v8))
	  (e58 (edge-new v5 v8))
	  
	  (e15 (edge-new v1 v5))
	  (e26 (edge-new v2 v6))
	  (e37 (edge-new v3 v7))
	  (e48 (edge-new v4 v8))
	  
	  (e13 (edge-new v1 v3))
	  (e16 (edge-new v1 v6))
	  (e18 (edge-new v1 v8))
	  (e27 (edge-new v2 v7))
	  (e47 (edge-new v4 v7))
	  (e57 (edge-new v5 v7))
	  
	  )

      (let ((faces     
	     (list (face-new e12 e23 e13) 
		   (face-new e13 e34 e14)
		   (face-new e12 e26 e16) 
		   (face-new e15 e56 e16)
		   (face-new e15 e58 e18) 
		   (face-new e14 e48 e18)
		   (face-new e58 e78 e57) 
		   (face-new e56 e67 e57)
		   (face-new e26 e67 e27) 
		   (face-new e37 e23 e27)
		   (face-new e37 e47 e34) 
		   (face-new e78 e48 e47) ))
	    (s (surface-new))
	    )

	(for-each (lambda (f) 
		    (if (not (face-is-compatible f s)) (triangle-revert f))
		    (surface-add-face s f))
		  faces)
	
	s
	))
    ))

;;
;; Returns a tetrahedron of side length 2*sqrt(2) centered at origin.
;; The edges of the tetrahedron are perpendicular to the cardinal
;; directions.
;;
;;        v4
;;         +
;;         | \ e6
;;   e5   '|e4 \ 
;;    v1 . +-e3-+ v3
;;        /   . 
;;      ./e1. e2
;;      / .
;;     +
;;    v2

(define (tetrahedron)

  (let (;; the vertices
	(v1 (vertex-new  1.0  1.0  1.0))
	(v2 (vertex-new -1.0 -1.0  1.0))
	(v3 (vertex-new -1.0  1.0 -1.0))
	(v4 (vertex-new  1.0 -1.0 -1.0))
	)

    ;; the edges
    (let ((e1 (edge-new v1 v2))
	  (e2 (edge-new v2 v3))
	  (e3 (edge-new v3 v1))
	  (e4 (edge-new v1 v4))
	  (e5 (edge-new v4 v2))
	  (e6 (edge-new v4 v3)))

      (let ((faces     
	     (list (face-new e1 e2 e3) 
		   (face-new e1 e4 e5) 
		   (face-new e2 e5 e6) 
		   (face-new e3 e4 e6)))
	    (s (surface-new))
	    )

	(for-each (lambda (f) 
		    (if (not (face-is-compatible f s)) (triangle-revert f))
		    (surface-add-face s f))
		  faces)
	s
	))
    ))

(define (sphere geodesation-order)
  (define surface_face_revert (foreign-value "&gts_triangle_revert" c-pointer))
  (let ((s (surface-new)))
    (surface-generate-sphere s (or (and (positive? geodesation-order) geodesation-order) 4))
    (surface-foreach-face s surface_face_revert s)
    s))



(define (isosurface f iso nx ny nz #!key (dx (/ 20 (- nx 1))) 
		                         (dy (/ 20 (- ny 1))) 
     		                         (dz (/ 20 (- nz 1)))
					 (method 'cartesian))
  (define surface_face_revert (foreign-value "&gts_triangle_revert" c-pointer))
  (let ((s (surface-new))
	(x -10.0) (y -10.0)  (z -10.0))
    (case method
	   ((cartesian)          (let ((s (isosurface-cartesian s iso nx ny nz dx dy dz x y z f)))
				   (surface-foreach-face s surface_face_revert s)
				   s))
	   ((tetra)              (isosurface-tetra s iso nx ny nz dx dy dz x y z f))
	   ((tetra-bounded)      (let ((s (isosurface-tetra-bounded s iso nx ny nz dx dy dz x y z f)))
				   (surface-foreach-face s surface_face_revert s)
				   s))
	   ((tetra-bcl)          (isosurface-tetra-bcl s iso nx ny nz dx dy dz x y z f)))
    ))


(define pi 3.1415926535)

(define inv2pi (/ 1.0 (* 2.0 pi)))

(define (atan-turns x z) ;; result between 0 and 1
  (let ((r (fp* (atan x z) inv2pi)))
    (if (fp< r 0.0) (fp+ 1.0 r) r)))


(define (cylinder-function b e r)
  (lambda (x y z) 
    (let ((d  (fp+ (fp* x x) (fp* y y))))
      (cond ((fp= z b) (abs z))
	    ((fp< (abs (fp- z e)) 1e-10) (abs z))
	    (else (fp- d (fp* r r)))))))

  
(define (cylinder radius iso)  
  (let ((fp (object-evict (cylinder-function -10.0 10.0 radius))))
    (let ((s (isosurface fp iso (* 10 iso) (* 10 iso) (* 10 iso) method: 'cartesian)))
    (object-release fp)
    s)))


(define (cone-function b e r)
  (let* ((h (abs (- e b)))
	 (c (/ (* h h) (* r r) )))
    (lambda (x y z) 
      (let ((d (fp+ (fp* x x) (fp* y y))))
	(cond ((fp= z b) (abs z))
	      ((fp< (abs (fp- z e)) 1e-10) (abs z))
	      (else (let ((dz (fp- z h)))
		      (fp- (fp* d c) (fp* dz dz))))
	      )))))
  
(define (cone radius iso)  
  (let ((fp (object-evict (cone-function -10.0 10.0 radius))))
    (let ((s (isosurface fp iso (* 10 iso) (* 10 iso) (* 10 iso) method: 'cartesian)))
    (object-release fp)
    s)))