N-Body Utilities Coming Real Soon Now (TM)

It's been a while since I posted here, but never fear: I've been busy. I've been working on an updated set of N-Body utilities for another paper examining the efficiency of our algorithm (see astro-ph/0611416, to be published in ApJ 663 (2007) 1420) relative to drift-kick-drift leapfrog in a treecode. I've got code up and running, and am beginning simulations to compare the two algorithms. The code is a mix of C (for anything O(N) and larger) and Scheme (to patch the various C algorithms together into actual simulations).

It looks, preliminarily, like our algorithm significantly outperforms leapfrog (i.e. tree-computed force errors dominate at larger timesteps with our algorithm) in this low-accuracy regime, which is very exciting. Right now Ed has an idea for quantifying the errors introduced into this process in terms of the change in the velocity diffusion coefficient in the Fokker-Planck equation when the accuracy parameter of the treecode is varied, so I'm writing some code to run simulations and compute the velocity diffusion coefficient.

Once I've got the code into a workable state (i.e. once it easily runs all the simulations I've got set up for it), I'll write up some documentation and post the whole mess on PlaneT. N-Body simulations for all!

I just read on this PLT Scheme mailing list thread about a nifty way to post syntax-highlighted scheme code to your blog, so I have to try it out. Here's some code which runs a 10K-body Plummer model up to core collapse (at which point, you'll have to use a code which models binary systems with small separations carefully, which this code does not).

(module up-to-core-collapse mzscheme
  (require (planet "bh-tree.ss" ("wmfarr" "nbody.plt" 1 0))
           (planet "ics.ss" ("wmfarr" "nbody.plt" 1 0))
           (planet "timestep.ss" ("wmfarr" "nbody.plt" 1 0))
           (planet "evolve.ss" ("wmfarr" "nbody.plt" 1 0))
           (planet "system.ss" ("wmfarr" "nbody.plt" 1 0));
           (planet "gl3pt.ss" ("wmfarr" "nbody.plt" 1 0)))
  
  
  ;; The bodies.
  (define bs (make-plummer* 10000))
  
  ;; Various parameters.
  (define T (* 20 (median-relaxation-time bs))) ; Core collapse within 20 median relax times.
  (define eps 1e-3) ; Force softening: don't resolve scales < 1e-3
  (define energy-beta 0.1)
  (define advance-beta 1.0)
  (define output-interval 1.0)
  
  ;; Total initial energy
  (define e (bh-tree-total-energy eps energy-beta bs))
  
  ;; Timestep function.
  (define frac 0.5)
  (define timestep (max-force-timestep eps frac))
  
  ;; Single-step advancer.
  (define adv (lambda (bs dt) (gl3pt-advance eps advance-beta bs dt)))
  
  ;; Filter between steps---output roughly at output-interval.
  (define last-output (make-parameter 0))
  (define t (make-parameter 0))
  (define (filter bs dt)
    (t (+ (t) dt))
    (when (> (t) (+ (last-output) output-interval))
      (last-output (t))
      (let ((e-new (bh-tree-total-energy eps energy-beta bs)))
        (printf (string-append "At time ~a, energy is ~a, 10% lagrangian "
                               "radius is ~a, energy error is ~a~%")
                (t) e-new (lagrangian-radius bs 0.1) (abs (/ (- e e-new) e))))))
  
  ;; Do it.
  ((evolve adv timestep filter) bs T))