docs: 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