;;;; Two-dimensional bouncing balls
;; Originally by Piotr Mitros
;; Modified by gjs to be a little bit cleaner, more Scheme-like. 
;; Distributed under GPL. 

;;; Kinetic Energy of a massive particle

(define ((KE m) v)
  (* 1/2 m (square v)))

;;; Potential energy of interaction: 
;;;   gravitational attraction 
;;;   and short-range repulsion.

(define ((PE G m1 m2 k) xy12)
  (let ((xy1 (ref xy12 0))
	(xy2 (ref xy12 1)))
    (let ((rsq (square (- xy1 xy2))))
      (- (/ k rsq)			; repulsion
	 (/ (* G m1 m2) (sqrt rsq))	; attraction
	 ))))

;;; xs = #( #(x1 y1) #(x2 y2) ... )
;;; vs = #( #(vx1 vy1) #(vx2 vy2) ... )

(define ((Lag G k) local)
  (let ((xs (s:->list (coordinate local)))
	(vs (s:->list (velocity local))))
    (- (reduce + 0
	       (map (KE 1) vs))
       (reduce + 0
	       (map (PE G 1 1 k)
		    (distinct-pairs xs))))))

;;; A useful utility.
(define s:->list
  (compose vector->list s:->vector))

#|
(set! *divide-out-terms* #f)

;;; Horrible expressions!
(pec ((Lag 1 'k)
      (up 't
	  (up (up 'x1 'y1) (up 'x2 'y2))
	  (up (up 'vx1 'vy1) (up 'vx2 'vy2)))))

(pec ((Lagrange-explicit (Lag 1 'k))
      (up 't
	  (up (up 'x1 'y1) (up 'x2 'y2))
	  (up (up 'vx1 'vy1) (up 'vx2 'vy2)))))
|#

;;; Numerical integration

(define (sysder G k) 
  (Lagrangian->state-derivative (Lag G k)))

;;; try G=1 k=1
;;; scale=+-10

(define win (frame -10 +10 -10 +10))

(define (clear) 
  (graphics-clear win))

(define (close) 
  (graphics-close win))

(define ((monitor win) state)
  (let ((xs (s:->list (coordinate state))))
    ;(clear)
    (for-each (lambda (xy)
		(plot-point win (ref xy 0) (ref xy 1)))
	      xs))
  state)

(begin
  (clear)
  ((evolve sysder 1 1)
   (up 0
       (up (up +7 0) (up -7 0))
       (up (up 0 -.003) (up 0 +.003)))
   (monitor win)
   .1
   50000
   1e-10))

;;; (close)