A couple of days ago, I mentioned calculating the sunrise and sunset times. This 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:09Astronomical 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), 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 the Julian Day fraction part when the 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 the northern hemisphere.
This code prints the first day after the change.
For example sun_is_winter(jdn) function will return:
-1– it is not winter- `` – 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
Download sun_example.tar.gz file and play with the source code.