The day the Earth stood still

README.md

The day the Earth stood still

The source code in this gist go together with this mind bytes post.

The two additional files are

• suncrash.py: Ultimate implementation of the problem in python
• main.rs: Same implementation as suncrash.py in Rust

main.rs

use std::time::Instant;

const GRAV: f64 = 6.67408e-11;
const M_SUN: f64 = 1.989e30;
const R_START: f64 = 1.495979e11;

fn main() {
    let ds = 1000.;

    let now = Instant::now();
    let t_total = calc_loop(ds);
    let elapsed = now.elapsed();
    println!(
        "Time until Earth crashes into the Sun: {:.2?} days",
        t_total / 24. / 3600.
    );
    println!("Computation time: {:.2?}.", elapsed);
}

fn calc_loop(ds: f64) -> f64 {
    let mut t_total: f64 = 0.; // total time
    let mut v_curr: f64 = 0.; // initial velocity
    let mut r_curr = R_START;

    let mut a_curr: f64;
    let mut t_curr: f64;

    while r_curr > 0. {
        a_curr = GRAV * M_SUN / r_curr.powf(2.);
        t_curr = (-v_curr + (v_curr.powf(2.) + 2. * a_curr * ds).sqrt()) / a_curr;

        v_curr += a_curr * t_curr;
        r_curr -= ds;
        t_total += t_curr;
    }
    t_total
}

suncrash.py

import time

from numba import njit
import numpy as np

# Constants
GRAV = 6.67408e-11 # m^3 kg^-1 s^-2
M_SUN = 1.98847e30 # kg
R_START = 1.495979e11 # m

# Step size
dr = 1e3  # m

@njit
def total_time(dr):
    # Initial conditions
    t_total = 0  # s: total time
    v_curr = 0  # m/s: initial velocity
    r_curr = R_START  # m: initial distance

    # Calculate time until Earth crashes into the Sun
    while r_curr > 0:
        # Calculate acceleration
        a_curr = GRAV * M_SUN / r_curr**2

        # Calculate time it takes to travel dr
        t_curr = (-v_curr + np.sqrt(v_curr**2 + 2 * a_curr * dr)) / a_curr

        # Update total time
        t_total += t_curr

        # Update velocity
        v_curr += a_curr * t_curr

        # Update distance
        r_curr -= dr
    
    return t_total

# Print result

_ = total_time(dr)  # Compile function
tic = time.time()
t_total = total_time(dr)
toc = time.time()
print(f"Time until Earth crashes into the Sun: {t_total / 3600 / 24:.2f} days")
print(f"Computation time: {toc - tic:.2f} seconds")