Kepler's Second Law - reaching an unsolvable equation?

[Upsilon]

I have been working on modeling the orbit of a satellite using Kepler's second law of planetary motion, and I have gotten to a point that is really quite bothersome to me. Essentially, my problem boils down to solving this equation for θ (angular position of the planet from the focus of the orbit in radians) in terms of t (time elapsed since perigee in seconds):

(b²sin θ)/(a²cos θ - ca) = tan (2πt/P)

Everything else is a constant; P is the orbital period in seconds, a is the semi-major axis, b is the semi-minor axis, c is the focal distance of the orbit, and π is, well, pi. I can produce my work leading up to this point, although I'm fairly certain it is correct and not terribly relevant to the problem at hand. Is anyone able to solve this?

 

[jfizzix]

I have been working on modeling the orbit of a satellite using Kepler's second law of planetary motion, and I have gotten to a point that is really quite bothersome to me. Essentially, my problem boils down to solving this equation for θ (angular position of the planet from the focus of the orbit in radians) in terms of t (time elapsed since perigee in seconds):

(b²sin θ)/(a²cos θ - ca) = tan (2πt/P)

Everything else is a constant; P is the orbital period in seconds, a is the semi-major axis, b is the semi-minor axis, c is the focal distance of the orbit, and π is, well, pi. I can produce my work leading up to this point, although I'm fairly certain it is correct and not terribly relevant to the problem at hand. Is anyone able to solve this?

Putting this into another form, the equation is:

\frac{W \sin(\theta)}{X \cos(\theta) - Y}=Z
for some W,X,Y, and Z.

Since \sin(\theta)=\sqrt{1-cos^{2}(\theta)}

you can use this substitution, square both sides of the equation, and then solve for Cos(\theta) with the quadratic formula.

 

[Upsilon]

Putting this into another form, the equation is:

\frac{W \sin(\theta)}{X \cos(\theta) - Y}=Z
for some W,X,Y, and Z.

Since \sin(\theta)=\sqrt{1-cos^{2}(\theta)}

you can use this substitution, square both sides of the equation, and then solve for Cos(\theta) with the quadratic formula.

Wow, thanks, I can't believe I missed that. I have tried something similar before, but I didn't use any substitutions and got all the constants mixed in with each other, eventually leading me to a nasty quartic. At that point I just dropped it.

Anyway, this is the work I went through:

(\frac{W\sqrt{1-\cos^2\theta}}{X\cos\theta - Y})^2 = Z^2

\frac{W^2(1 - \cos^2\theta)}{X^2\cos^2\theta - 2XY\cos\theta + Y^2} = Z^2

W^2 - W^2\cos^2\theta = Z^2X^2\cos^2\theta - 2XYZ^2\cos\theta + Z^2Y^2

0 = (Z^2X^2 + W^2)\cos^2\theta - 2XYZ^2\cos\theta + Z^2Y^2 - W^2

\cos\theta = \frac{2XYZ^2 \pm \sqrt{4X^2Y^2Z^4 - 16X^2Y^2 + 8X^2W^2 - 8Y^2W^2 + 4W^4}}{4X^2 + 2W^2}

\theta = \arccos(\frac{2XYZ^2 \pm \sqrt{4X^2Y^2Z^4 - 16X^2Y^2 + 8X^2W^2 - 8Y^2W^2 + 4W^4}}{4X^2 + 2W^2})

This will give me 2 answers for theta, no?

 

[jfizzix]

perhaps, but will only one make physical sense? What happens if t is near zero?

 

[Upsilon]

Ok, so I spent some time trying to simulate this, and it simply wasn't working correctly. I went back to my other work leading up to the equation above and realized I had made a critical mistake. So let me start from the beginning:

Kepler's Second Law: \frac{A}{t} = \frac{2\pi t}{P} where P is the orbital period of the object, t is time elapsed, and A is total area swept out after t has elapsed.

Area of an ellipse: \pi a b

