/*
  hora - An offline planetary hours calculator.

  Written by cidoku <cidoku.net>

  To the extent possible under law, the author(s) have dedicated all
  copyright and related and neighboring rights to this software to the public
  domain worldwide. This software is distributed without any warranty.

  You should have received a copy of the CC0 Public Domain Dedication along
  with this software. If not, see
  <http://creativecommons.org/publicdomain/zero/1.0/>.
*/

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

#define LATITUDE  40.712778
#define LONGITUDE -75.006111
#define ZENITH    90.833333

const char *planets[] = {"Sun\t", "Venus\t", "Mercury", "Moon\t",
  "Saturn", "Jupiter", "Mars\t"};
const char *symbols[] = {"☉", "♀", "☿", "☽", "♄", "♃", "♂"};
const char help[] =
  "Usage: %s [-option] [latitude] [longitude] [time_zone]\n"
  "Options:\n"
  "  -:  Displays the full table but allows you to set latitude, etc."
  "  -r: Display only the current day and hour rulers.\n"
  "  -R: Like -r, but the day ruler is replaced with the night ruler after"
  " sunset.\n"
  "  -d: Like -r, but the number of the hour is shown alongside its ruler.\n"
  "  -D: Like -d, but the day ruler is replaced with the night ruler after\n"
  "      sunset.\n"
  "  -h: Show this help.\n"
  ;

#define DEGTORAD(x) x * M_PI / 180.0
#define RADTODEG(x) x * 180.0 / M_PI
#define COSD(x)   cos(DEGTORAD(x))
#define SIND(x)   sin(DEGTORAD(x))
#define TAND(x)   tan(DEGTORAD(x))
#define ASIND(x)  RADTODEG(asin(x))
#define ACOSD(x)  RADTODEG(acos(x))
#define ATAND(x)  RADTODEG(atan(x))
#define ANOMALY(x)  (0.9856 * x) - 3.289
#define REALLONG(x) x + (1.916 * SIND(x)) + (0.020 * SIND(2 * x)) + 282.634
#define RIGHTASC(x) ATAND(0.91764 * TAND(x))
#define QUAD(x)   (floor(x / 90)) * 90
#define SINDEC(x) 0.39782 * SIND(x)
#define COSASIN(x)  COSD(ASIND(x))
#define COSH(lat, sinD, cosD)   (COSD(ZENITH) - (sinD * SIND(lat))) / (cosD * COSD(lat))
#define SUNCALC(H, RA, t, lng)    H + RA - (0.06571 * t) - 6.622 - lng

double adjust_to_range(double value, double l_limit, double r_limit) {
  if (value >= r_limit)
    return value - r_limit;
  else if (value < l_limit)
    return value + r_limit;
  return value;
}

// Almanac for Computers, 1990
// published by Nautical Almanac Office
// United States Naval Observatory
// Washington, DC 20392
// https://edwilliams.org/sunrise_sunset_algorithm.htm
void calculate_sunrise_sunset(double latitude, double longitude,
  int day, int month, int year, double *sunrise, double *sunset) {

  //first calculate the day of the year
  double N1 = floor(275 * month / 9);
  double N2 = floor((month + 9) / 12);
  double N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
  double N = N1 - (N2 * N3) + day - 30;

  //convert the longitude to hour value and calculate an approximate time
  double lng_hour = longitude / 15;
  double t_R = N + ((6 - lng_hour) / 24);
  double t_S = N + ((18 - lng_hour) / 24);

  //calculate the Sun's mean anomaly
  double M_R = ANOMALY(t_R);
  double M_S = ANOMALY(t_S);

  //calculate the Sun's true longitude
  double L_R = adjust_to_range(REALLONG(M_R), 0, 360);
  double L_S = adjust_to_range(REALLONG(M_S), 0, 360);

  //calculate the Sun's right ascension
  double RA_R = adjust_to_range(RIGHTASC(L_R), 0, 360);
  double RA_S = adjust_to_range(RIGHTASC(L_S), 0, 360);

  //right ascension value needs to be in the same quadrant as L
  double Lquadrant_R  = QUAD(L_R);
  double Lquadrant_S  = QUAD(L_S);
  double RAquadrant_R = QUAD(RA_R);
  double RAquadrant_S = QUAD(RA_S);
  RA_R += (Lquadrant_R - RAquadrant_R);
  RA_S += (Lquadrant_S - RAquadrant_S);

  //right ascension value needs to be converted into hours
  RA_R /= 15;
  RA_S /= 15;

  //calculate the Sun's declination
  double sinDec_R = SINDEC(L_R);
  double sinDec_S = SINDEC(L_S);
  double cosDec_R = COSASIN(sinDec_R);
  double cosDec_S = COSASIN(sinDec_S);

  //calculate the Sun's local hour angle
  double cosH_R = COSH(latitude, sinDec_R, cosDec_R);
  double cosH_S = COSH(latitude, sinDec_S, cosDec_S);

  //finish calculating H and convert into hours
  double H_R = 360 - ACOSD(cosH_R);
  double H_S = ACOSD(cosH_S);
  H_R /= 15;
  H_S /= 15;

  //calculate local mean time of rising/setting and adjust back to UTC
  *sunrise = adjust_to_range(SUNCALC(H_R, RA_R, t_R, lng_hour), 0, 24);
  *sunset = adjust_to_range(SUNCALC(H_S, RA_S, t_S, lng_hour), 0, 24);
}

