;;;; Simple 3D Viewer for state space trajectories and wire
;;;; frames

;;;; Wire frame representation

;;; The program displays a object represented as a list of
;;; lines.  A line is simply a cons of two vectors, each of
;;; length three.  The car and cdr are the start and end
;;; points of the line in three dimensional space. The
;;; program projects these lines into two dimensional space
;;; under control of the mouse moving in the angle window

;;;; Organisation of X11 code

;;; CLX code needs alot of set up. The data structures that
;;; get built are genuinely global so I use the property
;;; lists of symbols to hold them. 

;;; For example, this program has two windows, the view
;;; window that displays the object and the angle window, in
;;; which mouse motion controls the viewing angle.  For each
;;; window, the program must store an X11 window object and
;;; a graphics context. So I store these as the Window and
;;; Grackon properties of the View and Angle symbols.

;;; The function Setup separates out these preliminary operations
;;; from the code that runs interactively.

;;; I don't "use" the Xlib package, so my code is littered
;;; with package qualifiers. X11's client side refresh means
;;; that the drawing functions are all one level of
;;; abstraction down from where one wants to code. Once you
;;; have decided how to implement client side refresh you
;;; need to wrap the Xlib drawing routines in versions that
;;; make a note of what you are drawing ready for responding
;;; to an expose event. Even when, as here, the
;;; interactivity makes client side refresh a non-issue, Xlib is general
;;; enough for all applications, and any particular application will
;;; want to wrap the Xlib routines in ways that take advantage of
;;; its peculiar regularities. For example my draw-line is a thin wrapper
;;; on Xlib:draw-line

;;;; (gripe-p) => true

