def average(A_matrices, B_vectors, timesteps):
	""" 
	Average over matrices.

	Arguments:
	A_matrices: a dictionary of NxN matrices, indexed with their timestep.
	B_vectors: a dictionary of vectors with len N, indexed with their timestep.
	timesteps: a list or tuple of timesteps.

	Returns:
	A: the average of all A_matrices.
	B: the average of all B_matrices.
	"""

	# Initialize empty
	time = list(B_vectors.keys())[0]
	A = A_matrices[time] * 0
	B = B_vectors[time] * 0

	# Average by adding all objects and dividing by their number
	for timestep in timesteps:
		A += A_matrices[timestep]
		B += B_vectors[timestep]

	# Divide
	number_snapshots = len(timesteps)
	A /= number_snapshots
	B /= number_snapshots

	return A, B