How to determine the sunrise and sunset times

Couple days ago I mentioned about calculating the sunrise and sunset times. It is an interesting problem so after reading Astronomy Answers - Position of the Sun I quickly implemented it using C language. It is only an experiment so results are approximate.

What can be easily calculated?

Sunset/Sunrise

Sunrise and sunset calculations for 5/2013, Gdansk/Poland (GMT=1, DST=1)

  Day        JDN   Sunrise   Transit    Sunset
    1    2456414      5:08     12:44     20:19
    2    2456415      5:06     12:44     20:21
    3    2456416      5:04     12:43     20:23
    4    2456417      5:02     12:43     20:25
    5    2456418      5:00     12:43     20:26
    6    2456419      4:58     12:43     20:28
    7    2456420      4:56     12:43     20:30
    8    2456421      4:54     12:43     20:32
    9    2456422      4:52     12:43     20:34
   10    2456423      4:51     12:43     20:36
   11    2456424      4:49     12:43     20:37
   12    2456425      4:47     12:43     20:39
   13    2456426      4:45     12:43     20:41
   14    2456427      4:43     12:43     20:43
   15    2456428      4:42     12:43     20:45
   16    2456429      4:40     12:43     20:46
   17    2456430      4:38     12:43     20:48
   18    2456431      4:37     12:43     20:50
   19    2456432      4:35     12:43     20:51
   20    2456433      4:34     12:43     20:53
   21    2456434      4:32     12:43     20:55
   22    2456435      4:31     12:43     20:56
   23    2456436      4:29     12:44     20:58
   24    2456437      4:28     12:44     20:59
   25    2456438      4:27     12:44     21:01
   26    2456439      4:26     12:44     21:02
   27    2456440      4:24     12:44     21:04
   28    2456441      4:23     12:44     21:05
   29    2456442      4:22     12:44     21:07
   30    2456443      4:21     12:44     21:08
   31    2456444      4:20     12:45     21:09

Astronomical seasons

Astronomical seasons 2014, Gdansk/Poland (northern hemisphere)

  Day        JDN
   81    2456739  first day of spring
  174    2456832  first day of summer
  267    2456925  first day of autumn
  357    2457015  first day of winter

Day and night

Day and night, Gdansk/Poland, UTC

      Date & time                    Day and night (one character = 6 degrees)
 01/01/2013 00:00 +++++++++-------------------------------------------++++++++
 01/01/2013 06:00 -------------------------------------+++++++++++++++++------
 01/01/2013 12:00 ----------------------+++++++++++++++++---------------------
 01/01/2013 18:00 -------+++++++++++++++++------------------------------------

 21/05/2013 00:00 ++++++++++++++++++++--------------------++++++++++++++++++++
 21/05/2013 06:00 +++++--------------------+++++++++++++++++++++++++++++++++++
 21/05/2013 12:00 ----------++++++++++++++++++++++++++++++++++++++++----------
 21/05/2013 18:00 +++++++++++++++++++++++++++++++++++--------------------+++++

First step - JDN

At first you need to calculate Julian Day Number. It is the first and most important step but very easy to implement.

/* source: http://aa.quae.nl/en/reken/zonpositie.html */

/* convert Gregorian calendar date to Julian Day Number */
/* source: http://en.wikipedia.org/wiki/Julian_day      */
int julian_jdn(int day, int month, int year) {
  int a, y, m;

  a = (14 - month) / 12;
  y = year + 4800 - a;
  me = month + 12 * a - 3;
  return (day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045);
}

If you want to convert JDN to the Gregorian date (d/m/y) then read my previous blog post.

Second step - Sun calculations

In the next step I will use functions defined below.

/* source: http://aa.quae.nl/en/reken/zonpositie.html
 * milosz at sleeplessbeastie.eu
 */

/* convert degrees to radians */
double sun_to_rad(double angle)
{
    return (M_PI * angle / 180);
}

/* convert radians to degrees */
double sun_to_deg(double angle)
{
    return (180 * angle / M_PI);
}

