Calculating the RA, Dec, Alt and Azimuth of Sun

I think this website is useful: http://www.jgiesen.de/SME/details/basics/index.htm

I've created a program on RA, Dec, Alt and Azimuth of Sun. Please go to http://david.sam-siu.com/sunrise4.exe to have a copy.

The calculation of those parameters involves 4 major items, including
1. Conversion of date and time
2. Astronomical algorithms:
3. Conversion to RA, Declin.
4. Conversion to Alt, Azimuth

Some details are given here:
Julian Day (valid from 1900/3/1 to 2100/2/28)

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


k = 2*PI/360;

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

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

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);

L = L0 + DL; // true longitude, degree



convert ecliptic longitude L to right ascension RA and declination delta

X = cos(L); Y = cos(eps)*sin(L); Z = sin(eps)*sin(L); R = Math.sqrt(1.0-Z*Z);
eps = 23.43999; // obliquity of ecliptic

delta = (180/PI)*arctan(Z/R); // in degrees

RA = (24/PI)*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 to altitude h, azimuth az

sin h = sin beta sin delta + cos beta cos delta cos tau
tan az = (- sin tau) / (cos beta tan delta - sin beta cos tau)