#| 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))))))