factorize eccentricity computation
We might need this in different places:
Triplet r = planet.position();
Triplet v = planet.velocity();
Triplet h = cross(r,v);
real m = p.mass()+XMSTAR;
real d = r.norm();
real ener = 0.5*v.norm()*v.norm - m/d;
real a = -0.5*m/ener; // Semi-major axis
real e = sqrt(1- h.norm()*h.norm()/(m*a)); // eccentricity