Area of a segment of an ellipse: A = \frac{1}{2} ab\theta where θ is measured from the center of the ellipse counter-clockwise from the positive x-axis, and A is the area of the segment encompassed within θ.

Polar equation of an ellipse with a focus at the origin: r = \frac{b^2}{a - c\cos\theta}

I derived my initial equation from these. It required that I manipulate the area of a segment equation to where θ was measured from the focus instead of the center of the ellipse. Let's say θ₁ is the angle measured from the center of ellipse, and θ₂ is the angle measured from the focus, and r is the length of the segment from the focus to the border of the ellipse measured at θ₂. I developed this relation, unknowing at the time that it only works situationally:

\theta_1= \arctan(\frac{r \sin\theta_2}{r \cos\theta_2 - c})

It is correct in the 1st quadrant but nowhere else. How can I develop a correct equation modeling the relationship between angle from the center of the ellipse and the angle from its focus? Or do I need to use a different situational equation for each quadrant?

 

[Janus]

This problem ( position as a function of time) has to be done by iteration.

For how to do this, check here:

http://www.slideshare.net/IngesAerospace/orbital-mechanic-4-position-and-velocity-as-a-function-of-time

 

[Upsilon]

Yes, I'm starting to see that now. It is because of the inability to solve for the definite area of an ellipsoidal segment where θ is measured from the focus of the ellipse, no? In this illustration, that means finding the total area right of r and above the x-axis. That would mean finding the area of the segment θ₁ and adding the area of triangle rck to it. The relationship between θ₁ and θ₂ is:

θ_1 = \pi - \arccos(\frac{c - r\cos θ_2}{\sqrt{r^2 + c^2 - 2rc\cos θ_2}})

Regardless of that equation, the area that I want would be:

A = \frac{abθ_1 + cr\sin θ_2}{2}

Which I believe is unsolvable for θ₂ even with the correct relation between θ₁ and θ₂, right?

 

[Upsilon]

