How to compute time-dependent planetary motions

brandon-rhodes / pyephem

Scientific-grade astronomy routines for Python

Other

766 stars 121 forks source link

How to compute time-dependent planetary motions #216

Closed kabannister closed 2 years ago

kabannister commented 2 years ago

I understand that pyephem can be used to compute time-dependent planetary motions. A question I have encountered is" On average, what percent of the time is Mercury the closest planet (of the inner eight) to Pluto. For simplicity (if it can be described that way!), all orbits can be assumed to be circular and coplanar. I saw that Dr. Glen Miller describes how to access JPL data files and extract time-dependent XYZ coordinate data for the planets, but it is a very complicated process that requires detailed knowledge of how the JPL files are formatted. If there was a simpler way to download the [t, X(t), Y(t), Z(t)] data so that a MATLAB script could be written to read them and do calculations, that would be great.

brandon-rhodes commented 2 years ago

PyEphem is a Python library, so you will have to look elsewhere for resources that might help with MATLAB, which doesn’t use the Python language but an expression engine of its own. The Python library jplephem illustrates how Python reads JPL data files; maybe you could do something similar in MATLAB:

https://pypi.org/project/jplephem/

kabannister commented 2 years ago

Thank you, sir, that's a tremendous help!- I will take a look at jplephem. As I say, if I can generate output files that MATLAB can read, then I can get it to do the calcs. I forgot to include to include the following in my message: This fellow did a very nice job showing how he did a similar thing for Jupiter and Mercury a few years ago. No idea how he did it, and he only considered the inner 8 planets in his analysis. I guess because Pluto has fallen into disfavor lately as a "planet," it was not considered. https://www.ideatovalue.com/curi/nickskillicorn/2019/10/this-blew-my-mind-mercury-is-the-closest-planet-to-jupiter/

brandon-rhodes commented 2 years ago

I hope that jplephem lets you print the x,y,z positions in a format that Matlab can accept.

Thanks for the link to the blog post, which explains your interest in the calculation. I am not sure why the blog writer acts surprised? Given n planets going around a star, they will inevitably spend time half on one side of the star and the other half on the other side, and sometimes it will happen that the innermost planet alone winds up on the side of one of the outer planets. It will be a rarer and rarer event for each successive outermost planet, but it should still eventually be guaranteed to occur; so I don’t think he should be so surprised that it does?

drbitboy commented 2 years ago

The Geometry Finder (GF) subsystem JPL's NAIF/SPICE toolkit can do this. With a little bit of C programming*, you could get a set of windows when that was true, and then summing those windows and dividing by the total time observed, should yield that percentage you are looking for.

i.e. given a time, calculating whether Mercury is closest, which is about a dozen lines of code or even fewer if a loop is used. Plus some other code for which there are extensive examples.

kabannister commented 2 years ago

Thanking Dr. Carich for his great suggestion and the tip about the Geometry Finder (GF). Unfortunately, I don't know much about C programming and it would take me a long time to get proficient with it to work out how to do the programming task he outlined.. Again, it's a case of knowing what I want to do, but not knowing how to get it done. - KB

On Wed, Nov 3, 2021 at 4:09 PM Brian Carcich @.***> wrote:

The Geometry Finder (GF) subsystem JPL's NAIF/SPICE toolkit can do this. With a little bit of C programming*, you could get a set of windows when that was true, and then summing those windows and dividing by the total time observed, should yield that percentage you are looking for.

i.e. given a time, calculating whether Mercury is closest, which is about a dozen lines of code or even fewer if a loop is used. Plus some other code for which there are extensive examples.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-959935563, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS662IDSZ2CD2CCLB4F3UKGQIPANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

drbitboy commented 2 years ago

it's a case of knowing what I want to do, but not knowing how to get it done

You could pay someone.

kabannister commented 2 years ago

Thank you. I have in fact shelled out some bucks - to Mathworks, for their MATLAB Aerospace Toolbox, which may have exactly what I need. This Toolbox evidently has functions which will access the JPL database for all the planets and moons and extract what I need in the correct formats. I hope it works as advertized! Will keep you posted on how it goes. In the meantime, I downloaded Python 3.10, which is always a good thing to have in the arsenal. Keep up the good work! V/R, Ken Bannister

On Thu, Nov 4, 2021 at 10:11 AM Brian Carcich @.***> wrote:

it's a case of knowing what I want to do, but not knowing how to get it done

You could pay someone.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-961017711, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS62UDKEMQ5KNLINGOETUKKPDRANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

drbitboy commented 2 years ago

I just re-read the problem definition; we don't need ephemerides if we make enough first-order assumptions.

This code provides an approximation.

$ cat pm.py

