/* Program sunpos.c  (V. 1.0.1)
   Written 1997 by Thomas Spahni
   Calculation of the apparent position of the Sun.
   This is a "low precision" calculation of the coordinates.
   Results are presented as:
   - Apparent Ecliptical Longitude
   - Local Hour Angle, westward from the South
   - Apparent Right Ascension
   - Apparent Declination
   - Azimuth from North via East
   - Elevation

   Additional data (Nutation and Apparent Greenwich Sidereal Time)
   can be displayed by uncommenting the block from program
   lines 457 to 466.

   Even if the output conversion routines format the display
   to fractions of seconds and arc seconds this is not the
   computational accuracy of the program.
   
   The program does not take into account any perturbations of the
   Orbit. The latitude is assumed to be zero, because it is always
   less than 1.2 seconds of arc.

   The precision of the other modules used to calculate the apparent
   position (nutation, obliquity of ecliptic, sidereal time) is much
   better, but this will not increase overall accuracy.
   
   For a way to achieve (much) higher precision see: 
   ftp://cdsarc.u-strasbg.fr/pub/cats/VI/81/
   
   To compile: use your editor to save the Makefile code at
   the end as a separate file;
   then type 'make'.
   
   For direct compilation with the GNU C-compiler enter:
   cc sunpos.c -o sunpos -lm

   */
   
/* --- Default Values; may be changed by user ---*/

/* The program requires Longitude and Latitude to
   calculate Azimuth and Elevation.
   These values will be used as a default when user input is == 0
   Eastern longitudes are negative!
   */
#define DEFAULT_LONGITUDE -8.551250
#define DEFAULT_LATITUDE  47.377306

/* --- END of user configurable default values ---*/


#include <stdio.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>


#define MY_PI   3.14159265358979323846
#define MINUS   45
#define BLANK   32

double torad = MY_PI / 180.0;

void getdata(int *year, int *month, int *day, int *hour, int *minutes,
             double *lng, double *lat)
{
   printf("Enter Date and Time of the day in UTC:\n");
   printf("Year ......... [ 4-digits YYYY ] ..............: ");
   scanf("%d", year);
   printf("Month ........ [ 1- or 2-digits 1-12 ] ........: ");
   scanf("%d", month);
   printf("Day .......... [ 1- or 2-digits 1-31 ] ........: ");
   scanf("%d", day);
   printf("Hour (UTC) ... [ 1- or 2-digits 0-23 ] ........: ");
   scanf("%d", hour);
   printf("Minutes ...... [ 1- or 2-digits 0-59 ] ........: ");
   scanf("%d", minutes);
   
   printf("\nEnter Geographic Position in Degrees and fractions:\n");
   printf("Geogr. Long. ..[ -ddd.nnnn (East=neg.!) ] .....: ");
   scanf("%lf", lng); 
   printf("Geogr. Lat.    [   dd.nnnn ] ..................: ");
   scanf("%lf", lat);
}


double fix(double n)
{
   if (n >= 0.0) {
      return floor(n);
      }
   else {
      return ceil(n);
      }
}
            

double normal(double angle) {
/* --- reduce to the range 0 ... 360 degrees --- */
   double rev;

   rev = angle / 360.0;
   rev = rev - fix(rev);
   if (rev < 0.0)
      rev += 1.0;
   return (rev*360.0);
}


double set_time_0(double jd) {
   return ( fix(jd + 0.5) - 0.5);
}
                           

double date_jd(int year, int month, int day) {

   int y, m, a, b;
   double jd, d;
         
   if (month > 2) {
      y = year;
      m = month;
      }
   else {
      y = year - 1;
      m = month + 12;
      }

   if ((year > 0) && ((year+month/100+day/10000) >= 1582.1015)) {
      a = y / 100;
      b = 2 - a + (a/4);
      }
   else {
      b = 0;
      }

   if (y < 0) {
      d = ceil(365.25 * (double) y - 0.75);
      }
   else {
      d = floor(365.25 * (double) y);
      }

   jd = d+floor(30.6001*(double)(m+1))+(double)day+1720994.5+(double)b;
   return jd;
}
                                                                                                                                    