Ok, so in a last-ditch effort to find a closed-form solution to the problem (I'm stubborn like that) I tried to integrate the polar equation of the ellipse to get the area swept out at θ. I have gotten to a point where I must integrate this:

A = \int {\frac{1}{2} (\frac{b^2}{a - c\cos\theta})^2 d\theta}

And I honestly have no idea how to do so, and to my knowledge has not been done before. And even if I do, I'm not sure I will be able to derive a closed-form solution for θ in terms of A.

 

[Upsilon]

Alright, I've done some work and I've found the integral I needed and thus a closed-form solution for A as a function of θ. I then plug this in for A in Kepler's law:

\frac{A}{t} = \frac{\pi ab}{P_o}

And attempt to solve for θ. The integral is really quite nasty, so I made some substitutions to make the equation as a whole to begin with. In this image, line 1 is the raw form of Kepler's equation where I substituted the integral in for A. After that I made simple substitutions and simplifications and got to the last line; at this point I am pretty sure there is no closed-form solution for θ:

http://gyazo.com/66665ac0d6b0821c0177da6d8d2fdc28

If nobody has any secret methods or anything that would help solve this, then I guess I'm moving on to Newton's method.

 

[tfr000]

You are in good company in trying to solve this... google "Kepler's problem".
There are some series expansions that are nice, but they don't converge for all eccentricities.

 

[Upsilon]

You are in good company in trying to solve this... google "Kepler's problem".
There are some series expansions that are nice, but they don't converge for all eccentricities.

Yeah, I would have done that had I wanted to. Unfortunately like I said I'm stubborn and I was trying my best to solve this non-iteratively. However after much frustration I have given up that notion. Unfortunately again, I'm still stubborn and want to solve it myself given Kepler's equations. It seems quite simple and straightforward to me, so I feel like I'm doing something wrong. I'm going to start with my destination; the angular position of the planet θ. Wikipedia gives the following solution for θ as a function of the eccentric anomaly E (e is the eccentricity of the orbit):
\theta = \arccos (\frac{\cos E - e}{1 - e \cos E})

That means I must find E. Ok, Kepler has his nifty little equation relating E to the mean anomaly M:
M = E - e\sin E

which unfortunately is not solvable for E. Before I can get into iterative solving, I must find M. Again, another nifty equation:
M = \sqrt{\frac{G(M_1 + M_2)}{a^3}}t

where G is the gravitational constant, M₁ and M₂ are the masses of the two bodies in the system in kilograms, a is the semi-major axis of the orbit in meters, and t is time elapsed since perihelion in seconds. So I find M, then I plug it back into Kepler's equation, set the equation equal to 0, and begin iterative solving. After I'm done and I have an approximate solution for E, I then simply plug this solution into the first equation I referenced and suddenly I have the planet's angular position relative to the sun as a function of time. Is it really that easy?

 

[tfr000]

Is it really that easy?

Pretty much. There are literally hundreds of ways of solving it, from most of the great mathematicians from Kepler until now. Most celestial mechanics books have a chapter devoted to its solution. What you wrote just above is how most people do it these days.

 

[Upsilon]

One last question, what do you call the angle between the semi-major axis and the x-axis of an arbitrary coordinate system? I already know that this is the orbital inclination when viewing the system from the side, but I'm talking about when viewing the system top-down. In other words, the angle that defines how the orbit's ellipse is rotated around the sun (since not all planets have their perihelion perfectly lined up with the x-axis).

 

[D H]

One last question, what do you call the angle between the semi-major axis and the x-axis of an arbitrary coordinate system?.

As far as I know, there is no name for that angle. There is a somewhat related concept, the longitude of periapsis. This is the sum of the longitude of ascending node and the argument of periapsis. These two angles are measured on two planes, so the longitude of periapsis is not an "angle"; it's instead a dogleg or compound angle.

 

[Upsilon]

As far as I know, there is no name for that angle. There is a somewhat related concept, the longitude of periapsis. This is the sum of the longitude of ascending node and the argument of periapsis. These two angles are measured on two planes, so the longitude of periapsis is not an "angle"; it's instead a dogleg or compound angle.

I was looking at this image on Wikipedia, and its definition for the argument of periapsis is the angle on the orbital plane between the periapsis and where the orbital plane intersects the reference plane. Is this not what I want? If you rotate the orbital plane counter-clockwise (in a relative way, I'm talking about looking at the orbital plane face-on and "spinning" it around its center), that increases the argument of periapsis, which sounds like what I'm looking for. Am I missing something?

 

[D H]

I was looking at this image on Wikipedia, and its definition for the argument of periapsis is the angle on the orbital plane between the periapsis and where the orbital plane intersects the reference plane. Is this not what I want?

No, it's not what you want (at least it's not what you asked about). You asked for "the angle between the semi-major axis and the x-axis of an arbitrary coordinate system". That would be the angle through the origin from the reference direction to the periapsis in that wikipedia diagram you referenced. That angle has no name. The argument of periapsis is the angle on the orbital plane from the ascending node to periapsis.

 

[Upsilon]

Alright, I'm going to try my best to explain what exactly I'm looking for. Modeling the orbit of a planet with the solution I posted above will give me the orbit flat against the plane of reference with the perihelion lined up with the reference direction (the positive x-axis). I use the orbital inclination given to then elevate the orbital plane some angle relative to the reference plane. The last step in modeling the orbit correctly is finding how the new orbital plane is rotated about its own focus, which is what I'm asking about.

Like I said, it seems to me the argument of periapsis is what I want based on the image. I first elevate the orbital plane i degrees relative to the reference plane to get the true orbital plane. Afterwards, I rotate this true orbital plane ω degrees to find how this inclined orbital plane is rotated relative to the sun.

EDIT: Just thought of a simpler way of explaining it; if the orbital plane and the reference plane happen to be the same, how much do I need to rotate the orbital plane counterclockwise about its focus to get the true orbit?

 

[D H]

In that case, yes, you are asking about the argument of periapsis.