#|
nbody-ics-test.scm: Test suite for nbody-ics.scm.
Copyright (C) 2007 Will M. Farr <farr@mit.edu>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|#

(module nbody-ics-test mzscheme
  (require (planet "test.ss" ("schematics" "schemeunit.plt"))
           "nbody-ics.scm"
           (lib "42.ss" "srfi")
           (lib "list.ss"))
  
  (provide nbody-ics-test-suite)
  
  (define (norm-squared v)
    (sum-ec (:vector x v) (* x x)))
  
  (define (norm v) (sqrt (norm-squared v)))
  
  (define (distance-squared v w)
    (sum-ec (:parallel (:vector x v) (:vector y w))
            (let ((dx (- x y)))
              (* dx dx))))
  
  (define (distance v w) (sqrt (distance-squared v w)))
  
  (define (body-kinetic-energy b)
    (/ (norm-squared (body-p b))
       (* 2.0 (body-m b))))
  
  (define (body-potential-energy b1 b2)
    (let ((r (distance (body-q b1) (body-q b2))))
      (- (/ (* (body-m b1) (body-m b2))
            r))))
  
  (define (kinetic-energy bs)
    (sum-ec (:vector b bs) (body-kinetic-energy b)))
  
  (define (potential-energy bs)
    (sum-ec (:vector b1 (index i) bs)
            (sum-ec (:range j (+ i 1) (vector-length bs))
                    (let ((b2 (vector-ref bs j)))
                      (body-potential-energy b1 b2)))))
  
  (define (energy bs) (+ (potential-energy bs) (kinetic-energy bs)))
  
  (define ((close? eps) a b)
    (< (abs (- a b)) (abs eps)))
  
  (define-simple-check (check-close? eps a b)
    ((close? eps) a b))
  
  (define *eps* 1e-10)
  
  (define nbody-ics-test-suite
    (test-suite
     "nbody-ics.scm test suite"
     (test-case
      "adjust-frame! test"
      (let ((bs (adjust-frame! (make-hot-spherical 1000))))
        (check-close? *eps* (norm (center-of-mass bs)) 0.0)
        (check-close? *eps* (norm (total-momentum bs)) 0.0)))
     (test-case
      "make-cold-spherical test"
      (let ((bs (make-cold-spherical 100)))
        (check-close? *eps* (kinetic-energy bs) 0.0)
        (check-close? 5e-2 (potential-energy bs) -0.25)))
     (test-case
      "make-hot-spherical test"
      (let ((bs (make-hot-spherical 100)))
        (check-close? 5e-2 (kinetic-energy bs) 0.25)
        (check-close? 5e-2 (potential-energy bs) -0.5))
      (let ((Q (random)))
        (let ((bs (make-hot-spherical 100 Q)))
          (check-close? 1e-1 (energy bs) -0.25)
          (check-close? 5e-2 (/ (kinetic-energy bs)
                                (abs (potential-energy bs)))
                        Q))))
     (test-case
      "make-plummer test"
      (let ((bs (list->vector
                 (sort (vector->list (make-plummer 100))
                       (lambda (b1 b2)
                         (let ((r1 (norm-squared (body-q b1)))
                               (r2 (norm-squared (body-q b2))))
                           (< r1 r2)))))))
        (check-close? 5e-2 (kinetic-energy bs) 0.25)
        (check-close? 1e-1 (potential-energy bs) -0.5)
        (check-close? 1.5e-1 (norm (body-q (vector-ref bs 24))) 0.48)
        (check-close? 1.5e-1 (norm (body-q (vector-ref bs 49))) 0.77)
        (check-close? 2.5e-1 (norm (body-q (vector-ref bs 74))) 1.28))))))