/* draw a full circle (370' to (370'- 360') = 10') */
double sun_full_circle(double angle)
{
    return (angle - (int)(angle / 360) * 360);
}

/* Earth - mean anomaly */
double sun_earth_mean_anomaly(double jdn)
{
    return sun_full_circle(357.5291 + 0.98560028 * (jdn - 2451545));
}

/* Earth - the equation of center */
double sun_earth_equation_of_center(double mean_anomaly)
{
    return 1.9148 * sin(sun_to_rad(mean_anomaly)) +
           0.0200 * sin(2 * sun_to_rad(mean_anomaly)) +
           0.0003 * sin(3 * sun_to_rad(mean_anomaly));
}

/* Earth - true anomaly */
double sun_earth_true_anomaly(double jdn)
{
    return sun_earth_mean_anomaly(jdn) +
           sun_earth_equation_of_center(sun_earth_mean_anomaly(jdn));
}

/* Ecliptic longitude of the sun */
double sun_sun_ecliptic_longitude(double true_anomaly)
{
    return sun_full_circle(true_anomaly + 102.9372 + 180);
}

/* check for specific season using ecliptic longitude of the sun*/
/* for northern hemisphere */
/* 0   - 90    spring */
/* 90  - 180   summer */
/* 180 - 270   autumn */
/* 270 - 360   winter */
int sun_check_season(int jdn, int a, int b)
{
    double el_curr = sun_sun_ecliptic_longitude(sun_earth_true_anomaly(jdn));
    double el_prev = sun_sun_ecliptic_longitude(sun_earth_true_anomaly(jdn - 1));


    if ((el_prev <= b && el_curr >= b) || (b == 360 && el_prev >= a && el_prev <= b && el_curr >= 0 && el_curr < a))
        return 0;
    else if ((el_prev >= a && el_prev < b) && (el_prev >= a && el_prev < b))
        return 1;
    else
        return -1;
}

/* check if JDN is spring */
/*  1 - true              */
/*  0 - true, last day    */
/* -1 - false             */
int sun_is_spring(int jdn)
{
    return sun_check_season(jdn, 0, 90);
}

/* check if JDN is summer */
/*  1 - true              */
/*  0 - true, last day    */
/* -1 - false             */
int sun_is_summer(int jdn)
{
    return sun_check_season(jdn, 90, 180);
}

/* check if JDN is autumn */
/*  1 - true              */
/*  0 - true, last day    */
/* -1 - false             */
int sun_is_autumn(int jdn)
{
    return sun_check_season(jdn, 180, 270);
}

/* check if JDN is spring */
/*  1 - true              */
/*  0 - true, last day    */
/* -1 - false             */
int sun_is_winter(int jdn)
{
    return sun_check_season(jdn, 270, 360);
}

/* Earth - right ascension */
double sun_earth_right_ascension(double ecliptic_longitude)
{
    return sun_to_deg(
               atan2(
                   sin(sun_to_rad(ecliptic_longitude)) * cos(sun_to_rad(23.45)),
                   cos(sun_to_rad(ecliptic_longitude))
               )
           );
}

/* Earth - declination */
double sun_earth_declination(double ecliptic_longitude)
{
    return sun_to_deg(
               asin(
                   sin(sun_to_rad(ecliptic_longitude)) * sin(sun_to_rad(23.45))
               )
           );
}

/* Earth - sidereal time */
double sun_earth_sideral_time(double jdn, double lw)
{
    return sun_full_circle(
               280.1600 + 360.9856235 * (jdn - 2451545) - lw
           );
}

/* Earth - hour angle */
double sun_earth_hour_angle(double sideral_time, double right_ascension)
{
    return sideral_time - right_ascension;
}

/* Earth - azimuth */
double sun_earth_azimuth(double hour_angle, double declination,  double ln)
{
    return sun_full_circle(
               sun_to_deg(
                   atan2(
                       sin(sun_to_rad(hour_angle)),
                       cos(sun_to_rad(hour_angle)) * sin(sun_to_rad(ln)) -
                       tan(sun_to_rad(declination)) * cos(sun_to_rad(ln))
                   )
               )
           );
}

