/*--------------------------------------------------------------------------*/
/* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <float.h>

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