Calculating elliptic orbits in Cartesian coordinates

[AbsentMinded]

I have a function to plot the orbits of planets based on their orbital elements (Semi-major Axis, Eccentricity, Argument of periapsis, Inclination, and longitude of ascending node). I have the x and y coordinates working great using only the semi-major axis, eccentricity, and argument of periapsis. Now I'm having trouble taking things into the third dimension with the inclination and longitude of ascending node. I know I need to get the distance from the longitude of the ascending node, and multiply it by sin(inclination), but I'm having trouble. Here's what I have so far:


apoapsis = (1+eccentricity)*SemiMajor;
periapsis = (1-eccentricity)*SemiMajor;
Semi-Minor = sqrt(apoapsis*periapsis);
SunFocus = sqrt((SemiMajor)^2-(SemiMinor)^2);
xc = SunFocus*cos(ArgOfPeriapsis);
yc = SunFocus*sin(ArgOfPeriapsis);
x = xc + SemiMajor*cos(time)*cos(ArgOfPeriapsis) - SemiMinor*sin(time)*sin(ArgOfPeriapsis);
y = yc + SemiMajor*cos(time)*sin(ArgOfPeriapsis) + SemiMinor*sin(time)*cos(ArgOfPeriapsis);
z = ?

This is psuedocode of a Matlab function. The inputs are the 5 orbital elements listed above (no mean anomaly for now). The orbits plot perfectly in two dimensions for every object I've thrown at it (plantes, dwarfs, comets, etc.), but I really want to get the inclination in there for 3 dimensions. Any help would be greatly appreciated.

 

[Philosophaie]

This works for Heliocentric positions where:

R,X,Y,Z-Heliocentric Distances
TA - True Anomaly
N - Longitude of the Ascending Node
w - Argument of the Perihelion

R = a * (1 - e ^ 2) / (1 + e * Cos(TA))
X = R * (Cos(N) * Cos(TA + w) - Sin(N) * Sin(TA+w)*Cos(i)
Y = R * (Sin(N) * Cos(TA+w) + Cos(N) * Sin(TA+w)) * Cos(i))
Z = R * Sin(TA+w) * Sin(i)

TA can be found using Mean Anomaly and eccentricity.

Correct approximations for Keplerian Elements can be found at the bottom of the site:

JPL Keplerian Elements

 

[Philosophaie]

This works for Heliocentric positions where:

R,X,Y,Z-Heliocentric Distances
TA - True Anomaly
N - Longitude of the Ascending Node
w - Argument of the Perihelion

R = a * (1 - e ^ 2) / (1 + e * Cos(TA))
X = R * (Cos(N) * Cos(TA + w) - Sin(N) * Sin(TA+w)*Cos(i)
Y = R * (Sin(N) * Cos(TA+w) + Cos(N) * Sin(TA+w)) * Cos(i))
Z = R * Sin(TA+w) * Sin(i)

TA can be found using Mean Anomaly and eccentricity.

Correct approximations for Keplerian Elements can be found at the bottom of the site:

http://ssd.jpl.nasa.gov/?planet_pos

 

[AbsentMinded]

Awesome, I'll give this a try, thank you.

 

[Philosophaie]

One other way is to find the Eccentric Anomaly, EA:

https://www.physicsforums.com/showthread.php?t=712636


Then converting it to True Anomaly:
$$\tan{\frac{TA}{2}} = \sqrt{{{1+e} \over {1-e}}} \tan{EA \over 2}.$$

 

[AbsentMinded]

This works for Heliocentric positions where:

R,X,Y,Z-Heliocentric Distances
TA - True Anomaly
N - Longitude of the Ascending Node
w - Argument of the Perihelion

R = a * (1 - e ^ 2) / (1 + e * Cos(TA))
X = R * (Cos(N) * Cos(TA + w) - Sin(N) * Sin(TA+w)*Cos(i)
Y = R * (Sin(N) * Cos(TA+w) + Cos(N) * Sin(TA+w)) * Cos(i))
Z = R * Sin(TA+w) * Sin(i)

TA can be found using Mean Anomaly and eccentricity.

Correct approximations for Keplerian Elements can be found at the bottom of the site:

http://ssd.jpl.nasa.gov/?planet_pos

That worked perfectly, thank you!http://imgur.com/5Zafkj1
http://imgur.com/5Zafkj1

 

[iontom]

That worked perfectly, thank you!http://imgur.com/5Zafkj1
http://imgur.com/5Zafkj1

How did you solve for the TA? Did you use Bessel Functions? Is there any easy way to approximate it?