Orbit in space (extrasolar planets)

[Silly Lung]

I'm working on this project and I've kind of come to a roadblock.

Here's a little background: we have a 1 solar mass star and a planet orbiting it. The planet is defined by 5 orbital elements (period, mass, eccentricity, longitude of periastron, and mean anomaly). I wrote a code that will plot the radial velocity of the planet vs. time. I'm using the systemic console (www.oklo.org). I don't know how familiar you guys are with it; this is my first post here.

Now the radial velocity of the planet tells us the radial velocity of the star. The only difference is the semi-major axis is reduced by 1/(1 + Ms/Mp) and the path is flipped by 180°.

Here's the problem: I plot the radial velocity vs. time and I get an ok looking graph. Then I go to do a fitting test. The amplitude is dead on. So that's correct. But every time, the phase is off! Sometimes by ~30°, sometimes by ~180°. What's going on?

Here's my math (skip to the bottom to see the phase): I'm using the book Solar System Dynamics by Murray and Dermott. There's a neat matrix and I'm taking the radial velocity to be along the X axis. The equation for that is...

X = r(cosΩ*cos(ω + f) - sinΩ*sin(ω + f)*cosI)

where r is the radius, Ω is the longitude of ascending node, ω is the longitude of periastron, I is the inclination angle, and f is the true anomaly. Ω = I = 0, so our equation is reduced to...

X = r*cos(ω + f)

Take a derivative to get radial velocity. r and f change with time so you have to product rule.

There's useful formulas in the book for r dot and r*(f dot), but my final solution for the phase of the planet is...

-[e*sinω + sin(ω +f)]

e is eccentricity.

Just switch the sign for the radial velocity of the star. Now f can't be explained simply in terms of the mean anomaly, but it can with the eccentric anomaly, E.

tan(f/2) = sqrt[(1+e)/(1-e)]*tan(M/2)

It's messy but my new phase for the star looks like...

e*sinω + sin(ω + [2 * arctan(sqrt[(1+e)/(1-e)]*tan(E/2))])

If we consider a perfectly circular orbit (e = 0), then arctan-tan eliminate, 2s cancel and we have ω + M on the inside, as I'd expect.

I related E to M using Newton-Raphson iterations, and checked for convergence using Kepler's Equation. And it works. I think I'm overlooking something really obvious. I've been staring at this for a few days.

Anyone brave enough to conquer my wall of text and lend me some insight? It would be greatly appreciated! Thank you!

 

[tony873004]

On the systemic console, I think you have to tinker with the mean anomoly to get it in phase. But it sounds like you're trying to compute it, so no tinkering is required.

Your data file should be a long list of radial velocities. Is it possible that it begins at the wrong spot?

Also, be careful of the trig functions. They often provide answers in the wrong quadrant. For example, according to my calculator, tan(135) does not equal arctan(tan(135)). Rather than returning me to the 2nd quadrant, it gives me a 4th quadrant answer, so a little piecewise logic is required to get the desired quadrant.

 

[Silly Lung]

You are right. I'm generating my own mean anomalies so the value I get should match the value on the mean anomaly slider.

Interesting you mentioned the wrong spot thing. I hadn't thought about that until yesterday so I quickly fixed it. The phase is still off but not as bad...

Thanks for bringing the trig functions to my attention. That could very well be the problem. I'll take a look at them.

 

[Silly Lung]

It was the arctan-tan! Also making sin(ω +f) into a cosine made it work. I don't know why, but that's what happened.