Category Archives: Gravity

Planetary Orbits in Javascript

3 Replies

I have developed an enthusiasm (aka weird obsession) for celestial mechanics and developed several games and the Unity Asset Gravity Engine . In developing Gravity Engine I learned a great deal about high pedigree N-body simulations, but in some cases (e.g. a model of the solar system) it is not really necessary to simulate the system but rather just evolve it in the correct way. In this case Gravity Engine offers the option to use Kepler’s equation and move bodies in their elliptical orbit as a function of time. This uses far less CPU than doing the 9*8 mutual gravitational interactions (10*9, if you add in the “dwarf planet” Pluto).

[If you don’t see an animation, click the post title to see ONLY this post. There is WordPress JS bug when this animation and a YouTube link are present]

Creating code to move a planet in an elliptical orbit with the correct velocity is surprisingly tricky. This is one of those cases where you might expect you could just grab a formula from Wikipedia and bang out some code. This bias comes from all the examples in physics class where the goal is to “find/derive the formula” and get a tidy equation.

If you dig around on physics pages for an equation of an elliptical orbit you will generally encounter the equation for the shape of the orbit with eccentricity e, and semi-major axis, a:

r = \frac{a(1-e^2)}{1+\cos \theta_F}

I have attached the subscript F to the angle to indicate this is the angle from the focus of the ellipse between the position of the body and the long axis of the ellipse. Historically, this angle is named the True Anomaly.

Where is time in this equation?

Nowhere. This equation doesn’t tell us anything about time.

To get an equation for how the object moves as a function of time, we’ll need Kepler’s equation.  Kepler constructed his equation without Calculus (which came along about 60 years after Kepler did this work) using geometric arguments and the assumption that an objects speed in an elliptical orbit was inversely proportional to it’s distance from the focus. Kepler’s equation is:

M = E - e \sin E

Here M is the position of a body moving in a uniform circle at a constant rate (that we will relate to time) and E is the angle to the position on the ellipse from the origin, called the eccentric anomaly. Here a picture will help:

Eccentric_and_True

The eccentric anomaly E is NOT the same angle as \theta_F (f in the picture) however they can be related with a bit of geometry:

\cos \theta_F = \frac{\cos E - e}{1 - e \cos E}

\sin \theta_F = \frac{\sqrt{1-e^e} \sin E}{1 - e \cos E}.

If we have a specific time we want a position for we need to convert this into a value of M. This is done by dividing the time by the time per orbit, T. Kepler can also help us here with his third law that relates the size and eccentricity of the orbit to the mass of the bodies:

T = \frac{\sqrt{a^3}}{m}

where m is the combined mass of the central object and orbiting body. (Kepler did not know the proportionality constant was the mass, that came later).

Given M, Solve for E

Ok, we’ll just isolate E…hmmm. E appears by itself and inside the sine function. That sinks our chance of getting a tidy mathematical formula. This equation is legendary in Mathematics, since it is an early example of a transcendental equation with an import application. It has been studied extensively and the approaches are well summarized in the book “Solving Kepler’s Equation Over Three Centuries” by Peter Colwell.

There are some series approximations, but they are not valid for all eccentricities. The most common approach is to iterate the equation until we converge on a value that is “good enough”.

The “recipe” for tying this all together is:

  1. Determine the orbital period T
  2. For the time t we’re interested in, divide by the orbital period and use the remainder to find the angle if the body were moving in a uniform circle, M.
  3. Using iteration solve Kepler’s equation and find E
  4. Using E, determine \theta_F
  5. Find the corresponding r using the orbit equation.

In Javascript this becomes:

var orbitPeriod = 2.0 * Math.PI * Math.sqrt(a*a*a/(m*m)); // G=1

function orbitBody() {

 // hide last position drawn
 context.clearRect( last_x -10, last_y - 10, 20, 20);

 // 1) find the relative time in the orbit and convert to Radians
 var M = 2.0 * Math.PI * time/orbitPeriod;

 // 2) Seed with mean anomaly and solve Kepler's eqn for E
 var u = M; // seed with mean anomoly
 var u_next = 0;
 var loopCount = 0;
 // iterate until within 10-6
 while(loopCount++ < LOOP_LIMIT) {
 // this should always converge in a small number of iterations - but be paranoid
 u_next = u + (M - (u - e * Math.sin(u)))/(1 - e * Math.cos(u));
 if (Math.abs(u_next - u) < 1E-6)
 break;
 u = u_next;
 }

 // 2) eccentric anomoly is angle from center of ellipse, not focus (where centerObject is). Convert
 // to true anomoly, f - the angle measured from the focus. (see Fig 3.2 in Gravity) 
 var cos_f = (Math.cos(u) - e)/(1 - e * Math.cos(u));
 var sin_f = (Math.sqrt(1 - e*e) * Math.sin (u))/(1 - e * Math.cos(u));
 var r = a * (1 - e*e)/(1 + e * cos_f);

 time = time + 1;
 // animate
 last_x = focus_x + r*cos_f;
 last_y = focus_y + r*sin_f;
 drawBody( r* cos_f, r*sin_f, "blue");
 setTimeout(orbitBody, 30);
}