import math
###      Merc Venus Eart Mars   Jupite Saturn Uranus Neptun Pluto
radii = (.39, .723, 1.0, 1.524, 5.203, 9.539, 19.18, 30.06, 39.53,)
Rpluto = radii[-1]
twopi = 2.0 * math.pi
P = 100.0
for R in radii[1:8]:
     costheta = R / (Rpluto * 2.0)
     twotheta = 2.0 * math.acos(costheta)
     frac = (1.0 - (twotheta / twopi))
     P *= frac
print(P)

$ python pm.py 1.2984404882776617

brandon-rhodes commented 2 years ago

Wow, that’s a surprisingly large number! I had assumed the answer would be a couple of magnitudes less.

drbitboy commented 2 years ago

I thought so too, then I realized that Mercury is at the Sun wrt Pluto to zeroth order, and also every other planet is on the other side of the Sun from Pluto at least half the time to zeroth order, so 2**(-7) is a lower limit to zeroth order.

Sidebar: so you agree that the algorithm is valid (to zeroth order, LOL)?

drbitboy commented 2 years ago

I realized the geometry finder is overkill: SpiceyPy can do this in 15 lines (pm2.py); a more general approach, with command-line arguments like sample frequency, takes 40 (pm_full.py); see the attached .ZIP archive.

$ python pm2.py de430.bsp 
{'mercury_closest_percent': 1.6871913826059255}   <== 11 centuries of daily samples

 python pm2.py de432s.bsp 
{'mercury_closest_percent': 0.5199211908931699}   <== 1 century of daily samples

So the previous 1.3% geometric result is within 25% of that 1.7% over 11 centuries, and with orbital periods of some bodies measured in a century or more, even 11-centuries of samples may have some bias; I'd call that a win.

SPKs are available here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/

And here are summaries of those SPKs (Spacecraft-Planetary-Kernels, generically a.k.a. SPICE kernels)


BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066

Summary for: de430.bsp

Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
        VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
        EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
        MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
        JUPITER BARYCENTER (5)  SUN (10)
        Start of Interval (ET)              End of Interval (ET)
        -----------------------------       -----------------------------
        1549 DEC 31 00:00:00.000            2650 JAN 25 00:00:00.000

Summary for: de432s.bsp

Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
        VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
        EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
        MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
        JUPITER BARYCENTER (5)  SUN (10)
        Start of Interval (ET)              End of Interval (ET)
        -----------------------------       -----------------------------
        1949 DEC 14 00:00:00.000            2050 JAN 02 00:00:00.000

pm_update.zip

kabannister commented 2 years ago

Gentlemen: I see that you have been burning the midnight oil doing further research on this problem! Awesome!

Now, when I look at this listing by Dr. Carcich, I can only guess it is written in Python. Maybe in C? Is it that laborious time-dependent calculations can be avoided by a clever technique? That is to say, Mercury only spends a certain amount of time alone looking at Pluto, and this code narrows that time to a percentage?
I see that the radii of the planets (in AUs?) are stored in a vector named 'radii'?
No idea what 'Rpluto = radii[-1]" means.
I am guessing that "P = 100.0" is a multiplicative factor to be used later.
'for R in radii[1:8]' appears to be the beginning of a loop over the vector of radii values. How the loop ends, I don't know. I'm used to seeing 'end' in other languages. Some trig calculations involving Pluto's orbit radius are being done in the loop, but can only guess what is being done.
Evidently at the end of the loop, the value P is used somehow to compute a fraction (perhaps a percentage?).
A value of P is printed out. No idea what this value represents. Is it a percentage?

$ cat pm.py

import math

Merc Venus Eart Mars Jupite Saturn Uranus Neptun Pluto

radii = (.39, .723, 1.0, 1.524, 5.203, 9.539, 19.18, 30.06, 39.53,) Rpluto = radii[-1] twopi = 2.0 math.pi P = 100.0 for R in radii[1:8]: costheta = R / (Rpluto 2.0) twotheta = 2.0 math.acos(costheta) frac = (1.0 - (twotheta / twopi)) P = frac print(P) $ python pm.py 1.2984404882776617

V/R,

Ken Bannister

On Thu, Nov 4, 2021 at 4:00 PM Brian Carcich @.***> wrote:

I realized the geometry finder is overkill: SpiceyPy can do this in 15 lines (pm2.py); a more general approach, with command-line arguments like sample frequency, takes 40 (pm_full.py); see the attached .ZIP archive.

$ python pm2.py de430.bsp {'mercury_closest_percent': 1.6871913826059255} <== 11 centuries of daily samples

python pm2.py de432s.bsp {'mercury_closest_percent': 0.5199211908931699} <== 1 century of daily samples

SPKs are available here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/

And here are summaries of those SPKs (Spacecraft-Planetary-Kernels, generically a.k.a. SPICE kernels)

BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066

