C and Rust Leapfrog Comparison

leapfrog.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>


typedef struct {
    double x, y, z;
} Point;


typedef struct {
    double mass;
    Point position, velocity, acceleration;
} Particle;


const double G = 6.6742367e-11;

void integrator_leapfrog_part1(int n_particles, Particle* particles, double half_time_step){
    for (int i = 0; i < n_particles; i++) {
        particles[i].position.x += half_time_step * particles[i].velocity.x;
        particles[i].position.y += half_time_step * particles[i].velocity.y;
        particles[i].position.z += half_time_step * particles[i].velocity.z;
    }
}


void integrator_leapfrog_part2(int n_particles, Particle* particles, double time_step, double half_time_step){
    for (int i = 0; i < n_particles; i++) {
        particles[i].velocity.x += time_step * particles[i].acceleration.x;
        particles[i].velocity.y += time_step * particles[i].acceleration.y;
        particles[i].velocity.z += time_step * particles[i].acceleration.z;

        particles[i].position.x += half_time_step * particles[i].velocity.x;
        particles[i].position.y += half_time_step * particles[i].velocity.y;
        particles[i].position.z += half_time_step * particles[i].velocity.z;
    }
}


void gravity_calculate_acceleration(int n_particles, Particle* particles) {
    for (int i = 0; i < n_particles; i++) {
        particles[i].acceleration.x = 0.0;
        particles[i].acceleration.y = 0.0;
        particles[i].acceleration.z = 0.0;
        for (int j = 0; j < n_particles; j++) {
            if (j == i) {
                continue;
            }
            double dx = particles[i].position.x - particles[j].position.x;
            double dy = particles[i].position.y - particles[j].position.y;
            double dz = particles[i].position.z - particles[j].position.z;
            double r = sqrt(dx*dx + dy*dy + dz*dz);
            double prefact = -G/(r*r*r) * particles[j].mass;
            particles[i].acceleration.x += prefact * dx;
            particles[i].acceleration.y += prefact * dy;
            particles[i].acceleration.z += prefact * dz;
        }
    }
}


int main() {

    double time = 0.0;
    double time_step = 0.08;
    double half_time_step = time_step/2;
    double time_limit = 365.25 * 1e4;

    Particle star = {
            .mass = 0.08,
            .position = {.x = 0.0, .y = 0.0, .z = 0.0},
            .velocity = {.x =  0.0, .y = 0.0, .z = 0.0},
            .acceleration = {.x = 0.0, .y = 0.0, .z = 0.0}
    };
    Particle planet = {
            .mass = 3.0e-6,
            .position = {.x = 0.0162, .y = 6.57192058353e-15, .z = 5.74968548652e-16},
            .velocity = {.x = -1.48427302304e-14, .y = 0.0399408809121, .z = 0.00349437429104},
            .acceleration = {.x = 0.0, .y = 0.0, .z = 0.0},
    };

    int num_particles = 2;
    Particle* particles = malloc(2*sizeof(Particle));
    particles[0] = star;
    particles[1] = planet;

    while (time < time_limit) {
        integrator_leapfrog_part1(num_particles, particles, half_time_step);
        time += half_time_step;
        gravity_calculate_acceleration(num_particles, particles);
        integrator_leapfrog_part2(num_particles, particles, time_step, half_time_step);
        time += half_time_step;
    }

    printf("%f", particles[1].position.x);
    return 0;
}

leapfrog.rs

// Originally 1.24 seconds, using compile-time arrays
// New is 1.38 seconds


struct Point {
    x: f64,
    y: f64,
    z: f64,
}


struct Particle {
    mass: f64,
    position: Point,
    velocity: Point,
    acceleration: Point,
}


fn main() {
    let mut time: f64 = 0.;
    let time_step: f64 = 0.08;
    let half_time_step: f64 = time_step/2.;
    let time_limit: f64 = 365.25 * 1e4;

    let star = Particle {
        mass: 0.08,
        position: Point {x: 0.0, y: 0.0, z: 0.0},
        velocity: Point {x: 0.0, y: 0.0, z: 0.0},
        acceleration: Point {x: 0.0, y: 0.0, z: 0.0}
    };
    let planet = Particle {
        mass: 3.0e-6,
        position: Point {x: 0.0162, y: 6.57192058353e-15, z: 5.74968548652e-16},
        velocity: Point {x: -1.48427302304e-14, y: 0.0399408809121, z: 0.00349437429104},
        acceleration: Point {x: 0.0, y: 0.0, z: 0.0},
    };
    let mut particles = vec![star, planet];

    while time < time_limit {
        integrator_leapfrog_part1(&mut particles, half_time_step);
        time += half_time_step;
        gravity_calculate_acceleration(&mut particles);
        integrator_leapfrog_part2(&mut particles, time_step, half_time_step);
        time += half_time_step;
    }
    println!("{:?}", particles[1].position.x);
}


fn integrator_leapfrog_part1(particles: &mut Vec<Particle>, half_time_step: f64) {
    for ref mut particle in particles {
        particle.position.x += half_time_step * particle.velocity.x;
        particle.position.y += half_time_step * particle.velocity.y;
        particle.position.z += half_time_step * particle.velocity.z;
    }
}


fn integrator_leapfrog_part2(particles: &mut Vec<Particle>, time_step: f64, half_time_step: f64) {
    for ref mut particle in particles {
        particle.velocity.x += time_step * particle.acceleration.x;
        particle.velocity.y += time_step * particle.acceleration.y;
        particle.velocity.z += time_step * particle.acceleration.z;

        particle.position.x += half_time_step * particle.velocity.x;
        particle.position.y += half_time_step * particle.velocity.y;
        particle.position.z += half_time_step * particle.velocity.z;
    }
}


fn gravity_calculate_acceleration(particles: &mut Vec<Particle>) {
    let g = 6.6742367e-11; // m^3.kg^-1.s^-2

    for i in 0..particles.len() {
        let mut acceleration = Point {x: 0.0, y: 0.0, z: 0.0};
        for (j, other) in particles.iter().enumerate() {
            if i == j {
                continue;
            }
            let dx = particles[i].position.x - other.position.x;
            let dy = particles[i].position.y - other.position.y;
            let dz = particles[i].position.z - other.position.z;
            let r = (dx*dx + dy*dy + dz*dz).sqrt();
            let prefact = -g/(r*r*r) * other.mass;

            acceleration.x += prefact * dx;
            acceleration.y += prefact * dy;
            acceleration.z += prefact * dz;
        }
        particles[i].acceleration = acceleration;
    }
}