/* Earth - altitude */
double sun_earth_altitude(double hour_angle, double declination, double ln)
{
    return sun_full_circle(
               sun_to_deg(
                   asin(
                       sin(sun_to_rad(ln)) * sin(sun_to_rad(declination)) +
                       cos(sun_to_rad(ln)) * cos(sun_to_rad(declination)) *
                       cos(sun_to_rad(hour_angle))
                   )
               )
           );
}

/* Sun -  transtit */
double sun_sun_transit(double jdn, double lw, double mean_anomaly)
{
    return  2451545 + 0.0009 + (lw / 360) +
            (int)((jdn - 2451545 - 0.0009) / 1 - lw / 360) +
            0.0053 * sin(sun_to_rad(mean_anomaly)) -
            0.0069 * sin(sun_to_rad(2 * sun_full_circle(mean_anomaly + 102.9372 + 180)));
}

/* Sun - sunrise */
double sun_sun_sunrise(double ln, double declination, double jdn, double lw, double mean_anomaly)
{
    double h;

    h = -sun_full_circle(
            sun_to_deg(
                acos(
                    (
                        sin(sun_to_rad(-0.83)) -
                        sin(sun_to_rad(ln)) * sin(sun_to_rad(declination))
                    ) / (
                        cos(sun_to_rad(ln)) * cos(sun_to_rad(declination))
                    )
                )
            )
        );

    return 2451545 + 0.0009 + (h + lw) / 360 +
           (int)((jdn - 2451545 - 0.0009) / 1 - (h + lw) / 360) +
           0.0053 * sin(sun_to_rad(mean_anomaly)) -
           0.0069 * sin(sun_to_rad(2 * sun_full_circle(mean_anomaly + 102.9372 + 180)));

}

/* Sun - sunset */
double sun_sun_sunset(double ln, double declination, double jdn, double lw, double mean_anomaly)
{
    double h;

    h = sun_full_circle(
            sun_to_deg(
                acos(
                    (
                        sin(sun_to_rad(-0.83)) -
                        sin(sun_to_rad(ln)) * sin(sun_to_rad(declination))
                    ) / (
                        cos(sun_to_rad(ln)) * cos(sun_to_rad(declination))
                    )
                )
            )
        );

    return 2 + 2451545 + 0.0009 + (h + lw) / 360 +
           (int)((jdn - 2451545 - 0.0009) / 1 - (h + lw) / 360) +
           0.0053 * sin(sun_to_rad(mean_anomaly)) -
           0.0069 * sin(sun_to_rad(2 * sun_full_circle(mean_anomaly + 102.9372 + 180)));
}

Third step - Example

Sunrise/Sunset

Source code:

    double jdn; /* Julian Day Number */

    double ln, lw; /* position */
    ln =  54.4;    /* north latitude */
    lw = -18.5;    /* west longitude */

    /* time modifier */
    int gmt = +1;
    int dst = +1; /* it is very simplified as this is only an example */
    int time_mod = gmt + dst;

    int year = 2013, month = 6;
    double ta, l, delta;
    double transit, sunrise, sunset;
    double temp;
    int i;

    printf("\nSunrise and sunset calculations for %d/%d (GMT=%d, DST=%d)\n\n", month, year, gmt, dst);
    printf("%5s %10s %9s %9s %9s\n", "Day", "JDN", "Sunrise", "Transit", "Sunset");

    for (i = julian_jdn(1, month, 2013); i < julian_jdn(1, month + 1, 2013); i++) {
        ta    = sun_earth_true_anomaly(i);
        l     = sun_sun_ecliptic_longitude(ta);
        delta = sun_earth_declination(l);

        printf("%5d %10d ", i + 1 - julian_jdn(1, month, 2013), i);

        sunrise = sun_sun_sunrise(ln, delta, i, lw, ta);
        temp    = (sunrise - (int) sunrise) > 0.5 ? (sunrise - (int) sunrise) - 0.5 : (sunrise - (int) sunrise) + 0.5;
        printf("%6d:%02d ", (int)(temp * 24) + time_mod, (int)(temp * 24 * 60) % 60);

        transit = sun_sun_transit(i, lw, ta);
        temp    = (transit - (int) transit) > 0.5 ? (transit - (int) transit) - 0.5 : (transit - (int) transit) + 0.5;

        printf("%6d:%02d ", (int)(temp * 24) + time_mod, (int)(temp * 24 * 60) % 60);

        sunset  = sun_sun_sunset(ln, delta, i, lw, ta);
        temp    = (sunset - (int) sunset) > 0.5 ? (sunset - (int) sunset) - 0.5 : (sunset - (int) sunset) + 0.5;
        printf("%6d:%02d\n", (int)(temp * 24) + time_mod , (int)(temp * 24 * 60) % 60);
    }

