As a simple project to introduce myself to orbital mechanics I'm attempting to calculate impact times of an object close to Earth. Right now I am using only Keplerian mechanics (no air resistance or other perturbing forces). From Fundamentals of Astrodynamics, given an object with eccentric anomaly $E_0$ at $t_0$ and $E_1$ at $t_1$:
$$t_1 = t_0 + \sqrt{\frac{a^3}{\mu}}\left( E_1 - e\sin E_1 - (E_0 - e\sin E_0)\right)$$
where $a$ is the semi-major axis, $e$ is the eccentricity, and $\mu$ is the gravitational parameter. We can then calculate the radius of the periapsis. Assuming a spherical Earth, if this is less than the radius of the Earth the object will impact. If this is the case, letting $r_1 = R = $ Earth radius, we can calculate the eccentric anomalies at $t_0$ and $t_1$ using:
$$ \cos\nu_i =\frac{a(1 - e^2) - r_i}{e r_i}$$ $$ \cos E_i = \frac{e + cos\nu_i}{1 + e\cos\nu_i}$$
where $\nu_i$ is the true anomaly, and $r_i$ is the distance from the centre of the Earth at $t_i$. This method works well for non-parabolic trajectories (I'm open to modifications for parabolic trajectories).
The next step is to loosen the restriction of a spherical Earth and allow Earth to be an ellipsoid, in other words let $R = R(z)$ (symmetric around rotational axis). Specifically I'd like to allow for a bulge at the equator by prescribing an equatorial radius, a polar radius and an eccentricity. Is it possible to solve this problem analytically or does it have to be done numerically?
This is my first post on this stack exchange; I felt it was better suited here than scicomp for example. If that's not the case please let me know.