(require (planet "simulation.ss" ("williams" "simulation.plt")))
(require (planet "random-distributions.ss" ("williams" "science.plt")))
(require (planet "math.ss" ("williams" "science.plt")))

;;; simulation parameters
(define F 100.0)                        ; force constant
(define p-life 5000.0)                  ; maximum particle lifetime
(define-values (WIDTH HEIGHT) (values 600 600))
(define-values (n-min n-max) (values 10 30))
(define-values (x-min x-max) (values -1000.0 1000.0))
(define-values (y-min y-max) (values -1000.0 1000.0))
(define-values (v-min v-max) (values -1.0 1.0))

;;; global variables - set in run-simulation
(define n 0)                            ; total number of particles
(define p #f)                           ; particle vector
(define x #f)                           ; particle position x vector
(define y #f)                           ; particle position y vector
(define dx/dt #f)                       ; particle position dx/dt vector
(define dy/dt #f)                       ; particle position dy/dt vector

;;; convert model coordinates to screen coordinates
(define (model-coords->screen-coords x y)
  (values (* WIDTH (/ (- x x-min) (- x-max x-min)))
          (* HEIGHT (/ (- y y-min) (- y-max y-min)))))

;;; particle model
(define-process (particle i)
  (let ((death-time (+ (current-simulation-time) (random-flat 0.0 p-life))))
    ;; let others know we exist
    (vector-set! p i self)
    ;; initialize velocities
    (vector-set! dx/dt i (make-continuous-variable (random-flat v-min v-max)))
    (vector-set! dy/dt i (make-continuous-variable (random-flat v-min v-max)))
    ;; initialize positions
    (vector-set! x i (make-continuous-variable (random-flat x-min x-max)))
    (vector-set! y i (make-continuous-variable (random-flat y-min y-max)))
    ;; continuous model
    (work/continuously
      until (>= (current-simulation-time) death-time)
      (let ((xi (variable-value (vector-ref x i)))
            (yi (variable-value (vector-ref y i)))
            (x-dot 0.0)
            (y-dot 0.0))
        ;; compute accelerations
        (do ((ii 0 (+ ii 1)))
            ((= ii n) (void))
          (when (and (not (= ii i))
                     (vector-ref p ii))
            (let* ((xii (variable-value (vector-ref x ii)))
                   (yii (variable-value (vector-ref y ii)))
                   (r2 (max 100.0 
                            (+ (* (- xii xi) (- xii xi))
                               (* (- yii yi) (- yii yi))))))
              (set! x-dot (+ x-dot (* (sign (- xii xi)) F (/ r2))))
              (set! y-dot (+ y-dot (* (sign (- yii yi)) F (/ r2)))))))
        ;; set velocity derivatives (i.e. accelerations)
        (set-variable-dt! (vector-ref dx/dt i) x-dot)
        (set-variable-dt! (vector-ref dy/dt i) y-dot)
        ;; set position derivatives (i.e. velocities)
        (set-variable-dt! (vector-ref x i) (variable-value (vector-ref dx/dt i)))
        (set-variable-dt! (vector-ref y i) (variable-value (vector-ref dy/dt i)))))
    ;; let others know we no longer exist
    (vector-set! p i #f)))

;; simulation model
(define (run-simulation)
  (with-new-simulation-environment
    ;; randomize
    (random-source-randomize! (current-random-source))
    ;; fixed step size
    (current-simulation-step-size 1.0)
    (current-simulation-control #f)
    ;; set number of particles randomly
    (set! n (+ n-min (random-uniform-int (- n-max n-min))))
    ;; create global variable vectors
    (set! p (make-vector n #f))
    (set! x (make-vector n))
    (set! y (make-vector n))
    (set! dx/dt (make-vector n))
    (set! dy/dt (make-vector n))
    ;; create the particles
    (do ((i 0 (+ i 1)))
        ((= i n) (void))
      (schedule now (particle i)))
    ;; set simulation monitor
    (current-simulation-monitor
     (lambda ()
       (let ((dc (send canvas get-dc)))
         (do ((i 0 (+ i 1)))
             ((= i n) (void))
           (when (vector-ref p i)
             (let-values (((sx sy) (model-coords->screen-coords
                                    (variable-value (vector-ref x i))
                                    (variable-value (vector-ref y i)))))
               (send dc draw-point sx sy)))))))
  ;; start the simulation
  (start-simulation)))

;;; simulation graphics

(define frame (instantiate frame% ("Interacting Particles") (x 0) (y 0)))
(send frame show #t)

(define my-canvas%
  (class canvas%
    (inherit get-dc)
    (define/override (on-paint)
      (run-simulation))
    (super-instantiate ())))

(define canvas (instantiate my-canvas% (frame)
                 (min-width WIDTH)
                 (min-height HEIGHT)))