daytcalc2.c

/* ------------------------------------------------------------ *
 * file:        daytcalc2.c                                     *
 *                                                              *
 * purpose:     Calculate local sunrise and sunset times.       *
 *              Returns 1 for nighttime, 0 for daytime, or      *
 *              -1 for any errors.                              *
 *                                                              *
 * author:      05/11/2017 Frank4DD                             *
 *                                                              *
 * compile:     gcc daytcalc2.c -o daytcalc2 -lm                *
 *                                                              *
 * example run: fm@susie:~$ ./daytcalc2 1486784589 -v           *
 * 2017-02-05 DST: 0 Sunrise: 6:38 Sunset: 17:11 Duration: 10:33*
 * SunriseTS: 1486244280 SunsetTS: 1486282260 DurationSec: 37980*
 *                                                              *
 * example run: fm@susie:~$ ./daytcalc 1486784589 -f            *
 * date=2017-02-10 sunrise=6:33  sunset=17:18 daytime=10:33     *
 * ------------------------------------------------------------ */
/* http://stackoverflow.com/questions/7064531/sunrise-sunset-times-in-c */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#include <ctype.h>
#include <stdbool.h>

#define PI 3.1415926
#define ZENITH -.83

/* ------------------------------------------------------------ *
 * print_usage() prints the programs commandline instructions.  *
 * ------------------------------------------------------------ */
void usage() {
   static char const usage[] = "Usage: daytcalc -t timestamp -x longitude -y latitude -z offset -d -f\n\n\
Command line parameters have the following format:\n\
   -t   Unix timestamp, example: 1486784589, optional, defaults to now\n\
   -x   longitude, example: 12.45277778\n\
   -y   latitude, example: 51.340277778\n\
   -z   timezone, example: 9, optional, defaults to 0\n\
   -d   daylight savings time flag\n\
   -f   output text for redirect into file\n\
   -v   verbose output flag\n\
   -h   print usage flag\n\n\
Usage example:\n\
./daytcalc -t 1486784589 -x 12.45277778 -y 51.340277778 -z 1 -d -f\n";
   printf(usage);
}

bool is_float(const char *s, float *dest) {
  if (s == NULL) return false;
  char *endptr;
  *dest = (float) strtod(s, &endptr);
  if (s == endptr) return false; // no conversion

  // Look at trailing text
  while (isspace((unsigned char ) *endptr))
    endptr++;
  return *endptr == '\0';
}

/* ------------------------------------------------------------ *
 * Global variables and defaults                                *
 * ------------------------------------------------------------ */
int verbose = 1;
time_t calc_t = 0;
float latitude = 0;
float longitude = 0;
int tzoffset = 0;
int dst = 0;
int txt = 0;

extern char *optarg;
extern int optind, opterr, optopt;

/* ------------------------------------------------------------ *
 * parseargs() checks the commandline arguments with C getopt   *
 * ------------------------------------------------------------ */
void parseargs(int argc, char* argv[]) {
   int arg, index;
   opterr = 0;

   if(argc == 1) { usage(); exit(-1); }

   while ((arg = (int) getopt (argc, argv, "t:x:y:z:dvhf")) != -1)
      switch (arg) {
         // arg -t timestamp, type: time_t, example: 1486784589
         // optional, defaults to now
         case 't':
            if(verbose == 1) printf("arg t, value %s\n", optarg);
            calc_t = (time_t) atoll(optarg);
            if (calc_t < 1) {
               printf("Error: Cannot get valid -t timestamp argument.\n");
               exit(-1);
            }
            break;

         // arg -n longitude, type: float, example: 12.45277777777777
         // optional, defaults to prime meridian 0�0'5.3?W (0.001389)
         // One meter resolution can be represented using 5 decimal places
         case 'x':
            //longitude = atof(optarg);
            if(is_float(optarg, &longitude)) {
               if(verbose == 1) printf("arg x, string %s, value (float) %e\n", optarg, longitude);
            } else {
               printf("Error: Cannot get valid -y latitude argument.\n");
               exit(-1);
            }
            if(longitude < -180.0f || longitude > 180.0f) {
               printf("Error: longitude value %e is out of range (< -180 or > 180).\n", longitude);
               exit(-1);
            }
            break;

         // arg -n latitude, type: float, example: 51.34027777777778
         // optional, defaults to prime meridian 51�28'40.1?N (51.477778)
         // One meter resolution can be represented using 5 decimal places
         case 'y':
            //latitude = atof(optarg);
            if(is_float(optarg, &latitude)) {
               if(verbose == 1) printf("arg y, string %s, value (float) %e\n", optarg, latitude);
            } else {
               printf("Error: Cannot get valid -y latitude argument.\n");
               exit(-1);
            }
            if(latitude < -90.0f || latitude > 90.0f) {
               printf("Error: latitude value %e is out of range (< -90 or > 90).\n", latitude);
               exit(-1);
            }
            break;

         // arg -z timezone, type: integer, example: 9
         // optional, defaults to 0
         case 'z':
            if(verbose == 1) printf("arg z, value %s\n", optarg);
            tzoffset = atoi(optarg);
            if (tzoffset < -11 || tzoffset > 11) {
               printf("Error: Cannot get valid -z timezone offset argument.\n");
               exit(-1);
            }
            break;

         // arg -d daylight savings time, type: integer, example: 1
         // optional, defaults to 0
         case 'd':
            if(verbose == 1) printf("arg d, setting DST\n");
            dst = 1; break;

         // arg -v verbose, type: flag, optional
         case 'v':
            verbose = 1; break;

         // arg -h usage, type: flag, optional
         case 'h':
            usage(); break;

         // arg -f outfile, type: char[255], example: /path/to/file
         // optional; defaults to /tmp/daytcalc.txt
         case 'f':
            if(verbose == 1) printf("arg f, creating text\n");
            txt = 1; break;

         case '?':
            if (isprint (optopt))
               fprintf (stderr, "Unknown option `-%c'.\n", optopt);
            else
               fprintf (stderr, "Unknown option character `\\x%x'.\n", optopt);

         default:
            usage();
    }
}

