back to "Positional Astronomy"   "Sun, Moon & Earth Applet"

Astronomical Algorithms

The motions of Earth and planets are usually computed in ecliptic coordinates, based on the plane of the ecliptic. The position of an object is defined by the ecliptic latitude (=0 for the sun), the ecliptic longitude, and the distance.

The figure represents the elliptical orbit of a body K, the Sun situated in the focus S:

We consider a fictitious body K' describing a circular orbit around S with constant velocity, with the same period as the real body K, and situated at P' at the instance when the real body is at the perihelion P. The angle PSK' is called mean anomaly M, increasing linearly with time. The problem consists in finding the true anomaly (angle PSK) at a given instant, when the mean anomaly M and the eccentricity of the ellipse are known.


Julian Day (valid from 1900/3/1 to 2100/2/28)

Julian day: 86400 s, Julian year: 365.25 d, Julian Century: 36525 d

double JulianDay (int date, int month, int year, double UT) {

if (month<=2) {month=month+12; year=year-1;}
return (int)(365.25*year) + (int)(30.6001*(month+1)) - 15 + 1720996.5 + date + UT/24.0

}


Solar Coordinates (according to: Jean Meeus: Astronomical Algorithms), accuracy of 0.01 degree

k = 2*PI/360;

//mean anomaly, degree
M = 357.52910 + 35999.05030*T - 0.0001559*T*T - 0.00000048*T*T*T

// mean longitude, degree
L0 = 280.46645 + 36000.76983*T + 0.0003032*T*T 

// Sun's equation of center
DL = (1.914600 - 0.004817*T - 0.000014*T*T)*sin(k*M)
+ (0.019993 - 0.000101*T)*sin(k*2*M) + 0.000290*sin(k*3*M)

// true longitude, degree
L = L0 + DL


convert ecliptic longitude L to right ascension RA and declination delta
(the ecliptic latitude of the Sun is assumed to be zero):

// number of Julian centuries since Jan 1, 2000, 12 UT
T = (JD-2451545.0) / 36525

// obliquity eps of ecliptic:
eps = 23.0 + 26.0/60.0 + 21.448/3600.0 - (46.8150*T + 0.00059*T*T - 0.001813*T*T*T)/3600

X = cos(L)
Y = cos(eps)*sin(L)
Z = sin(eps)*sin(L)
R = sqrt(1.0-Z*Z)

delta = arctan(Z/R) // in degrees
RA = (24/180)*arctan(Y/(X+R)) // in hours


compute sidereal time at Greenwich (according to: Jean Meeus: Astronomical Algorithms)

T = (JD - 2451545.0 ) / 36525

theta0 = 280.46061837 + 360.98564736629*(JD-2451545.0) + 0.000387933*T*T - T*T*T/38710000.0


convert tau, delta to horizon coordinates of the observer (altitude h, azimuth az):

sin h = sin(latitude)*sin(delta) + cos(latitude)*cos(delta)*cos(tau)

tan(az) = - sin(tau) / [cos(latitude)*tan(delta) - sin(latitude)*cos(tau)]


Home    back to "Positional Astronomy"   "Sun, Moon & Earth Applet"

 

Last update: 2019, Apr 25