nbody.rs

/* The Computer Language Benchmarks Game

   http://benchmarksgame.alioth.debian.org/



   contributed by Mark C. Lewis

   modified slightly by Chad Whipkey

   converted from java to rust by Timofey Prodanov

*/

use std::iter::AdditiveIterator;
use std::os;

struct Body {
   pub x:      f64,
   pub y:      f64,
   pub z:      f64,
   pub vx:      f64,
   pub vy:      f64,
   pub vz:      f64,
   pub mass:   f64
}

static pi: f64 = 3.141592653589793;
static solar_mass: f64 = 4. * pi * pi;
static days_per_year: f64 = 365.24;
static planets_count: uint = 5;

impl Body {

   pub fn bodies() -> [Body, ..planets_count] {
      [
         // sun

         Body {
            x: 0.,
            y: 0.,
            z: 0.,
            vx: 0.,
            vy: 0.,
            vz: 0.,
            mass: solar_mass
         },
         // jupiter

         Body {
            x: 4.84143144246472090e+00,
            y: -1.16032004402742839e+00,
            z: -1.03622044471123109e-01,
            vx: 1.66007664274403694e-03 * days_per_year,
            vy: 7.69901118419740425e-03 * days_per_year,
            vz: -6.90460016972063023e-05 * days_per_year,
            mass: 9.54791938424326609e-04 * solar_mass
         },
         // saturn

         Body {
            x: 8.34336671824457987e+00,
            y: 4.12479856412430479e+00,
            z: -4.03523417114321381e-01,
            vx: -2.76742510726862411e-03 * days_per_year,
            vy: 4.99852801234917238e-03 * days_per_year,
            vz: 2.30417297573763929e-05 * days_per_year,
            mass: 2.85885980666130812e-04 * solar_mass
         },
         // uranus

         Body {
            x: 1.28943695621391310e+01,
            y: -1.51111514016986312e+01,
            z: -2.23307578892655734e-01,
            vx: 2.96460137564761618e-03 * days_per_year,
            vy: 2.37847173959480950e-03 * days_per_year,
            vz: -2.96589568540237556e-05 * days_per_year,
            mass: 4.36624404335156298e-05 * solar_mass
         },
         // neptune

         Body {
            x: 1.53796971148509165e+01,
            y: -2.59193146099879641e+01,
            z: 1.79258772950371181e-01,
            vx: 2.68067772490389322e-03 * days_per_year,
            vy: 1.62824170038242295e-03 * days_per_year,
            vz: -9.51592254519715870e-05 * days_per_year,
            mass: 5.15138902046611451e-05 * solar_mass
         }
      ]
   }
   
   fn offset_momentum(&mut self, px: f64, py: f64, pz: f64) {
      self.vx = -px / solar_mass;
      self.vy = -py / solar_mass;
      self.vz = -pz / solar_mass;
   }

}

struct NBodySystem {
   bodies: [Body, ..planets_count]
}

impl NBodySystem {
   
   pub fn new() -> NBodySystem {
      let mut s = NBodySystem {
         bodies: Body::bodies()
      };
      
      let px = s.bodies.iter().map(|&x| x.vx * x.mass).sum();
      let py = s.bodies.iter().map(|&x| x.vy * x.mass).sum();
      let pz = s.bodies.iter().map(|&x| x.vz * x.mass).sum();
      s.bodies[0].offset_momentum(px, py, pz);
      s
   }
   
   pub fn advance(&mut self, dt: f64) {
      for i in range(0, planets_count - 1) {
         for j in range(i + 1, planets_count) {
            let dx = self.bodies[i].x - self.bodies[j].x;
            let dy = self.bodies[i].y - self.bodies[j].y;
            let dz = self.bodies[i].z - self.bodies[j].z;
            
            let dSquared = dx * dx + dy * dy + dz * dz;
            let dist = dSquared.sqrt();
            let mag = dt / (dSquared * dist);
            
            self.bodies[i].vx -= dx * self.bodies[j].mass * mag;
            self.bodies[i].vy -= dy * self.bodies[j].mass * mag;
            self.bodies[i].vz -= dz * self.bodies[j].mass * mag;
            
            self.bodies[j].vx += dx * self.bodies[i].mass * mag;
            self.bodies[j].vy += dy * self.bodies[i].mass * mag;
            self.bodies[j].vz += dz * self.bodies[i].mass * mag;
         }
      }
      
      for b in self.bodies.mut_iter() {
         b.x += dt * b.vx;
         b.y += dt * b.vy;
         b.z += dt * b.vz;
      }
   }
   
   pub fn energy(&self) -> f64 {
      let mut e = 0f64;
      for i in range(0, planets_count) {
         e += 0.5 * self.bodies[i].mass * (
            self.bodies[i].vx * self.bodies[i].vx +
            self.bodies[i].vy * self.bodies[i].vy +
            self.bodies[i].vz * self.bodies[i].vz);
         
         for j in range(i + 1, planets_count) {
            let dx = self.bodies[i].x - self.bodies[j].x;
            let dy = self.bodies[i].y - self.bodies[j].y;
            let dz = self.bodies[i].z - self.bodies[j].z;
            
            let dist = (dx * dx + dy * dy + dz * dz).sqrt();
            e -= (self.bodies[i].mass * self.bodies[j].mass) / dist;
         }
      }
      e
   }
   
}

fn main() {
   let n: int = from_str(os::args().get(1).as_slice()).unwrap();
   let mut nbs = NBodySystem::new();
   println!("{:.9f}", nbs.energy());
   for _ in range(0, n) {
      nbs.advance(0.01);
   }
   println!("{:.9f}", nbs.energy());
}