float calculateSunrise(int year,int month,int day,float lat, float lng,int localOffset, int daylightSavings) {
   //1. first calculate the day of the year
   float N1 = floor(275 * month / 9);
   float N2 = floor((month + 9) / 12);
   float N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
   float N = N1 - (N2 * N3) + day - 30;

   //2. convert the longitude to hour value and calculate an approximate time
   float lngHour = lng / 15.0;
   float t = N + ((6 - lngHour) / 24);   //if rising time is desired:

   //3. calculate the Sun's mean anomaly
   float M = (0.9856 * t) - 3.289;

   //4. calculate the Sun's true longitude
   float L = fmod(M + (1.916 * sin((PI/180)*M)) + (0.020 * sin(2 *(PI/180) * M)) + 282.634,360.0);

   //5a. calculate the Sun's right ascension
   float RA = fmod(180/PI*atan(0.91764 * tan((PI/180)*L)),360.0);

   //5b. right ascension value needs to be in the same quadrant as L
   float Lquadrant  = floor( L/90) * 90;
   float RAquadrant = floor(RA/90) * 90;
   RA = RA + (Lquadrant - RAquadrant);

   //5c. right ascension value needs to be converted into hours
   RA = RA / 15;

   //6. calculate the Sun's declination
   float sinDec = 0.39782 * sin((PI/180)*L);
   float cosDec = cos(asin(sinDec));

   //7a. calculate the Sun's local hour angle
   float cosH = (sin((PI/180)*ZENITH) - (sinDec * sin((PI/180)*lat))) / (cosDec * cos((PI/180)*lat));
   /*
   if (cosH >  1)
   the sun never rises on this location (on the specified date)
   if (cosH < -1)
   the sun never sets on this location (on the specified date)
   */

   //7b. finish calculating H and convert into hours
   float H = 360 - (180/PI)*acos(cosH);   //   if if rising time is desired:
   H = H / 15;

   //8. calculate local mean time of rising/setting
   float T = H + RA - (0.06571 * t) - 6.622;

   //9. adjust back to UTC
   float UT = fmod(T - lngHour,24.0);

   //10. convert UT value to local time zone of latitude/longitude
   return UT + localOffset + daylightSavings;
}

float calculateSunset(int year,int month,int day,float lat, float lng,int localOffset, int daylightSavings) {
   //1. first calculate the day of the year
   float N1 = floor(275 * month / 9);
   float N2 = floor((month + 9) / 12);
   float N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
   float N = N1 - (N2 * N3) + day - 30;

   //2. convert the longitude to hour value and calculate an approximate time
   float lngHour = lng / 15.0;
   float t = N + ((18 - lngHour) / 24);   // if setting time is desired:

   //3. calculate the Sun's mean anomaly
   float M = (0.9856 * t) - 3.289;

   //4. calculate the Sun's true longitude
   float L = fmod(M + (1.916 * sin((PI/180)*M)) + (0.020 * sin(2 *(PI/180) * M)) + 282.634,360.0);

   //5a. calculate the Sun's right ascension
   float RA = fmod(180/PI*atan(0.91764 * tan((PI/180)*L)),360.0);

   //5b. right ascension value needs to be in the same quadrant as L
   float Lquadrant  = floor( L/90) * 90;
   float RAquadrant = floor(RA/90) * 90;
   RA = RA + (Lquadrant - RAquadrant);

   //5c. right ascension value needs to be converted into hours
   RA = RA / 15;

   //6. calculate the Sun's declination
   float sinDec = 0.39782 * sin((PI/180)*L);
   float cosDec = cos(asin(sinDec));

   //7a. calculate the Sun's local hour angle
   float cosH = (sin((PI/180)*ZENITH) - (sinDec * sin((PI/180)*lat))) / (cosDec * cos((PI/180)*lat));
   /*
   if (cosH >  1)
   the sun never rises on this location (on the specified date)
   if (cosH < -1)
   the sun never sets on this location (on the specified date)
   */

   //7b. finish calculating H and convert into hours
   float H = (180/PI)*acos(cosH);    // if setting time is desired:
   H = H / 15;

   //8. calculate local mean time of rising/setting
   float T = H + RA - (0.06571 * t) - 6.622;

   //9. adjust back to UTC
   float UT = fmod(T - lngHour,24.0);

   //10. convert UT value to local time zone of latitude/longitude
   return UT + localOffset + daylightSavings;
}