/* Convert Julian Date to Julian Centuries */
double jd_t2000(double jd) {
   return ( (jd-2451545.0) / 36525.0 );
   }


double ecliptic( double jdtdt ) {
/* calculate mean obliquity of ecliptic */
   int i;
   double ecl;
   static double u[10];
   static double eclarg[] = {
      - 1.300258,   - 0.000430,   + 0.555347,   - 0.014272,
      - 0.069353,   - 0.010847,   + 0.001978,   + 0.007742,
      + 0.001608,   + 0.000681 };
                              
   /* make all powers of u */
   u[0] = jd_t2000(jdtdt) / 100.0;
   for(i=1; i<=9; i++) {
      u[i] = u[i-1] * u[0];
      }

   ecl = 23.439291;
   for(i=0; i<=9; i++) {
      ecl += eclarg[i] * u[i];
      }
   return ecl;
}
                                                                     
/* --- DATA FOR NUTATION ---*/
/* --- number of terms for Nutation: --- */
#define TERMS 63

static char argtbl[5][TERMS] = {
 { 0, -2,  0,  0,  0,  0, -2,  0,  0, -2, -2, -2,  0,  2,  0,  2,  0,  0,
  -2,  0,  2,  0,  0, -2,  0, -2,  0,  0,  2, -2,  0, -2,  0,  0,  2,  2,
   0, -2,  0,  2,  2, -2, -2,  2,  2,  0, -2, -2,  0, -2, -2,  0, -1, -2,
   1,  0,  0, -1,  0,  0,  2,  0,  2 },

 { 0,  0,  0,  0,  1,  0,  1,  0,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0,
   0,  0,  0,  0,  0,  0,  0,  0,  0,  2,  0,  2,  1,  0, -1,  0,  0,  0,
   1,  1, -1,  0,  0,  0,  0,  0,  0, -1, -1,  0,  0,  0,  1,  0,  0,  1,
   0,  0,  0, -1,  1, -1, -1,  0, -1 },

 { 0,  0,  0,  0,  0,  1,  0,  0,  1,  0,  1,  0, -1,  0,  1, -1, -1,  1,
   2, -2,  0,  2,  2,  1,  0,  0, -1,  0, -1,  0,  0,  1,  0,  2, -1,  1,
   0,  1,  0,  0,  1,  2,  1, -2,  0,  1,  0,  0,  2,  2,  0,  1,  1,  0,
   0,  1, -2,  1,  1,  1, -1,  3,  0 },

 { 0,  2,  2,  0,  0,  0,  2,  2,  2,  2,  0,  2,  2,  0,  0,  2,  0,  2,
   0,  2,  2,  2,  0,  2,  2,  2,  2,  0,  0,  2,  0,  0,  0, -2,  2,  2,
   2,  0,  2,  2,  0,  2,  2,  0,  0,  0,  2,  0,  2,  0,  2, -2,  0,  0,
   0,  2,  2,  0,  0,  2,  2,  2,  2 },

 { 1,  2,  2,  2,  0,  0,  2,  1,  2,  2,  0,  1,  2,  0,  1,  2,  1,  1,
   0,  1,  2,  2,  0,  2,  0,  0,  1,  0,  1,  2,  1,  1,  1,  0,  1,  2,
   2,  0,  2,  1,  0,  2,  1,  1,  1,  0,  1,  1,  1,  1,  1,  0,  0,  0,
   0,  0,  2,  0,  0,  2,  2,  2,  2 }
};

