gistfile1.txt

user=> (time (nbody/main 100000)) ;;one above
(time (nbody/main 100000))
-0.16907516382852444
-0.16907985939166476
"Elapsed time: 2124.75375 msecs"
nil

user=> (time (main 100000)) ;;; nbody_v2
(time (main 100000))
-0.16907516382852444 
-0.16907985939168302 
"Elapsed time: 2784.951194 msecs"

;; You'll notice I didn't do any of the fancy loop unrolling I was talking about.
;; Should also note that I changed my clojure core to inline all versions of the 'nth' function
;; and that it inlines all of the math ops. (so that will effect performance on other machines)

my-nbody

(ns nbody)

(let [pi 3.141592653589793,
      +days-per-year+ 365.24,
      +solar-mass+ (* 4 pi pi)
      make-body (fn [x y z vx vy vz mass]
		  [x y z 
		   (* +days-per-year+ vx) 
		   (* +days-per-year+ vy) 
		   (* +days-per-year+ vz) 
		   (* +solar-mass+ mass)])

      apply-forces (fn 
		     [dt a [jx jy jz jvx jvy jvz jm]]
		     (let [[ix iy iz ivx ivy ivz im] @a
			   dx (- ix jx)
			   dy (- iy jy)
			   dz (- iz jz)
			   distance-squared (+ (* dx dx) (* dy dy) (* dz dz))
			   distance (. java.lang.Math (sqrt distance-squared))
			   mag (/ dt (* distance-squared distance))
			   dxm (* dx mag)
			   dym (* dy mag) 
			   dzm (* dz mag)]
      
		       (reset! a [ix iy iz 
				  (- ivx (* dxm jm))
				  (- ivy (* dym jm))
				  (- ivz (* dzm jm)) im])
		       [jx jy jz 
			(+ jvx (* dxm im))
			(+ jvy (* dym im))
			(+ jvz (* dzm im)) jm]))

      advance (fn [bodies dt]						
					;(println bodies)
		(let [bodies2 (reverse (loop 
					   [a (atom (first bodies)) bs (rest bodies) li (list)]
					 (let [new-bs (map apply-forces (repeat dt) (repeat a) bs)]
					   (if (== 0 (count new-bs))
					     (conj li @a)
					     (recur (atom (first new-bs)) (rest new-bs) (conj li @a))))))]
					;(println bodies2)
		  (let  [bodies3
			 (map (fn [[x y z vx vy vz m]] 
				[(+ x (* dt vx)) (+ y (* dt vy)) (+ z (* dt vz)) vx vy vz m]) bodies2)]
					;(println bodies3)
		    bodies3)))
 
      offset-momentum (fn [[jup sat nep ura [sx sy sz _ _ _ sm] :as bodies]]
			(let [[px py pz] (reduce (fn [[p1 p2 p3] [p4 p5 p6]]
						   [(+ p1 p4) (+ p2 p5) (+ p3 p6)])
						 (map (fn [[_ _ _ vx vy vz m]]
							[(* vx m) (* vy m) (* vz m)]) bodies))]
			  [jup sat nep ura [sx sy sz 
					    (/  (- px) +solar-mass+) 
					    (/ (- py) +solar-mass+) 
					    (/ (- pz) +solar-mass+) sm]]))

      energy (fn [system]
	       (loop [[xa ya za vxa vya vza mass-a] (first system) bs (rest system) e0 0.0]
		 (let [e1 (+ e0 (* 0.5
				   mass-a
				   (+ (* vxa vxa)
				      (* vya vya)
				      (* vza vza))))]

		   (if (first bs)
		     (recur (first bs) (rest bs) 
			    (- e1 (reduce + (map (fn [[xb yb zb _ _ _ mass-b]]
						   (let [dx (- xa xb)
							 dy (- ya yb)
							 dz (- za zb)
							 distance (. java.lang.Math 
								     (sqrt (+ (* dx dx) (* dy dy) (* dz dz))))]
						     (/ (* mass-a mass-b) distance)))
						 bs))))
		     (- e1 (reduce + (map (fn [[xb yb zb _ _ _ mass-b]]
					    (let [dx (- xa xb)
						  dy (- ya yb)
						  dz (- za zb)
						  distance (. java.lang.Math 
							      (sqrt (+ (* dx dx) (* dy dy) (* dz dz))))]
					      (/ (* mass-a mass-b) distance)))
					  bs)))))))]


  (let [bodies (offset-momentum
		[(make-body 4.84143144246472090e+00
			    -1.16032004402742839e+00
			    -1.03622044471123109e-01
			    1.66007664274403694e-03
			    7.69901118419740425e-03 
			    -6.90460016972063023e-05 
			    9.54791938424326609e-04),
		 (make-body 8.34336671824457987e+00
			    4.12479856412430479e+00
			    -4.03523417114321381e-01
			    -2.76742510726862411e-03
			    4.99852801234917238e-03
			    2.30417297573763929e-05
			    2.85885980666130812e-04)
		 (make-body 1.28943695621391310e+01
			    -1.51111514016986312e+01
			    -2.23307578892655734e-01
			    2.96460137564761618e-03
			    2.37847173959480950e-03
			    -2.96589568540237556e-05
			    4.36624404335156298e-05)
		 (make-body 1.53796971148509165e+01
			    -2.59193146099879641e+01
			    1.79258772950371181e-01
			    2.68067772490389322e-03
			    1.62824170038242295e-03
			    -9.51592254519715870e-05
			    5.15138902046611451e-05)
		 (make-body 0.0 0.0 0.0 0.0 0.0 0.0 1.0)])]
    (defn main [n]      
      (do
	(println (energy bodies))
	(println (energy (let [max n]
			   (loop [system bodies c 0]
			     (if (< c max)
			       (recur (advance system 0.01) (+ c 1))
			       system)))))))))