The position of the Sun

The position of the Sun is computed just like the position of any other planet, but since the Sun always is moving in the ecliptic, and since the eccentricity of the orbit is quite small, a few simplifications can be made. Therefore, a separate presentation for the Sun is given.

Of course, we're here really computing the position of the Earth in its orbit around the Sun, but since we're viewing the sky from an Earth-centered perspective, we'll pretend that the Sun is in orbit around the Earth instead.

First, compute the eccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):

E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )

or (if E and M are expressed in radians):

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

Note that the formulae for computing E are not exact; however they're accurate enough here.

Then compute the Sun's distance r and its true anomaly v from:

xv = r * cos(v) = cos(E) - e
yv = r * sin(v) = sqrt(1.0 - e*e) * sin(E)

v = atan2( yv, xv )
r = sqrt( xv*xv + yv*yv )

(note that the r computed here is later used as rs)

atan2() is a function that converts an x,y coordinate pair to the correct angle in all four quadrants. It is available as a library function in Fortran, C and C++. In other languages, one has to write one's own atan2() function. It's not that difficult:

atan2( y, x ) = atan(y/x) if x positive
atan2( y, x ) = atan(y/x) +- 180 degrees if x negative
atan2( y, x ) = sign(y) * 90 degrees if x zero
Now, compute the Sun's true longitude:

lonsun = v + w

Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys:

xs = r * cos(lonsun)
ys = r * sin(lonsun)

(since the Sun always is in the ecliptic plane, zs is of course zero). xs,ys is the Sun's position in a coordinate system in the plane of the ecliptic. To convert this to equatorial, rectangular, geocentric coordinates, compute:

xe = xs
ye = ys * cos(ecl)
ze = ys * sin(ecl)

Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):

RA = atan2( ye, xe )
Dec = atan2( ze, sqrt(xe*xe+ye*ye) )