int main(int argc, char* argv[]) {

   parseargs(argc, argv);

   /* We convert the timestamp string to a local time based struct */
   struct tm *calc_tm = gmtime(&calc_t);
   int year = calc_tm->tm_year+1900;
   int month = calc_tm->tm_mon+1;
   int day = calc_tm->tm_mday;

   if(verbose == 1) {
      printf("Calc TS UTC: %ld Date: %d-%d-%d Time: %d:%.2d\n",
        (long) calc_t, year, month, day,  calc_tm->tm_hour, calc_tm->tm_min);
   }

   /* ------------------------------------------------------------ *
    * calculation argument list:                                   *
    * --------------------------                                   *
    * int year, int month, int day, float lat, loat lng,           *
    * int localOffset, int dst                                     *
    *                                                              *
    * argument description:                                        *
    * ---------------------                                        *
    * year, month, day: The date to calculate sunrise/sunset for.  *
    * lat, lng: the locations latitude and longitude, in decimal.  *
    * localOffset: the timezone offset. It is <0 for the western   *
    * hemisphere; and >0 for eastern hemisphere locations.         *
    * dst: daylightSavings, 1 if it is in effect, otherwise set 0. *
    * ------------------------------------------------------------ */
   float localsunriseT = calculateSunrise(year, month, day, latitude, longitude, tzoffset, dst);
   double sunrise_hr = fmod(24 + localsunriseT,24.0);
   double sunrise_min = modf(fmod(24+localsunriseT,24.0),&sunrise_hr)*60;

   float localsunsetT = calculateSunset(year, month, day, latitude, longitude, tzoffset, dst);
   double sunset_hr = fmod(24 + localsunsetT,24.0);
   double sunset_min = modf(fmod(24+localsunsetT,24.0),&sunset_hr)*60;

   struct tm sunrise_tm;
   sunrise_tm.tm_year = year-1900;
   sunrise_tm.tm_mon = month-1;
   sunrise_tm.tm_mday = day;
   sunrise_tm.tm_hour = sunrise_hr;
   sunrise_tm.tm_min = sunrise_min;
   sunrise_tm.tm_sec = 0;
   sunrise_tm.tm_isdst = dst;

   time_t sunrise = mktime(&sunrise_tm);
   if(sunrise == -1) printf("Error");

   struct tm sunset_tm;
   sunset_tm.tm_year = year-1900;
   sunset_tm.tm_mon = month-1;
   sunset_tm.tm_mday = day;
   sunset_tm.tm_hour = sunset_hr;
   sunset_tm.tm_min = sunset_min+1;
   sunset_tm.tm_sec = 0;
   sunset_tm.tm_isdst = dst;

   time_t sunset = mktime(&sunset_tm);
   if(sunset == -1) printf("Error");

   int daytimeflag;
   if(calc_t < sunrise) daytimeflag = 1;
   else if(calc_t > sunset) daytimeflag = 1;
   else daytimeflag = 0;

   long daytime = sunset-sunrise;
   int daytime_hr = daytime / 3600;
   int daytime_min = (daytime % 3600) / 60;


   /* ------------------------------------------------------------ *
    * The calculation for sunrise_hr/min came up with 5:60, so we  *
    * use the timestamp-converted time data which fixes the issue. *
    * ------------------------------------------------------------ */
   char rise[6];
   strftime(rise, sizeof(rise), "%H:%M", &sunrise_tm);
   char sset[6];
   strftime(sset, sizeof(sset), "%H:%M", &sunset_tm);

   if(verbose == 1) {
      printf("SunriseT: %.0f:%.0f SunsetT: %.0f:%.0f\n", sunrise_hr, sunrise_min, sunset_hr, sunset_min);
      printf("Sunrise: %s", ctime(&sunrise));
       printf("Sunset: %s", ctime(&sunset));
      printf("Calc TS: %ld SunriseTS: %ld SunsetTS: %ld\n", calc_t, sunrise, sunset);
      printf("RRD return value: %d\n", daytimeflag);
   }

   if(txt == 1) {
      printf("date=%d-%.2d-%.2d sunrise=%s sunset=%s daytime=%.2d:%.2d\n",
              year, month, day, rise, sset, daytime_hr, daytime_min);
   }
   exit(daytimeflag);
}