def calculate_velocities(cartesians, timestep=1):
    """
    Calculate velocities at each timestep given Cartesian coordinates. Velocities at the first and last point are
    extrapolated.
    :param cartesians: Cartesian coordinates along trajectory
    :param timestep: time step between frames in units of fs, default=1
    :return: velocities
    """

    velocities = []
    for i in range(0, len(cartesians)):
        if i == 0:
            velocity = (cartesians[i + 1] - cartesians[i]) / timestep
        elif i == len(cartesians) - 1:
            velocity = (cartesians[i] - cartesians[i - 1]) / timestep
        else:
            velocity = (cartesians[i + 1] - cartesians[i - 1]) / 2 * timestep

        velocities.append(velocity)

    return velocities