6
$\begingroup$

I am currently trying to perform the same interpolation that NASA uses to generate their Development Ephemeris files (Chebyshev coefficients). I have a table with position and velocity of the major planets in equal-spaced intervals. And I need to get this interpolation done, and get the Nth degree Chebyshev coefficients.

This process is well described at Newhall 1988 Numerical Representation of Planetary Ephemerides Celestial Mechanics, vol 45, 1–3, March 1988, pp 305–310 (viewable here as well). But my question is, do I need to implement this (as it will take lots of time) or someone has programmed it before? I mean, is there some MATLAB (could be any other program actually but I'm currently using MATLAB) routine which could make this for me? Maybe some SPICE subroutine?

I'm finding very hard to get this done due to the lack of information available online. If someone knows some explanation/example different from the Newhall one and can share it would be very appreciated.

Improve this question
$\endgroup$
6
  • $\begingroup$ Chebyshev interpolation uses unequal intervals. You need samples at specific locations (densely near endpoints, sparsely in the middle) to minimize interpolation error. It’s just like regular interpolation: you have a basis set of polynomials (chebyshev polynomials) with unknown coefficients and you use the values at the unevenly sampled intervals to solve for the coefficients. Have a look at the MATLAB function called Chebfun.m. $\endgroup$
    – Paul
    Apr 27, 2019 at 19:34
  • $\begingroup$ @Paul no, in this case the points will be evenly spaced. Have a look at the linked paper in the question (I added a viewable source). $\endgroup$
    – uhoh
    Apr 28, 2019 at 0:28
  • $\begingroup$ Great question! In an answer somewhere here (I'll look for it) I think it was explained that the attractiveness of the Chebyshev polynomials was that they provide a clear handle on how large interpolation errors are, and not that they are necessarily best for fitting to orbits. You could use other interpolators for personal work. I'm guessing you want to make Chebyshev polynomial-based ephemerides in order to use the same interpolators that work for the Development Ephemerides, is that right? If so, it seems you'd like to get a hold of a copy of subroutine PVCH as described in Newhall 1988. $\endgroup$
    – uhoh
    Apr 28, 2019 at 0:28
  • $\begingroup$ Answers to Is there a way to extract the Chebyshev coefficients for a body from a SPICE kernel? give some more insight into the construction of the kernels, but not on their generation. This is linked from "roll your own" in this answer but alas doesn't fit the bill. I was wondering if Project Pluto would have something to generate kernels but I don't see it. Section IV in adsabs.harvard.edu/full/1983A%26A...125..150N is interesting, but no help $\endgroup$
    – uhoh
    Apr 28, 2019 at 0:48
  • 2
    $\begingroup$ @uhoh - As an aside, "Doctor Mohawk" and I briefly worked for the same employer. He asked the referenced question when both of us worked for that employer. I didn't know who "Doctor Mohawk" was at the time I answered his question. He knew who I was; my ID is my real name. That I happened to answer in the middle of the night was what prompted his remark "Dave, do ever sleep?" Apparently not. $\endgroup$
    – David Hammen
    Apr 28, 2019 at 7:54

2 Answers 2

Reset to default
4
$\begingroup$

The standard (non-custom) tools for generating Chebyshev polynomial approximations involve evaluating the function to be approximated at the non-uniformly spaced Chebyshev nodes and computing the Chebyshev inner product.

What JPL did in 1988 was highly custom. It was not anything close to the standard Chebyshev fit. They used uniformly spaced data rather than the Chebyshev nodes, and they used an ad hoc least squares approach with exact matches at the end points rather the Chebyshev inner product, and with velocity (in units of au/day) arbitrarily weighted as 2/5 as important as position (in units of au). You will not find an off-the-shelf function that does this.

Moreover, this is what JPL did in 1988. It's been 30+ years since that paper was published. In the interim, vehicles have been sent to Mercury and Pluto (and beyond), and places in between. I highly doubt that the algorithm JPL used in 1988 is what they use now.

Improve this answer
$\endgroup$
0
3
$\begingroup$

It's possible that additional work has been done on this topic since the Newhall paper you linked, but if so, it does not cite that paper.

There is a non-peer-reviewed paper that also allows fitting acceleration data, as well, but acceleration is not used in the JPL ephemerides.

There are probably pieces of software in various places for producing Chebyshev fits of functions. For example, Numerical Recipes includes an algorithm which works when the original function is available. However, for fitting Chebyshev polynomials to coordinates and velocities at discrete times, the Newhall algorithm is my recommendation for several reasons:

  • If you only fit positions, but not velocities, the Chebyshev representation will not be faithful to actual spacecraft or planetary trajectories; it will only be correct at the times used for the fit and only for position.

  • Even at the points used for the fit, your error will be high.

  • Your terminal points (first and last time) will not be sufficiently constrained, so if you represent a trajectory using multiple Chebyshev polynomials, there will be discontinuities.

The Newhall algorithm addresses these problems. I don't see anything better out there, and I've looked around a bit.

As Dave Hammen mentions, the Newhall paper does include a velocity weighting of 0.4, which they say they arrived at experimentally, though they provide no data.

It is also beneficial, but not necessary, to use points which are arbitrarily spaced in specific circumstances. That is, if we fit a $n$-degree polynomial at the $n+1$ roots of the $n+1$ degree Chebyshev, we minimize the upper bound on the error; see Theorem 16.10 at the link. The roots of a degree $n+1$ Chebyshev polynomial are given by

$$t_i = \cos\left(\frac{2i+1}{2n+2}\pi\right)\,\text{ for } i = 0,\dots,n\,\text{.}$$

(Of course, $t$ needs to be scaled from $[-1,1]$ to your full ephemeris time interval $[a,b]$.)

The Newhall paper seems to use points which are equally spaced. According to the equation above, equally spaced points are not optimal, however. Indeed, if you're using the Newhall method, you should fit to the roots (the equation above) and, I think, to the position extrema as well (I have not proved this to my satisfaction yet). Including the extrema ought to also give you the terminal points.

The good news is it ought to only be two hundred or so lines of Python code (based on my own implementation). Might take you about a day to implement.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ @uhoh Fixed it :) $\endgroup$
    – Translunar
    Aug 1, 2019 at 21:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.