;; I hate getting caught out by the fact that CL:T is a constant
;; I want to use t for time
(defconstant true 'cl:t)
(shadow 't)
;; Ofcourse, now I get caught out when (format t ...) doesn't work.
;; It needs to be (format true ...) so that format sees CL:T
;; Format doesn't accept CL-USER:T 

;; CLX is mostly functions, so I can almost leave loading CLX until
;; load time. Almost. Xlib:event-case is a macro, so it had better 
;; be available at compile time
(eval-when (:compile-toplevel :load-toplevel :execute)
	   (require :clx))

;;;; MACRO repeat-interactively
;;;
;;; Allows the programmer to work on running code
;;; 
;;; Fakes multi-tasking with Sleep and Listen.
;;; 
;;; Your code executes in a dynamic environment with
;;; a 3 special variables 
;;; 
;;; PROMPT containing the format string for the prompt
;;; TICK containing the delay in the loop
;;; /your own choice/ a count of how many times round the loop
;;;
;;; and a catch-tag, STOP, for exiting

(defmacro repeat-interactively ((&optional (tick 1)
					   (count-var 'count)
					   (prompt "~&>>-~A-> "))
				&body code)
  `(flet ((prompt-user()
	    (declare (special prompt ,count-var ))
	    (format true prompt ,count-var)
	    (force-output)))
      (let ((,count-var 0)
	    (prompt ,prompt)
	    (tick ,tick))
	(declare (special ,count-var prompt tick))
	(prompt-user)
	(catch 'stop
	  (loop
	   (restart-case
	     (progn ,@code (incf ,count-var)(sleep tick))
	     (fixed () :report "Resume processing."))
	   (when (listen)
	     (tagbody top
	       (restart-case
		 (progn
		   (print (eval (read)))
		   (prompt-user))
		 (shrug () :report "Continue processing."
			(prompt-user))
		 (retry ()
			:report "Evaluate another form before continuing."
			(prompt-user)
			(go top))))))))))

;;;; MACRO define vector funtion

;;; A common idea in this kind of code is to have a function
;;; that maps a state space to itself, often to a related space
;;; of the same dimension. The function usually has some
;;; parameters that control the mapping. The awkwardness of
;;; coding the function arise from the inhomogeneity of the
;;; state space. It is natural think in terms of named
;;; components, not numbered components, which suggests
;;; using a structure.  It is natural elsewhere to treat the
;;; state space as R^n, which suggests using a vector.

;;; A typical example is the Lorenz model. Your text book will
;;; probably notate the equations thus
#|
.
x = p(y-x)
.
y = -xy+rx-y
.
z = xy-bz
|#
;;; We code it
;;; (defvecfun lorenz (p r b)
;;;   (x (* p (- y x)))
;;;   (y (- (* r x)
;;;         (* x y)
;;;         y)
;;;   (z (- (* x y)
;;;         (* b z))))

;;; generating a function called lorenz, which takes four
;;; arguments, a vector and the three parameters p,r, and b,
;;; and which returns a vector.  The macro uses the names of
;;; the components to set up the accessors so that the code works.

(defmacro defvecfun (name parameter-list &rest defining-equations)
  (let ((vector (make-symbol "V"))
        (variables (mapcar #'car defining-equations)))
    `(defun ,name (,vector ,@parameter-list)
       (let ,(mapcar
              (lambda (var)
                (list var
                      (list 'aref
                            vector
                            (position var
                                      variables))))
              variables)
         (vector ,@(mapcar #'second defining-equations))))))

(defvar *scale* 100)
(defvar *theta* 0 "elevation in radians")
(defvar *phi* 0 "azimuth in radians")
(defvar *rotating* nil)
(defvar *old-x* 180)
(defvar *old-y* 90)
(defvar *lines* '() "list of lines in the wire frame")



;;;; Entry point
(defun start (&optional (host ""))
  (unwind-protect
      (interactive-event-case host)
    (xlib:close-display (get 'display 'x))))

(defun interactive-event-case (host)
  (setup host)
  (map-windows '(view angle))
  (repeat-interactively (0) ; no delay, we use the time out in the event loop
			;; We call a function rather than including the actual code
			;; This lets us redefine the event-loop while the code is running
			(event-loop)))

(defun setup (host)
  (setf (get 'display 'x) (xlib:open-display host)
	(get 'screen 'x) (first (xlib:display-roots
				 (get 'display 'x)))
	(get 'root 'window) (xlib:screen-root (get 'screen 'x))
	(get 'black 'pixel) (xlib:screen-black-pixel (get 'screen 'x))
	(get 'white 'pixel) (xlib:screen-white-pixel (get 'screen 'x))
	(get 'angle 'window) (xlib:create-window
				 :parent (get 'root 'window)
				 :x 0 :y 0
				 :width 360
				 :height 180
				 :background (get 'white 'pixel)
				 :event-mask '(:exposure
					       :pointer-motion
					       :button-press
					       :enter-window))
	(get 'angle 'grackon) (xlib:create-gcontext
				  :drawable (get 'angle 'window)
				  :foreground (get 'black 'pixel)
				  :background (get 'white 'pixel)
				  :font "-*-lucida-medium-r-*-*-12-*-*-*-*-*-*")
	(get 'view 'window) (xlib:create-window
				:parent (get 'root 'window)
				:x 0 :y 0
				:width 400
				:height 300
				:event-mask '(:exposure)
				:background (get 'white 'pixel))
	(get 'view 'mid-width) 200
	(get 'view 'mid-height) 150
	(get 'view 'grackon) (xlib:create-gcontext
				 :drawable (get 'view 'window)
				 :foreground (get 'black 'pixel)
				 :background (get 'white 'pixel))))

(defun map-windows (window-list)
  ;; A flight of technique Since I only have two windows it
  ;; would have been simpler just to call xlib:map-window
  ;; twice.  On the other hand this code shows off the idiom
  ;; of using a symbol to name a window and storing the X11
  ;; window object as a symbol's window property. The code
  ;; is the comment and this is the meta-comment.
  (dolist (window window-list)
    (xlib:map-window (get window 'window))))


(defun event-loop ()
  (xlib:event-case ((get 'display 'x)
		    :force-output-p true
		    :discard-p true
		    :timeout 0.2)
		   (:exposure (window count)
			      (when (slot-exists-p window 'alan)
				(funcall (slot-value window 'alan)))
			      (when (zerop count)
				(refresh window))
			      nil)
		   (:motion-notify (x y)
				   ;; We know this is in the angle window
				   ;; because the view window doesn't accept
				   ;; pointer motion events
				   (process-angle-input x y)
				   nil)
		   (:button-press (x y)
				  ;; We know this is in the angle window
				  ;; because the view window doesn't accept
				  ;; button press events
				  ;; 
				  ;; The user wants to leave the angle window
				  ;; so stop rotating the object and remember
				  ;; the angles at which the user stopped
				  ;; twiddling
				  (setf *rotating* nil 
					*old-x* x
					*old-y* y)
				  ;; The user will need a cross on the screen
				  ;; to find his way back
				  (mark-old-position x y))))

(defun mark-old-position (x y)
  (let ((size 10))
    (draw-line 'angle
	       (- x size)
	       y
	       (+ x size)
	       y)
    (draw-line 'angle
	       x
	       (+ y size)
	       x
	       (- y size))))

(defun process-angle-input (x y)
  ;; The user is either rotating the object
  ;; or he is trying to pick up the cross with the mouse
  ;; or maybe the mouse is just traversing the angle window
  ;; as it journeys to a window beyond.
  (cond (*rotating*
	 (setf *theta* (/ (* 2 pi (- x 180))
			360)
	       *phi* (/ (* 2 pi (- y 90))
		      360))
	 (refresh-view-window))
	((and (= x *old-x*)
	      (= y *old-y*))
	 ;; if pixel precise is too fiddley on a big screen
	 ;; this is where to put in a tolerance
	 (xlib:clear-area (get 'angle 'window))
	 (setf *rotating* true))))

(defun refresh (window)
  "Dispatch according to the identity of the X11 window object."
  (cond ((eql window (get 'angle 'window)) (refresh-angle-window))
	((eql window (get 'view 'window)) (refresh-view-window))))

(defun refresh-angle-window ()
  (unless *rotating*
    (mark-old-position *old-x* *old-y*)))

(defun refresh-view-window ()
    (xlib:clear-area (get 'view 'window))
    (dolist (line *lines*)
      ;; The inner loop
      ;; Building up my argument list like this
      ;; has got to be a waste of cons cells
      (apply #'draw-line
	     'view
	     (append (project (car line))
		     (project (cdr line))))))


(defun draw-line (paper x1 y1 x2 y2 &optional relative-p)
  (xlib:draw-line (get paper 'window)
		  (get paper 'grackon)
		  x1 y1 x2 y2 relative-p))

(defun project (xyz)
  (let ((transformed-vector (elevation (azimuth xyz *theta*) *phi*)))
    (list (round (+ (get 'view 'mid-width)
		    (* *scale* (aref transformed-vector 1)))) ;to right
	  (round (- (get 'view 'mid-height)
		    (* *scale* (aref transformed-vector 2)))))))

;;; I start with a right handed co-ordinate system
;;; x towards me
;;; y to the right
;;; z going up

;;; I apply an azimuth transformation

(defvecfun azimuth (theta)
  (x (- (* x (cos theta))
	(* y (sin theta))))
  (y (+ (* y (cos theta))
	(* x (sin theta))))
  (z z))

;;; Now my co-ordinates are towards-me, to-right, and z

(defvecfun elevation (phi)
  (towards-me (+ (* towards-me (cos phi))
		 (* z (sin phi))))
  (to-right to-right)
  (z (- (* z (cos phi))
	(* towards-me (sin phi)))))

;;;; THE END

;;; This is the end of the 3D Viewer code, but it is no fun with
;;; nothing to look at

;;;; Make a cube
     
;;; I succeed in getting the computer to do all the work
;;; by writing a six deep nested loop to generate lines
;;; between every corner of the cube, and then filtering to
;;; pick out the edges, which are shorter than the diagonals.

(defun list-dist (u v)
  (if u (+ (expt (- (car u)
		    (car v)) 2)
	   (list-dist (cdr u) (cdr v)))
    0))

(defun make-cube ()
  (let ((lines '())
	(k (list -1 1)))
    (dolist (a k)
      (dolist (b k)
	(dolist (c k)
	  (dolist (x k)
	    (dolist (y k)
	      (dolist (z k)
		(when (= 4 (list-dist (list a b c)
				      (list x y z)))
		  (push (cons (vector a b c)
			      (vector x y z))
			lines))))))))
    (remove-duplicates lines
		       :test (lambda (x y)(or (equalp x y)
					      (and (equalp (car x)(cdr y))
						   (equalp (cdr x)(car y))))))))

(setf *lines* (make-cube))

;;;; Make icosahedron

;;; The top of an icosahedron is a cap of five triangles
;;; r is the radius of the circle around the pentagonal base
;;; of the cap
(defparameter r (/ (* 2 (sin (/ pi 5)))))
;;; So we can use pythagoras to find the height h of the cap
(defparameter h (sqrt (- 1 (expt r 2))))
;;; The bottom of an icosahedron is a similar cap of five
;;; triangles, but upside down (ofcourse) and 36 degrees out of phase
;;; The horizontal displacement, d, of the corners between upper and lower
;;; rings is
(defparameter d (* 2 r (sin (/ pi 10))))
;;; Another application of Pytharoras give the height of each ring
;;; above and below the x=0 y=0 plane
(defparameter w (/ (sqrt (- 1 (expt d 2))) 2))

;;; Now we can write down the x,y,z coordinates of
;;; the 6 corners of the top cap, t0 the apex, and t1-5 the ring of 5
(defparameter t0 (vector 0 0 (+ h w)))
(defparameter t1 (vector r 0 w))
(defparameter t2 (vector (* r (cos (* 2/5 pi)))
			 (* r (sin (* 2/5 pi)))
			 w))
(defparameter t3 (vector (* r (cos (* 4/5 pi)))
			 (* r (sin (* 4/5 pi)))
			 w))
(defparameter t4 (vector (* r (cos (* 6/5 pi)))
			 (* r (sin (* 6/5 pi)))
			 w))
(defparameter t5 (vector (* r (cos (* 8/5 pi)))
			 (* r (sin (* 8/5 pi)))
			 w))

;;; the bottom is just a negation of the top
(defvecfun neg ()
  (x (- x))
  (y (- y))
  (z (- z)))

(defparameter b0 (neg t0))
(defparameter b1 (neg t1))
(defparameter b2 (neg t2))
(defparameter b3 (neg t3))
(defparameter b4 (neg t4))
(defparameter b5 (neg t5))

;;; This time I'm just going to type in a list of the edges
;;; and make the list of lines from that
(defun make-from-symbol-pairs (pairs)
  (loop for (from to) in pairs collect (cons (symbol-value from)
					     (symbol-value to))))

;;; The edges of an icosahedron
(defparameter icosahedron-pairs
  '((t0 t1)
    (t0 t2)
    (t0 t3)
    (t0 t4)
    (t0 t5)
    (t1 t2)
    (t2 t3)
    (t3 t4)
    (t4 t5)
    (t5 t1)
    (b0 b1)
    (b0 b2)
    (b0 b3)
    (b0 b4)
    (b0 b5)
    (b1 b2)
    (b2 b3)
    (b3 b4)
    (b4 b5)
    (b5 b1)
    (t1 b3)
    (t1 b4)
    (t2 b4)
    (t2 b5)
    (t3 b5)
    (t3 b1)
    (t4 b1)
    (t4 b2)
    (t5 b2)
    (t5 b3)))

(defun make-icosahedron ()
  (make-from-symbol-pairs icosahedron-pairs))


(defun make-globe (longitude latitude)
  (flet ((global-position (lat long)
			  (let ((lat (/ (* pi lat) latitude))
				(long (/ (* 2 pi long)
					 longitude)))
			    (vector (* (sin lat)(cos long))
				    (* (sin lat)(sin long))
				    (cos lat)))))
	(let ((lines '()))
	  (dotimes (lat latitude)
	    (dotimes (long longitude)
	      (push (cons (global-position lat long)
			  (global-position (+ lat 1) long))
		    lines)
	      (push (cons (global-position lat long)
			  (global-position lat (+ 1 long)))
		    lines)))
	  lines)))


(defun torus (i n wind loop)
  "Junction i of n in winding and looping round a torus"
  (let ((R 1)
	(\r 0.20)
	(theta (/ (* 2 pi i loop) n))
	(phi (/ (* 2 pi i wind) n)))
  (vector (+ (* R (cos theta))
	     (* \r (cos theta) (cos phi)))
	  (+ (* R (sin theta))
	     (* \r (sin theta) (cos phi)))
	  (* \r (sin phi)))))

(defun make-torus (n wind loop)
  "n line segments wind round a torus in the course of looping it"
  (let ((lines '()))
    (dotimes (i n lines)
      (push (cons (torus i n wind loop)
		  (torus (+ i 1) n wind loop))
	    lines))))