; This benchmark was obtained from Andrew Wright.
; 970215 / wdc Changed {box, unbox, set-box!} to {list, car, set-car!},
;              flushed #% prefixes, defined void,
;              and added nbody-benchmark.
; 981116 / wdc Replaced nbody-benchmark by main, added apply:+.

(library (rnrs-benchmarks nbody)
  (export main)
  (import
    (rnrs)
    (rnrs arithmetic flonums) 
    (rnrs r5rs)
    (rnrs mutable-pairs) 
    (rnrs-benchmarks))
  
  (define void
    (let ((invisible (string->symbol "")))
      (lambda () invisible)))
  
  (define (apply:+ xs)
    (do ((result 0.0 (fl+ result (car xs)))
         (xs xs (cdr xs)))
        ((null? xs) result)))
  
  ;; code is slightly broken...
  
  (define vect (vector))
  ; minimal standard random number generator
  ; 32 bit integer version
  ; cacm 31 10, oct 88
  ;
  
  (define *seed* (list 1))
  
  (define (srand seed)
    (set-car! *seed* seed))
  
  (define (rand)
    (let* ((hi (quotient (car *seed*) 127773))
           (lo (modulo (car *seed*) 127773))
           (test (- (* 16807 lo) (* 2836 hi))))
      (if (> test 0)
          (set-car! *seed* test)
          (set-car! *seed* (+ test 2147483647)))
      (car *seed*)))
  
  ;; return a random number in the interval [0,n)
  (define random
    (lambda (n)
      (modulo (abs (rand)) n)))
      
  
  (define (array-ref a . indices)
    (let loop ((a a) (indices indices))
      (if (null? indices)
          a
  	(loop (vector-ref a (car indices)) (cdr indices)))))
  
  (define (atan0 y x)
    (if (and (= x y) (= x 0)) 0 (atan y x)))
  
  
  ;change this to desired precision
  ;measured in order of expansions calculated
  (define precision 10)
  
  ;;; =========================================================
  ;;;                        Algorithm
  ;;; =========================================================
  
  (define (cartesian-algorithm tree)
    (up! tree
         cartesian-make-multipole-expansion
         cartesian-multipole-shift
         cartesian-expansion-sum)
    (down! tree
           cartesian-multipole-to-local-convert
           cartesian-local-shift
           cartesian-eval-local-expansion
           cartesian-expansion-sum
           cartesian-zero-expansion))
  
  (define (spherical-algorithm tree)
    (up! tree
         spherical-make-multipole-expansion
         spherical-multipole-shift
         spherical-expansion-sum)
    (down! tree
           spherical-multipole-to-local-convert
           spherical-local-shift
           spherical-eval-local-expansion
           spherical-expansion-sum
           spherical-zero-expansion))
  
  ;;; The upward path in the algorithm calculates
  ;;; multipole expansions at every node.
  
  (define (up! tree
               make-multipole-expansion
               multipole-shift
               expansion-sum)
    (let loop ((node (tree-body tree)))
      (define center (node-center node))
  
      (if (leaf-node? node)
  
  	(let ((multipole-expansion
  	       (expansion-sum
  		(map (lambda (particle)
  		       (make-multipole-expansion
  			(pt- (particle-position particle) center)
  			(particle-strength particle)))
  		     (node-particles node)))))
  	  (set-node-multipole-expansion! node multipole-expansion)
  	  (cons center multipole-expansion))
  
  	(let ((multipole-expansion
  	       (expansion-sum
  		(map (lambda (child)
  		       (define center-and-expansion (loop child))
  		       (multipole-shift (pt- (car center-and-expansion)
  					     center)
  					(cdr center-and-expansion)))
  		     (node-children node)))))
  	  (set-node-multipole-expansion! node multipole-expansion)
  	  (cons center multipole-expansion)))))
  
  ;;; Downward path of the algorithm which calculates local expansionss
  ;;; at every node and accelerations and potentials at each particle.
  
  (define (down! tree
                 multipole-to-local-convert
                 local-shift
                 eval-local-expansion
                 expansion-sum
                 zero-expansion)
    (let loop ((node (tree-body tree))
  	     (parent-local-expansion (zero-expansion))
  	     (parent-center (node-center (tree-body tree))))
      (let* ((center (node-center node))
  	   (interactive-sum
              (expansion-sum
  	     (map (lambda (interactive)
  		    (multipole-to-local-convert
  		     (pt- (node-center interactive) center)
  		     (node-multipole-expansion interactive)))
  		  (node-interactive-field node))))
  	   (local-expansion
  	    (expansion-sum (list (local-shift (pt- center parent-center)
  					      parent-local-expansion)
  				 interactive-sum))))
        (if (leaf-node? node)
  	  (eval-potentials-and-accelerations
  	   node
  	   local-expansion
  	   eval-local-expansion)
  	  (for-each (lambda (child) (loop child local-expansion center))
  		    (node-children node))))))
  
  (define (eval-potentials-and-accelerations
           node
           local-expansion
           eval-local-expansion)
    (let ((center (node-center node))
          (near-field (apply append (map node-particles
  					(node-near-field node)))))
      (for-each (lambda (particle)
  		(let* ((pos (particle-position particle))
  		       (far-field-accel-and-poten
  			(eval-local-expansion (pt- pos center)
  					      local-expansion)))
  		  (set-particle-acceleration!
  		   particle
  		   (pt+ (car far-field-accel-and-poten)
  			(sum-vectors
  			 (map (lambda (near)
  				(direct-accel (pt- (particle-position near)
  						   pos)
  					      (particle-strength near)))
  			      (nfilter near-field
  				      (lambda (near)
  					(not (eq? near particle))))))))
  		  (set-particle-potential!
  		   particle
  		   (+ (cdr far-field-accel-and-poten)
  		      (apply:+ (map (lambda (near)
  				      (direct-poten
  				       (pt- (particle-position near) pos)
  				       (particle-strength near)))
  				    (nfilter near-field
  					    (lambda (near)
  					      (not
  					       (eq? near particle))))))))))
  	      (node-particles node))))
  
  ;;; ================================================================
  ;;;                    Expansion Theorems
  ;;; ================================================================
  
  (define (cartesian-make-multipole-expansion pt strength)
    (let ((min-x (- (pt-x pt)))
          (min-y (- (pt-y pt)))
          (min-z (- (pt-z pt))))
      (make-cartesian-expansion
       (lambda (i j k) (* strength
                          (expt min-x i)
                          (expt min-y j)
                          (expt min-z k)
                          (n1/prod-fac i j k))))))
  
  (define (spherical-make-multipole-expansion pt strength)
    (let ((r (pt-r pt))
          (theta (pt-theta pt))
          (phi (pt-phi pt)))
      (make-spherical-expansion
       (lambda (l m)
         (* strength
            (expt r l)
            (eval-spherical-harmonic l (- m) theta phi))))))
  
  ;;; ==================================================================
  ;;;                    Shifting lemmas
  ;;; ==================================================================
  
  ;;; Shift multipole expansion
  
  (define (cartesian-multipole-shift pt multipole-expansion)
    (let ((pt-expansion (cartesian-make-multipole-expansion pt 1)))
      (make-cartesian-expansion
       (lambda (i j k)
         (sum-3d i j k
  	       (lambda (l m n)
  		 (* (array-ref multipole-expansion l m n)
  		    (array-ref pt-expansion (- i l) (- j m) (- k n)))))))))
  
  (define (spherical-multipole-shift pt multipole-expansion)
    (let ((r (pt-r pt))
          (theta (pt-theta pt))
          (phi (pt-phi pt)))
      (letrec ((foo (lambda (a b) 
  		    (if (< (* a b) 0)
  			(expt -1 (min (abs a) (abs b)))
  			1)))
  	     (bar (lambda (a b) 
  		    (/ (expt -1 a)
  		       (sqrt (* (fac (- a b)) (fac (+ a b))))))))
        (let ((pt-expansion (make-spherical-expansion
  			   (lambda (l m)
  			     (* (eval-spherical-harmonic l m theta phi)
  				(bar l m)
  				(expt r l))))))
  	(make-spherical-expansion
  	 (lambda (j k)
  	   (sum-2d j
  		   (lambda (l m)
  		     (if (> (abs (- k m)) (- j l))
  			 0
  			 (* (spherical-ref multipole-expansion (- j l) (- k m))
  			    (foo m (- k m))
  			    (bar (- j l) (- k m))
  			    (spherical-ref pt-expansion l (- m))
  			    (/ (bar j k))))))))))))
  
  ;;; Convert multipole to local
  
  (define (cartesian-multipole-to-local-convert pt multipole-expansion)
    (define pt-expansion
      (let* ((n1/radius (/ (pt-r vect)))
  	   (n2cosines (pt-scalar* (* 2 n1/radius) vect))
  	   (x (pt-x n2cosines))
  	   (y (pt-y n2cosines))
  	   (z (pt-z n2cosines)))
        (make-cartesian-expansion
         (lambda (i j k)
  	 (define ijk (+ i j k))
  	 (* (expt -1 ijk)
  	    (expt n1/radius (+ 1 ijk))
  	    (prod-fac i j k)
  	    (sum-3d (/ i 2) (/ j 2) (/ k 2)
  		    (lambda (l m n)
  		      (* (fac-1 (- ijk l m n))
  			 (n1/prod-fac l m n)
  			 (n1/prod-fac (- i (* 2 l))
  				     (- j (* 2 m))
  				     (- k (* 2 n)))
  			 (expt x (- i (* 2 l)))
  			 (expt y (- j (* 2 m)))
  			 (expt z (- k (* 2 n)))))))))))
    (make-cartesian-expansion
     (lambda (i j k)
       (sum2-3d (- precision i j k)
  	      (lambda (l m n)
  		(* (array-ref multipole-expansion l m n)
  		   (array-ref pt-expansion (+ i l) (+ j m) (+ k n))))))))
  
  (define (spherical-multipole-to-local-convert pt multipole-expansion)
    (let* ((r (pt-r pt))
  	 (theta (pt-theta pt))
  	 (phi (pt-phi pt))
  	 (foo (lambda (a b c) 
  		(* (expt -1 b)
  		   (if (> (* a c) 0)
  		       (expt -1 (min (abs a) (abs c)))
  		       1))))
  	 (bar (lambda (a b) 
  		(/ (expt -1 a)
  		   (sqrt (* (fac (- a b)) (fac (+ a b)))))))
  	 (pt-expansion (make-spherical-expansion
  			(lambda (l m)
  			  (/ (eval-spherical-harmonic l m theta phi)
  			     (bar l m)
  			     (expt r (+ 1 l)))))))
      (make-spherical-expansion
       (lambda (j k)
         (* (bar j k)
  	  (sum-2d (- precision j 1)
  		  (lambda (l m)
  		    (* (spherical-ref multipole-expansion l m)
  		       (foo k l m)
  		       (bar l m)
  		       (spherical-ref pt-expansion (+ j l) (- m k))))))))))
  
  ;;; Shift local expansion
  
  (define (cartesian-local-shift pt local-expansion)
    (let* ((x (pt-x pt))
  	 (y (pt-y pt))
  	 (z (pt-z pt))
  	 (expts (make-cartesian-expansion
  		 (lambda (l m n) (* (expt x l)
  				    (expt y m)
  				    (expt z n))))))
      (make-cartesian-expansion
       (lambda (i j k)
         (sum2-3d (- precision i j k)
  		(lambda (l m n)
  		  (* (array-ref local-expansion (+ i l) (+ j m) (+ k n))
  		     (array-ref expts l m n)
  		     (n1/prod-fac l m n))))))))
  
  (define (spherical-local-shift pt local-expansion)
    (let* ((pt (pt- (make-pt 0 0 0) pt))
  	 (r (pt-r pt))
  	 (theta (pt-theta pt))
  	 (phi (pt-phi pt))
  	 (foo (lambda (a b c) 
  		(* (expt -1 a)
  		   (expt -1 (/ (+ (abs (- b c))
  				  (abs b)
  				  (- (abs c)))
  			       2)))))
  	 (bar (lambda (a b) 
  		(/ (expt -1 a)
  		   (sqrt (* (fac (- a b))
  			    (fac (+ a b)))))))
  	 (stuff (make-spherical-expansion
  		   (lambda (l m)
  		     (* (eval-spherical-harmonic l m theta phi)
  			(expt r l))))))
      (make-spherical-expansion
       (lambda (j k)
         (sum2-2d j (lambda (l m)
  		    (if (> (abs (- m k)) (- l j))
  			0
  			(* (spherical-ref local-expansion l m)
  			   (bar (- l j) (- m k))
  			   (bar j k)
  			   (spherical-ref stuff (- l j) (- m k))
  			   (/ (foo (- l j) (- m k) m) (bar l m))))))))))
  
  ;;; Evaluate the resulting local expansion at point pt.
  
  (define (cartesian-eval-local-expansion pt local-expansion)
    (let* ((x (pt-x pt))
  	 (y (pt-y pt))
  	 (z (pt-z pt))
  	 (local-expansion
  	  (make-cartesian-expansion
  	   (lambda (i j k)
  	     (* (array-ref local-expansion i j k)
  		(n1/prod-fac i j k)))))
  	 (expts (make-cartesian-expansion
  		 (lambda (i j k)
  		   (* (expt x i) (expt y j) (expt z k))))))
      (cons (make-pt
             (sum2-3d (- precision 1)
                      (lambda (l m n)
                        (* (array-ref expts l m n)
                           (+ 1 l)
                           (array-ref local-expansion (+ 1 l) m n))))
             (sum2-3d (- precision 1)
                      (lambda (l m n)
                        (* (array-ref expts l m n)
                           (+ 1 m)
                           (array-ref local-expansion l (+ 1 m) n))))
             (sum2-3d (- precision 1)
                      (lambda (l m n)
                        (* (array-ref expts l m n)
                           (+ 1 n)
                           (array-ref local-expansion l m (+ 1 n))))))
            (sum2-3d precision
                     (lambda (l m n)
                       (* (array-ref expts l m n)
                          (array-ref local-expansion l m n)))))))
  
  (define (spherical-eval-local-expansion pt local-expansion)
    (let* ((r (pt-r pt))
  	 (x (pt-x pt))
  	 (y (pt-y pt))
  	 (z (pt-z pt))
  	 (rho-sq (+ (* x x) (* y y)))
  	 (theta (pt-theta pt))
  	 (phi (pt-phi pt))
  	 (r-deriv
  	  (real-part
  	   (sum-2d (- precision 1)
  		   (lambda (l m)
  		     (* (spherical-ref local-expansion l m)
  			(expt r (- l 1))
  			l
  			(eval-spherical-harmonic l m theta phi))))))
  	 (theta-deriv
  	  (real-part
  	   (sum-2d (- precision 1)
  		   (lambda (l m)
  		     (* (spherical-ref local-expansion l m)
  			(expt r l)
  			(eval-spher-harm-theta-deriv l m theta phi))))))
  	 (phi-deriv
  	  (real-part
  	   (sum-2d (- precision 1)
  		   (lambda (l m)
  		     (* (spherical-ref local-expansion l m)
  			(expt r l)
  			(eval-spher-harm-phi-deriv l m theta phi)))))))
      (cons (make-pt
  	   (+ (* r-deriv (/ x r))
  	      (* theta-deriv (/ x (sqrt rho-sq) (+ z (/ rho-sq z))))
  	      (* phi-deriv -1 (/ y rho-sq)))
  	   (+ (* r-deriv (/ y r))
  	      (* theta-deriv (/ y (sqrt rho-sq) (+ z (/ rho-sq z))))
  	      (/ phi-deriv x (+ 1 (/ (* y y) x x))))
  	   (+ (* r-deriv (/ z r))
  	      (* theta-deriv (/ -1 r r) (sqrt rho-sq))))
  	  (real-part (sum-2d (- precision 1)
  			     (lambda (l m)
  			       (* (spherical-ref local-expansion l m)
  				  (eval-spherical-harmonic l m theta phi)
  				  (expt r l))))))))
  
  ;;; Direct calculation of acceleration and potential
  
  (define (direct-accel pt strength)
    (pt-scalar* (/ strength (expt (pt-r pt) 3))
                pt))
  
  (define (direct-poten pt strength)
    (/ strength (pt-r pt)))
  
  ;;; =================================================================
  ;;;                 TREES NODES PARTICLES and POINTS
  ;;; =================================================================
  
  (begin (begin (begin (define make-raw-tree
                         (lambda (tree-1 tree-2 tree-3)
                           (vector '<tree> tree-1 tree-2 tree-3)))
                       (define tree?
                         (lambda (obj)
                           (if (vector? obj)
                               (if (= (vector-length obj) 4)
                                   (eq? (vector-ref obj 0) '<tree>)
                                   #f)
                               #f)))
                       (define tree-1
                         (lambda (obj)
                           (if (tree? obj)
                               (void)
                               (error
                                 'tree-1
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-ref obj 1)))
                       (define tree-2
                         (lambda (obj)
                           (if (tree? obj)
                               (void)
                               (error
                                 'tree-2
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-ref obj 2)))
                       (define tree-3
                         (lambda (obj)
                           (if (tree? obj)
                               (void)
                               (error
                                 'tree-3
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-ref obj 3)))
                       (define set-tree-1!
                         (lambda (obj newval)
                           (if (tree? obj)
                               (void)
                               (error
                                 'set-tree-1!
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-set! obj 1 newval)))
                       (define set-tree-2!
                         (lambda (obj newval)
                           (if (tree? obj)
                               (void)
                               (error
                                 'set-tree-2!
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-set! obj 2 newval)))
                       (define set-tree-3!
                         (lambda (obj newval)
                           (if (tree? obj)
                               (void)
                               (error
                                 'set-tree-3!
                                 "~s is not a ~s"
                                 obj
                                 '<tree>))
                           (vector-set! obj 3 newval))))
                (define make-tree
                  (lambda (body low-left-front-vertex up-right-back-vertex)
                    ((lambda ()
                       (make-raw-tree
                         body
                         low-left-front-vertex
                         up-right-back-vertex)))))
                (define tree-body tree-1)
                (define tree-low-left-front-vertex tree-2)
                (define tree-up-right-back-vertex tree-3)
                (define set-tree-body! set-tree-1!)
                (define set-tree-low-left-front-vertex! set-tree-2!)
                (define set-tree-up-right-back-vertex! set-tree-3!))
         (begin (begin (define make-raw-node
                         (lambda (node-1
                                  node-2
                                  node-3
                                  node-4
                                  node-5
                                  node-6
                                  node-7
                                  node-8)
                           (vector
                             '<node>
                             node-1
                             node-2
                             node-3
                             node-4
                             node-5
                             node-6
                             node-7
                             node-8)))
                       (define node?
                         (lambda (obj)
                           (if (vector? obj)
                               (if (= (vector-length obj) 9)
                                   (eq? (vector-ref obj 0) '<node>)
                                   #f)
                               #f)))
                       (define node-1
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-1
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 1)))
                       (define node-2
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-2
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 2)))
                       (define node-3
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-3
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 3)))
                       (define node-4
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-4
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 4)))
                       (define node-5
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-5
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 5)))
                       (define node-6
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-6
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 6)))
                       (define node-7
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-7
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 7)))
                       (define node-8
                         (lambda (obj)
                           (if (node? obj)
                               (void)
                               (error
                                 'node-8
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-ref obj 8)))
                       (define set-node-1!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-1!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 1 newval)))
                       (define set-node-2!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-2!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 2 newval)))
                       (define set-node-3!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-3!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 3 newval)))
                       (define set-node-4!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-4!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 4 newval)))
                       (define set-node-5!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-5!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 5 newval)))
                       (define set-node-6!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-6!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 6 newval)))
                       (define set-node-7!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-7!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 7 newval)))
                       (define set-node-8!
                         (lambda (obj newval)
                           (if (node? obj)
                               (void)
                               (error
                                 'set-node-8!
                                 "~s is not a ~s"
                                 obj
                                 '<node>))
                           (vector-set! obj 8 newval))))
                (define make-node
                  (lambda (center
                           low-left-front-vertex
                           up-right-back-vertex
                           children
                           particles
                           multipole-expansion
                           near-field
                           interactive-field)
                    ((lambda ()
                       (make-raw-node
                         center
                         low-left-front-vertex
                         up-right-back-vertex
                         children
                         particles
                         multipole-expansion
                         near-field
                         interactive-field)))))
                (define node-center node-1)
                (define node-low-left-front-vertex node-2)
                (define node-up-right-back-vertex node-3)
                (define node-children node-4)
                (define node-particles node-5)
                (define node-multipole-expansion node-6)
                (define node-near-field node-7)
                (define node-interactive-field node-8)
                (define set-node-center! set-node-1!)
                (define set-node-low-left-front-vertex! set-node-2!)
                (define set-node-up-right-back-vertex! set-node-3!)
                (define set-node-children! set-node-4!)
                (define set-node-particles! set-node-5!)
                (define set-node-multipole-expansion! set-node-6!)
                (define set-node-near-field! set-node-7!)
                (define set-node-interactive-field! set-node-8!))
         (define leaf-node?
           (lambda (node) (null? (node-children node))))
         (begin (begin (define make-raw-particle
                         (lambda (particle-1
                                  particle-2
                                  particle-3
                                  particle-4
                                  particle-5
                                  particle-6)
                           (vector
                             '<particle>
                             particle-1
                             particle-2
                             particle-3
                             particle-4
                             particle-5
                             particle-6)))
                       (define particle?
                         (lambda (obj)
                           (if (vector? obj)
                               (if (= (vector-length obj) 7)
                                   (eq? (vector-ref obj 0) '<particle>)
                                   #f)
                               #f)))
                       (define particle-1
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-1
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 1)))
                       (define particle-2
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-2
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 2)))
                       (define particle-3
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-3
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 3)))
                       (define particle-4
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-4
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 4)))
                       (define particle-5
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-5
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 5)))
                       (define particle-6
                         (lambda (obj)
                           (if (particle? obj)
                               (void)
                               (error
                                 'particle-6
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-ref obj 6)))
                       (define set-particle-1!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-1!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 1 newval)))
                       (define set-particle-2!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-2!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 2 newval)))
                       (define set-particle-3!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-3!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 3 newval)))
                       (define set-particle-4!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-4!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 4 newval)))
                       (define set-particle-5!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-5!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 5 newval)))
                       (define set-particle-6!
                         (lambda (obj newval)
                           (if (particle? obj)
                               (void)
                               (error
                                 'set-particle-6!
                                 "~s is not a ~s"
                                 obj
                                 '<particle>))
                           (vector-set! obj 6 newval))))
                (define make-particle
                  (lambda (position
                           acceleration
                           d-acceleration
                           potential
                           d-potential
                           strength)
                    ((lambda ()
                       (make-raw-particle
                         position
                         acceleration
                         d-acceleration
                         potential
                         d-potential
                         strength)))))
                (define particle-position particle-1)
                (define particle-acceleration particle-2)
                (define particle-d-acceleration particle-3)
                (define particle-potential particle-4)
                (define particle-d-potential particle-5)
                (define particle-strength particle-6)
                (define set-particle-position! set-particle-1!)
                (define set-particle-acceleration! set-particle-2!)
                (define set-particle-d-acceleration! set-particle-3!)
                (define set-particle-potential! set-particle-4!)
                (define set-particle-d-potential! set-particle-5!)
                (define set-particle-strength! set-particle-6!))
         (begin (begin (define make-raw-pt
                         (lambda (pt-1 pt-2 pt-3)
                           (vector '<pt> pt-1 pt-2 pt-3)))
                       (define pt?
                         (lambda (obj)
                           (if (vector? obj)
                               (if (= (vector-length obj) 4)
                                   (eq? (vector-ref obj 0) '<pt>)
                                   #f)
                               #f)))
                       (define pt-1
                         (lambda (obj)
                           (if (pt? obj)
                               (void)
                               (error 'pt-1 "~s is not a ~s" obj '<pt>))
                           (vector-ref obj 1)))
                       (define pt-2
                         (lambda (obj)
                           (if (pt? obj)
                               (void)
                               (error 'pt-2 "~s is not a ~s" obj '<pt>))
                           (vector-ref obj 2)))
                       (define pt-3
                         (lambda (obj)
                           (if (pt? obj)
                               (void)
                               (error 'pt-3 "~s is not a ~s" obj '<pt>))
                           (vector-ref obj 3)))
                       (define set-pt-1!
                         (lambda (obj newval)
                           (if (pt? obj)
                               (void)
                               (error
                                 'set-pt-1!
                                 "~s is not a ~s"
                                 obj
                                 '<pt>))
                           (vector-set! obj 1 newval)))
                       (define set-pt-2!
                         (lambda (obj newval)
                           (if (pt? obj)
                               (void)
                               (error
                                 'set-pt-2!
                                 "~s is not a ~s"
                                 obj
                                 '<pt>))
                           (vector-set! obj 2 newval)))
                       (define set-pt-3!
                         (lambda (obj newval)
                           (if (pt? obj)
                               (void)
                               (error
                                 'set-pt-3!
                                 "~s is not a ~s"
                                 obj
                                 '<pt>))
                           (vector-set! obj 3 newval))))
                (define make-pt
                  (lambda (x y z) ((lambda () (make-raw-pt x y z)))))
                (define pt-x pt-1)
                (define pt-y pt-2)
                (define pt-z pt-3)
                (define set-pt-x! set-pt-1!)
                (define set-pt-y! set-pt-2!)
                (define set-pt-z! set-pt-3!)))
  
  ;(define-structure (tree
  ;		   body
  ;		   low-left-front-vertex
  ;		   up-right-back-vertex))
  ;
  ;(define-structure (node
  ;		   center
  ;		   low-left-front-vertex
  ;		   up-right-back-vertex
  ;		   children
  ;		   particles
  ;		   multipole-expansion
  ;		   near-field
  ;		   interactive-field))
  ;
  ;(define (leaf-node? node)
  ;  (null? (node-children node)))
  ;
  ;(define-structure (particle
  ;		   position
  ;		   acceleration
  ;		   d-acceleration
  ;		   potential
  ;		   d-potential
  ;		   strength))
  ;
  ;(define-structure (pt x y z))
  ;
  (define (pt-r pt)
    (sqrt (+ (* (pt-x pt) (pt-x pt))
             (* (pt-y pt) (pt-y pt))
             (* (pt-z pt) (pt-z pt)))))
  
  (define (pt-theta pt)
    (let ((x (pt-x pt))
          (y (pt-y pt))
          (z (pt-z pt)))
      (atan0 (sqrt (+ (* x x) (* y y))) z)))
  
  (define (pt-phi pt)
    (let ((x (pt-x pt)) (y (pt-y pt)))
      (atan0 y x)))
  
  (define (pt+ pt1 pt2)
    (make-pt (+ (pt-x pt1) (pt-x pt2))
             (+ (pt-y pt1) (pt-y pt2))
             (+ (pt-z pt1) (pt-z pt2))))
  
  (define (sum-vectors vectors)
    (make-pt (apply:+ (map pt-x vectors))
             (apply:+ (map pt-y vectors))
             (apply:+ (map pt-z vectors))))
  
  (define (pt- pt1 pt2)
    (make-pt (- (pt-x pt1) (pt-x pt2))
             (- (pt-y pt1) (pt-y pt2))
             (- (pt-z pt1) (pt-z pt2))))
  
  (define (pt-average pt1 pt2)
    (pt-scalar* .5 (pt+ pt1 pt2)))
  
  (define (pt-scalar* scalar pt)
    (make-pt (* scalar (pt-x pt))
             (* scalar (pt-y pt))
             (* scalar (pt-z pt))))
  
  (define (within-box? pt pt1 pt2)
    (and (<= (pt-x pt) (pt-x pt2)) (> (pt-x pt) (pt-x pt1))
         (<= (pt-y pt) (pt-y pt2)) (> (pt-y pt) (pt-y pt1))
         (<= (pt-z pt) (pt-z pt2)) (> (pt-z pt) (pt-z pt1))))
  
  ;;; ==========================================================
  ;;;                     Useful Things
  ;;; ==========================================================
  
  (define (nfilter list predicate)
    (let loop ((list list))
      (cond ((null? list) '())
            ((predicate (car list)) (cons (car list) (loop (cdr list))))
            (else (loop (cdr list))))))
  
  ;;; array in the shape of a pyramid with each
  ;;; element a function of the indices
  
  (define (make-cartesian-expansion func)
    (let ((expansion (make-vector precision 0)))
      (let loop1 ((i 0))
        (if (= i precision)
            expansion
            (let ((foo (make-vector (- precision i) 0)))
              (vector-set! expansion i foo)
              (let loop2 ((j 0))
                (if (= j (- precision i))
                    (loop1 (+ 1 i))
                    (let ((bar (make-vector (- precision i j) 0)))
                      (vector-set! foo j bar)
                      (let loop3 ((k 0))
                        (if (= k (- precision i j))
                            (loop2 (+ 1 j))
                            (begin (vector-set! bar k (func i j k))
                                   (loop3 (+ 1 k)))))))))))))
  
  ;;; array in the shape of a triangle with each
  ;;; element a function of the indices
  
  (define (make-spherical-expansion func)
    (let ((expansion (make-vector precision 0)))
      (let loop1 ((l 0))
        (if (= l precision)
            expansion
            (let ((foo (make-vector (+ 1 l) 0)))
              (vector-set! expansion l foo)
              (let loop2 ((m 0))
                (if (= m (+ 1 l))
                    (loop1 (+ 1 l))
                    (begin (vector-set! foo m (func l m))
                           (loop2 (+ 1 m))))))))))
  
  (define (spherical-ref expansion l m)
    (let ((conj (lambda (z) (make-rectangular (real-part z) (- (imag-part z))))))
      (if (negative? m)
  	(conj (array-ref expansion l (- m)))
  	(array-ref expansion l m))))
  
  (define (cartesian-expansion-sum expansions)
    (make-cartesian-expansion
     (lambda (i j k)
       (apply:+ (map (lambda (expansion)
                       (array-ref expansion i j k))
                     expansions)))))
  
  (define (spherical-expansion-sum expansions)
    (make-spherical-expansion
     (lambda (l m)
       (apply:+ (map (lambda (expansion)
                       (spherical-ref expansion l m))
                     expansions)))))
  
  (define (cartesian-zero-expansion)
    (make-cartesian-expansion (lambda (i j k) 0)))
  
  (define (spherical-zero-expansion)
    (make-spherical-expansion (lambda (l m) 0)))
  
  (define (sum-3d end1 end2 end3 func)
    (let loop1 ((l 0) (sum 0))
      (if (> l end1)
          sum
          (loop1
           (+ 1 l)
           (+ sum
              (let loop2 ((m 0) (sum 0))
                (if (> m end2)
                    sum
                    (loop2
                     (+ 1 m)
                     (+ sum
                        (let loop3 ((n 0) (sum 0))
                          (if (> n end3)
                              sum
                              (loop3 (+ 1 n)
                                     (+ sum (func l m n))))))))))))))
  
  (define (sum2-3d end func)
    (let loop1 ((l 0) (sum 0))
      (if (= l end)
          sum
          (loop1
           (+ 1 l)
           (+ sum
              (let loop2 ((m 0) (sum 0))
                (if (= (+ l m) end)
                    sum
                    (loop2
                     (+ 1 m)
                     (+ sum
                        (let loop3 ((n 0) (sum 0))
                          (if (= (+ l m n) end)
                              sum
                              (loop3 (+ 1 n)
                                     (+ sum (func l m n))))))))))))))
  
  (define (sum-2d end func)
    (let loop1 ((l 0) (sum 0))
      (if (> l end)
          sum
          (loop1 (+ 1 l)
                 (+ sum (let loop2 ((m (- l)) (sum 0))
                          (if (> m l)
                              sum
                              (loop2 (+ 1 m)
                                     (+ sum (func l m))))))))))
  
  (define (sum2-2d init func)
    (let loop1 ((l init) (sum 0))
      (if (= l precision)
          sum
          (loop1 (+ 1 l)
                 (+ sum (let loop2 ((m (- l)) (sum 0))
                          (if (> m l)
                              sum
                              (loop2 (+ 1 m)
                                     (+ sum (func l m))))))))))
  
  ;;; Storing factorials in a table
  
  (define fac
    (let ((table (make-vector (* 4 precision) 0)))
      (vector-set! table 0 1)
      (let loop ((n 1))
        (if (= n (* 4 precision))
            (lambda (x) (vector-ref table x))
            (begin (vector-set! table
                                n
                                (* n (vector-ref table (- n 1))))
                   (loop (+ 1 n)))))))
  
  ;;; The table for (* (-0.5) (-1.5) (-2.5) ... (+ -0.5 -i 1))
  
  (define fac-1
    (let ((table (make-vector precision 0)))
      (vector-set! table 0 1)
      (let loop ((n 1))
        (if (= n precision)
            (lambda (x) (vector-ref table x))
            (begin (vector-set! table
                                n
                                (* (- .5 n)
                                   (vector-ref table (- n 1))))
                   (loop (+ 1 n)))))))
  
  (define fac-2
    (let ((table (make-vector (* 4 precision) 0)))
      (vector-set! table 0 1)
      (let loop ((n 1))
        (if (= n (* 4 precision))
            (lambda (n) (if (< n 0)
                            1
                            (vector-ref table n)))
            (begin (vector-set! table n (* (if (even? n) 1 n)
                                           (vector-ref table (- n 1))))
                   (loop (+ 1 n)))))))
  
  ;;; Storing the products of factorials in a table.
  
  (define prod-fac
    (let ((table (make-cartesian-expansion
                  (lambda (i j k) (* (fac i) (fac j) (fac k))))))
      (lambda (i j k) (array-ref table i j k))))
  
  (define n1/prod-fac
    (let ((table (make-cartesian-expansion
                  (lambda (i j k) (/ (prod-fac i j k))))))
      (lambda (i j k) (array-ref table i j k))))
  
  (define (assoc-legendre l m x)
    (cond ((= l m) (* (expt -1 m)
                      (fac-2 (- (* 2 m) 1))
                      (expt (- 1 (* x x)) (/ m 2))))
          ((= l (+ 1 m)) (* x (+ 1 (* 2 m)) (assoc-legendre m m x)))
          (else (/ (- (* x (- (* 2 l) 1) (assoc-legendre (- l 1) m x))
                      (* (+ l m -1) (assoc-legendre (- l 2) m x)))
                   (- l m)))))
  
  (define (eval-spherical-harmonic l m theta phi)
    (let ((mm (abs m)))
      (* (sqrt (/ (fac (- l mm)) (fac (+ l mm))))
         (assoc-legendre l mm (cos theta))
         (make-polar 1 (* m phi)))))
  
  (define (eval-spher-harm-phi-deriv l m theta phi)
    (* (eval-spherical-harmonic l m theta phi)
       m
       (make-rectangular 0 1)))
  
  (define (eval-spher-harm-theta-deriv l m theta phi)
    (let ((mm (abs m)))
      (* (sqrt (/ (fac (- l mm)) (fac (+ l mm))))
         (make-polar 1 (* m phi))
         (- (sin theta))
         (assoc-legendre-deriv l mm (cos theta)))))
  
  (define (assoc-legendre-deriv l m x)
    (cond ((= l m) (* (expt -1 (+ 1 m))
                      (fac-2 (- (* 2 m) 1))
                      m
                      (expt (- 1 (* x x)) (- (/ m 2) 1))
                      x))
          ((= l (+ 1 m)) (* (+ 1 (* 2 m))
                            (+ (assoc-legendre m m x)
                               (* x (assoc-legendre-deriv m m x)))))
          (else (/ (- (* (- (* 2 l) 1)
                         (+ (assoc-legendre (- l 1) m x)
                            (* x (assoc-legendre-deriv (- l 1) m x))))
                      (* (+ l m -1) (assoc-legendre-deriv (- l 2) m x)))
                   (- l m)))))
  
  ;;; ================================================================
  ;;;                          TREE CODE
  ;;; ================================================================
  
  (define (build-tree height near-size)
    (let* ((vertex1 (make-pt -10 -10 -10))
           (vertex2 (make-pt 10 10 10))
  	 (tree (make-tree '() vertex1 vertex2)))
      (let loop ((level 0) (pt1 vertex1) (pt2 vertex2))
        (let* ((half-diagonal (pt-scalar* .5 (pt- pt2 pt1)))
  	     (diag-length/2 (pt-x half-diagonal)))
  	(insert-node tree level pt1 pt2)
  	(if (< level height)
  	    (let ((child-pt1s
  		   (map (lambda (offset)
  			  (pt+ pt1 (pt-scalar* diag-length/2 offset)))
  			(list (make-pt 0 0 0)
  			      (make-pt 0 0 1)
  			      (make-pt 0 1 0)
  			      (make-pt 1 0 0)
  			      (make-pt 0 1 1)
  			      (make-pt 1 0 1)
  			      (make-pt 1 1 0)
  			      (make-pt 1 1 1)))))
  	      (for-each (lambda (child-pt1)
  			  (loop (+ 1 level)
  				child-pt1
  				(pt+ child-pt1 half-diagonal)))
  			child-pt1s)))))
      (calc-near-and-interaction tree near-size)
      tree))
  
  (define (insert-node tree level pt1 pt2)
    (let* ((center (pt-average pt1 pt2))
           (new-node (make-node center pt1 pt2 '() '() '() '() '())))
      (letrec ((insert-internal (lambda (node depth)
  				(if (= level depth)
  				    (set-node-children! node 
  							(cons new-node 
  							      (node-children node)))
  				    (insert-internal (find-child node center) 
  						     (+ 1 depth))))))
        (if (= level 0)
  	  (set-tree-body! tree new-node)
  	  (insert-internal (tree-body tree) 1)))))
  
  (define (find-child node pos)
    (let loop ((children (node-children node)))
      (let ((child (car children)))
        (if (within-box? pos
  		       (node-low-left-front-vertex child)
  		       (node-up-right-back-vertex child))
  	  child
  	  (loop (cdr children))))))
  
  (define (insert-particle tree particle)
    (let* ((pos (particle-position particle)))
      (letrec ((insert-internal (lambda (node)
  				(if (leaf-node? node)
  				    (set-node-particles!
  				     node
  				     (cons particle (node-particles node)))
  				    (insert-internal (find-child node pos))))))
        (if (within-box? pos
  		       (tree-low-left-front-vertex tree)
  		       (tree-up-right-back-vertex tree))
  	  (insert-internal (tree-body tree))
  	  (error 'insert-particle "particle not within boundaries of tree" particle)))))
  
  ;;; This function finds the near and
  ;;; interaction fields for every node in the tree.
  
  (define (calc-near-and-interaction tree near-size)
    (set-node-near-field! (tree-body tree) (list (tree-body tree)))
    (let loop ((node (tree-body tree)) (parent #f))
      (if parent
  	(let* ((center (node-center node))
  	       (dist (* near-size
  			(abs (- (pt-x center)
  				(pt-x (node-center parent)))))))
  	  (for-each
  	   (lambda (parent-near)
  	     (let ((interactives (list '())))
  	       (for-each
  		(lambda (child)
  		  (if (> (pt-r (pt- center (node-center child))) dist)
  		      (set-car! interactives (cons child (car interactives)))
  		      (set-node-near-field!
  		       node
  		       (cons child (node-near-field node)))))
  		(node-children parent-near))
  	       (set-node-interactive-field!
  		node
  		(append (car interactives) (node-interactive-field node)))))
  	   (node-near-field parent))))
      (for-each (lambda (child) (loop child node)) (node-children node))))
  
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; GO
  
  (define (initial-particle x y z m)
      (make-particle (make-pt x y z)
  		   (make-pt 0 0 0)
  		   (make-pt 0 0 0)
  		   0 0 m))
  
  (define (random-float bot top)
    (+ (* (- top bot)
          (/ (* (random 1000000) 1.0) 1000000.0))
       bot))
  
  (define (random-particle)
    (make-particle (make-pt (random-float -10.0 10.0)
  			  (random-float -10.0 10.0)
  			  (random-float -10.0 10.0))
  		 (make-pt 0 0 0)
  		 (make-pt 0 0 0)
  		 0 0 1.0))
  
  (define *particles* (list '()))
  
  (define (go depth precision n-particles)
    (let ((tree (build-tree depth 27))
  	(particles (let next ((i 0) (ps '()))
  		     (if (<= i n-particles)
  		       (next (+ i 1) (cons (random-particle) ps))
  		       ps))))
      (for-each (lambda (p) (insert-particle tree p)) particles)
      (cartesian-algorithm tree)
      (set-car! *particles* particles)))
  
  ;;; virtual time for cartesian-algorithm step
  ;;; (go 1 3 10)      0.31 seconds
  ;;; (go 3 5 128)  1397.31
  ;;; (go 3 5 256)  1625.29
  ;;; (go 3 5 512)  2380.35
  ;;; (go 2 5 128)    27.44 seconds
  
  (define (main . args)
    (run-benchmark
      "nbody"
      nbody-iters
      (lambda () #t)
      (lambda (i j k) (go i j k))
      2 5 128)))