I have left out some of the init code for clarity. If you view the source for this page you can find all this.

(Eccentric anomaly image created by CheCheDaWaff/Creative Commons License).

ThreeBody – Can Gravity “keep it all together”?

Leave a reply

I have returned to N-body, my first mobile effort , with the goal of making it cross -platform and adding new ideas. Given all my recent positive experiences with the Unity game engine the idea is to port as much as possible and avoid “the rewrite”. (This has a very tempestuous history in software and is one of the things Joel Spolsky thinks you should never do.)

The good news for the small community who have enjoyed N-body is that moving it to Unity is going quite well, and I am FULL of cool ideas I want to add: gravitational fields for non-spherical objects, galactic potentials, dust accretion – it’s almost endless. The bad news is that I read a great Sci-Fi novel “The Three-Body Problem” by one of China’s top sci-fi writers, and this has caused me to set aside Nbody for awhile. “The Three Body Problem” reminded me of the non Sci-Fi book with the same title (you could see how that would happen). One of the things I made a note of when I read the more academic of these books was the static three body problem. The question is: can three bodies be placed in initial conditions with zero velocity such that they will stay in a bound configuration? Would you expect that since gravity is attractive they would “keep together”, perhaps with some cool triple orbits?

As I was reading the sci-fi 3BP I was feeling like I needed a break from NBody and figured I’d try a quick mock up of the static three body problem in Unity. As I get more kung-fu with Unity (and using the leap frog integrator I had ported for Nbody) this kind of side project becomes very do-able. I spent an afternoon doing the basics and found that there was a nice little physics “time waster” game here.

Shortly after, during my usual lunch time perusal of arxiv, I found a paper about using GPUs for N-body simulations which mentioned replacing the three bodies with three binaries. Another addition for the ThreeBody app.

This was all about a month ago and since then I have added the minimum of “gamification” so others can take a swing at this problem on their mobile devices. It is weirdly addictive and a good way to see how sensitive the problem is to initial conditions. Hope someone out there also finds it fun.

You can find ThreeBody for Android and Blackberry (coming soon to  Amazon and iOS). There is some further technical detail.

[social4i size=”large” align=”float-right”]

N-body Choreographies

Leave a reply

You would think that after three hundred years the consequences of the law of gravity for “simple” systems would have been pretty well studied and nothing surprising would emerge. Amazingly there is such richness and complexity in the physics of three masses that new discoveries are still being made.

The opening of the door was the paper “Braids in Classical Gravity” by Christopher Moore 1. This did open the door but the paper seems to have missed its audience a bit since one of its results was rediscovered a few years later 2. What these papers discovered was new families of solutions to the gravitational N-body problem not just for N=3 but for much more complex systems. They involve a mathematical assumption that the N bodies follow the same path and were termed choreographies. The gem at the center of all these solutions is a stable figure-8 solution for three equal masses.

These solutions are hypnotizing and as soon as I learned of them I had to get them into the N-body app. The solutions have very special initial conditions and when I first put them in they did not evolve as expected. The version of N-body at the time had a very simple numerical integrator, known as Leapfrog with a fixed time step 3. The choreographies required a more careful approach. Even though N-body was close to being ready to go I could just not give up on these choreographies. So I put in a 4th order Hermite integrator with some adaptive time step capabilities. This was enough to get the simple chain choreographies running. (I then found some more cool choreographies on arXiv, in fact new classes of choreographies. These will need more integrator tuning but look for them in a future update to N-body. In addition the Java gallery by Vanderbei has some solutions I’d like to include in the next release).

I’d love to understand how these solutions are found, but I’m just not there yet. Scholarpedia N-Body Choreographies  is great place to start. A useful expository paper that also finds even more solutions (including some new non-choreography solutions) is Vanderbei’s New orbits for the N-body problem. His paper provides code to calculate choreographies in the specialized language AMPL. This makes the details of the search quite clear, even if you (like me) have never seen AMPL before.

