/*--------------------------------------------------------------------------*/ /* Ephemerides of the Sun */ /* Algorithm based on */ /* T.C. van Flandern, K.F. Pulkkinen: Low-Precision Formulae */ /* for Planetary Positions, */ /* Astrophysical Journal Supplement, vol 41, #2, 391-411, November 1979 */ /* Accuracy about 1' for 1700-2200 */ /* In Linux compile with */ /* gcc -o sun sun.c -lm */ /* Usage: */ /* sun year month day hour minute second */ /*--------------------------------------------------------------------------*/ #include #include #include #include #include #include #define two_pi 6.28318531 double julian (int y, int m, int d) ; void COORDS(int OBJ, double V, double U, double W, double XL, double *RA, double *DEC) ; void SUN (double T, double *ra, double *dec) ; double A(double B, double C, double D) ; void INITS(double JS) ; double C[35] ; main(int argc, char **argv) { double ra,dec,jd,T, dist ; int i,j, vuosi, kuu, paiva, tunti, minuutti, sekunti, d,m,s,sig ; if (argc != 7) { fprintf(stderr,"usage: sun year month day hour minute second \n") ; exit(1) ; } vuosi = atoi(argv[1]) ; kuu = atoi(argv[2]) ; paiva = atoi(argv[3]) ; tunti = atoi(argv[4]) ; minuutti = atoi(argv[5]) ; sekunti = atoi(argv[6]) ; jd=julian (vuosi, kuu, paiva) ; jd += (tunti + minuutti/60.0 + sekunti/3600.0)/24.0 ; INITS (jd) ; T=jd/36525.0 + 1.0 ; SUN (T, &ra, &dec) ; d=(int) ra ; ra = 60.0*(ra - d) ; m = (int) ra ; s = (int)(60.0*(ra-m)) ; printf("ra = %3d h %2d min %2d s\n", d,m,s) ; if (dec < 0.0) {sig=-1; dec = -dec;} else sig=1; d=(int) dec ; dec = 60.0*(dec - d) ; m = (int) dec ; s = (int)(60.0*(dec-m)) ; printf("d = %3d deg %2d' %2d''\n", sig*d,m,s) ; } /*--------------------------------------------------------------------------*/ /* julian computes modified julian date */ /* zero point on January 1.5, 2000 */ /*--------------------------------------------------------------------------*/ double julian (int y, int m, int d) { if (m <= 2) { y-- ; m+=12 ; } return(floor(365.25*y) + floor(30.6001*(m+1)) + d + (1720994.5 - 2451545.0) +2.0 - y / 100 + y / 400) ; } /*--------------------------------------------------------------------------*/ /* COORDS converts the quantities U,V,W to equatorial coordinates */ /*--------------------------------------------------------------------------*/ void COORDS(int OBJ, double V, double U, double W, double XL, double *RA, double *DEC) { double x,y,z ; x=V ; y=U; z=W ; W=W/sqrt(U-V*V) ; V=V/sqrt(U) ; if (abs(W) > 1.0) { printf (" ===> W/sqrt(U-V*V)= %f\n",W) ; if (W<0.0) W=-1.0 ; else W=1.0 ; } if (abs(V) > 1.0) { printf (" ===>V/sqrt(U)= %f\n",V) ; if (V<0.0) V=-1.0 ; else V=1.0 ; } *RA=XL+asin(W) ; *DEC=asin(V) ; if (*RA < 0.0) (*RA)+=two_pi ; *RA=(*RA)*24.0/two_pi ; *DEC=(*DEC)*360.0/two_pi ; } /*--------------------------------------------------------------------------*/ /* Sun */ /*--------------------------------------------------------------------------*/ void SUN(double T, double *RA, double *DEC) { double U,V,W ; V= 0.39785*sin(C[7]) -0.01000*sin(C[7]-C[8]) +0.00333*sin(C[7]+C[8]) -0.00021*sin(C[7])*T +0.00004*sin(C[7]+2.0*C[8]) -0.00004*cos(C[7]) -0.00004*sin(C[5]-C[7]) +0.00003*sin(C[7]-C[8])*T ; U= 1.0000 -0.03349*cos(C[8]) -0.00014*cos(C[8]*2.0) +0.00008*cos(C[8])*T -0.00003*sin(C[8]-C[19]) ; W= 0.03211*sin(C[8]) -0.04129*sin(C[7]*2.0) +0.00104*sin(2.0*C[7]-C[8]) -0.00035*sin(2.0*C[7]+C[8]) -0.00010 -0.00008*sin(C[8])*T -0.00008*sin(C[5]) +0.00007*sin(2.0*C[8]) +0.00005*sin(2.0*C[7])*T +0.00003*sin(C[1]-C[7]) -0.00002*cos(C[8]-C[19]) +0.00002*sin(4.0*C[8]-8.0*C[16]+3.0*C[19]) -0.00002*sin(C[8]-C[13]) -0.00002*cos(2.0*C[8]-2.0*C[13]) ; COORDS(0,V,U,W,C[7],RA,DEC) ; } /*--------------------------------------------------------------------------*/ /* INITS initializes tabel C for the given date */ /* JS=modified julian date, JD-2451545.0 */ /*--------------------------------------------------------------------------*/ void INITS(double JS) { C[1]=A(0.606434,0.03660110129,JS) ; C[5]=A(0.347343,-0.00014709391,JS)-two_pi; C[7]=A(0.779072,0.00273790931,JS) ; C[8]=A(0.993126,0.00273777850,JS) ; C[13]=A(0.140023,0.00445036173,JS); C[16]=A(0.053856,0.00145561327,JS); C[19]=A(0.056531,0.00023080893,JS); } /*--------------------------------------------------------------------------*/ double A(double B, double C, double D) { double E,sign,i ; E=B+C*D ; sign=(E >= 0.0 ? 1.0 : -1.0) ; return(sign*modf(fabs(E),&i)*6.283185307) ; }