Simulate the earth-sun system from t0 to t1 using Leapfrog Integration.
    INPUTS:
    config_func: function taking a date or date range and returning position and velocity of bodies
    accel_func:  function taking positions of the bodies and returning their accelerations
    t0: start date of the simulation; a python date
    t1: end date of the simulation (exclusive); a python date
    dt: time step in days.
    num_bodies: the number of celestial bodies in the simulation

def simulate_leapfrog(config_func: Callable, accel_func: Callable, 
                      t0: date, t1: date, steps_per_day: int):
    """
    Simulate the earth-sun system from t0 to t1 using Leapfrog Integration.
    INPUTS:
    config_func: function taking a date or date range and returning position and velocity of bodies
    accel_func:  function taking positions of the bodies and returning their accelerations
    t0: start date of the simulation; a python date
    t1: end date of the simulation (exclusive); a python date
    dt: time step in days.
    num_bodies: the number of celestial bodies in the simulation
    """
    
    # Length of the simulation (number of steps)
    N: int = (t1 - t0).days * steps_per_day

    # Get the initial conditions
    q0, v0 = config_func(t0)

    # Infer the number of dimensions from the shape of q0
    dims: int = q0.shape[1]
    
    # The time step in seconds
    dt = float(day2sec) / float(steps_per_day)
    # Square of the time step
    dt2: float = dt * dt
    
    # Initialize arrays to store computed positions and velocities
    q: np.ndarray = np.zeros((N, dims))
    v: np.ndarray = np.zeros((N, dims))
    
    # Initialize the first row with the initial conditions from the JPL ephemerides
    q[0, :] = q0
    v[0, :] = v0
    
    # Initialize an array to store the acceleration at each time step
    a: np.ndarray = np.zeros((N, dims))
    # First row of accelerations
    a[0, :] = accel_func(q[0])
    
    # Perform leapfrog integration simulation
    # https://en.wikipedia.org/wiki/Leapfrog_integration
    print(f'Performing leapfrog integration with {N} steps...')
    for i in tqdm(range(N-1)):
        # Positions at the next time step
        q[i+1,:] = q[i,:] + v[i,:] * dt + 0.5 * a[i,:] * dt2
        # Accelerations of each body in the system at the next time step
        a[i+1,:] = accel_func(q[i+1])        
        # Velocities of each body at the next time step
        v[i+1,:] = v[i,:] + 0.5 * (a[i,:] + a[i+1,:]) * dt
    return q, v