(= pi (* 4 (atan 1))
days-per-year* 365.24d0
solar-mass* (* 4d0 pi pi))
(def make-body (x y z vx vy vz mass)
(obj x x y y z z vx vx vy vy vz vz mass mass))
(= jupiter* (make-body 4.84143144246472090d0
-1.16032004402742839d0
-1.03622044471123109d-1
(* 1.66007664274403694d-3 days-per-year*)
(* 7.69901118419740425d-3 days-per-year*)
(* -6.90460016972063023d-5 days-per-year*)
(* 9.54791938424326609d-4 solar-mass*)))
(= saturn* (make-body 8.34336671824457987d0
4.12479856412430479d0
-4.03523417114321381d-1
(* -2.76742510726862411d-3 days-per-year*)
(* 4.99852801234917238d-3 days-per-year*)
(* 2.30417297573763929d-5 days-per-year*)
(* 2.85885980666130812d-4 solar-mass*)))
(= uranus* (make-body 1.28943695621391310d1
-1.51111514016986312d1
-2.23307578892655734d-1
(* 2.96460137564761618d-03 days-per-year*)
(* 2.37847173959480950d-03 days-per-year*)
(* -2.96589568540237556d-05 days-per-year*)
(* 4.36624404335156298d-05 solar-mass*)))
(= neptune* (make-body 1.53796971148509165d+01
-2.59193146099879641d+01
1.79258772950371181d-01
(* 2.68067772490389322d-03 days-per-year*)
(* 1.62824170038242295d-03 days-per-year*)
(* -9.51592254519715870d-05 days-per-year*)
(* 5.15138902046611451d-05 solar-mass*)))
(= sun* (make-body 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 solar-mass*))
(def apply-forces (a b dt)
(withs (dx (- (a 'x) (b 'x))
dy (- (a 'y) (b 'y))
dz (- (a 'z) (b 'z))
d (sqrt (+ (expt dx 2) (expt dy 2) (expt dz 2)))
mag (/ dt (expt d 3))
dxmag (* dx mag)
dymag (* dy mag)
dzmag (* dz mag))
(-- (a 'vx) (* dxmag (b 'mass)))
(-- (a 'vy) (* dymag (b 'mass)))
(-- (a 'vz) (* dzmag (b 'mass)))
(++ (b 'vx) (* dxmag (a 'mass)))
(++ (b 'vy) (* dymag (a 'mass)))
(++ (b 'vz) (* dzmag (a 'mass)))))
(def advance (solar-system dt)
(on a solar-system
(for i (+ index 1) (- (len solar-system) 1)
(apply-forces a (solar-system i) dt)))
(each b solar-system
(++ (b 'x) (* dt (b 'vx)))
(++ (b 'y) (* dt (b 'vy)))
(++ (b 'z) (* dt (b 'vz)))))
(def energy (solar-system)
(let e 0.0d0
(on a solar-system
(++ e (* 0.5d0
(a 'mass)
(+ (expt (a 'vx) 2)
(expt (a 'vy) 2)
(expt (a 'vz) 2))))
(for i (+ index 1) (- (len solar-system) 1)
(withs (b (solar-system i)
dx (- (a 'x) (b 'x))
dy (- (a 'y) (b 'y))
dz (- (a 'z) (b 'z))
d (sqrt (+ (expt dx 2) (expt dy 2) (expt dz 2))))
(-- e (/ (* (a 'mass) (b 'mass)) d)))))
e))
(def offset-momentum (solar-system)
(with (px 0.0d0
py 0.0d0
pz 0.0d0)
(each p solar-system
(++ px (* (p 'vx) (p 'mass)))
(++ py (* (p 'vy) (p 'mass)))
(++ pz (* (p 'vz) (p 'mass))))
(= ((car solar-system) 'vx) (/ (- px) solar-mass*)
((car solar-system) 'vy) (/ (- py) solar-mass*)
((car solar-system) 'vz) (/ (- pz) solar-mass*))))
(def nbody (n)
(let solar-system (list sun* jupiter* saturn* uranus* neptune*)
(offset-momentum solar-system)
(prn (energy solar-system))
(repeat n (advance solar-system 0.01d0))
(prn (energy solar-system))
nil))
I didn't have the patience to wait for Arc to finish the benchmark's n = 50,000,000 (gives you an idea of how well it performs!), but on smaller values it seems correct.
The reference implementation, run from a GNU clisp REPL (though I suppose I should use SBCL for a fairer time comparison, I'm just interested in the accuracy of the Arc results with this run):
> (time (nbody 50000))
-0.169075164
-0.169078071
Real time: 37.575706 sec.
Run time: 32.862053 sec.
Space: 374412048 Bytes
GC: 477, GC time: 1.836092 sec.
NIL
