Permalink
Browse files

Particles initialized to realistic values

  • Loading branch information...
marblestation committed Feb 14, 2017
1 parent 7c641cd commit 055d752eae75c3d64b6d75629c6da28c1aa5b7d8
Showing with 52 additions and 5 deletions.
  1. +17 −0 README.md
  2. +9 −3 c/leapfrog.c
  3. +8 −0 fortran/leapfrog.f90
  4. +8 −0 go/leapfrog.go
  5. +10 −2 rust/leapfrog.rs
View
@@ -2,6 +2,23 @@
Implementation in C, Fortran, Go and Rust of a very simple N-Body simulator with 3 particles using a LeapFrog integrator. Presented in [What can the programming language Rust do for astrophysics?](https://arxiv.org/abs/1702.02951), to appear in the Proceedings of the IAU Symposium 325 on Astroinformatics.
The code has evolved since publication, implementing several suggestion made by the [Hacker News](https://news.ycombinator.com/item?id=13632894) and [Rust subreddit](https://www.reddit.com/r/rust/comments/5trref/what_can_rust_do_for_astrophysics/) communities. The current times on a 1,6 GHz Intel Core i5 machine:
- C: 2m53.504s
- Fortran: 3m14.648s
- Rust: 2m33.082s
- Go: 4m6.579s
Output positions for the two particles after a one million year simulation:
```
[
[0.00011259641961937496, 0.00011235312238407964, 0.00000982962683070781],
[-2.986377308237493, 14588403.613911422, 1276320.0322206458]
]
```
The original article's conclusions are still valid, Rust can be as fast and precise as C or Fortran. Additionally, Rust characteristics ensures that scientific results are not affected by memory management issues and this is a great advantage for reliable scientific computation.
## Compilation & execution
View
@@ -65,6 +65,14 @@ int main(int argc, char* argv[]) {
a[i][1] = 0;
a[i][2] = 0;
}
m[0] = 0.08; // M_SUN
m[1] = 3.0e-6; // M_SUN
x[1][0] = 0.0162; // AU
x[1][1] = 6.57192058353e-15; // AU
x[1][2] = 5.74968548652e-16; // AU
v[1][0] = -1.48427302304e-14;
v[1][1] = 0.0399408809121;
v[1][2] = 0.00349437429104;
while(time <= time_limit) {
integrator_leapfrog_part1(n_particles, x, v, half_time_step);
@@ -73,8 +81,6 @@ int main(int argc, char* argv[]) {
integrator_leapfrog_part2(n_particles, x, v, a, time_step, half_time_step);
time += half_time_step;
}
printf("Position particle 0: %f %f %f\n", x[0][0], x[0][1], x[0][2]);
printf("Position particle 1: %f %f %f\n", x[1][0], x[1][1], x[1][2]);
printf("Position particle 2: %f %f %f\n", x[2][0], x[2][1], x[2][2]);
printf("Position 1: %f %f %f | Position 2: %f %f %f\n", x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2]);
}
View
@@ -27,6 +27,14 @@ program leapfrog
a(1,1:n_particles) = 0.0d0
a(2,1:n_particles) = 0.0d0
a(3,1:n_particles) = 0.0d0
m(1) = 0.08 ! M_SUN
m(2) = 3.0e-6 ! M_SUN
x(1, 2) = 0.0162 ! AU
x(2, 2) = 6.57192058353e-15 ! AU
x(3, 2) = 5.74968548652e-16 ! AU
v(1, 2) = -1.48427302304e-14
v(2, 2) = 0.0399408809121
v(3, 2) = 0.00349437429104
half_time_step = 0.5d0*time_step
do while (time.le.time_limit)
View
@@ -23,6 +23,14 @@ func main() {
v := make([][3]float64, n_particles)
a := make([][3]float64, n_particles)
m := make([]float64, n_particles)
m[0] = 0.08 // M_SUN
m[1] = 3.0e-6 // M_SUN
x[1][0] = 0.0162 // AU
x[1][1] = 6.57192058353e-15 // AU
x[1][2] = 5.74968548652e-16 // AU
v[1][0] = -1.48427302304e-14
v[1][1] = 0.0399408809121
v[1][2] = 0.00349437429104
for time <= time_limit {
integrator_leapfrog_part1(n_particles, x, v, half_time_step)
View
@@ -8,7 +8,15 @@ fn main() {
let mut x: [[f64; 3]; N_PARTICLES] = [[0.; 3]; N_PARTICLES];
let mut v: [[f64; 3]; N_PARTICLES] = [[0.; 3]; N_PARTICLES];
let mut a: [[f64; 3]; N_PARTICLES] = [[0.; 3]; N_PARTICLES];
let m: [f64; 3] = [0.; 3];
let mut m: [f64; N_PARTICLES] = [0.; N_PARTICLES];
m[0] = 0.08; // M_SUN
m[1] = 3.0e-6; // M_SUN
x[1][0] = 0.0162; // AU
x[1][1] = 6.57192058353e-15; // AU
x[1][2] = 5.74968548652e-16; // AU
v[1][0] = -1.48427302304e-14;
v[1][1] = 0.0399408809121;
v[1][2] = 0.00349437429104;
while time <= time_limit {
integrator_leapfrog_part1(N_PARTICLES, &mut x, &v, half_time_step);
@@ -40,7 +48,7 @@ fn integrator_leapfrog_part2(n_particles: usize, x: &mut [[f64; 3]; N_PARTICLES]
}
}
fn gravity_calculate_acceleration(n_particles: usize, m: &[f64; 3], x: &[[f64; 3]; N_PARTICLES], a: &mut [[f64; 3]; N_PARTICLES]) {
fn gravity_calculate_acceleration(n_particles: usize, m: &[f64; N_PARTICLES], x: &[[f64; 3]; N_PARTICLES], a: &mut [[f64; 3]; N_PARTICLES]) {
let g = 6.6742367e-11; // m^3.kg^-1.s^-2
for i in 0..n_particles {
a[i][0] = 0.;

0 comments on commit 055d752

Please sign in to comment.