SUN RA DEC

From Tauwiki

void SUN_RA_DEC(double t_ut, double day, double month, double year, double *ra, double *dec) { double d_0, t_0, s_0, d, t, L_0, M_0, C, L_S, K, tan_r_s, R_S, sin_dec, X, Y; /*Number of days and centuries from J2000*/ d_0 = 367 *year - floor((7./4.)*(year + floor((month + 9)/12))) + floor(275 * month/9) + day - 730531.5; t_0 = d_0/36525; /*Sidereal time*/ s_0 = 6.6974 + 2400.0513*t_0; /*Solar Coordinates*/ d = d_0 +t_ut/24;/*Number of days from J2000*/ t = d/36525; /*Number of centruries from J2000*/ L_0 = 280.466 + 36000.770 * t;/*Sun's mean longitude*/ /*Bring L_0 back down to the proper range*/ if (L_0 < 0) while (L_0 < 0) L_0 += 360; else while (L_0 > 360) L_0 -=360; M_0 = 357.529 + 35999.050 * t;/*Mean anomaly*/ /*Bring M_0 back down to the proper range*/ if (M_0 < 0) while (M_0 < 0) M_0 += 360; else while (M_0 > 360) M_0 -=360; C = (1.915 - 0.005*t)*sin(M_0 * RAD_PER_DEG) + 0.020*sin(2*M_0 * RAD_PER_DEG);/*Equation of Center*/ L_S = L_0 + C;/*True ecliptic longitude of Sun*/ K = 23.439 - 0.013*t;/*Obliquity of ecliptic*/ X = cos(L_S * RAD_PER_DEG); Y = cos(K * RAD_PER_DEG) * sin(L_S * RAD_PER_DEG); tan_r_s = Y/X; /*Sun RA*/ R_S = atan(tan_r_s)/RAD_PER_DEG;/*RA of Sun*/ if (X < 0) /*Make sure we're in the right quadrant*/ *ra = R_S + 180; else if ((X > 0) && (Y < 0)) *ra = R_S + 360; else *ra = 360 + R_S; sin_dec = sin(L_S * RAD_PER_DEG) * sin(K * RAD_PER_DEG); *dec = asin(sin_dec)/RAD_PER_DEG; } /*End calculate solar coordinates*/