A few weeks ago, the ESA Nearth Earth Object Coordination Center started a series of NEOCC riddles about Near Earth Object orbits and related topics. The first riddle was about orbits with a peculiar characteristic: they spend 50% of the time inside some fixed radius from the Sun (1.3au in the riddle), and the remaining 50% of the time outside this radius. It was published on June 4. Shortly after that I submitted my solution. The deadline for sending solutions ended yesterday, so today NEOCC has published their solution together with the list of people that solved the riddle correctly. In this post I publish my solution and make some additional comments.
First, something I find worthy of mention is that only 6 people have submitted a correct solution. It would be interesting to know the total number of people that submitted a solution, but I guess it is not much larger than that. The riddle was perhaps quite difficult depending on the approach to the problem (more on that later), but it didn’t require very advanced knowledge. Just basic knowledge of the mathematical description of Keplerian orbits (so you can do as me and look up the formulas in Wikipedia). Therefore, I would like to encourage more people to have a go with the next riddles. The next one will be published on June 30, which is Asteroid day.
My solution is now published in this Jupyter notebook. If your read the introduction of this solution, or the introduction of the solution published by NEOCC, you’ll see that after writing everything in mathematical terms, the riddle boils down to solving a constrained optimization problem.
At first, I tried to solve this problem analytically. This seems quite difficult if at all possible, and indeed I filled many pages with different approaches using Lagrange multipliers before giving up and doing it numerically with the computer.
My numerical approach actually boiled down to computing the maximum of the function\[f(E_0) = \frac{\sin E_0 + E_0 – \frac{\pi}{2}}{\sin E_0 ( 1 – (E_0 – \frac{\pi}{2}) \tan E_0))},\]where \(E_0\) denotes the eccentric anomaly at the moment when the mean anomaly is \(\pi/2\). I was quite lazy and found a numerical approximation to the maximum by evaluating \(f\) over a find grid of points (which you have to do anyway if you want to plot things to gain some insight into the problem), and picking the maximum over this grid.
Looking now at NEOCC’s solution, I see that they have made more involved calculations, since they have also introduced the true anomaly \(\theta\) and they approach the problem by using the fact that the derivative of the function to maximize must vanish at the maximum point. After some calculations, they arrive at\[\tag{4}0 = 1 + \cos \theta_0 – \sqrt{1-\varepsilon^2}\varepsilon\sin \theta_0 \frac{\sin E_0}{1 – \varepsilon \cos E_0},\]where \(\theta_0\) is the true anomaly at the moment when the mean anomaly is \(\pi/2\), and \(\varepsilon\) is the eccentricity (NEOCC uses a slightly different notation; here I’m keeping to the notation I used in my solution). This equation must be understood as an equation for \(\varepsilon\), since \(E_0\) is defined implicitly in terms of \(\varepsilon\) by Kepler’s equation\[E_0 – \varepsilon \sin E_0 = \frac{\pi}{2},\]and \(\theta_0\) can be written as an explicit function of \(\varepsilon\) and \(E_0\), since\[\tan^2 \frac{\theta_0}{2} = \frac{1+\varepsilon}{1-\varepsilon} \tan^2 \frac{E}{2}.\]
Now, NEOCC’s solution involves finding numerically a value of \(\varepsilon\) that satisfies (4) by solving for \(E_0\) and \(\theta_0\) as indicated above. This involves solving Kepler’s equation by iteration.
So drawing upon NEOCC’s idea to find the point where the derivative of the function to maximize vanishes, and trying to improve my laziness, we can follow my approach and compute and solve \[f'(E_0) = 0.\]The derivative \(f'(E_0)\) has a rather ugly expression but is straightforward to compute explicitly and only involves trigonometric functions in \(E_0\). Afterwards, your favourite root-finding numerical method can be used to find the point where the derivative vanishes. So maybe this will give an algorithm that converges to the maximum rather quickly, which is useful if we want to solve to many decimal places (actually several numerical methods for maximum finding are disguised root-finding methods for the derivative). Note that this doesn’t involve solving Kepler’s equation.
So in the end I like my approach more than NEOCC’s. By using \(E_0\) instead of \(\varepsilon\) as the independent variable I don’t need to solve Kepler’s equation, and I think that introducing the true anomaly into the problem only complicates calculations more. One thing is clear: it seems that it’s not easy to solve this without doing some numerical approximation in one way or another, so I would be gladly surprised if someone can show us another approach.
By the way, reviewing my solution I’ve found a silly mistake. When scanning the NEOCC database for objects that spend approximately 50% of their time inside 1.3 au and have an aphelion as large as possible, I limited the aphelion search to those between 1.5 and 1.6 au, while actually the maximum aphelion can go up to 1.603 au, as I showed in the first part of the riddle.
In this respect, I really like NEOCC’s simple and tidy solution for the second part of the riddle. But they don’t tell us how to compute the percentage of time below 1.3 au using the values in the database! Well, this is simple. If we denote by \(C\) the quotient between the aphelion and 1.3au, we have\[C = \frac{1+\varepsilon}{1 – \varepsilon \cos E},\]where \(E\) denotes the eccentric anomaly at the moment when the 1.3au radius is crossed. Since the eccentricity \(\varepsilon\) is also in the database, we can get \(E\) from this equation and then find the corresponding mean anomaly as\[M = E – \varepsilon \sin E.\]Finally, \(M/\pi\) gives us the fraction of time that the object spends inside the 1.3au radius.