Output:

Sunrise and sunset calculations for 6/2013 (GMT=1, DST=1)

  Day        JDN   Sunrise   Transit    Sunset
    1    2456445      4:19     12:45     21:10
    2    2456446      4:18     12:45     21:12
    3    2456447      4:17     12:45     21:13
    4    2456448      4:16     12:45     21:14
    5    2456449      4:16     12:45     21:15
    6    2456450      4:15     12:46     21:16
    7    2456451      4:14     12:46     21:17
    8    2456452      4:14     12:46     21:18
    9    2456453      4:13     12:46     21:19
   10    2456454      4:13     12:46     21:20
   11    2456455      4:12     12:46     21:21
   12    2456456      4:12     12:47     21:21
   13    2456457      4:12     12:47     21:22
   14    2456458      4:11     12:47     21:23
   15    2456459      4:11     12:47     21:23
   16    2456460      4:11     12:47     21:24
   17    2456461      4:11     12:48     21:24
   18    2456462      4:11     12:48     21:25
   19    2456463      4:11     12:48     21:25
   20    2456464      4:11     12:48     21:25
   21    2456465      4:11     12:49     21:26
   22    2456466      4:12     12:49     21:26
   23    2456467      4:12     12:49     21:26
   24    2456468      4:12     12:49     21:26
   25    2456469      4:13     12:49     21:26
   26    2456470      4:13     12:50     21:26
   27    2456471      4:14     12:50     21:26
   28    2456472      4:14     12:50     21:26
   29    2456473      4:15     12:50     21:25
   30    2456474      4:16     12:50     21:25

Notes

You need to note that Julian Day 2456474 starts at 2456473.5 (00:00) and ends at 2456474.5. 2456474.0 is equal to 30/06/2013 12:00.

temp  = (sunrise - (int) sunrise) > 0.5 ? (sunrise - (int) sunrise) - 0.5 : (sunrise - (int) sunrise) + 0.5;
printf("%6d:%02d ", (int)(temp * 24) + time_mod, (int)(temp * 24 * 60) % 60);

It is easier to determine time using Julian Day fraction part when above range is shifted.

You need to apply GMT and DST as by default sunrise/transit/sunset/... functions return values using Universal timezone.

Astronomical seasons

Source code:

    year = 2013;
    printf("\nAstronomical seasons %d\n\n", year);
    printf("%5s %10s\n", "Day", "JDN");

    for (i = julian_jdn(1, 1, year); i < julian_jdn(1, 1, year + 1); i++) {
        if (sun_is_winter(i) == 0) { /* -1 - no, 0 - last day, 1 - yes */
            printf("%5d %10d %20s\n", i + 1 - julian_jdn(1, 1, year) + 1, i + 1, "first day of spring");
        } else if (sun_is_spring(i) == 0) {
            printf("%5d %10d %20s\n", i + 1 - julian_jdn(1, 1, year) + 1, i + 1, "first day of summer");
        } else if (sun_is_summer(i) == 0) {
            printf("%5d %10d %20s\n", i + 1 - julian_jdn(1, 1, year) + 1, i + 1, "first day of autumn");
        } else if (sun_is_autumn(i) == 0) {
            printf("%5d %10d %20s\n", i + 1 - julian_jdn(1, 1, year) + 1, i + 1, "first day of winter");
        }
    }