static float ks[2][TERMS] = {
 { -171996.0, -13187.0, -2274.0, +2062.0, +1426.0, +712.0, -517.0, -386.0,
      -301.0,   +217.0,  -158.0,  +129.0,  +123.0,  +63.0,  +63.0,  -59.0,
       -58.0,    -51.0,   +48.0,   +46.0,   -38.0,  -31.0,  +29.0,  +29.0,
       +26.0,    -22.0,   +21.0,   +17.0,   +16.0,  -16.0,  -15.0,  -13.0,
       -12.0,    +11.0,   -10.0,    -8.0,    +7.0,   -7.0,   -7.0,   -7.0,
        +6.0,     +6.0,    +6.0,    -6.0,    -6.0,   +5.0,   -5.0,   -5.0,
        -5.0,     +4.0,    +4.0,    +4.0,    -4.0,   -4.0,   -4.0,   +3.0,
        -3.0,     -3.0,    -3.0,    -3.0,    -3.0,   -3.0,   -3.0 },

 {    -174.2,     -1.6,    -0.2,    +0.2,    -3.4,   +0.1,   +1.2,   -0.4,
         0.0,     -0.5,     0.0,    +0.1,     0.0,    0.0,   +0.1,    0.0,
        -0.1,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,    -0.1,     0.0,   +0.1,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0 }
};

static float kc[2][TERMS] = {
 {  +92025.0,  +5736.0,  +977.0,  -895.0,   +54.0,   -7.0, +224.0, +200.0,
      +129.0,    -95.0,     0.0,   -70.0,   -53.0,    0.0,  -33.0,  +26.0,
       +32.0,    +27.0,     0.0,   -24.0,   +16.0,  +13.0,    0.0,  -12.0,
         0.0,      0.0,   -10.0,     0.0,    -8.0,   +7.0,   +9.0,   +7.0,
        +6.0,      0.0,    +5.0,    +3.0,    -3.0,    0.0,   +3.0,   +3.0,
         0.0,     -3.0,    -3.0,    +3.0,    +3.0,    0.0,   +3.0,   +3.0,
        +3.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0 },

 {      +8.9,     -3.1,    -0.5,    +0.5,    -0.1,    0.0,   -0.6,    0.0,
        -0.1,     +0.3,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0,    0.0,
         0.0,      0.0,     0.0,     0.0,     0.0,    0.0,    0.0 }
};


void nutation(double jdtdt, double *nutl, double *nuto) {

   int i, j;
   double t, t2, t3;
   double d, m, ms, f, k;
   double *w[5];
   double dphi, deps;
   double arg, wk;
   double torad;


   t  = jd_t2000(jdtdt);
   t2 = t * t;
   t3 = t2 * t;
   torad = MY_PI / 180.0;

   d  = normal(297.85036 + 445267.111480*t - 0.0019142*t2 + t3/189474.0);
   m  = normal(357.52772 +  35999.050340*t - 0.0001603*t2 - t3/300000.0);
   ms = normal(134.96298 + 477198.867398*t + 0.0086972*t2 + t3/ 56250.0);
   f  = normal( 93.27191 + 483202.017538*t - 0.0036825*t2 + t3/327270.0);
   k  = normal(125.04452 -   1934.136261*t + 0.0020708*t2 + t3/450000.0);

   w[0] = &d;  w[1]=&m; w[2]=&ms; w[3]=&f; w[4]=&k;
   dphi = 0.0;
   deps = 0.0;
   for (i=0; i < TERMS; i++) {
      arg = 0.0;
      for (j=0; j<5; j++) {
         if (argtbl[j][i] != 0) {
            arg += argtbl[j][i] * *w[j];
         }
      }
      wk    = normal(arg) * torad;
      dphi += sin(wk) * ( ks[0][i] + ks[1][i] * t);
      deps += cos(wk) * ( kc[0][i] + kc[1][i] * t);
   }
   *nutl = dphi / 36000000.0;
   *nuto = deps / 36000000.0;
}


