Found a solution thanks to  http://mcngraphics.com/thelab/blender/connect/

def FindRotation(ep1, ep2):
ep1 = numpy.array(ep1)
ep2 = numpy.array(ep2)
center = ep1 + 0.5 * (ep2 - ep1)

Position of ep2 if an imaginary sphere encompassing

both objects was transformed to 0, 0, 0

ep2new = ep2 - center

--------------------------------------------

Spherical coordinates (radius r, inclination theta, azimuth phi) from

the Cartesian coordinates (x, y, z) of o1 on our imaginary sphere

r = numpy.linalg.norm(ep2new) # radius of sphere
theta = math.acos(ep2new[2]/r)
phi = math.atan2( ep2new[1], ep2new[0])
return (0,theta, phi)