Closed bryanwweber closed 2 years ago
bryanwweber
commented
2 years ago By the way, an easy workaround is to just fix the center of the mars object when it's created:
mars = planets["Mars Barycenter"]
mars.center = 10This carries over when .at() is called.
brandon-rhodes
commented
2 years ago We might need to call in @JoshPaterson to help answer this one!
The mars.at(t) call is performing correctly. It is supposed to return a position relative to 0 the Solar System Barycenter (SSB), because without that a later .observe() call could not properly apply relativity (that’s why JPL ephemerides report positions relative to the SSB, and not the Sun; because the Sun isn’t inertial, it’s in its own little orbit around the SSB with Jupiter and Saturn as the major partners).
It looks like the code, given an orbit around the SSB, assumes the SSB has no mass:
mu = GM_dict.get(position.center, 0) + GM_dict.get(position.target, 0)This is a bit confusing. Shouldn't it raise an error demanding that the user provide the mass of the center, if the center doesn't have a mass in the table? @JoshPaterson, this might need an exception instead of a .get(), at least for the center.
@bryanwweber, what you might want to do here is ask for the position of Mars relative to the Sun, then set the orbital-elements logic loose on the resulting position. Something like:
sun = planets['Sun']
b = (mars - sun).at(t)Try that out and let us know the result!
bryanwweber
commented
2 years ago It looks like the code, given an orbit around the SSB, assumes the SSB has no mass
This makes sense, and I had identified the same line as the culprit. Would the fix be to set the "mass" of the SSB equal to the sum of the GM parameters for all the solar system bodies? Horizons has no trouble calculating the orbital elements relative to the SSB:
The Keplerian GM from the HORIZONS output is 1.3289039033282275E+11 km^3/s^2 and
print(f"{sum(GM_dict.values()):.5e}")gives
'1.33069e+11'So it's close but not quite the same.
Edit: It just occurred to me that Horizons might not use all the masses, but only the masses of the major planetary systems, so
print(f"{sum(GM_dict[k] for k in range(1, 11))}")gives
'1.328905186666e+11'which is much closer. Let me check with the update DE 441 values from https://ssd.jpl.nasa.gov/ftp/xfr/gm_Horizons.pck
Edit 2: The DE441 values give '1.328905186713e+11' as the sum of the 10 major-body elements. So essentially the same as the existing values, and still slightly larger than the HORIZONS output.
bryanwweber
commented
2 years ago what you might want to do here is ask for the position of Mars relative to the Sun, then set the orbital-elements logic loose on the resulting position. Something like:
sun = planets['Sun'] b = (mars - sun).at(t)Try that out and let us know the result!
To actually answer this request, this works just fine for my purposes, so this is what I'll go with for now!
brandon-rhodes
commented
2 years ago Great, I'm glad that a sun-centered position is working for you! That lets you move forward while we wait and see if Josh responds with any thoughts on his default-to-zero code and on how to best handle the SSB.
JoshPaterson
commented
2 years ago Here are answers to some questions!
Why doesn't GM_dict have a value for the SSB, with a key of 0?
I used data from Horizons (ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck) to create that dictionary. I didn't notice that the file doesn't have a value for the SSB, so the dictionary doesn't either. I agree that it should be added. I wonder why it's not in that file?
If an object is not in the dictionary why does its mass default to 0? The value for GM is calculated as the sum of the masses of the center and the target, for example if you want the elements of mars orbiting the sun, GM would be the sum of the mass of both objects. The code assumes that if an object isn't in the dictionary it's mass is negligible. The goal here was to include the mass of the target in the calculation of GM if it is significant. I believe this is how Horizons works; if I'm wrong about that it should be considered a bug.
Why would the GM values used by Horizons be slightly different from GM_dict?
The values used by HORIZONS are newer and slightly different than those in GM_dict. In #420 we discussed updating the data in skyfield but decided against it.
What if the primary object is a barycenter? I think this is the problem that needs fixing. Currently, if the primary is a barycenter the secondary's mass is double counted. As discussed in #420 it probably would have been better to always have the user provide their own GM, but for now I'll look into how Horizons calculates GM if the center or target is a barycenter. I think it would be best to do the same thing they do. What do you all think?
Should an error be raised in osculating_elements_of if GM is 0?
Yes, I think so. I also think osculating_elements_of should have a keyword argument added so users can specify their own GM value, which would be passed though to the OsculatingElements object. Then the error can say something like: "unable to calculate GM for these objects. Use keyword argument _ to specify your own." Does that sound OK?
bryanwweber
commented
2 years ago As discussed in #420 it probably would have been better to always have the user provide their own GM, but for now I'll look into how Horizons calculates GM if the center or target is a barycenter. I think it would be best to do the same thing they do. What do you all think?
@JoshPaterson This sounds like a great idea. I think if Skyfield could grow the ability to load the PCK file, which could be passed to osculating_elements_of, that'd be super handy, and would seem to somewhat resolve the concerns in #420.
On a somewhat unrelated note, although in the same function, the use of a vector of times and a reference frame at the same time doesn't work at the moment. The position comes in as a 3xN array, and the reference frame is a 3x3xN array, so the .dot() product doesn't broadcast correctly. I got it working with this einsum:
np.einsum("mnr,nr->mr", reference_frame, position.position.au)and similar for the velocity. I'm not sure how to reconcile that with the 3x3 array used for a single time though. A more complete example showing the error:
planets = load_file("de440.bsp")
ts = load.timescale()
start_date = datetime(2021, 10, 28, tzinfo=tz.UTC)
end_date = datetime(2021, 12, 3, tzinfo=tz.UTC)
dates: list[datetime] = list(rrule(DAILY, dtstart=start_date, until=end_date))
t = ts.from_datetimes(dates)
ecliptic_frame = build_ecliptic_matrix(t)
sun = planets["Sun"]
mars = planets["Mars Barycenter"]
position = (mars - sun).at(t)
osculating_elements_of(position, ecliptic_frame)
Traceback (most recent call last):
File "/Users/bryan/GitHub/orbital-mechanics/reference/scripts/planet-orbital-elements.py", line 31, in <module>
osculating_elements_of(position, ecliptic_frame)
File "/Users/bryan/.pyenv/versions/3.10.0/envs/orbital-notes/lib/python3.10/site-packages/skyfield/elementslib.py", line 25, in osculating_elements_of
position_vec = Distance(reference_frame.dot(position.position.au))
ValueError: shapes (3,3,37) and (3,37) not aligned: 37 (dim 2) != 3 (dim 0)It makes sense to be that there would be no mass for the SSB, since — so far as I know? — nothing orbits around it. Small objects like comets and asteroids orbit the Sun. (As I always remember when I try asking a comet for its position, and it turns out that it needs to be combined with a Sun vector!) If someone tries to compute an orbit around a non-object like the SSB, we should probably raise a polite error explaining what to do instead?
The default-to-zero rule can probably be modified here. Something like:
Check if the target has a mass — and we could give asteroid and comet orbits a default mass of zero — and if it has one, then use it.
Else check if the target is a simple integer, in which case look it up in the GM_dict (for backward compatibility; but as discussed in #420, we really want folks to provide a GM_dict of their own selection, maybe from a group of dated GM_dict’s that we provide, or from an explicitly loaded .pck file).
If a mass hasn’t been found, raise an error explaining how the user might find and provide one.
The two optional parameters to osculating_elements_of() could be something like center_GM_km_s= and target_GM_km_s=. Oh — except that we seem to use mu instead of GM? Then, mu, I suppose.
While I’ll leave it to Josh to put together a more complete answer about the time vector and reference frame, I’ll at least offer an initial guess that if a vector time has been offered, then the routine should quietly add an additional dimension so the rotation matrix so that it has an extra dimension to match that of the time vector. (And thanks for the example showing the error, @bryanwweber!)
brandon-rhodes
commented
2 years ago Oh — and I should answer an earlier comment from @bryanwweber:
This makes sense, and I had identified the same line as the culprit. Would the fix be to set the "mass" of the SSB equal to the sum of the GM parameters for all the solar system bodies?
If I remember Particle Dynamics class, that would only make sense for a mass outside of the orbits of all of those planets, maybe very far outside, like an Oort Cloud object — but I don’t think you can include, say, the outer planets in the SSB mass if the object were an asteroid in the asteroid belt?
bryanwweber
commented
2 years ago Thanks for the replies @brandon-rhodes!
- It makes sense to be that there would be no mass for the SSB, since — so far as I know? — nothing orbits around it.
My particle dynamics is also rusty, but I thought in an n-body problem, the bodies orbited the common center of gravity? In the case of the solar system, the SSB is sufficiently close to the Sun that the error in assuming the Sun as the center is minimal for the orbits of the planets, but it's not strictly true.
I suppose the fact that HORIZONS also includes a GM value for the SSB/allows you to calculate elements relative to the SSB means that folks (like me 🤪) will expect Skyfield to do something with respect to the SSB. In Park et al.'s paper for DE 440/441, they include a plot of the SSB relative to the Sun (Fig. 2) and a plot of the Earth/Moon Barycenter relative to the SSB. So I hope raising an error isn't where this ends up 🙃
brandon-rhodes
commented
2 years ago @JoshPaterson — I spent some time this morning with Skyfield’s PCK text file parser, since it wasn't quite able to read gm_Horizons.pck. (It expected a space between variable name and equals sign, whereas the GM file starts with BODY000_GMLIST= without a space between them.) I’ve split it out from planetarylib.py, improved it, and will now go add a bit of documentation.
You can try reading the raw file by installing the dev version of Skyfield:
pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip— and then:
f = open('gm_Horizons.pck', 'rb')
d = {}
g = text_pck.load(f, d)
print(d)Do you have a guess as to whether reading this file should be the job of a GM-constants class, or become a new feature of the existing PlanetaryConstants manager?
It looks like what you’re missing is the ability to add together masses to form the SSB, since the file doesn’t give a mass for the SSB, but instead says:
BODY000_GMLIST= ( 1 2 3 4 5 6 7 8 9 10
199 299
301 399
401 402 499
501 502 503 504 505 506 514 515 516 599
601 602 603 604 605 606 607 608 609 610 611 612 615 616 699
701 702 703 704 705 799
801 803 804 805 806 807 808 899
901 902 903 904 905 999
2000001 2000002 2000003 2000004 2000007 2000010 2000015 2000016
2000031 2000052 2000065 2000087 2000088 2000107 2000433
2000511 2000704 )
Hi! Thanks for the wonderful package. I'm using the DE440 kernels, and trying to get the osculating elements. It seems like the
data.gravitational_parameters.GM_dictis missing a0key for the SSB, so the results end up being for Mars orbiting around itself.prints
Continuing with
prints
If the fix is as simple as adding a new key to
GM_dictwith the value for the sun (currently at key10), I'm happy to make a pull request! Otherwise, I'm a little bit out of my depth here 😄