Summary for: de430.bsp

Bodies: MERCURY BARYCENTER (1) SATURN BARYCENTER (6) MERCURY (199) VENUS BARYCENTER (2) URANUS BARYCENTER (7) VENUS (299) EARTH BARYCENTER (3) NEPTUNE BARYCENTER (8) MOON (301) MARS BARYCENTER (4) PLUTO BARYCENTER (9) EARTH (399) JUPITER BARYCENTER (5) SUN (10) Start of Interval (ET) End of Interval (ET)
    1549 DEC 31 00:00:00.000            2650 JAN 25 00:00:00.000
Summary for: de432s.bsp

Bodies: MERCURY BARYCENTER (1) SATURN BARYCENTER (6) MERCURY (199) VENUS BARYCENTER (2) URANUS BARYCENTER (7) VENUS (299) EARTH BARYCENTER (3) NEPTUNE BARYCENTER (8) MOON (301) MARS BARYCENTER (4) PLUTO BARYCENTER (9) EARTH (399) JUPITER BARYCENTER (5) SUN (10) Start of Interval (ET) End of Interval (ET)
    1949 DEC 14 00:00:00.000            2050 JAN 02 00:00:00.000
pm.zip https://github.com/brandon-rhodes/pyephem/files/7478130/pm.zip

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-961372944, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS64QXR34JG7SBV2UJO3UKLX6XANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

drbitboy commented 2 years ago

Burning the midnight oil? Not really; it took me probably less than half an hour to work out the geometric formulae and then a few minutes to code the prototype. The code to sample the data in the SPICE kernels took about the same amount of time (I am into my fourth decade of using SPICE).

Yes, Python

Clever technique? No, by using the already-existing libraries and tools e.g. the SPICE toolkit, SpiceyPy and list comprehension in Python. Because good coders borrow, but great coders steal.

Yes, radii is a tuple (Python sequence-of-objects object) of orbital radii, ordered by increasing value, of all nine bodies in AU. I originally retrieved the values from a random source, then got slight different values from from Wikipedia. The code uses them as radii of circular orbits (eccentricity = e = 1) per the original problem statement, even though they are actually semi-major axes of elliptical orbits (0 < e < 1).

radii[-1] is last element of the tuple i.e. the radius of Pluto

P is initialized to 100.0 because the problem statement asks for the fraction of time in percent. P represents the probability that Mercury is closer to Pluto than all other planets tested so far, so initially, before any planets were tested, that probably is 100% = 1.0.

As you inferred, for R in radii[1:8]: starts a loop over each radius, placed in R, except the first (radii[0], Mercury) and last (radii[8], Pluto). It could also be for R in radii[1:-1].

The extent of Python loops are determined by the indentation level of the statements that follow; there is no end or endloop or } statement. The final print(P) statement is at the same indentation as the for statement, which means it is the first statement after the end of the loop.

For first three statements in the loop, costheta = ...; twotheta = ...; frac = ..., estimate geometrically the fraction (frac; range will be [0:1]) of time that radius' (R's) body is not closer to Pluto than Mercury. Caveats:

Again, per the problem statement, assumes circular orbits which means constant angular velocity so time at any range of angles is proportional to that range
Assumes Pluto's position is static, or at least that its motion does not affect the value of frac significantly
Assumes Mercury in this context is coincident with the Sun, or more pedantically with the Solar System Barycenter. So the algorithm actually estimates the fraction of time that the Sun is closer to Pluto than any of the other seven besides Mercury. However, since (1) Mercury is closest to the Sun, (2) Mercury much closer than any of the other bodies, and (3) Mercury orbits the Sun so it is sometimes further and sometimes closer than the Sun, the algorithm considers Mercury's position as a "noisy" Sun and simply uses the Sun's position itself as a proxy for Mercury's position.

See the image below, pm_orrery.png; (I did the formula on the fly while coding, not noticing that the two cancel):

pm_orrery

The last statement in the loop, P *= frac would be more commonly understood as P = P * frac. Using this approach, the loop answers a restatement of the problem: what is the average percent of time Mercury is closer to Pluto than Venus AND Earth AND Mars AND ... AND Neptune? with A geometric estimate of the probability that Mercury is closer to Pluto than Venus TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Earth TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Mars TIMES ... TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Neptune Caveats:

Assumes those individual probabilities are uncorrelated

kabannister commented 2 years ago

Thanks so much for the detailed, line-by-line explanation. I was able to download the MATLAB Aerospace Toolbox and a necessary database, and was at least able to get, for a specific time, a read-out of one planet's XYZ coordinates in a Sun-based coordinate system. I will push on and try to get time-dependent coordinates for all the planets and then try to figure out an answer to my question that way. Will leep you posted. Thanks again! - Ken Bannister

