;;; -*- Mode:Common-Lisp; Package:TV; Base:10 -*-
;;;
;;; three-body
;;;
;;; This displays a three-dimensional view of three bodies in quasi-gravitational attraction.
;;;
;;; jsp 28-November-88
;;;
;;; (c) John S. Pezaris 1988.  All rights reserved.


(defvar body-list nil)
(defvar body-list-length nil)
(defvar body-font fonts:balls-16)

(defstruct body
  (acc )
  (vel )
  (pos )
  (mass 0)
  (size 0))

(defun three-body (&optional (stream *terminal-io*) (length 5) (times nil)
		   (max-size 16))
  
  (macrolet
    ((x-acc (b) `(first  (body-acc ,b)))
     (y-acc (b) `(second (body-acc ,b)))
     (z-acc (b) `(third  (body-acc ,b)))
     (x-vel (b) `(first  (body-vel ,b)))
     (y-vel (b) `(second (body-vel ,b)))
     (z-vel (b) `(third  (body-vel ,b)))
     (x-pos (b) `(first  (body-pos ,b)))
     (y-pos (b) `(second (body-pos ,b)))
     (z-pos (b) `(third  (body-pos ,b)))
     (x-pix (b xlim) `(min ,xlim (max 0 (floor (+ (/ ,xlim 2.) (x-pos ,b))))))
     (y-pix (b ylim) `(min ,ylim (max 0 (floor (+ (/ ,ylim 2.) (y-pos ,b))))))
     (z-pix (b)      `(max 1 (floor (/ (z-pos ,b) 10.0))))
     )
    
    
    (multiple-value-bind (xlim ylim)
	(send stream :inside-size)
      
      (let* ((xlim/2 (/ xlim 2))
	     (ylim/2 (/ ylim 2))
	     (-xlim/2 (* -1.0 xlim/2))
	     (-ylim/2 (* -1.0 ylim/2))
	     (diameter (- (min xlim/2 ylim/2) 50))
	     (list (progn
		     
		     (if (not (and body-list (= length body-list-length)))
			 (progn
			   (setq body-list-length length
				 body-list (loop repeat length collecting (make-body)))))
		     
		     (dotimes (n length)
		       (let ((new-size (min 16 (+ 8 (random (- max-size 9)))))
			     (b (nth n body-list)))
			 (setf (body-acc b) (list 0.0 0.0 0.0))
			 (setf (body-vel b) (list (/ (- 6.0 (random 11)) 8.0)
						  (/ (- 6.0 (random 11)) 8.0)
						  (/ (- 6.0 (random 11)) 8.0)))
			 (setf (body-pos b) (list (* diameter (cos (* n (/ (* 2.0 pi) length))))
						  (* diameter (sin (* n (/ (* 2.0 pi) length))))
						  (* diameter (sin (* n (/ (* 2.0 pi) length))))))
			 (setf (body-mass b) (* 10.0 new-size new-size))
			 (setf (body-size b) new-size)
			 ))
		     body-list)))
	
	(send stream :clear-screen)
	
	(loop until (if times (= (setf times (1- times)) 0) NIL)
	      do
	      
	      (dotimes (l length)
		;; calculate attraction of this body to the other bodies, this will set the x-acc and y-acc.
		(let ((new-x-acc 0.0)
		      (new-y-acc 0.0)
		      (new-z-acc 0.0)
		      (body (nth l body-list)))
		  (dotimes (ll length)
		    (let ((other-body (nth ll body-list)))
		      (if (not (eq body other-body))	; don't do it to yourself !!!
			  (let* ((x-dist (- (x-pos other-body) (x-pos body)))
				 (y-dist (- (y-pos other-body) (y-pos body)))
				 (z-dist (- (z-pos other-body) (z-pos body)))
				 (dist^2 (+ (* x-dist x-dist)
					    (* y-dist y-dist)
					    (* z-dist z-dist)))
				 (dist (sqrt dist^2))
				 (new-acc 0.0)
				 (new-acc/dist 0.0))

			    (if (> dist 0.1)
				(progn
				  (setq new-acc (* (/ (body-mass other-body) dist^2)
						   (cond ((< dist 100.0) -1.0)
;							 ((< dist 3.0)   -10.0)
;							 ((< dist 1.0)   -100.0)
							 (t 1.0))))
				  (setq new-acc/dist (/ new-acc dist))
				  (incf new-x-acc (* new-acc/dist x-dist))
				  (incf new-y-acc (* new-acc/dist y-dist))
				  (incf new-z-acc (* new-acc/dist z-dist))
				  )
				(progn
				  (incf new-x-acc (- 5.0 (random 10.0)))
				  (incf new-y-acc (- 5.0 (random 10.0)))
				  (incf new-z-acc (- 5.0 (random 10.0)))
				  )
				)
			    ))
		      ))
		  (setf (x-acc body) new-x-acc)
		  (setf (y-acc body) new-y-acc)
		  (setf (z-acc body) new-z-acc)
		  )
		)
	      
	      (dotimes (l length)
		(let* ((body (nth l body-list))
		       (old-x-pix (x-pix body xlim))
		       (old-y-pix (y-pix body ylim))
		       (old-z-pix (z-pix body))
		       (new-x-pix 0)
		       (new-y-pix 0)
		       (new-z-pix 0)
		       (body-size (body-size body)))
		  
		  ;; set the new velocities
		  (incf (x-vel body) (x-acc body))
		  (incf (y-vel body) (y-acc body))
		  (incf (z-vel body) (z-acc body))
		  
		  ;; make sure that if things get too fast, they slow down!
		  (if (< 10.0 (abs (x-vel body)))
		      (progn
			(setf (x-vel body) (* (x-vel body) 0.9))
			(setf (x-acc body) 0.0)))
		  
		  (if (< 10.0 (abs (y-vel body)))
		      (progn
			(setf (y-vel body) (* (y-vel body) 0.9))
			(setf (y-acc body) 0.0)))
		  
		  (if (< 10.0 (abs (z-vel body)))
		      (progn
			(setf (z-vel body) (* (z-vel body) 0.9))
			(setf (z-acc body) 0.0)))
		  
		  ;; set the new positions
		  (incf (x-pos body) (x-vel body))
		  (incf (y-pos body) (y-vel body))
		  (incf (z-pos body) (z-vel body))
		  
		  ;; check for maximum position in x direction
		  (if (<= (- xlim/2 max-size 4) (x-pos body))
		      (progn
			(setf (x-pos body) (- xlim/2 max-size 5))
			(setf (x-vel body) (* -1.0 (max 0.1 (abs (x-vel body))))))
		      (if (>= (+ -xlim/2 1) (x-pos body))
			  (progn
			    (setf (x-pos body) (+ -xlim/2 2))
			    (setf (x-vel body) (* 1.0 (max 0.1 (abs (x-vel body))))))))
		  
		  ;; check for maximum position in y direction
		  (if (<= (- ylim/2 max-size 4) (y-pos body))
		      (progn
			(setf (y-pos body) (- ylim/2 max-size 5))
			(setf (y-vel body) (* -1.0 (max 0.1 (abs (y-vel body))))))
		      (if (>= (+ -ylim/2 1) (y-pos body))
			  (progn
			    (setf (y-pos body) (+ -ylim/2 2))
			    (setf (y-vel body) (* 1.0 (max 0.1 (abs (y-vel body))))))))
		  
		  ;; check for maximum position in z direction
		  (if (<= 100.0 (z-pos body))
		      (progn
			(setf (z-pos body) 100.0)
			(setf (z-vel body) (* -1.0 (max 0.1 (abs (z-vel body))))))
		      (if (>= 0.0 (z-pos body))
			  (progn
			    (setf (z-pos body) 0.0)
			    (setf (z-vel body) (* 1.0 (max 0.1 (abs (z-vel body))))))))
		  
		  ;; calculate the new pixel position
		  (setq new-x-pix (x-pix body xlim))
		  (setq new-y-pix (y-pix body ylim))
		  (setq new-z-pix (z-pix body))
		  
		  (tv:prepare-sheet (stream)
		    ;; erase old image
		    (sys:%draw-rectangle old-z-pix old-z-pix
					 old-x-pix old-y-pix
					 tv:alu-setz stream)
		    ;; draw new one
		    (sys:%draw-rectangle new-z-pix new-z-pix
					 new-x-pix new-y-pix
					 tv:alu-seta stream)
		    )
		  
		  )
		)
	      ))
      )
    ))