/* --- calculate Greenwich mean sidereal time
 *
 * warning:
 * Julian day argument must express any instant of UNIVERSAL TIME
 */
    
 double jd_gmst(double jd) {
 double jd0, jdf;
 double t, t2, t3;
 double gt;
             
 jd0 = set_time_0(jd);
 jdf = jd - jd0;
 t   = jd_t2000(jd0);
 t2  = t * t;
 t3  = t * t2;
 gt  = 0.27905727326389+100.0021390378009*t+1.0775926e-6*t2 - 7.176e-11*t3;
 gt  += jdf * 1.00273790935;
 gt  = gt - fix(gt);
 if (gt < 0)
    gt += 1.0;
 return (gt * 360.0);
}

void dsphms(double time)
{
/* --- display a time in the format hh mm ss.fff
 *     input is in the form ddd.fffff range 0.0 ... 360.0
 */
  
   int hrs, min, sec, frc;
     
   if ( time < 0.0 )   {
      putchar(MINUS);
   } else   {
      putchar(BLANK);
   }
                          
   /* take |time|; round to units of 0.001 seconds of time: */
   time  = ( fix(fabs(time)*240000.0 + 0.5) + 0.5) / 3600000.0;
                                
   hrs = fix ( time );
   min = (int) fix( (time        - fix(time)        ) *  60.0);
   sec = (int) fix( (time*60.0   - fix(time*60.0)   ) *  60.0);
   frc = (int) fix( (time*3600.0 - fix(time*3600.0) ) *1000.0);
                                            
   printf("%02.2dh %0.2dm %0.2ds.%0.3d", hrs, min, sec, frc);
}


void dspdms(double time)
{
/* --- display a time in the format +/-ddd mm ss.ff
 *     input is in the form ddd.fffff range 0.0 ... +/-360.0
 */

   int deg, min, sec, frc;
     
   if ( time < 0.0 )  {
      putchar(MINUS);
   } else  {
      putchar(BLANK);
   }
                          
   /* take |time|; round to units of 0.01 seconds of arc: */
   time  = ( fix(fabs(time)*360000.0 + 0.5) + 0.5) / 360000.0;
                                
   deg = (int) fix(time);
   min = (int) fix( (time        - fix(time)        ) *  60.0);
   sec = (int) fix( (time*60.0   - fix(time*60.0)   ) *  60.0);
   frc = (int) fix( (time*3600.0 - fix(time*3600.0) ) * 100.0);

   printf("%03.3d\xb0 %02.2d\' %02.2d\".%02.2d", deg, min, sec, frc);
}