Output:

  Day        JDN
   81    2456374  first day of spring
  173    2456466  first day of summer
  267    2456560  first day of autumn
  357    2456650  first day of winter

This is true only for northern hemisphere.

This code prints first day after change.

For example sun_is_winter(jdn) function will return:

  • -1 - it is not winter
  • 0 - it is the last day of winter (time of change)
  • 1 - it is winter

Day and night

Source code:

   printf("\nDay and night\n\n");

    int date[2][3] = {
        {21, 5, 2013},
        {22, 5, 2013}
    };

    double time[4][2] = {
        { -0.5 , 0  },  /* 00:00 hour */
        { -0.25, 6  },  /* 06:00 hour */
        {  0   , 12 },  /* 12:00 hour */
        {  0.25, 18  }   /* 18:00 hour */
    };

    int day_and_night[360];

    int k, j;
    int h, r, s;

    printf("%17s %12s  %60s\n", "Date & time", "JDN", "Day and night");

    for (j = 0; j < (sizeof(date) / sizeof(date[0])); j++) {
        for (k = 0; k < (sizeof(time) / sizeof(time[0])); k++) {
            jdn = julian_jdn(date[j][0], date[j][1], date[j][2]);
            jdn += time[k][0];

            ta = sun_earth_true_anomaly(jdn);

            printf(" %02d/%02d/%4d %02.0f:00 %12.2f  ", date[j][0], date[j][1], date[j][2], time[k][1], jdn);

            for (i = -180; i <= 180; i++) {
                l     = sun_sun_ecliptic_longitude(ta);
                delta = sun_earth_declination(l);

                s = sun_earth_sideral_time(jdn, -i);
                r = sun_earth_right_ascension(l);
                h = sun_earth_hour_angle(s, r);

                day_and_night[i + 180] =  sun_earth_altitude(h, delta, ln);
            }

            for (i = 0; i < 60; i++) {
                temp = (day_and_night[6 * i] +
                        day_and_night[6 * i + 1] +
                        day_and_night[6 * i + 2] +
                        day_and_night[6 * i + 3] +
                        day_and_night[6 * i + 4] +
                        day_and_night[6 * i + 5]) / 6;
                temp > 0 ? printf("+") : printf("-");
            }

            printf("\n");
        }
    }

Output:

Day and night

      Date & time          JDN                                                 Day and night
 21/05/2013 00:00   2456433.50  ++++++++++++++++++++--------------------++++++++++++++++++++
 21/05/2013 06:00   2456433.75  +++++--------------------+++++++++++++++++++++++++++++++++++
 21/05/2013 12:00   2456434.00  ----------++++++++++++++++++++++++++++++++++++++++----------
 21/05/2013 18:00   2456434.25  +++++++++++++++++++++++++++++++++++--------------------+++++
 22/05/2013 00:00   2456434.50  ++++++++++++++++++++--------------------++++++++++++++++++++
 22/05/2013 06:00   2456434.75  +++++--------------------+++++++++++++++++++++++++++++++++++
 22/05/2013 12:00   2456435.00  ----------++++++++++++++++++++++++++++++++++++++++----------
 22/05/2013 18:00   2456435.25  +++++++++++++++++++++++++++++++++++--------------------+++++

Notes

Each +/- character is equal to six degrees as there is no need to use full resolution in this example (it wouldn't be readable).

This code does not take timezone and daylight saving time into consideration.

Additional notes

If you want to play with the source code then download sun_example.tar.gz file.

Milosz Galazka's Picture

About Milosz Galazka

Milosz is a Linux Foundation Certified Engineer working for a successful Polish company as a system administrator and a long time supporter of Free Software Foundation and Debian operating system.

Gdansk, Poland https://sleeplessbeastie.eu