I Position and Velocity in Heliocentric Ecliptic Coordinates

There have been many questions on this forum about celestial mechanics in general, and concerning position and velocity in an orbit in particular. So I offer this post as a summary and reference.

Here's a method for finding heliocentric position and sun-relative velocity in ecliptic coordinates from the Keplerian orbital elements (a, e, i, Ω, ω, T) at a time of interest, t, for both elliptical and hyperbolic orbits relative to the sun.

a : semimajor axis
* here usually in astronomical units

e : eccentricity

i : inclination to the ecliptic
* entered in degrees, usually converted to radians

Ω : longitude of the ascending node
* entered in degrees, usually converted to radians
* Ω is an angle, subtended at the sun, measured in the ecliptic plane in the counterclockwise direction from the perspective of someone positioned far along the North Ecliptic Pole, from a ray extending toward the Vernal Equinox to another ray extending toward the place where the orbit crosses the ecliptic going from south to north.

ω : argument of the perihelion
* entered in degrees, usually converted to radians
* ω is an angle, subtended at the sun, measured in the plane of the orbit, in the angular direction of motion, from a ray extending toward the ascending node to another ray extending toward the perihelion of the orbit.

T : time of perihelion passage
* here usually in Julian date, with differences in time being usually in days

.
For ELLIPTICAL ORBITS.

An elliptical orbit has a positive semimajor axis and an eccentricity strictly between zero and one. It is a bound orbit with a period. The planets, the asteroids, and some comets are in elliptical orbits.
The period of the orbit, P.

The eccentric anomaly, u, will be in radians.
* u is an angle, subtended at the geometric center of the elliptical orbit (not at the sun), measured in the plane of the orbit in the direction of motion, from a ray extending through the orbit's perihelion, to another ray extending through a point on a circle that circumscribes the orbital ellipse which is in the +y direction of the position of the body in its orbit at time t.

The canonical heliocentric position vector [x''',y''',z'''].

x''' = a (cos u − e)
y''' = a sin u √(1− e²)
z''' = 0

The canonical coordinate system has the sun at the origin, the XY plane in the plane of the orbit, with the orbit's perihelion being on the +x axis, and the +z axis such that an observer somewhere along it would observe a counterclockwise direction of motion.

The true anomaly, θ.

θ = arctan( y''' , x''' )

* θ is an angle, subtended at the sun, measured in the plane of the orbit, in the angular direction of motion, from a ray extending toward the perihelion of the orbit to another ray extending toward the orbiting body's position at time t.

Beginning with a position vector and a velocity vector in heliocentric ecliptic coordinates, at a time t, for an object in orbit around the sun, you can derive its Keplerian orbital elements (a, e, i, Ω, ω, T). The combination of the vectors of position and of velocity in an orbit is sometimes called the "state vector."

If you are given the position in spherical ecliptic coordinates (r = distance, λ = ecliptic longitude, β = ecliptic latitude), then do this to get the position in rectangular ecliptic coordinates:

x = r cos λ cos β
y = r sin λ cos β
z = r sin β

Let's assume that x, y, and z are in meters, and that Vx, Vy, Vz are the components of the sun-relative velocity in ecliptic coordinates in meters per second.

Note: the math to this point works for either elliptical or hyperbolic orbits. However, some of the stuff to follow works only for elliptical orbits. You'll know that the orbit is a hyperbola if (1) the semimajor axis is negative and (2) the eccentricity is greater than one.

Given the heliocentric positions of Earth [x₀,y₀,z₀] and of another planet in ecliptic coordinates [x,y,z] for the same moment of time, t, expressed as a Julian Date, we will find the geocentric position of the planet in celestial coordinates; i.e., we want the planet's right ascension and declination.

dx = x − x₀
dy = y − y₀
dz = z − z₀

We find the obliquity of Earth, ε, at time t.

T = t − 2451545.0

ε = 23.4392911° − 3.562266e-7 T − 1.22848e-16 T² + 1.03353e-20 T³

This equation gives the obliquity in degrees.

dx' = dx
dy' = dy cos ε − dz sin ε
dz' = dy sin ε + dz cos ε

The distance between Earth and the planet, dr, at time t is found from

dr = √{ (dx)² + (dy)² + (dz)² }

The planet's geocentric right ascension in decimal hours, α, at time t is found from

α = (12/π) Arctan( dy' , dx' )

The planet's geocentric declination in decimal degrees, δ, at time t is found from