int main(int argc, char *argv[]) {
  double latitude, longitude;
  int time_zone = 0, gmt_off;
  int opt = argc >= 2 && argv[1][0] == '-' ;

  if (argc == 1 + opt) {
    latitude = LATITUDE;
    longitude = LONGITUDE;
  } else if (argc == 3 + opt || argc == 4 + opt) {
    latitude = atof(argv[1 + opt]);
    longitude = atof(argv[2 + opt]);
    if (argc == 4 + opt)
      time_zone = atoi(argv[3 + opt]);
  } else {
    printf(help, argv[0]);
    return 1;
  }

  // Get current time and date
  time_t current_time = time(NULL);
  struct tm *local_time = localtime(&current_time);
  int hour = local_time->tm_hour;
  double current_hour = hour + (double)local_time->tm_min / 60;
  int day = local_time->tm_mday;
  int month = local_time->tm_mon + 1;
  int year = local_time->tm_year + 1900;
  int planetary_day = local_time->tm_wday;
  gmt_off = hour - gmtime(&current_time)->tm_hour;
  if (argc != 4 + opt)
    time_zone = gmt_off;
  else
    current_hour += time_zone - gmt_off;

  // Get times of sunrises and sunsets and adjust to local time zone
  double sunrise, sunset, next_sunrise, next_sunset;
  calculate_sunrise_sunset(latitude, longitude, day, month, year,
    &sunrise, &sunset);
  calculate_sunrise_sunset(latitude, longitude, day + 1, month, year,
    &next_sunrise, &next_sunset);
  sunrise = adjust_to_range(sunrise + time_zone, 0, 24);
  sunset = adjust_to_range(sunset + time_zone, 0, 24);
  next_sunrise = adjust_to_range(next_sunrise + time_zone, 0, 24);

  int day_ruler, planetary_hour, hour_ruler;
  double day_length, night_length, day_hour, night_hour, hour_start,
    hour_end, prev_sunrise, prev_sunset;

  // Figure out day and night lengths
  if (current_hour < sunrise) {
    calculate_sunrise_sunset(latitude, longitude, day - 1, month, year,
      &prev_sunrise, &prev_sunset);
    prev_sunrise = adjust_to_range(prev_sunrise + time_zone, 0, 24);
    prev_sunset = adjust_to_range(prev_sunset + time_zone, 0, 24);

    day_length = prev_sunset - prev_sunrise;
    night_length = sunrise + 24  - prev_sunset;
  } else {
    day_length = sunset - sunrise;
    night_length = next_sunrise + 24 - sunset;
  }
  day_hour = day_length / 12;
  night_hour = night_length / 12;

  // Figure out rulerships depending on current hour
  if (current_hour < sunset && current_hour >= sunrise) {
    current_hour = current_hour - sunrise;
    planetary_hour = (int)floor(current_hour / day_hour);
    hour_start = sunrise + day_hour * planetary_hour;
    hour_end = hour_start + day_hour;
  } else {
    if (current_hour < sunrise) {
      day_ruler = (planetary_day + 6) * 24 % 7;
      planetary_day = (planetary_day + 6) % 7;
      current_hour = current_hour + 24 - prev_sunset;
      hour_start = prev_sunset;
    } else {
      current_hour = current_hour - sunset;
      hour_start = sunset;
    }

    planetary_hour = (int)floor(current_hour / night_hour) + 12;
    hour_start += night_hour * (planetary_hour - 12);
    hour_end = hour_start + night_hour;
  }
  day_ruler = planetary_day * 24 % 7;
  hour_ruler = (planetary_day * 24 + planetary_hour) % 7;

  // Show rulerships
  if (opt) {
    switch (argv[1][1]) {
      case 'r':
        printf("%s %s\n", symbols[day_ruler], symbols[hour_ruler]);
        break;
      case 'R':
        printf("%s %s\n", symbols[planetary_hour < 12 ? day_ruler
            : (day_ruler + 12) % 7],
            symbols[hour_ruler]);
        break;
      case 'd':
        printf("%s %d%s\n", symbols[day_ruler], planetary_hour + 1,
            symbols[hour_ruler]);
        break;
      case 'D':
        printf("%s %d%s\n", symbols[planetary_hour < 12 ? day_ruler
            : (day_ruler + 12) % 7], planetary_hour + 1,
            symbols[hour_ruler]);
        break;
      case 'h':
        printf(help, argv[0]);
        break;
      default:
        goto table;
    }
    return 0;
  }

table:
  ;

  int loop = 0;
  double f_h_end, f_h_start;
  for (; planetary_hour < 24; planetary_hour++) {
    if (loop == 1) {
      printf("Next hours:\n");
      loop = -1;
    }
    f_h_start = floor(hour_start);
    f_h_end = floor(hour_end);
    printf("Hour %02d:\t%s %s\t(%02d:%02d - %02d:%02d)\n",
      planetary_hour + 1, symbols[hour_ruler], planets[hour_ruler],
      (int)fmod(f_h_start, 24),
      (int)((hour_start - f_h_start) * 60),
      (int)fmod(f_h_end, 24),
      (int)((hour_end - f_h_end) * 60));
    if (loop == 0) {
      int night_ruler = (day_ruler + 12) % 7;
      printf("Day:\t\t%s %s\t%s\n",
        symbols[day_ruler], planets[day_ruler],
        day_ruler == (int)hour_ruler ? "(*)" : "");
      if (planetary_hour >= 12)
        printf("Night:\t\t%s %s\t%s\n\n",
          symbols[night_ruler], planets[night_ruler],
          night_ruler == (int)hour_ruler ? "(*)" : "");
      else
        printf("\n");
      loop = 1;
    }
    hour_ruler = (hour_ruler + 1) % 7;
    hour_start += planetary_hour < 12 ? day_hour : night_hour;
    hour_end += planetary_hour < 11 ? day_hour : night_hour;
  }
  printf("\n");

  // Display other data
  double f_sunr, f_suns, f_d_len, f_n_len, f_d_hour, f_n_hour;
  f_sunr = floor(sunrise);
  f_suns = floor(sunset);
  f_d_len = floor(day_length);
  f_n_len = floor(night_length);
  f_d_hour = floor(day_hour);
  f_n_hour = floor(night_hour);
  printf("Sunrise:\t%02d:%02d\t",
  (int)f_sunr,
  (int)((sunrise - f_sunr) * 60));
  printf("Sunset:\t\t%02d:%02d\t\n",
  (int)f_suns,
  (int)((sunset - f_suns) * 60));
  printf("Day length:\t%02d:%02d\t",
  (int)f_d_len,
  (int)((day_length - f_d_len) * 60));
  printf("Night length:\t%02d:%02d\n",
  (int)f_n_len,
  (int)((night_length - f_n_len) * 60));
  printf("Day hour:\t%02d:%02d\t",
  (int)f_d_hour,
  (int)((day_hour - f_d_hour) * 60));
  printf("Night hour:\t%02d:%02d\n",
  (int)f_n_hour,
  (int)((night_hour - f_n_hour) * 60));

  return 0;
}