Get the positions and velocities of the sun and eight planets
    Returned as a tuple q, v
    q: Nx3 array of positions (x, y, z) in the J2000.0 coordinate frame.

def configuration(t0: date, t1: Optional[date] = None,  
                  steps_per_day: int = None) -> Tuple[np.ndarray, np.ndarray]:
    """
    Get the positions and velocities of the sun and eight planets
    Returned as a tuple q, v
    q: Nx3 array of positions (x, y, z) in the J2000.0 coordinate frame.
    """
    # Default steps_per_day = 1
    if steps_per_day is None:
        steps_per_day = 1
        
    # Time step dt is 1.0 over steps per day
    dt: float = 1.0 / float(steps_per_day)

    # Default t1 to one day after t0
    if t1 is not None:
        # Convert t to a julian day
        jd0: int = julian_day(t0)
        jd1: int = julian_day(t1)
    else:
        jd0: int = julian_day(t0)
        jd1: int = jd0 + dt

    # Pass the times as an array of julian days
    jd: np.ndarray = np.arange(jd0, jd1, dt)

    # Number of time steps
    N: int = len(jd)
    
    # bodies is a list of the celestial bodies considered; should be in an enclosing scope
    # Number of bodies
    B: int = len(bodies)
    # Number of dimensions
    dims: int = B * 3
    
    # Initialize empty arrays for position q and velocity v
    q: np.ndarray = np.zeros((N, dims))
    v: np.ndarray = np.zeros((N, dims))

    # Position and velocity of the sun as arrays of length 3
    body_ids: List[int] = [jpl_body_id[body] for body in bodies]

    # Fill in the position and velocity for each body in order
    for i, body_id in enumerate(body_ids):
        # The slice of columns for this body (same in q and v)
        slice_i = slice(3*i, 3*(i+1))
        # Extract the position and velocity from jpl
        qi, vi = jpl_kernel[0, body_id].compute_and_differentiate(jd)
        # Convert positions from km to meters (multiply by km2m)
        q[:, slice_i] = qi.T * km2m
        # Convert velocities from km / day to meters / sec (multiply by km2m, divide by day2sec)
        v[:, slice_i] = vi.T * (km2m / day2sec)

    # Return tuple of Tx6 arrays for position q and velocity v
    return q, v