/* https://gml.noaa.gov/grad/solcalc/solareqns.PDF
*
* Bad time offset equation - corrections sent to NOAA
*
* gcc -Wall -Wextra -pedantic -Wshadow -Werror -std=c11 -Ofast -o sunrise-sunset sunrise-sunset.c -lm
*
* provide lattitude longitude as first two arguments to program
* (uses Nacogdoches lat/lon if none provided)
* Nacogdoches, Addison Airport (KADS) and San Gabriel in Austin, shown below.
* Or use iPhone Maps and look at the pin info to get exact current lat/lon.
*
* NOTE: longitude west of the Prime Meridian is *Negative* (that's the US)
*
*/
#define _POSIX_C_SOURCE 2
#include <stdio.h>
#include <time.h>
#include <math.h>
/* Nacogdoches lat/lon */
#define NACLAT 31.62965
#define NACLON -94.62127
/* KADS 32.97 -96.83
* San Gabriel Austin 30.28 -97.74
*/
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif
typedef struct {
double lat, /* lattitude (deg) */
lon, /* longitude (deg) */
fracyr, /* fractional year (rad) */
eqtime, /* equation time */
decl, /* declination angle */
toff, /* time offset (min) */
tst, /* true solar time (min) */
ha, /* solar hour angle (deg) */
sza, /* solar zenith angle */
sa, /* solar azimuth */
sr, /* sunrise (min) */
ss, /* sunset (min) */
sn; /* solar noon */
int tzone; /* local timezone offset (+ East of Prime Meridian) */
} solar_position;
typedef struct {
int h, m ,s;
} h_m_s;
/* get nowtm in broken-down time GMT */
struct tm *gmt_tm (struct tm *nowtm)
{
time_t now;
time (&now);
return gmtime_r (&now, nowtm);
}
/* convert radians to degrees */
double r2deg (const double r) {
return r * 180.0 / M_PI;
}
/* convert degrees to radians */
double deg2r (const double d) {
return d * M_PI / 180.0;;
}
/* is year leapyear */
int is_leap_year (int year)
{
return (year % 4 == 0 && year % 100 != 0) || (year % 400 == 0);
}
/* convert seconds to h, m, s */
h_m_s sec2hms (double s)
{
h_m_s hms = { .h = s / 3600 }; /* initialize hours */
s -= hms.h * 3600; /* subtract hours no. of secs */
hms.m = s / 60; /* compute minutes */
hms.s = (int)s - hms.m * 60; /* subtract min no. of secs for seconds */
return hms; /* return struct */
}
void test (struct tm ntm, solar_position *sp)
{
printf ("Current Date/Time (UTC) : %02d/%02d/%d - %02d:%02d:%02d\n"
"Day of week : %3d (0 = Sunday)\n"
"Day of year : %3d (0 - 365, 366 for leapyear)\n"
"Daylight Savings Time : %3d (0 - not in effect)\n"
"Timezone offset : %3d (hours)\n\n",
ntm.tm_mon + 1, ntm.tm_mday, ntm.tm_year + 1900,
ntm.tm_hour, ntm.tm_min, ntm.tm_sec,
ntm.tm_wday, ntm.tm_yday + 1, ntm.tm_isdst, sp->tzone);
printf ("Lattitude: % g\nLongitude: % g\n\n", sp->lat, sp->lon);
printf (" sp->fy : % 12.6f fractional year\n"
" sp->et : % 12.6f equation of time (min)\n"
" sp->dcl : % 12.6f declination (radians)\n"
" sp->ha : % 12.6f hour angle (radians)\n\n"
" sp->sr : % 12.6f (UTC min)\n"
" sp->ss : % 12.6f (UTC min)\n"
" sp->sn : % 12.6f (UTC min)\n\n"
" len day : % 12.6f (daylight min)\n\n",
sp->fracyr, sp->eqtime, r2deg(sp->decl), sp->ha,
sp->sr, sp->ss, sp->sn, sp->ss - sp->sr);
}
int main (int argc, char **argv) {
/* general solar position calculations */
solar_position sp = { .lat = argc > 1 ? 0. : NACLAT,
.lon = argc > 2 ? 0. : NACLON };
struct tm ntm = { .tm_sec = 0 };
sp.tzone = -6; /* default Central Standard Time zone offset */
if (argc > 1 && sscanf (argv[1], "%lf", &sp.lat) != 1) {
fputs ("error: invalid lattitude provided as first-argument.\n", stderr);
return 1;
}
if (argc > 2 && sscanf (argv[2], "%lf", &sp.lon) != 1) {
fputs ("error: invalid longitude provided as second-argument.\n", stderr);
return 1;
}
if (argc > 3 && sscanf (argv[3], "%d", &sp.tzone) != 1) {
fputs ("error: invalid timezone offset provided as third-argument.\n",
stderr);
return 1;
}
if (gmt_tm (&ntm) == NULL) {
fputs ("error: gmtime_r failed.\n", stderr);
return 1;
}
ntm.tm_isdst = 0;
/* = 2 * M_PI / 365 * (day_of_year - 1 + (hour - 12) / 24
* so for struct tm, ntm.tm_mon + 1 - 1 => ntm.tm_mon
*/
sp.fracyr = 2 * M_PI / (365. + is_leap_year (ntm.tm_year + 1900)) *
(ntm.tm_yday + (ntm.tm_hour - 12) / 24.);
/* estimate the equation of time (in minutes) */
sp.eqtime = 229.18 * (0.000075 + 0.001868 * cos(sp.fracyr) -
0.032077 * sin(sp.fracyr) - 0.014615 * cos(2 * sp.fracyr) -
0.040849 * sin(2 * sp.fracyr));
/* estimate the solar declination angle (in radians) */
sp.decl = 0.006918 - 0.399912 * cos(sp.fracyr) + 0.070257 * sin(sp.fracyr) -
0.006758 * cos(2 * sp.fracyr) + 0.000907 * sin(2 * sp.fracyr) -
0.002697 * cos(3 * sp.fracyr) + 0.00148 * sin(3 * sp.fracyr);
/* ha for special case of sunrise/sunset.
* zenith is set to 90.833 (the approximate correction for atmospheric
* refraction at sunrise and sunset, and the size of the solar disk)
*/
sp.ha = acos (cos(deg2r (90.833)) / cos(deg2r (sp.lat)) * cos(sp.decl) -
tan(deg2r (sp.lat)) * tan(sp.decl));
/* UTC time of sunrise [+] (or sunset [-]) in minutes */
sp.sr = 720 - 4 * (sp.lon + r2deg (sp.ha)) - sp.eqtime;
sp.ss = 720 - 4 * (sp.lon - r2deg (sp.ha)) - sp.eqtime;
/* UTC time for solar noon in minutes */
sp.sn = 720 - 4 * sp.lon - sp.eqtime;
test (ntm, &sp);
printf (" sp.ha : % 12.6f (deg)\n\n", r2deg(sp.ha));
puts ("Localtime Minutes:\n");
/* local time (minutes) */
double sr = sp.sr + 60 * (sp.tzone + ntm.tm_isdst),
ss = sp.ss + 60 * (sp.tzone + ntm.tm_isdst),
sn = sp.sn + 60 * (sp.tzone + ntm.tm_isdst);
/* convert minutes to h, m, s */
h_m_s hms_sr = sec2hms (sr * 60),
hms_ss = sec2hms (ss * 60),
hms_sn = sec2hms (sn * 60);
printf (" sr.sp : % 12.6f (min)\n"
" sp.ss : % 12.6f (min)\n"
" sp.sn : % 12.6f (min)\n\n"
" sunrise : %02d:%02d:%02d\n"
" sunset : %02d:%02d:%02d\n"
" solor noon : %02d:%02d:%02d\n\n",
sr, ss, sn,
hms_sr.h, hms_sr.m, hms_sr.s,
hms_ss.h, hms_ss.m, hms_ss.s,
hms_sn.h, hms_sn.m, hms_sn.s);
}