1
#!/usr/bin/env gsi-script
;; gsc-nbody.scm
;; NBODY exercise. Gambit Scheme version
;; compile: $ gsc -exe gsc-nbody
;;     run: $ ./gsc-nbody <n-iterations>
;;

(declare (standard-bindings)
         #;(extended-bindings)
         (not safe)
         #;(not interrupts-enabled)
         )

(define PI (* 4.0 (atan 1)))
(define SOLAR-MASS (* 4 PI PI))
(define DAYS-PER-YEAR 365.24)
(define-syntax dpe (syntax-rules () ((_ v) (fl* v DAYS-PER-YEAR))))
(define-syntax slr (syntax-rules () ((_ m) (fl* m SOLAR-MASS))))

(define-syntax dotimes (syntax-rules () ((_ n body ...) (do ((i_ n (fx- i_ 1))) ((fxzero? i_)) body ...))))
(define-syntax subn! (syntax-rules () ((_ what withs) (set! what (fl- what withs)))))
(define-syntax addn! (syntax-rules () ((_ what withs) (set! what (fl+ what withs)))))

(define (take-energy n)
  (define bodies
    (f64vector
     SOLAR-MASS 0.0 0.0 0.0 0.0 0.0 0.0  ;; SUN
     (slr 9.54791938424326609e-04) 4.84143144246472090e+00 -1.16032004402742839e+00 -1.03622044471123109e-01 (dpe 1.66007664274403694e-03) (dpe 7.69901118419740425e-03) (dpe -6.90460016972063023e-05)
     (slr 2.85885980666130812e-04) 8.34336671824457987e+00 4.12479856412430479e+00 -4.03523417114321381e-01 (dpe -2.76742510726862411e-03) (dpe 4.99852801234917238e-03) (dpe 2.30417297573763929e-05)
     (slr 4.36624404335156298e-05) 1.28943695621391310e+01 -1.51111514016986312e+01 -2.23307578892655734e-01 (dpe 2.96460137564761618e-03) (dpe 2.37847173959480950e-03) (dpe -2.96589568540237556e-05)
     (slr 5.15138902046611451e-05) 1.53796971148509165e+01 -2.59193146099879641e+01 1.79258772950371181e-01 (dpe 2.68067772490389322e-03) (dpe 1.62824170038242295e-03) (dpe -9.51592254519715870e-05)
     ))
  (define body-size 7)
  (define bodies-len (f64vector-length bodies))    ;; (2)
  (define dt 0.01)

  (define-syntax let-attr
    (syntax-rules () ((_ ((a ix offset) ...) body ...)
                      (let ((a (f64vector-ref bodies (fx+ ix offset))) ...)
                        body ...))))
  (define (subnv! ix offset val deltaval) (f64vector-set! bodies (fx+ ix offset) (fl- val deltaval)))
  (define (addnv! ix offset val deltaval) (f64vector-set! bodies (fx+ ix offset) (fl+ val deltaval)))
  (define-syntax ^2 (syntax-rules () ((_ x) (flsquare x))))
  (define-syntax sum^2 (syntax-rules () ((_ x y z) (fl+ (^2 x) (^2 y) (^2 z)))))
  (define (dsqt v) (fl* v (flsqrt v)))
  (define E 0.0)

  (dotimes
   n
   (do ((i 0 (fx+ i body-size)))
       ((fx= i bodies-len))
     (do ((j (fx+ i body-size) (fx+ j body-size)))
         ((fx= j bodies-len))
       (let-attr
        ((massi i 0) (Xi i 1) (Yi i 2) (Zi i 3)
         (massj j 0) (Xj j 1) (Yj j 2) (Zj j 3)
         (VXi i 4) (VYi i 5) (VZi i 6)
         (VXj j 4) (VYj j 5) (VZj j 6))
        (let* ((dx (fl- Xi Xj))
               (dy (fl- Yi Yj))
               (dz (fl- Zi Zj))
               (mag (fl/ dt (dsqt (sum^2 dx dy dz)))))

          (subnv! i 4 VXi (fl* dx massj mag))
          (subnv! i 5 VYi (fl* dy massj mag))
          (subnv! i 6 VZi (fl* dz massj mag))
          (addnv! j 4 VXj (fl* dx massi mag))
          (addnv! j 5 VYj (fl* dy massi mag))
          (addnv! j 6 VZj (fl* dz massi mag))
        ))))
   (do ((i 0 (fx+ i body-size)))
       ((fx= i bodies-len))
     (let-attr
      ((X i 1) (Y i 2) (Z i 3)
       (VX i 4) (VY i 5) (VZ i 6))
      (addnv! i 1 X (fl* dt VX))
      (addnv! i 2 Y (fl* dt VY))
      (addnv! i 3 Z (fl* dt VZ)))))

  (do ((i 0 (fx+ i body-size)))
      ((fx= i bodies-len) E)
    (let-attr
     ((massi i 0)
      (VX i 4)
      (VY i 5)
      (VZ i 6))
     (addn! E (fl* 0.5 massi (sum^2 VX VY VZ)))
     (do ((j (fx+ i body-size) (fx+ body-size j)))
         ((fx= j bodies-len))
       (let-attr
        ((Xi i 1) (Yi i 2) (Zi i 3)
         (Xj j 1) (Yj j 2) (Zj j 3) (massj j 0))
        (subn! E (fl/ (fl* massi massj)
                      (flsqrt (sum^2
                               (fl- Xi Xj)
                               (fl- Yi Yj)
                               (fl- Zi Zj))))))
       )))
  )

(define (main . args)
  (cond ((pair? args)
         (let* ((n (string->number (car args)))
                (st (current-second))
                (e (take-energy n))
                (et (current-second)))
           (print e " per " (- et st) " sec on " n "\n")))
        (else (print "Usage: gsc-nbody <n-iterations>~%"))))

For immediate assistance, please email our customer support: [email protected]

Download RAW File