BenEngebreth.org

Implied Orbits and SAD Tracklets

HelioLinC and its versions and variants (Lambert tracklet IOD, RR phase space, Rubin's implementation) require you guess the heliocentric range of an object at an epoch. By guessing and checking a number of different initial range values you can find objects that occupy (approximately) the orbit described by that range hypothesis. Prior to now I've been modeling how that range evolves in time from the epoch as a simple quadratic: \(r(t) = r + \dot{r}t + \ddot{r}t^2\) where you guess the instantateous initial parameters \(r\), \(\dot{r}\), and \(\ddot{r}\) at the epoch \(t_{epoch}\). This works very well most of the time. It effectively works perfectly for all objects more distant than ~1.5AU (and even most < 1.5AU). But the quadratic model does not represent the actual \(r(t)\) as an object gets closer to the Sun. To give credit where credit is due, this reality was pointed out to me by Ari Heinze who I believe received this wisdom from Matthew Holman. So I started to pull the surprisingly productive thread of modeling \(r(t)\) from basic orbital mechanics. I'll jump right into the math and highlight some implications at the end.

Start with a heliocentric range assertion

We assert \(r\), \(\dot{r}\), and \(\ddot{r}\) for an observation at \(t_{epoch}\). The asserted radial acceleration \(\ddot{r}\) is the sum of the radial acceleration due to gravity and the radial acceleration due to the tangential component of the velocity.

\[\ddot{r} = \frac{-GM}{r^2} + \frac{v_{t}^2}{r}\]

We can solve for tangential velocity \(v_{t}\) by rearranging the equation above.

\[{v_t}\ = \sqrt{r(\ddot{r} + \frac{GM}{r^2}})\]

Now let me make the claim that we can select any plane perpendicular to \(\hat{r}\) to apply \({v_t}\) in to get a full velocity vector \(\vec{v}\) at \(t_{epoch}\) because the range \(r(t)\) does not depend on the plane the tangential velocity is in (more on this in a bit). For instance, we can cross \(\hat{r}\) with the unit vector pointing in the z-direction \(\hat{z}\) to get a unit vector \(\hat{v_{t}}\) perpendicular to \(\hat{r}\) to apply the tangential velocity in.

\[\hat{v_{t}} = \frac{\hat{r} \times \hat{z}}{|\hat{r} \times \hat{z}|} \quad where \quad \hat{z} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}\]

Then, combined with the \(\dot{r}\) we asserted (\(\dot{r}\) = \(v_{r}\) at \(t_{epoch}\)) and the unit vector \(\hat{r}\) pointing to the observation at \(t_{epoch}\), we can create a full velocity vector at the epoch.

\[\vec{v} = v_{r}\hat{r} + v_{t}\hat{v_{t}} \]

We also know \(\vec{r}\) at \(t_{epoch}\) from the range assertion \(r\) applied to our heliocentric projection of the observation position on the sky \(\hat{r}\) at \(t_{epoch}\). Thus we have a full state vector \(\vec{r}, \vec{v}, t_{epoch}\) and can propagate this state vector to all observation times to get \(\vec{r}(t)\) and \(\vec{v}(t)\). From \(\vec{r}(t)\) we can easily calculate the range as a function of time \(r(t)\). The \(\vec{r}(t)\) and \(\vec{v}(t)\) vectors are specific to the arbitrary orbital plane we've chosen, but again the range \(r(t)\) does not depend on the plane.

The values \({v_r}(t)\) and \({v_t}(t)\) also do not depend on the plane and we can calculate them for each observation time as below.

\[{v_r} = \frac{\vec{r} \cdot \vec{v}}{|\vec{r}|} \quad {v_t} = \frac{|\vec{r} \times \vec{v}|}{|\vec{r}|}\]

And there is one final piece of information we can glean. We can use the \(\vec{r}(t)\) we generated on our arbitrary plane to calculate the angular displacement in time \(\Delta\theta(t)\) from the observation at \(t_{epoch}\) for each observation time given this hypothesis. Once again, this information does not depend on the plane the object is moving in.

