src/test/bench/shootout-nbody.rs

   1 // Copyright 2014 The Rust Project Developers. See the COPYRIGHT
   2 // file at the top-level directory of this distribution and at
   3 // http://rust-lang.org/COPYRIGHT.
   4 //
   5 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
   6 // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
   7 // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
   8 // option. This file may not be copied, modified, or distributed
   9 // except according to those terms.
  10
  11 static PI: f64 = 3.141592653589793;
  12 static SOLAR_MASS: f64 = 4.0 * PI * PI;
  13 static YEAR: f64 = 365.24;
  14 static N_BODIES: uint = 5;
  15
  16 static BODIES: [Planet, ..N_BODIES] = [
  17     // Sun
  18     Planet {
  19         x: 0.0, y: 0.0, z: 0.0,
  20         vx: 0.0, vy: 0.0, vz: 0.0,
  21         mass: SOLAR_MASS,
  22     },
  23     // Jupiter
  24     Planet {
  25         x: 4.84143144246472090e+00,
  26         y: -1.16032004402742839e+00,
  27         z: -1.03622044471123109e-01,
  28         vx: 1.66007664274403694e-03 * YEAR,
  29         vy: 7.69901118419740425e-03 * YEAR,
  30         vz: -6.90460016972063023e-05 * YEAR,
  31         mass: 9.54791938424326609e-04 * SOLAR_MASS,
  32     },
  33     // Saturn
  34     Planet {
  35         x: 8.34336671824457987e+00,
  36         y: 4.12479856412430479e+00,
  37         z: -4.03523417114321381e-01,
  38         vx: -2.76742510726862411e-03 * YEAR,
  39         vy: 4.99852801234917238e-03 * YEAR,
  40         vz: 2.30417297573763929e-05 * YEAR,
  41         mass: 2.85885980666130812e-04 * SOLAR_MASS,
  42     },
  43     // Uranus
  44     Planet {
  45         x: 1.28943695621391310e+01,
  46         y: -1.51111514016986312e+01,
  47         z: -2.23307578892655734e-01,
  48         vx: 2.96460137564761618e-03 * YEAR,
  49         vy: 2.37847173959480950e-03 * YEAR,
  50         vz: -2.96589568540237556e-05 * YEAR,
  51         mass: 4.36624404335156298e-05 * SOLAR_MASS,
  52     },
  53     // Neptune
  54     Planet {
  55         x: 1.53796971148509165e+01,
  56         y: -2.59193146099879641e+01,
  57         z: 1.79258772950371181e-01,
  58         vx: 2.68067772490389322e-03 * YEAR,
  59         vy: 1.62824170038242295e-03 * YEAR,
  60         vz: -9.51592254519715870e-05 * YEAR,
  61         mass: 5.15138902046611451e-05 * SOLAR_MASS,
  62     },
  63 ];
  64
  65 struct Planet {
  66     x: f64, y: f64, z: f64,
  67     vx: f64, vy: f64, vz: f64,
  68     mass: f64,
  69 }
  70
  71 fn advance(bodies: &mut [Planet, ..N_BODIES], dt: f64, steps: int) {
  72     for _ in range(0, steps) {
  73         let mut b_slice = bodies.as_mut_slice();
  74         loop {
  75             let bi = match b_slice.mut_shift_ref() {
  76                 Some(bi) => bi,
  77                 None => break
  78             };
  79             for bj in b_slice.mut_iter() {
  80                 let dx = bi.x - bj.x;
  81                 let dy = bi.y - bj.y;
  82                 let dz = bi.z - bj.z;
  83
  84                 let d2 = dx * dx + dy * dy + dz * dz;
  85                 let mag = dt / (d2 * d2.sqrt());
  86
  87                 let massj_mag = bj.mass * mag;
  88                 bi.vx -= dx * massj_mag;
  89                 bi.vy -= dy * massj_mag;
  90                 bi.vz -= dz * massj_mag;
  91
  92                 let massi_mag = bi.mass * mag;
  93                 bj.vx += dx * massi_mag;
  94                 bj.vy += dy * massi_mag;
  95                 bj.vz += dz * massi_mag;
  96             }
  97             bi.x += dt * bi.vx;
  98             bi.y += dt * bi.vy;
  99             bi.z += dt * bi.vz;
 100         }
 101     }
 102 }
 103
 104 fn energy(bodies: &[Planet, ..N_BODIES]) -> f64 {
 105     let mut e = 0.0;
 106     let mut bodies = bodies.as_slice();
 107     loop {
 108         let bi = match bodies.shift_ref() {
 109             Some(bi) => bi,
 110             None => break
 111         };
 112         e += (bi.vx * bi.vx + bi.vy * bi.vy + bi.vz * bi.vz) * bi.mass / 2.0;
 113         for bj in bodies.iter() {
 114             let dx = bi.x - bj.x;
 115             let dy = bi.y - bj.y;
 116             let dz = bi.z - bj.z;
 117             let dist = (dx * dx + dy * dy + dz * dz).sqrt();
 118             e -= bi.mass * bj.mass / dist;
 119         }
 120     }
 121     e
 122 }
 123
 124 fn offset_momentum(bodies: &mut [Planet, ..N_BODIES]) {
 125     let mut px = 0.0;
 126     let mut py = 0.0;
 127     let mut pz = 0.0;
 128     for bi in bodies.iter() {
 129         px += bi.vx * bi.mass;
 130         py += bi.vy * bi.mass;
 131         pz += bi.vz * bi.mass;
 132     }
 133     let sun = &mut bodies[0];
 134     sun.vx = - px / SOLAR_MASS;
 135     sun.vy = - py / SOLAR_MASS;
 136     sun.vz = - pz / SOLAR_MASS;
 137 }
 138
 139 fn main() {
 140     let n = if std::os::getenv("RUST_BENCH").is_some() {
 141         5000000
 142     } else {
 143         std::os::args().as_slice().get(1)
 144             .and_then(|arg| from_str(arg.as_slice()))
 145             .unwrap_or(1000)
 146     };
 147     let mut bodies = BODIES;
 148
 149     offset_momentum(&mut bodies);
 150     println!("{:.9f}", energy(&bodies));
 151
 152     advance(&mut bodies, 0.01, n);
 153
 154     println!("{:.9f}", energy(&bodies));
 155 }