int main(void)
{
   double jd, t, t2, t3;
   double gml, mas, exc, eqc, wls, ecl, rek, dec, ome, lam;
   double mst, ast, hra, azi, ele;

   int day, month, year, hour, minutes;
   double lat, lng;
   double nutl, nuto, equ_of_equinox;
      
   printf("\nCalculation of the local position of the Sun:\n");
   printf(  "============================================\n\n");
   getdata(&year, &month, &day, &hour, &minutes, &lng, &lat);

   if ( lng == 0.0 ) {
      printf("Using default value ");
      printf("%f", DEFAULT_LONGITUDE);
      printf(" for longitude\n");
      lng = DEFAULT_LONGITUDE;
      }
   
   if ( lat == 0.0 ) {
      printf("Using default value ");
      printf("%f", DEFAULT_LATITUDE);
      printf(" for latitude\n");
      lat = DEFAULT_LATITUDE;
      }

   jd = date_jd(year, month, day);
   /* printf("%f\n", jd); */
   jd = jd + hour/24.0 + minutes/(24.0*60.0);

   t = jd_t2000(jd);
   /* printf("T = %.9f\n", t); */
   
   t2 = t * t;   t3 = t2 * t;
   
   /*--- Here begins the calculation of the position of the sun ---*/
   /* true geometric longitude referred to the mean equinox of the date */
   gml = 280.46645 + 36000.76983 * t + 0.0003032 * t2;
   /* printf("GML = %.5f\n", normal(gml)); */

   mas = 357.52910 + 35999.05030 * t - 0.0001559 * t2 - 0.00000048 * t3;
   /* printf("Mean Anomaly = %.5f\n", normal(mas)); */

   exc = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2;
   /* printf("Excentricity = %.9f\n", normal(exc)); */

   eqc =  (1.914600 - 0.004817 * t - 0.000014 * t2) * sin(normal(mas)*torad)
        + (0.019993 - 0.000101 * t) * sin(normal(2.0*mas) * torad)
        +  0.000290 * sin(normal(3.0*mas) * torad);

   wls = gml + eqc;
   /*--- Here ends the calculation of the longitude of the sun ---*/

   nutation(jd, &nutl, &nuto);

   /* Correction for Aberration and for Nutation: */
   lam = normal(wls - 0.00569 + nutl);
   

   /* ecl = ecliptic(jd) + 0.00256 * cos(ome * torad); */
   ecl = ecliptic(jd) + nuto;

   rek = atan2((cos(ecl*torad) * sin(lam*torad)), cos(lam*torad)) / torad;
   dec = asin(sin(ecl*torad) * sin(lam*torad));


   mst = jd_gmst(jd);
   equ_of_equinox = nutl * cos(ecl * torad);
   ast = normal(mst + equ_of_equinox);
   /* local hour angle of the object */
   hra = normal(mst - lng - rek);
      

   lat = lat * torad;
   azi = atan2( sin(hra*torad), 
               (cos(hra*torad) * sin(lat) - tan(dec) * cos(lat)));

   /* Navigators count Azimut from North via East */
   azi = normal(azi/torad+180.0);
   
   ele = asin(sin(lat) * sin(dec) + 
              cos(lat) * cos(dec) * cos(hra*torad)) / torad;

   printf("\n>>> Position of the Sun: <<<\n");

   /*
   printf("Nutation in Longitude ...............: ");
   dspdms(nutl); printf("\n");
   printf("Nutation in Latitude ................: ");
   dspdms(nuto); printf("\n");
   
   printf("Apparent Greenwich Sidereal Time ....: ");
   dsphms(ast); printf("\n");

   printf("Obliquity of the Ecliptic ...........: ");
   dspdms(ecl); printf("\n");
   */
      
   printf("Apparent Ecliptical Longitude .......: ");
   dspdms(lam); printf("\n");

   printf("Hour Angle, westward from the South .: ");
   dsphms(normal(hra)); printf("\n");

   printf("Apparent Right Ascension ............: ");
   dsphms(normal(rek)); printf("\n");

   printf("Apparent Declination ................: ");
   dspdms(dec/torad); printf("\n");

   printf("Azimuth from North via East .........: ");
   dspdms(azi); printf("\n");
   
   printf("Elevation ...........................: ");
   dspdms(ele); printf("\n");
   
   printf("\n");
   return (EXIT_SUCCESS);
}

/* --- The Code for the Makefile, extract and save as 'Makefile'
---> cut here <---
# Makefile for sunpos  
#
# to make the program just type 'make'
# to install in /usr/local/bin type 'make install'
# edit INSTALL_DIR as needed
WARNINGS         = 
COMPILE          = -DLINUX
DEBUG            =
LNFLAGS          = -lm
# - LINUX uses gcc
CC               = gcc

# Where to install the executable:
INSTALL_DIR      = /usr/local/bin


CFLAGS = $(COMPILE) $(DEBUG) $(WARNINGS)

OFILES           = sunpos.o

HFILES           = 

.c.o:
	$(CC) -c $(CFLAGS) $*.c

sunpos:         $(OFILES) $(HFILES)
		$(CC) $(CFLAGS) -o sunpos $(OFILES) $(LNFLAGS)


.PHONY: clean
clean :
	-rm -f sunpos $(OFILES)
	-rm -f *.?~
	-rm -f *~
	-rm -f ?akefile~

install: sunpos
	install ./sunpos $(INSTALL_DIR)

# */