On Fri, Nov 5, 2021 at 10:53 AM Brian Carcich @.***> wrote:

Yes, Python

Yes, radii is a tuple (vector-like Python object) of orbital radii, ordered by increasing value, of all nine bodies in AU. I originally retrieved the values from a random source, then got slight different values from from Wikipedia. The code uses them as radii of circular orbits (eccentricity = e = 1) per the original problem statement, even though they are actually semi-major axes of elliptical orbits (0 < e < 1).

radii[-1] is last element of the tuple i.e. the radius of Pluto

P is initialized to 100.0 because the problem statement asks for the fraction of time in percent. P represents the probability that Mercury is closer to Pluto than all other planets tested so far, so initially, before any planets were tested, that probably is 100% = 1.0.

for R in radii[1:8]: starts a loop over each radius, placed in R, except the first (radii[0], Mercury) and last (radii[8], Pluto). It should could also be for R in radii[1:-1].

The extent of Python loops are determined by the indentation level of the statements that follow; there is no end or endloop or } statement. The final print(P) statement is at the same indentation as the for statement, which means it is the first statement after the end of the loop.

For first three statements in the loop, costheta = ...; twotheta = ...; frac = ..., estimate geometrically the fraction of time (frac; range will be [0:1]) that radius' (R's) body is not closer to Pluto than Mercury. Caveats:

Again, per the problem statement, assumes circular orbits which means constant angular velocity so time at any range of angles is proportional to that range

Assumes Pluto's position is static, or at least that its motion does not affect the value of frac significantly

Assumes Mercury in this context is coincident with the Sun, or more pedantically with the Solar System Barycenter. So the algorithm actually estimates the fraction of time that the Sun is closer to Pluto than any of the other seven besides Mercury. However, since (1) Mercury is closest to the Sun, (2) Mercury much closer than any of the other bodies, and (3) Mercury orbits the Sun so it is sometimes further and sometimes closer than the Sun, the algorithm considers Mercury's position as a "noisy" Sun and simply uses the Sun's position itself as a proxy for Mercury's position.

The last statement in the loop, P = frac would be more commonly understood as P = P frac. Using this approach, the loop answers a restatement of the problem: what is the average percent of time Mercury is closer to Pluto than Venus AND Earth AND Mars AND ... AND Neptune? with A geometric estimate of the probability that Mercury is closer to Pluto than Venus TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Earth TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Mars TIMES ... TIMES A geometric estimate of the probability that Mercury is closer to Pluto than Neptune Caveats:

Assumes those individual probabilities are uncorrelated

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-961959992, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS62Z5MAWF4C5F5MUTVDUKP4YBANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

drbitboy commented 2 years ago

You can also duplicate my code in MATLAB and use SPICE kernels:

https://naif.jpl.nasa.gov/naif/toolkit_MATLAB.html

kabannister commented 2 years ago

Sir, many thanks again, I will follow up on your suggestion. V/R, Ken Bannister

On Sat, Nov 6, 2021 at 11:41 AM Brian Carcich @.***> wrote:

You can also duplicate my code in MATLAB and use SPICE kernels:

https://naif.jpl.nasa.gov/naif/toolkit_MATLAB.html

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-962469394, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS67FI5CEMI7CWRHFEH3UKVLEVANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

kabannister commented 2 years ago

Sir: I have not forgotten this topic - I have some home & hearth issues to attend to just now. Soon as I am done with those, in a week or so, I will get back to the Mercury-Pluto problem. I am just wondering if it is worth going through the agony of a time-dependent analysis, or does the product-of-probabilities approach described above knock the problem in the head? With a time-dependent analysis, it seems one would have to march along in time from say, 1/1/1900, to the present, get the daily (or hourly?) stand-off distances between Pluto and its eight companions, and then catalogue those time intervals when Mercury is the closest to Pluto, then sum up all those intervals and divide that sum by the total time and multiply the ratio by 100%. Can be done, but is it mad overkill?

V/R,

Ken Bannister

On Sun, Nov 7, 2021 at 7:42 AM Ken Bannister @.***> wrote:

Sir, many thanks again, I will follow up on your suggestion. V/R, Ken Bannister

On Sat, Nov 6, 2021 at 11:41 AM Brian Carcich @.***> wrote:

You can also duplicate my code in MATLAB and use SPICE kernels:

https://naif.jpl.nasa.gov/naif/toolkit_MATLAB.html

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brandon-rhodes/pyephem/issues/216#issuecomment-962469394, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWKUS67FI5CEMI7CWRHFEH3UKVLEVANCNFSM5HJLCUFQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

drbitboy commented 2 years ago

CPU and diskspace (the de440 SPKs total "only" 3GB) are cheap; it's hard to say if anything is overkill.

Best regards,

Brian Carcich