Based on what I’ve absorbed thus far it seems the general idea involves using a path optimization approach using variational principles. This comes from an idea which is in the very core of physics, the notion of finding the path along which a specific quantity (the action) is an extreme value. I learned how to do this when I was in graduate school but I have always been a bit baffled by how the mathematical machinery is searching over paths. There are such a huge number of possible paths! The key mathematics here is the calculus of variations. I’ve gotten as far a buying a copy of the classic staring point Calculus of Variations (aka Gelfand and Fomin) but it is still mocking me quietly from the shelf…

 

Notes:

  1. Phys.Rev.Lett 70:3675-3679. This was in turn anticipated by M. H´enon (1976) “A family of periodic solutions to the planar threebody problem, and their stability.” Celestial Mechanics 13 267–285. ↩
  2.  A. Chenciner and R. Montgomery (2001) “A remarkable periodic solution of the three body problem in the case of equal masses.” Annals of Mathematics 152 881–901. ↩
  3. This is a good place to acknowledge the places I read about N-body evolution codes.  Firstly, Moving Stars Around which is written in a quirky Socratic-dialog style. As part of the learning arc some of the code shown has deliberate errors that are then explained, so be careful. Everything you ever wanted to know about the subject is in the comprehensive Gravitational N-Body Simulations: Tools and Algorithms (Cambridge Monographs on Mathematical Physics). For a more concise summary of N-body codes and their pros and cons the section in Galactic Dynamics: (Second Edition) (Princeton Series in Astrophysics) is very to the point.  ↩

+1 Changes Everything

Leave a reply

There are times in physics when what seems like a simple extension to a problem causes the floor to drop out. The motion of three bodies due to gravity is one of those cases. For the two-body solution there is a modest differential equation and a solution that describes the motion as an ellipse, parabola or hyperbola. Adding a third body seems like “no big deal”. You might think that the math gets a bit messier so maybe the solution will not be quite as clean as a simple ellipse – but there should still be some kind of sensible solution. This is one of those cases where our optimism/intuition fails.

Newton’s work was published in the 1680s. It took other great mathematicians to get further. Euler found an unstable 3-body solution in the early 1700s. Lagrange found another unstable solution as well as some stable solutions in special cases for the 3-body problem. His work was published in 1772 – almost a hundred years later. [Lagrange also fundamentally changed physics with his Mechanique Analytique. Today huge parts of modern physics is still rooted in “the Lagrangian”].

Lagrange found a few special cases where three bodies could stay in a stable solution by making simplifications that result in what is called the “restricted three body problem”. In this case one body is “heavy”, the other has moderate mass and the third has no mass. As a concrete example think of the Sun, Jupiter and an asteroid in the same orbit around the sun as Jupiter. As a further simplification, assume Jupiter moves in a perfectly circular orbit. It is then convenient to think in terms of the position of the third massless body in relation to Jupiter i.e. in a co-ordinate frame that rotates around the sun at exactly the same rate as Jupiter. Lagrange found there were five points at which the third body experienced zero net force and could stay fixed in relation to Jupiter in this rotating frame.

This in itself was an interesting result. However these five Lagrange Points would only be really interesting if the massless body actually stayed in this relation to Jupiter even if it was bumped away from this zero force location i.e. if the motion around a Lagrange point is stable. It turns out that two of the five Lagrange points are stable but there’s a catch. The catch is that the Lagrange points L4 and L5 are stable only if Jupiter is not “too heavy” compared to the sun. Exactly how heavy? L4/L5 are stable as long as Jupiter is not more than 3.85 percent the mass of the Sun.

