;;;-*-lisp-*-

(eval-when (eval compile)
  (or (get 'graphm 'version)
      (load (list (car (namelist infile)) 'graphm))))

(graphs-module plot3)

(or (get 'plot 'version)
    (load '((graphs)plot)))

(or (get 'graph-3d 'version)
    (load '((graphs)graph3)))

(defvar 3d-perspective-stream
	(make-z-perspective-stream graphic-stream))
(CALL 3d-perspective-stream
      'set-z-screen 5.0)
(CALL 3d-perspective-stream
      'set-z-eye 10.0)

(defvar 3d-clip-stream
	(make-z-clip-stream 3d-perspective-stream))
(CALL 3d-clip-stream 'set-clip-z 10.0)
(defvar 3d-stream (make-orthogonal-3d-stream 3d-perspective-stream))
(CALL 3d-stream 'set-theta 4.188)
(CALL 3d-stream 'set-phi 0.1)
(CALL 3d-stream 'set-psi 4.188)

;;;X = (A-RHO*SIN(PHI/2))*COS(PHI);
;;;Y = (A-RHO*SIN(PHI/2))*SIN(PHI);
;;;Z = RHO*COS(PHI/2);
;;;
;;;" abs(rho)<H, minf<phi<inf, 0<h<a<inf "$

(defvar mobius-a 0.5)
(defvar mobius-h 0.2)

(declare (flonum (mobius-x flonum flonum)
		 (mobius-y flonum flonum)
		 (mobius-z flonum flonum)))

(defun mobius-x (rho phi)
       (*$ (cos phi) (-$ mobius-a (*$ rho (sin (*$ 0.5 phi))))))
(defun mobius-y (rho phi)
       (*$ (sin phi) (-$ mobius-a (*$ rho (sin (*$ 0.5 phi))))))
(defun mobius-z (rho phi)
       (*$ rho (cos (*$ phi 0.5))))

(defun mobius (n m)
       (do ((dphi (//$ (*$ 2.0 3.1416) (float n)))
	    (phi 0.0 (+$ phi dphi))
	    (rho (-$ mobius-h) (-$ mobius-h))
	    (drho (//$ (*$ 2.0 mobius-h) (float m)))
	    (j 0 (1+ j)))
	   ((> j n))
	   (CALL 3d-stream
		 'set-pen
		 (mobius-x rho phi)
		 (mobius-y rho phi)
		 (mobius-z rho phi))
	   (do ((k 0 (1+ k)))
	       ((> k m))
	       (setq rho (+$ rho drho))
	       (CALL 3d-stream
		 'move-pen
		 (mobius-x rho phi)
		 (mobius-y rho phi)
		 (mobius-z rho phi)))))

;;; toroidal spiral.

(defvar r1 1.0)
(defvar r2 0.3)


(defun torus-point (com theta phi)
       (let ((r2-projection (+$ r1 (*$ r2 (cos phi)))))
	    (CALL 3d-stream
		  com
		  (*$ r2-projection (cos theta))
		  (*$ r2-projection (sin theta))
		  (*$ r2 (sin phi)))))

(defun torus (n m
		&aux
		(dtheta (//$ (*$ 2.0 3.1416) (float n)))
		(dphi (//$ (*$ 2.0 3.1416) (float m))))
       ;; this is a slow torus because of all the sin and cos calculation,
       ;; which are not really needed.
       (do ((j 0 (1+ j))
	    (theta 0.0 (+$ theta dtheta)))
	   ((> j n))
	   (torus-point 'set-pen theta 0.0)
	   (do ((k 0 (1+ k))
		(phi 0.0 (+$ phi dphi)))
	       ((> k m))
	       (torus-point 'move-pen theta phi))))

(defmacro dotrig ((j n sin cos) . body)
	  `(let* ((n-temp ,n)
		  (dt (//$ #.(*$ 2.0 3.14159265) (float n-temp))))
		 (do ((,j 0 (1+ ,j))
		      (,cos 1.0 (-$ (*$ cdt ,cos) (*$ sdt ,sin)))
		      (,sin 0.0 (+$ (*$ cdt ,sin) (*$ sdt ,cos)))
		      (cdt (cos dt))
		      (sdt (sin dt)))
		     ((> ,j n-temp))
		     ,@body)))

(defun ftorus-point (com sin-theta cos-theta sin-phi cos-phi)
       (let ((r2-projection (+$ r1 (*$ r2 cos-phi))))
	    (CALL 3d-stream
		  com
		  (*$ r2-projection cos-theta)
		  (*$ r2-projection sin-theta)
		  (*$ r2 sin-phi))))

;; this function is shorter and much faster, the fence-post test
;; for SET-PEN is insignificant.

(defun ftorus (n m)
       (dotrig (j n sin-theta cos-theta)
	       (dotrig (k m sin-phi cos-phi)
		       (ftorus-point (if (zerop j) 'set-pen 'move-pen)
				     sin-theta cos-theta
				     sin-phi cos-phi))))

(comment
(defun fat-line (stream x1 y2 z1 x2 y2 z2 &optional (w% 0.05)
			(n 3))
       ; first get basis vectors for the planes perpendicular to the
       ; line.
))

       

(defun draw-segment-list (stream list &aux p)
       (cond ((null list))
	     (t
	      (setq p (pop list))
	      (CALL stream 'set-pen (car p) (cadr p) (caddr p))
	      (do ()
		  ((null list) t)
		  (setq p (pop list))
		  (CALL stream 'move-pen (car p) (cadr p) (caddr p))))))

(defun cube ()
       (draw-segment-list 3d-stream
			  '((0.0 0.0 0.0)
			    (0.0 0.0 1.0)
			    (0.0 1.0 1.0)
			    (0.0 1.0 0.0)
			    (0.0 0.0 0.0)
			    (1.0 0.0 0.0)
			    (1.0 1.0 0.0)
			    (1.0 1.0 1.0)
			    (1.0 0.0 1.0)
			    (0.0 0.0 1.0)))
       (draw-segment-list 3d-stream
			  '((1.0 1.0 0.0)(0.0 1.0 0.0)))
       (draw-segment-list 3d-stream
			  '((1.0 0.0 0.0) (1.0 0.0 1.0)))
       (draw-segment-list 3d-stream
			  '((0.0 1.0 1.0) (1.0 1.0 1.0))))


(defvar 2pie (*$ 2.0 3.14159265))

(defun rotate-theta (f &optional (n 10.))
  (CALL 3d-stream 'set-dtheta (//$ 2pie (float n)))
  (do ()
      ((zerop (Setq n (1- n))))
    (CALL 3d-stream 'translate-theta)
    (funcall f)))