So by asserting \(r\), \(\dot{r}\), and \(\ddot{r}\) at \(t_{epoch}\), we can calculate \(\vec{r}(t)\), \(\vec{v_r}(t)\), \({v_t}(t)\) and \(\Delta\theta(t)\) for all orbits with that asserted range hypothesis. We know the vector form \(\vec{v_r}(t)\) because \({v_r}\) always points in the radial direction \(\hat{r}\). The only thing our asserted range hypothesis can't tell us is the plane of the orbit which would determine \(\hat{v_t}(t)\).

What we have now is something much more informative than an \(r(t)\) model alone. We now have the description of an orbit that includes not only \(r(t)\), but \(\vec{v_r}(t)\), \({v_t}(t)\) and \(\Delta\theta(t)\) as well.

Plane independence

In the section above I kept stating that certain measures did not depend on the orbital plane. To support that claim, consider the image below. In it, the position vector \(\vec{r}\) is pointing from the Sun to an object for our hypothesized range \(r\). The radial velocity vector \(\vec{v_r}\) points in the \(\hat{r}\) direction (we'll assume \({v_r}\) is a negative value in this instance and thus \(\vec{v_r}\) points towards the Sun). \({v_t}\) is the magnitude of the tangential component of the velocity. Our assertion of \(r\), \(\dot{r}\), and \(\ddot{r}\) give us all of these values at \(t_{epoch}\). If we knew which way \(\hat{v_t}\) was pointing we'd know the full state vector. But we know that \(\hat{v_t}\) can only point in so many ways; it must be perpendicular to \(\hat{r}\). You can imagine all of its possible pointings then by "holding" the position vector \(\vec{r}\) in between your fingers and spinning it about its long axis like a pencil. As \({v_t}\) rotates through 3 dimensions, it draws a circle that is bisected by the "page" as we view it from above. These are the only planes that \(\hat{v_t}\) can be in.

The components of the \(\vec{r}\) and \(\vec{v}\) that we know at the epoch given our assertion


To test the claim that \(r(t)\), \(\vec{v_r}(t)\), \({v_t}(t)\) and \(\Delta\theta(t)\) do not depend on the plane \(\hat{v_t}\) is in (as long as \(\hat{v_t}\) is perpendicular to \(\hat{r}\)), we can simply test all of those \(\hat{r}\) perpendicular pointings of \(\hat{v_t}\) for a given initial range assertion to see if any of the time dependent values change based on the \(\hat{v_t}\) we choose. The answer is: they do not. Here is an outline of the process that describes this test.

  1. Select a random heliocentric unit vector \(\hat{r}\) (e.g. [1,0,0])
  2. Apply a range hypothesis to this vector (e.g. \(r\)=3, \(\dot{r}\)=0, \(\ddot{r}\)=-1e-5) to get \(\vec{r}\) at the hypothesis epoch (let's call it t=0)
  3. Calculate \(v_{t}\) at the hypothesis epoch from the equation above
  4. Cross \(\hat{r}\) with a random vector (e.g. [0,0,1]) and normalize to get a direction perpendicular to \(\hat{r}\) to apply \(v_{t}\) in - this is a \(\hat{v_{t}}\)
  5. Rotate \(\hat{v_{t}}\) through 360 degrees about the axis pointing towards \(\hat{r}\) and keep each of these unit vectors to representing all \(\hat{v_{t}}\) values
  6. Create the velocity vectors at the hypothesis epoch with the equation \(\vec{v} = v_{r}\hat{r} + v_{t}\hat{v_{t}}\) for each \(\hat{v_{t}}\)
  7. Propagate \(\vec{r}\) and each \(\vec{v}\) to N different times (e.g. t=[1,2,3]) for all velocity vectors
  8. Calculate \(r(t)\), \({v_r}(t)\), \({v_t}(t)\) and \(\Delta\theta(t)\) for all N times and all velocity vectors using equations above
  9. Repeat for more range hypotheses in step 2 then compare with the results from step 8 (they should all be the same)

An illustration of what \(\Delta\theta(t)\) is measuring


HelioLinC implications: Accurate \(r(t)\) and SAD tracklet rejection

What do we gain with this new information? For one, we can calculate a more accurate \(r(t)\) because we're modeling gravity now instead of using a quadratic approximation that assumes constant acceleration. The plots below use Rubin DP0.3 synthetic solar system object data and show the norm of the modeled \(r(t)\) minus the true \(r(t)\) as a function of heliocentric distance in AU (left) or geocentric distance in AU (right). You can see in the plot on the left below that the orbit model (red) of \(r(t)\) generally has a lower error than a quadratic (blue) to the data. The two instances where the orbit model does worse than the quadratic model are where the object in question is very close to the earth as seen in the plot of the data as a function of geocentric distance on the right. This is due to the two-body model I'm using for propagation that does not account for the effect of earth's gravity for objects that closely approach it.

But we've also gained three new bits of information. The measures \(\vec{v_r}(t)\), \({v_t}(t)\) and \(\Delta\theta(t)\) were not pieces of the puzzle I was using previously, at least. With HelioLinC we use the solution to Lambert's problem to determine the \(\vec{r}(t)\) and \(\vec{v}(t)\) that connect a tracklet with observations at times (say) \({t_m}\) and \({t_n}\) with an orbit for the hypothesized range. We can now test to see if the angle swept by the tracklet projected to heliocentric coordinates at the time varying range from our new \(r(t)\) model agrees with the angle sweep \(\Delta\theta({t_n}) - \Delta\theta({t_m})\) we'd expect for the observations along the orbit given by our range hypothesis at those same times. If there's a large angular discrepancy, the tracklet doesn't conform to the orbit our range hypothesis implies and we don't have to propagate that tracklet. The heliocentric projected positions of these tracklets have what I'll call swept angle discrepancies relative to the angular displacements represented by our initial range hypothesis implied orbit. They are SAD tracklets for this hypothesized range.

A tracklet (black dots) with a swept angle discrepancy to the orbit implied by the range hypothesis (gray dots)


There are two benefits from SAD tracklet rejection: 1) You spend less compute time propagating tracklets because you skip the tracklets that don't conform to the orbit implied by the initial range hypothesis. 2) The clustering phase space has fewer orbit non-conforming propagated tracklets in it that might generate mixed linkages. Fewer propagated tracklets in the phase space also makes clustering faster, of course. The speedup benefit seems to be the larger one in practice because the RR phase space already casts non-conforming tracklets to positions that wouldn't cluster anyway, but I do see a small increase in pure linkages and a small decrease in mixed linkages in my testing of SAD tracklet rejection compared to no tracklet rejection. But the really big win seems to be the amount of time you need to spend propagating and clustering. In my tests, the total time spent in the propagation and clustering steps of HelioLinC are reduced by ~40% for NEOs. And the benefit of this should be greater as you increase the number of hypotheses you test because you can use a tighter tolerance for SAD tracklet rejection. For this reason, SAD tracklet rejection should be most helpful for linking NEOs where you have to vary \(r\), \(\dot{r}\), and \(\ddot{r}\) in a fine-grained manner in order to accruately model \(r(t)\) and thus have to test a lot of hypotheses.

Conclusion

As I said at the beginning, modeling \(r(t)\) correctly was an unexpectedly productive thread to pull. Not only did we end up with a better \(r(t)\) model, but we gained this new understanding of the expected angular displacement as a function of time for a tracklet and given range hypothesis that is directly applicable to our HelioLinC implementation. We can now reject tracklets that don't have the expected angular sweep for the hypothesis being tested before we propagate and (inevitably fail to) cluster them, which significantly improves HelioLinC's runtime.

Published: 01/23/2025