/*
 *  SunRise.c - a program to calculate sunrise and sunset
 *  Amiga C version by Ralph Siemsen 30-Jun-92 @ Ottawa
 *
 *  This program calculates the time of sunrise, sunset, and twilight
 *  (astronomical, nautical, and civil) to within a minute or two
 *  during the second half of the twentieth century.
 *
 */

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

#define PI 3.1415926536
#define RISE 0
#define SET 1

double	CalcTime(), StoD(), DtoS(), Norm();
double	ApproxTime, SRA, SDec;
double	Latitude, Longtitude, TimeZone;
short	Month, Day, Year, DayOfYear, Mode;

main ()
{
  short i;  /* variable for local storage */

  /* Query the user for latitude,longtitude,etc. */
    GetInputs();

  /* Compute the day-of-year */
    i = (Month + 9) / 12;
    if (Year % 4) i*=2;		/* check for leap year */
    DayOfYear = (275 * Month)/9 + Day - i - 30;

  /* Rising Phenomena */
    Mode = RISE;
    Prep ();
    printf ("\n");
    printf ("  A. DAWN: %5.2f\n", CalcTime(-0.309017) );
    printf ("  N. DAWN: %5.2f\n", CalcTime(-0.207912) );
    printf ("  C. DAWN: %5.2f\n", CalcTime(-0.104528) );
    printf ("  SUNRISE: %5.2f\n", CalcTime(-0.0145439) );

  /* Setting Phenomena */
    Mode = SET;
    Prep ();
    printf ("\n");
    printf ("  SUNSET : %5.2f\n", CalcTime(-0.0145439) );
    printf ("  C. DUSK: %5.2f\n", CalcTime(-0.104528) );
    printf ("  N. DUSK: %5.2f\n", CalcTime(-0.207912) );
    printf ("  A. DUSK: %5.2f\n", CalcTime(-0.309017) );
}

GetInputs ()
{
    double temp;

    printf ("Latitude? ");
        scanf ("%lf",&temp);
        Latitude = 0.0174533 * StoD(temp);
    printf ("Longtitude? ");
        scanf ("%lf",&temp);
        Longtitude = 0.0174533 * StoD(temp);
    printf ("Time Zone? ");
        scanf ("%lf",&temp);
        TimeZone = 0.261799 * StoD(temp);
    printf ("Year? ");
        scanf ("%d",&Year);
    printf ("Month? ");
        scanf ("%d",&Month);
    printf ("Day? ");
        scanf ("%d",&Day);
}

Prep ()
{
    /* This routine initializes the global variables */
    /* ApproxTime, SRA, and SDec */

    double SMA;		/* Solar Mean Anomaly */
    double STL;		/* Solar True Longtitude */

    ApproxTime = DayOfYear + (((Mode ? 3*PI/2 : PI/2) + Longtitude) / (2*PI));

    SMA = ApproxTime * 0.017202 - 0.0574039;
    STL = SMA + 0.0334405*sin(SMA) + 0.000349066*sin(2*SMA) + 4.93289;
    STL = Norm(STL);

    SRA = atan (0.91746*tan(STL));
    /* Quadrant Adjustment: */
    if (STL > 3*PI/2) SRA = SRA + 2*PI;
    else if (STL > PI/2) SRA = SRA + PI;

    SDec = 0.39782*sin(STL);
    SDec = atan (SDec / sqrt (-SDec*SDec + 1));
}

double CalcTime (magic)
double magic;
{
    double S,T;
    
    S = magic - sin(SDec) * sin(Latitude);
    S = S / (cos(SDec) * cos(Latitude));
    
    if (fabs(S) <= 1) {
        S = PI/2 - atan(S / sqrt(-S*S+1));		/* Adjustment */
        if (Mode == RISE) S = 2*PI - S;
        T = S + SRA - 0.0172028*ApproxTime - 1.73364 + Longtitude - TimeZone;
	return (DtoS(3.81972*Norm(T)));
    }
    else return (0.0);
}

double Norm (value)
double (value);
{
    /* This routine adjusts its argument so that is in the range 0 - 2PI */

    while (value < 0) value = value + 2*PI;
    while (value >= 2*PI) value = value - 2*PI;
    return (value);
}

double StoD (value)
double value;
{
    /* convert the sexagesimal 'value' into decimal representation */
    /* ie: 12.30 (12 hours,30 minutes) becomes 12.5 hours */

    short hours,minutes,seconds;

    hours = value;			/* real->integer */
    value = (value - hours) * 100;
    minutes = value;			/* real->integer */
    value = (value - minutes) * 100;
    seconds = value;			/* real->integer */

    return (hours + minutes/60.0 + seconds/3600.0);
}

double DtoS (value)
double value;
{
    /* convert decimal into sexagesimal; opposite of StoD() */

    short hours,minutes;

    hours = value;			/* real->integer */
    value = (value - hours) * 60;
    minutes = value+0.5;		/* real->integer */

    return (hours + minutes/100.0);
}