# Thread: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel?

1. ## How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel?

I'm trying to plot the location of an eccentrically-orbiting body based on time, and to do so I need to solve Kepler's equation for the Eccentric Anomally.
https://www.math.ubc.ca/~cass/course...01a/orbits.pdf

The Mean Anomaly: M = 2 pi * (t - T)/(P)

t = time
T = time at perhelion (closest to sun)
P = period, or year

The Eccentric Anomaly, E, is a bit harder to come by, and is defined (in a roundabout way) by Kepler's Equation:

M = E - e * sin(E)

where e = the eccentricity. 0 = circle, 1 = parabola.

If you have E, you can then plot out the location of the body relative to the sun. But first you have to find E.

I'm trying to do this in Excel, so I'd like to be able to do this automatically, without having to run a macro or use Goal Seek.

Does anyone know any numerical methods I could use to solve for E, that could be relatively easily implemented in Excel?

Edit: I'm doing this for fun, not for a class, so no worries about me using you guys to cheat on a homework assignment. I'm actually trying to plot the position of The Messenger, a comet in Dark Sun, with a 45-year period. Where an Athasian year is 375 days precisely. If I can get it to work, I may be able to plot out other "odd bodies" as well.

Also, it's been longer than I care to admit since my college class on Differential Equations, and I was never really that good at them back then.   Reply With Quote

2. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel?

First, the slightly snarky answer - the best way to solve for anything in Excel is to use the Excel sheet to hold data for Matlab to manipulate, and export back to Excel.

For the less snarky point, this looks like somewhere where a bisection method is entirely usable. The solution space is bounded with E somewhere between 0 and 2 pi, so if you use 0 = E - e*sin(E) - M, you can reliably start with a positive and negative point, take the mid point, keep whichever side keeps you around 0, repeat until you get within a certain proximity threshold.

I also wouldn't be surprised if this was a built in numerical solver.  Reply With Quote

3. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by Lord Torath I'm doing this for fun
Well, in that case, I recommend finding a monumentally wealthy Dane and stealing their life's work upon their death, which you will feel pretty guilty about but c'mon you can put it to much greater use than their illegitimate child (because they married a commoner) who just isn't as scientifically inclined as you.  Reply With Quote

4. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by Knaight First, the slightly snarky answer - the best way to solve for anything in Excel is to use the Excel sheet to hold data for Matlab to manipulate, and export back to Excel. Originally Posted by Knaight For the less snarky point, this looks like somewhere where a bisection method is entirely usable. The solution space is bounded with E somewhere between 0 and 2 pi, so if you use 0 = E - e*sin(E) - M, you can reliably start with a positive and negative point, take the mid point, keep whichever side keeps you around 0, repeat until you get within a certain proximity threshold.

I also wouldn't be surprised if this was a built in numerical solver.
OK. I think I could do this, but I'd need a ton of IF statements to help me decide which side to pick. I think it might be easier to just make a table for values of E ranging from 0 to 2 pi (or perhaps tau), and calculate M based on that, then use Index-Match to snag the value of E that gives the closest value of M. It's not quite the bisecting method, but I think it should work. If I use 101 values for E that gives me increments of 0.063, which is probably precise enough for my purposes. Maybe. With a period of 16,875 days, that means it will only move once every 168 days. I think I want this a bit more precise than that. I could go to 1000 cells, but I'm not sure I want to dedicate that much real estate to that.

I think I can implement Newton's Method, using much less real estate (I've set the eccentricity of The Messenger is 0.9751, so methods that assume a small eccentricty will not work.). Cell A1 sets E = 0. Cell A2 calculates EA2=EA1[EA1e*sin(EA1)M]/[1e*cos(EA1)].

Cell A3 = EA3=EA2[EA2e*sin(EA2)M]/[1e*cos(EA2)] and so on. I suspect I can probably get this to converge within 10 or 20 steps. I like that a lot better than calculating 1000 different values for E, and picking the closest.

1 Combined with a period of 16875 days, this gives a semi-major axis of 657.9 million miles, and a semi-minor axis of 146.19 million miles. So the orbit is 1316 million miles long and 292 million miles wide.

Edit: OK, zero is a bad starting point for E. After 23 iterations E was diverging to values in the tens of thousands. Starting at tau/2 (okay, fine, pi) converged to 5 decimal points on the ninth iteration. At least for T= 115.  Reply With Quote

5. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by Lord Torath M = E - e * sin(E)

where e = the eccentricity. 0 = circle, 1 = parabola.

If you have E, you can then plot out the location of the body relative to the sun. But first you have to find E.
I'm assuming e and M are known?

E2=M+e*sin(E1) seems to converge (which then at least gets rid of the multiple ifs) and reasonably fast where I've tried it.
You know 0<M<6.3 and 0<e<1, so there shouldn't be too many pathological sizes.  Reply With Quote

6. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by jayem I'm assuming e and M are known?
Yes. M is determined by time elapsed, time at perihelion, and the period (or year). Eccentricity dances with the semi-major axis, semi minor axis, and the period. Pick two of those and the other two can be determined. Originally Posted by jayem E2=M+e*sin(E1) seems to converge (which then at least gets rid of the multiple ifs) and reasonably fast where I've tried it.
You know 0<M<6.3 and 0<e<1, so there shouldn't be too many pathological sizes.
You're right, that does converge much more quickly. And avoids some of the oscillation I was getting with Newton's Method.  Reply With Quote

7. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by Lord Torath Yes. M is determined by time elapsed, time at perihelion, and the period (or year). Eccentricity dances with the semi-major axis, semi minor axis, and the period. Pick two of those and the other two can be determined.

You're right, that does converge much more quickly. And avoids some of the oscillation I was getting with Newton's Method.
It goes the wrong way when e>1 (which I was going to say doesn't happen but is a hyperbola, and is a bit rubbish when e>0.9 (which will)
I'm sure there must be a nice way to invert it.

But if you average the previous 2 terms before you iterate, that seems to do the trick well enough as well to buy a bit more time.  Reply With Quote

8. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel?

I think I've got something wrong with my calcs. Perihelion occurs at Time = 2250 (full period is 16875), and has a sun-comet distance of 657.9 million miles. But at T=1500, the sun-comet distance is only 284.2 million miles. Time to plot the full orbit, and not just its current position.

Yup. I've got the sun in the center of the orbit, rather than at one of the foci.

Huh. I could have figured that out much earlier of I'd just double-checked my semi-major axis (a), which is, you guessed it, 657.9 Mmiles.

Edit:
Oh. I forgot to subtract the center-sun distance (a*e) from the x-coordinate. OK, that looks much better!  Reply With Quote

9. ## Re: How Do I Solve Kepler's Equation for the Eccentric Anomaly in Excel? Originally Posted by Lord Torath OK. I think I could do this, but I'd need a ton of IF statements to help me decide which side to pick.
Not really - there are useful properties of this particular equation (mostly being monotonically increasing* with E) that makes this really easy. We can be confident that at E=0 you undercount (or are just there), and your 0 sided equation is negative. At E= 2*pi your equation will be positive. So whenever you generate a new point if you get a positive result, use it as the new upper bound. If you get a negative result, use it as the new lower bound.

This is especially useful when dealing with a wobbly Newton method, which is often a worry for sinusoidal functions in particular. One of the major upsides to the bisection method is that if you have correct outer bounds for a single intersection the domain will converge reliably around it - which is true of closed methods in general, which I'd recommend for this particular problem. Most open methods struggle with regions of near 0 derivative, which we have.

*Technically just monotonically non decreasing if e = 1, but at the 0.975 you set it's monotonically increasing.  Reply With Quote

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•