These points are of more than academic interest. This solution motivated a search for asteroids near these locations in Jupiter’s orbit and as the technology of telescopes improved the first such asteroid, Achilles, was found in 1906. There are now thought to be as many of these “Trojan” asteroids as there are asteroids in the belt between Mars and Jupiter. Other planets have such asteroids as well. Recently, more asteroids were found at the L4 point of Mars (http://arxiv.org/pdf/1303.0124.pdf).

What if the test particle is not quite at the stable Lagrange points L4/L5 or does not have quite the right velocity? Then the resulting motion is very interesting. This leads to orbital paths relative to Jupiter that get names such as tadpole or horseshoe orbits, due to their shape – you can get a sense of why these names arise by playing the levels in the Lagrange group in N-body. The image below shows oscillations about the L5 point in the N-body iPad app.

Tadpole

Oscillations about Lagrange points also have real-world implications. For example, the detection of a dust ring around Uranus puzzled astronomers since location of the ring could not be stable. It was proposed that there were “sheparding satellites” 1. Subsequent discovery of two small moons, Ophelia and Cordillia by Voyager lend credence to this idea. 2

Amazingly (and a testament to just how opaque the physics becomes when N=3) another stable three body solution for equal masses was found in the 1990s more than two hundred years after Lagrange. It was a result of considering topologies of paths. We’ll look at this fascinating solution and some of it’s unstable cousins in the choreography group of N-body levels.

One aside: In general relativity the step of going from one to two makes the floor fall out. In the case of Einstein’s general theory of relativity (GR) finding a solution for ONE body is really hard. The spherically symmetric solution was found by Schwarzschild in 1916 about a year after GR was published. It took until the 1960s before the axially symmetric one-body solution was found. There are no two body solutions corresponding to a simple orbit that can be written in terms of modest mathematical functions. This is part because in GR the idea is that gravitational interactions do not operate via a force but by defining the shape of the space they are part of. As bodies move around each other there are ripples in space-time (waves) that need to be taken into account. This makes things very messy. Numerical modelling and simplifications can be employed and things can be learned – but there is not a tidy mathematical description.

Notes:

  1. Goldreich and Tremaine (1979) Nature 277 97-99 ↩
  2. See the discussion in 10.5.3 of Solar System Dynamics ↩

Did Gravity have to be an inverse square law?

3 Replies

How fundamental is the fact that gravity is proportional to 1/r^2 versus a more general 1/r^n or something more complicated like f(r) = e^{-r}? Is this solution unique? Does it have any special properties other forces would not have?

This was settled in a publication in to the Royal Society in 1684 by Newton. The idea that the force might be 1/r^2 was not new.  This had been proposed by Hooke, at the time the resident experimenter at the recently formed Royal Society. Hooke postulated an inverse square law but did not prove it, lacking the mathematical insight of Newton. He was always a bit bitter about not being recognized for his contribution in postulating the correct formula and was combative in his exchanges with Newton. It is possible that Newton’s “If I have seen further it is by standing on the shoulders of giants…” was partly to acknowledge Hooke. It’s also possible this was a bit of a dig at Hooke who was a tall, somewhat stooped gentleman. 1

To understand what was really done by Newton it helps to pause and think about what he had to work with. The most precious thing was accurate data – the positions of planets measured carefully over sustained periods of time. This data measured the shape of the planets orbits and their speed at different points of the orbit. Kepler had done a lot of data collection and in 1609 using this information he was able to see patterns that helped him determine relationships in the data. Kepler’s first law was the radical notion that planets moved in ellipses (not circles) and the second was the bodies moved such that the line from the body to the center swept out an equal area in equal times.

So what do you think? Is that in itself enough to determine what force gravity is? Enter Newton.

Newton proved a number of things about generic forces acting on a line between two bodies. In particular he established that for such a force the motion will be in a plane and there will be a constant of motion (angular momentum) that is conserved. He also demonstrated that the shape of the path of one body around the other is the same as each body around their common centre of mass. He showed that for a general force a body will sweep out equal areas in equal times. This fact alone does not distinguish between forces. It is true for any central force.

How about the path of the orbit? Does that uniquely determine the force? Halley (of Halley’s comet fame) had the same question and in 1684 he went to Cambridge to see Newton. He asked what shape a force of 1/r^2 would produce. Newton’s responded that it would be an ellipse and that he had proved it. Halley was amazed and asked for the proof, leaving Newton mumbling “Gee, I know that used to be around here somewhere…”. Newton promised to provide it and several months later he provided it and a lot more in his note The Motion of Gyrating Bodies (De Motu Corporum in Gyrum). This was read to the Royal Society on Dec. 10, 1684. This then sparked Newton into writing the Principia which was then published in 1687. (Halley had to foot the bill for the publication, since the Royal Society had “bet big” on a a previous publication The History of Fishes that turned out to be a flop and consumed all their funds. At one point the society paid Robert Hooke in copies of the book.)

Principia is simply amazing. Any one of the theorems in it would have given Newton a place in history but the book has dozens and dozens of ground breaking results. The way mathematics was written at the time and the terminology used makes a translation of this work pretty tough going, it’s an improvement on trying to read the original latin but not by much. Fortunately there is a marvellous book by Chadrasekhar Newton’s Principia for the Common Reader. It is a masterful book from a physicist whose ability to do difficult calculations is astounding. (His The Mathematical Theory of Black Holes (Oxford Classic Texts in the Physical Sciences)
is perhaps the most densely packed book of equations I ever studied.) The term “common reader” is perhaps a bit generous to the reader, since there are differential equation at about a second year undergrad level.

Newton did WAY more than just figure out the law of gravity. He looked at central forces in general, proved that a spiral would be the path due to a force of 1/r^3. (Depending on the specific type of spiral the force could be slightly different. For example the spiral r = k \theta is of the form f(r)\propto \frac{1}{r^3} + \frac{2 k^2}{r^5} 2)

Newtonian gravity does have one special characteristic. The orbits are closed. After one orbit the object returns to exactly where it started. That is the path is an ellipse and not a “spirograph” type pattern. This characteristic also occurs for a force of f(r)\propto r^2. No other forces have this characteristic, so in a way 1/r^2 is unique.

Note I was careful to say “Newtonian gravity”. Actual gravity is best described by Einstein’s general relativity. One of the drivers for revising our view of gravity was the fact the the orbit of Mercury was not closed. It’s point of closest approach to the sun shifts around very slightly over time (even after all the effects of the other planets are taken into account). As an approximation this can be modelled as a slightly different central force. The careful measurements of the orbits in the solar system (and increasingly the orbits of binaries and exoplanets) are still scrutinized for signs that our understanding of gravity is correct. (See e.g. http://arxiv.org/pdf/gr-qc/0302048.pdf) Another reason for re-thinking gravity is that a central classical force has a a kind of spooky action at a distance built in. If a mass doubles then that effect is immediately felt everywhere in the universe and that would be hard to reconcile with the idea that nothing moves faster than light. In general relativity the changes in the curvature propagate with finite speed so there is a time delay (although you can play games with the equations in some cases using negative energy density to make changes go at faster than the speed of light 3)

Notes:

  1. Hooke was an interesting figure who did a great deal and is probably a bit under-appreciated; see The Forgotten Genius: The Biography Of Robert Hooke 1635-1703
    . ↩
  2. This comes from a book I used to bridge the gap between my engineering undergrad education and “real” physics. Introduction to Classical Mechanics by Atam P. Arya (first edition, 1990). Example 7.1  ↩
  3. I first bumped into this in graduate school in a paper:M. Alcubierre, “The Warp Drive: Hyper-fast Travel within General Relativity”, Class. Quantum Grav., 11, pp.L73-L77, 1994. This is still an active area of investigation; see http://arxiv.org/pdf/1202.5708.pdf ↩

Why are orbits not spirals?

1 Reply

As an engineering undergrad I had to take a course in probability and statistics. This course made me crazy because I could seldom get the right answer to the assigned problems. These were inevitably questions about counting certain groups of events and I would come up with a counting scheme that seemed logical but was often wrong. I could look at the provided answer which counted in some other way and it made sense BUT I could not see why my answer was wrong.

I was reminded of this when I showed N-body to a friend and his first question was “Why don’t the orbits spiral in?”. I kept trying to repeat the answers from my physics education but we kept coming back to “I’m sure your answer is correct, but why is mine wrong?”. To deal with that it finally occurred to me to ask him why he though they were spirals – instead of continuing to blather on about ellipses.

“If I throw a ball (even really far) – it hits the ground. This would happen in a vacuum too. Is this not a spiral path?” was the response. This is actually a good observation – which when linked with other observations about decaying orbits of satellites could easily leave you with the idea that orbits are spirals.

The explanation I came up with was that if the earth was replaced by a single point of the same mass and you threw the ball then the path would be an ellipse and it would come back to where it was. When the real earth is in the way, the ball hits the earth. To a person on the earth it does appear that the object is getting closer to the centre – it is. If you looked really carefully you would discover the path you thought was a spiral was actually a segment of an ellipse.

Why an ellipse? That’s something I used to be able to derive. It turns out to be not as simple a process as you might hope. After all it took someone with the stature of Newton to figure it out. The way we now do that calculation using differential equations is quite different from the geometric approach Newton used. Ironically it draws far more on the Calculus notation of Leibnitz (who independently developed Calculus, something Newton never acknowledged and worked hard to undermine).

As part of developing N-body I tried to dust off that part of my brain and sat down with a pad of paper to see how far I could get. Embarrassingly, I was a little rusty since my PhD is from the mid-nineties and since then I have not done much math. Nevertheless, I peeked only as much as needed and did finally get to the answer. I’ll post about that as a separate topic.

The details do really matter. That’s what science IS. One of my favourite illustrations of this viewpoint is in the movie Insignificancewhich has the characters of Marilyn Monroe and Albert Einstein. Marylyn explains special relativity in essentially correct lay-person terms and says:

Marilyn: I understand the results and the premise. I guess that’s the main thing, huh?

Einstein: That’s nothing.